1 #define PETSCSNES_DLL 2 3 #include "../src/snes/impls/vi/viimpl.h" /*I "petscsnes.h" I*/ 4 #include "../include/private/kspimpl.h" 5 6 #undef __FUNCT__ 7 #define __FUNCT__ "SNESMonitorVI" 8 PetscErrorCode SNESMonitorVI(SNES snes,PetscInt its,PetscReal fgnorm,void *dummy) 9 { 10 PetscErrorCode ierr; 11 SNES_VI *vi = (SNES_VI*)snes->data; 12 PetscViewerASCIIMonitor viewer = (PetscViewerASCIIMonitor) dummy; 13 const PetscScalar *x,*xl,*xu,*f; 14 PetscInt i,n,act = 0; 15 PetscReal rnorm,fnorm; 16 17 PetscFunctionBegin; 18 if (!dummy) { 19 ierr = PetscViewerASCIIMonitorCreate(((PetscObject)snes)->comm,"stdout",0,&viewer);CHKERRQ(ierr); 20 } 21 ierr = VecGetLocalSize(snes->vec_sol,&n);CHKERRQ(ierr); 22 ierr = VecGetArrayRead(vi->xl,&xl);CHKERRQ(ierr); 23 ierr = VecGetArrayRead(vi->xu,&xu);CHKERRQ(ierr); 24 ierr = VecGetArrayRead(snes->vec_sol,&x);CHKERRQ(ierr); 25 ierr = VecGetArrayRead(snes->vec_func,&f);CHKERRQ(ierr); 26 27 rnorm = 0.0; 28 for (i=0; i<n; i++) { 29 if (((x[i] > xl[i] + 1.e-8 || (f[i] < 0.0)) && ((x[i] < xu[i] - 1.e-8) || f[i] > 0.0))) rnorm += f[i]*f[i]; 30 else act++; 31 } 32 ierr = VecRestoreArrayRead(snes->vec_func,&f);CHKERRQ(ierr); 33 ierr = VecRestoreArrayRead(vi->xl,&xl);CHKERRQ(ierr); 34 ierr = VecRestoreArrayRead(vi->xu,&xu);CHKERRQ(ierr); 35 ierr = VecRestoreArrayRead(snes->vec_sol,&x);CHKERRQ(ierr); 36 ierr = MPI_Allreduce(&rnorm,&fnorm,1,MPIU_REAL,MPI_SUM,((PetscObject)snes)->comm);CHKERRQ(ierr); 37 fnorm = sqrt(fnorm); 38 ierr = PetscViewerASCIIMonitorPrintf(viewer,"%3D SNES VI Function norm %14.12e Active constraints %D\n",its,fnorm,act);CHKERRQ(ierr); 39 if (!dummy) { 40 ierr = PetscViewerASCIIMonitorDestroy(viewer);CHKERRQ(ierr); 41 } 42 PetscFunctionReturn(0); 43 } 44 45 /* 46 Checks if J^T F = 0 which implies we've found a local minimum of the norm of the function, 47 || F(u) ||_2 but not a zero, F(u) = 0. In the case when one cannot compute J^T F we use the fact that 48 0 = (J^T F)^T W = F^T J W iff W not in the null space of J. Thanks for Jorge More 49 for this trick. One assumes that the probability that W is in the null space of J is very, very small. 50 */ 51 #undef __FUNCT__ 52 #define __FUNCT__ "SNESVICheckLocalMin_Private" 53 PetscErrorCode SNESVICheckLocalMin_Private(SNES snes,Mat A,Vec F,Vec W,PetscReal fnorm,PetscBool *ismin) 54 { 55 PetscReal a1; 56 PetscErrorCode ierr; 57 PetscBool hastranspose; 58 59 PetscFunctionBegin; 60 *ismin = PETSC_FALSE; 61 ierr = MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);CHKERRQ(ierr); 62 if (hastranspose) { 63 /* Compute || J^T F|| */ 64 ierr = MatMultTranspose(A,F,W);CHKERRQ(ierr); 65 ierr = VecNorm(W,NORM_2,&a1);CHKERRQ(ierr); 66 ierr = PetscInfo1(snes,"|| J^T F|| %G near zero implies found a local minimum\n",a1/fnorm);CHKERRQ(ierr); 67 if (a1/fnorm < 1.e-4) *ismin = PETSC_TRUE; 68 } else { 69 Vec work; 70 PetscScalar result; 71 PetscReal wnorm; 72 73 ierr = VecSetRandom(W,PETSC_NULL);CHKERRQ(ierr); 74 ierr = VecNorm(W,NORM_2,&wnorm);CHKERRQ(ierr); 75 ierr = VecDuplicate(W,&work);CHKERRQ(ierr); 76 ierr = MatMult(A,W,work);CHKERRQ(ierr); 77 ierr = VecDot(F,work,&result);CHKERRQ(ierr); 78 ierr = VecDestroy(work);CHKERRQ(ierr); 79 a1 = PetscAbsScalar(result)/(fnorm*wnorm); 80 ierr = PetscInfo1(snes,"(F^T J random)/(|| F ||*||random|| %G near zero implies found a local minimum\n",a1);CHKERRQ(ierr); 81 if (a1 < 1.e-4) *ismin = PETSC_TRUE; 82 } 83 PetscFunctionReturn(0); 84 } 85 86 /* 87 Checks if J^T(F - J*X) = 0 88 */ 89 #undef __FUNCT__ 90 #define __FUNCT__ "SNESVICheckResidual_Private" 91 PetscErrorCode SNESVICheckResidual_Private(SNES snes,Mat A,Vec F,Vec X,Vec W1,Vec W2) 92 { 93 PetscReal a1,a2; 94 PetscErrorCode ierr; 95 PetscBool hastranspose; 96 97 PetscFunctionBegin; 98 ierr = MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);CHKERRQ(ierr); 99 if (hastranspose) { 100 ierr = MatMult(A,X,W1);CHKERRQ(ierr); 101 ierr = VecAXPY(W1,-1.0,F);CHKERRQ(ierr); 102 103 /* Compute || J^T W|| */ 104 ierr = MatMultTranspose(A,W1,W2);CHKERRQ(ierr); 105 ierr = VecNorm(W1,NORM_2,&a1);CHKERRQ(ierr); 106 ierr = VecNorm(W2,NORM_2,&a2);CHKERRQ(ierr); 107 if (a1 != 0.0) { 108 ierr = PetscInfo1(snes,"||J^T(F-Ax)||/||F-AX|| %G near zero implies inconsistent rhs\n",a2/a1);CHKERRQ(ierr); 109 } 110 } 111 PetscFunctionReturn(0); 112 } 113 114 /* 115 SNESDefaultConverged_VI - Checks the convergence of the semismooth newton algorithm. 116 117 Notes: 118 The convergence criterion currently implemented is 119 merit < abstol 120 merit < rtol*merit_initial 121 */ 122 #undef __FUNCT__ 123 #define __FUNCT__ "SNESDefaultConverged_VI" 124 PetscErrorCode SNESDefaultConverged_VI(SNES snes,PetscInt it,PetscReal xnorm,PetscReal gradnorm,PetscReal merit,SNESConvergedReason *reason,void *dummy) 125 { 126 PetscErrorCode ierr; 127 128 PetscFunctionBegin; 129 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 130 PetscValidPointer(reason,6); 131 132 *reason = SNES_CONVERGED_ITERATING; 133 134 if (!it) { 135 /* set parameter for default relative tolerance convergence test */ 136 snes->ttol = merit*snes->rtol; 137 } 138 if (merit != merit) { 139 ierr = PetscInfo(snes,"Failed to converged, function norm is NaN\n");CHKERRQ(ierr); 140 *reason = SNES_DIVERGED_FNORM_NAN; 141 } else if (merit < snes->abstol) { 142 ierr = PetscInfo2(snes,"Converged due to merit function %G < %G\n",merit,snes->abstol);CHKERRQ(ierr); 143 *reason = SNES_CONVERGED_FNORM_ABS; 144 } else if (snes->nfuncs >= snes->max_funcs) { 145 ierr = PetscInfo2(snes,"Exceeded maximum number of function evaluations: %D > %D\n",snes->nfuncs,snes->max_funcs);CHKERRQ(ierr); 146 *reason = SNES_DIVERGED_FUNCTION_COUNT; 147 } 148 149 if (it && !*reason) { 150 if (merit < snes->ttol) { 151 ierr = PetscInfo2(snes,"Converged due to merit function %G < %G (relative tolerance)\n",merit,snes->ttol);CHKERRQ(ierr); 152 *reason = SNES_CONVERGED_FNORM_RELATIVE; 153 } 154 } 155 PetscFunctionReturn(0); 156 } 157 158 /* 159 SNESVIComputeMeritFunction - Evaluates the merit function for the mixed complementarity problem. 160 161 Input Parameter: 162 . phi - the semismooth function 163 164 Output Parameter: 165 . merit - the merit function 166 . phinorm - ||phi|| 167 168 Notes: 169 The merit function for the mixed complementarity problem is defined as 170 merit = 0.5*phi^T*phi 171 */ 172 #undef __FUNCT__ 173 #define __FUNCT__ "SNESVIComputeMeritFunction" 174 static PetscErrorCode SNESVIComputeMeritFunction(Vec phi, PetscReal* merit,PetscReal* phinorm) 175 { 176 PetscErrorCode ierr; 177 178 PetscFunctionBegin; 179 ierr = VecNormBegin(phi,NORM_2,phinorm); 180 ierr = VecNormEnd(phi,NORM_2,phinorm); 181 182 *merit = 0.5*(*phinorm)*(*phinorm); 183 PetscFunctionReturn(0); 184 } 185 186 static inline PetscScalar Phi(PetscScalar a,PetscScalar b) 187 { 188 return a + b - sqrt(a*a + b*b); 189 } 190 191 static inline PetscScalar DPhi(PetscScalar a,PetscScalar b) 192 { 193 if ((PetscAbsScalar(a) >= 1.e-6) || (PetscAbsScalar(b) >= 1.e-6)) return 1.0 - a/ sqrt(a*a + b*b); 194 else return .5; 195 } 196 197 /* 198 SNESVIComputeFunction - Reformulates a system of nonlinear equations in mixed complementarity form to a system of nonlinear equations in semismooth form. 199 200 Input Parameters: 201 . snes - the SNES context 202 . x - current iterate 203 . functx - user defined function context 204 205 Output Parameters: 206 . phi - Semismooth function 207 208 */ 209 #undef __FUNCT__ 210 #define __FUNCT__ "SNESVIComputeFunction" 211 static PetscErrorCode SNESVIComputeFunction(SNES snes,Vec X,Vec phi,void* functx) 212 { 213 PetscErrorCode ierr; 214 SNES_VI *vi = (SNES_VI*)snes->data; 215 Vec Xl = vi->xl,Xu = vi->xu,F = snes->vec_func; 216 PetscScalar *phi_arr,*x_arr,*f_arr,*l,*u; 217 PetscInt i,nlocal; 218 219 PetscFunctionBegin; 220 ierr = (*vi->computeuserfunction)(snes,X,F,functx);CHKERRQ(ierr); 221 ierr = VecGetLocalSize(X,&nlocal);CHKERRQ(ierr); 222 ierr = VecGetArray(X,&x_arr);CHKERRQ(ierr); 223 ierr = VecGetArray(F,&f_arr);CHKERRQ(ierr); 224 ierr = VecGetArray(Xl,&l);CHKERRQ(ierr); 225 ierr = VecGetArray(Xu,&u);CHKERRQ(ierr); 226 ierr = VecGetArray(phi,&phi_arr);CHKERRQ(ierr); 227 228 for (i=0;i < nlocal;i++) { 229 if ((l[i] <= PETSC_VI_NINF) && (u[i] >= PETSC_VI_INF)) { /* no constraints on variable */ 230 phi_arr[i] = f_arr[i]; 231 } else if (l[i] <= PETSC_VI_NINF) { /* upper bound on variable only */ 232 phi_arr[i] = -Phi(u[i] - x_arr[i],-f_arr[i]); 233 } else if (u[i] >= PETSC_VI_INF) { /* lower bound on variable only */ 234 phi_arr[i] = Phi(x_arr[i] - l[i],f_arr[i]); 235 } else if (l[i] == u[i]) { 236 phi_arr[i] = l[i] - x_arr[i]; 237 } else { /* both bounds on variable */ 238 phi_arr[i] = Phi(x_arr[i] - l[i],-Phi(u[i] - x_arr[i],-f_arr[i])); 239 } 240 } 241 242 ierr = VecRestoreArray(X,&x_arr);CHKERRQ(ierr); 243 ierr = VecRestoreArray(F,&f_arr);CHKERRQ(ierr); 244 ierr = VecRestoreArray(Xl,&l);CHKERRQ(ierr); 245 ierr = VecRestoreArray(Xu,&u);CHKERRQ(ierr); 246 ierr = VecRestoreArray(phi,&phi_arr);CHKERRQ(ierr); 247 PetscFunctionReturn(0); 248 } 249 250 /* 251 SNESVIComputeBsubdifferentialVectors - Computes the diagonal shift (Da) and row scaling (Db) vectors needed for the 252 the semismooth jacobian. 253 */ 254 #undef __FUNCT__ 255 #define __FUNCT__ "SNESVIComputeBsubdifferentialVectors" 256 PetscErrorCode SNESVIComputeBsubdifferentialVectors(SNES snes,Vec X,Vec F,Mat jac,Vec Da,Vec Db) 257 { 258 PetscErrorCode ierr; 259 SNES_VI *vi = (SNES_VI*)snes->data; 260 PetscScalar *l,*u,*x,*f,*da,*db,da1,da2,db1,db2; 261 PetscInt i,nlocal; 262 263 PetscFunctionBegin; 264 265 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 266 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 267 ierr = VecGetArray(vi->xl,&l);CHKERRQ(ierr); 268 ierr = VecGetArray(vi->xu,&u);CHKERRQ(ierr); 269 ierr = VecGetArray(Da,&da);CHKERRQ(ierr); 270 ierr = VecGetArray(Db,&db);CHKERRQ(ierr); 271 ierr = VecGetLocalSize(X,&nlocal);CHKERRQ(ierr); 272 273 for (i=0;i< nlocal;i++) { 274 if ((l[i] <= PETSC_VI_NINF) && (u[i] >= PETSC_VI_INF)) {/* no constraints on variable */ 275 da[i] = 0; 276 db[i] = 1; 277 } else if (l[i] <= PETSC_VI_NINF) { /* upper bound on variable only */ 278 da[i] = DPhi(u[i] - x[i], -f[i]); 279 db[i] = DPhi(-f[i],u[i] - x[i]); 280 } else if (u[i] >= PETSC_VI_INF) { /* lower bound on variable only */ 281 da[i] = DPhi(x[i] - l[i], f[i]); 282 db[i] = DPhi(f[i],x[i] - l[i]); 283 } else if (l[i] == u[i]) { /* fixed variable */ 284 da[i] = 1; 285 db[i] = 0; 286 } else { /* upper and lower bounds on variable */ 287 da1 = DPhi(x[i] - l[i], -Phi(u[i] - x[i], -f[i])); 288 db1 = DPhi(-Phi(u[i] - x[i], -f[i]),x[i] - l[i]); 289 da2 = DPhi(u[i] - x[i], -f[i]); 290 db2 = DPhi(-f[i],u[i] - x[i]); 291 da[i] = da1 + db1*da2; 292 db[i] = db1*db2; 293 } 294 } 295 296 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr); 297 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 298 ierr = VecRestoreArray(vi->xl,&l);CHKERRQ(ierr); 299 ierr = VecRestoreArray(vi->xu,&u);CHKERRQ(ierr); 300 ierr = VecRestoreArray(Da,&da);CHKERRQ(ierr); 301 ierr = VecRestoreArray(Db,&db);CHKERRQ(ierr); 302 PetscFunctionReturn(0); 303 } 304 305 /* 306 SNESVIComputeJacobian - Computes the jacobian of the semismooth function.The Jacobian for the semismooth function is an element of the B-subdifferential of the Fischer-Burmeister function for complementarity problems. 307 308 Input Parameters: 309 . Da - Diagonal shift vector for the semismooth jacobian. 310 . Db - Row scaling vector for the semismooth jacobian. 311 312 Output Parameters: 313 . jac - semismooth jacobian 314 . jac_pre - optional preconditioning matrix 315 316 Notes: 317 The semismooth jacobian matrix is given by 318 jac = Da + Db*jacfun 319 where Db is the row scaling matrix stored as a vector, 320 Da is the diagonal perturbation matrix stored as a vector 321 and jacfun is the jacobian of the original nonlinear function. 322 */ 323 #undef __FUNCT__ 324 #define __FUNCT__ "SNESVIComputeJacobian" 325 PetscErrorCode SNESVIComputeJacobian(Mat jac, Mat jac_pre,Vec Da, Vec Db) 326 { 327 PetscErrorCode ierr; 328 329 /* Do row scaling and add diagonal perturbation */ 330 ierr = MatDiagonalScale(jac,Db,PETSC_NULL);CHKERRQ(ierr); 331 ierr = MatDiagonalSet(jac,Da,ADD_VALUES);CHKERRQ(ierr); 332 if (jac != jac_pre) { /* If jac and jac_pre are different */ 333 ierr = MatDiagonalScale(jac_pre,Db,PETSC_NULL); 334 ierr = MatDiagonalSet(jac_pre,Da,ADD_VALUES);CHKERRQ(ierr); 335 } 336 PetscFunctionReturn(0); 337 } 338 339 /* 340 SNESVIComputeMeritFunctionGradient - Computes the gradient of the merit function psi. 341 342 Input Parameters: 343 phi - semismooth function. 344 H - semismooth jacobian 345 346 Output Parameters: 347 dpsi - merit function gradient 348 349 Notes: 350 The merit function gradient is computed as follows 351 dpsi = H^T*phi 352 */ 353 #undef __FUNCT__ 354 #define __FUNCT__ "SNESVIComputeMeritFunctionGradient" 355 PetscErrorCode SNESVIComputeMeritFunctionGradient(Mat H, Vec phi, Vec dpsi) 356 { 357 PetscErrorCode ierr; 358 359 PetscFunctionBegin; 360 ierr = MatMultTranspose(H,phi,dpsi); 361 PetscFunctionReturn(0); 362 } 363 364 /* -------------------------------------------------------------------------- */ 365 /* 366 SNESVIProjectOntoBounds - Projects X onto the feasible region so that Xl[i] <= X[i] <= Xu[i] for i = 1...n. 367 368 Input Parameters: 369 . SNES - nonlinear solver context 370 371 Output Parameters: 372 . X - Bound projected X 373 374 */ 375 376 #undef __FUNCT__ 377 #define __FUNCT__ "SNESVIProjectOntoBounds" 378 PetscErrorCode SNESVIProjectOntoBounds(SNES snes,Vec X) 379 { 380 PetscErrorCode ierr; 381 SNES_VI *vi = (SNES_VI*)snes->data; 382 const PetscScalar *xl,*xu; 383 PetscScalar *x; 384 PetscInt i,n; 385 386 PetscFunctionBegin; 387 ierr = VecGetLocalSize(X,&n);CHKERRQ(ierr); 388 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 389 ierr = VecGetArrayRead(vi->xl,&xl);CHKERRQ(ierr); 390 ierr = VecGetArrayRead(vi->xu,&xu);CHKERRQ(ierr); 391 392 for(i = 0;i<n;i++) { 393 if (x[i] < xl[i]) x[i] = xl[i]; 394 else if (x[i] > xu[i]) x[i] = xu[i]; 395 } 396 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr); 397 ierr = VecRestoreArrayRead(vi->xl,&xl);CHKERRQ(ierr); 398 ierr = VecRestoreArrayRead(vi->xu,&xu);CHKERRQ(ierr); 399 PetscFunctionReturn(0); 400 } 401 402 403 /* 404 SNESVIAdjustInitialGuess - Readjusts the initial guess to the SNES solver supplied by the user so that the initial guess lies inside the feasible region . 405 406 Input Parameters: 407 . lb - lower bound. 408 . ub - upper bound. 409 410 Output Parameters: 411 . X - the readjusted initial guess. 412 413 Notes: 414 The readjusted initial guess X[i] = max(max(min(l[i],X[i]),min(X[i],u[i])),min(l[i],u[i])) 415 */ 416 #undef __FUNCT__ 417 #define __FUNCT__ "SNESVIAdjustInitialGuess" 418 PetscErrorCode SNESVIAdjustInitialGuess(Vec X, Vec lb, Vec ub) 419 { 420 PetscErrorCode ierr; 421 PetscInt i,nlocal; 422 PetscScalar *x,*l,*u; 423 424 PetscFunctionBegin; 425 ierr = VecGetLocalSize(X,&nlocal);CHKERRQ(ierr); 426 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 427 ierr = VecGetArray(lb,&l);CHKERRQ(ierr); 428 ierr = VecGetArray(ub,&u);CHKERRQ(ierr); 429 430 for(i = 0; i < nlocal; i++) { 431 x[i] = PetscMax(PetscMax(PetscMin(x[i],l[i]),PetscMin(x[i],u[i])),PetscMin(l[i],u[i])); 432 } 433 434 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr); 435 ierr = VecRestoreArray(lb,&l);CHKERRQ(ierr); 436 ierr = VecRestoreArray(ub,&u);CHKERRQ(ierr); 437 PetscFunctionReturn(0); 438 } 439 440 /* -------------------------------------------------------------------- 441 442 This file implements a semismooth truncated Newton method with a line search, 443 for solving a system of nonlinear equations in complementarity form, using the KSP, Vec, 444 and Mat interfaces for linear solvers, vectors, and matrices, 445 respectively. 446 447 The following basic routines are required for each nonlinear solver: 448 SNESCreate_XXX() - Creates a nonlinear solver context 449 SNESSetFromOptions_XXX() - Sets runtime options 450 SNESSolve_XXX() - Solves the nonlinear system 451 SNESDestroy_XXX() - Destroys the nonlinear solver context 452 The suffix "_XXX" denotes a particular implementation, in this case 453 we use _VI (e.g., SNESCreate_VI, SNESSolve_VI) for solving 454 systems of nonlinear equations with a line search (LS) method. 455 These routines are actually called via the common user interface 456 routines SNESCreate(), SNESSetFromOptions(), SNESSolve(), and 457 SNESDestroy(), so the application code interface remains identical 458 for all nonlinear solvers. 459 460 Another key routine is: 461 SNESSetUp_XXX() - Prepares for the use of a nonlinear solver 462 by setting data structures and options. The interface routine SNESSetUp() 463 is not usually called directly by the user, but instead is called by 464 SNESSolve() if necessary. 465 466 Additional basic routines are: 467 SNESView_XXX() - Prints details of runtime options that 468 have actually been used. 469 These are called by application codes via the interface routines 470 SNESView(). 471 472 The various types of solvers (preconditioners, Krylov subspace methods, 473 nonlinear solvers, timesteppers) are all organized similarly, so the 474 above description applies to these categories also. 475 476 -------------------------------------------------------------------- */ 477 /* 478 SNESSolveVI_SS - Solves the complementarity problem with a semismooth Newton 479 method using a line search. 480 481 Input Parameters: 482 . snes - the SNES context 483 484 Output Parameter: 485 . outits - number of iterations until termination 486 487 Application Interface Routine: SNESSolve() 488 489 Notes: 490 This implements essentially a semismooth Newton method with a 491 line search. By default a cubic backtracking line search 492 is employed, as described in the text "Numerical Methods for 493 Unconstrained Optimization and Nonlinear Equations" by Dennis 494 and Schnabel. 495 */ 496 #undef __FUNCT__ 497 #define __FUNCT__ "SNESSolveVI_SS" 498 PetscErrorCode SNESSolveVI_SS(SNES snes) 499 { 500 SNES_VI *vi = (SNES_VI*)snes->data; 501 PetscErrorCode ierr; 502 PetscInt maxits,i,lits; 503 PetscBool lssucceed; 504 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 505 PetscReal gnorm,xnorm=0,ynorm; 506 Vec Y,X,F,G,W; 507 KSPConvergedReason kspreason; 508 509 PetscFunctionBegin; 510 vi->computeuserfunction = snes->ops->computefunction; 511 snes->ops->computefunction = SNESVIComputeFunction; 512 513 snes->numFailures = 0; 514 snes->numLinearSolveFailures = 0; 515 snes->reason = SNES_CONVERGED_ITERATING; 516 517 maxits = snes->max_its; /* maximum number of iterations */ 518 X = snes->vec_sol; /* solution vector */ 519 F = snes->vec_func; /* residual vector */ 520 Y = snes->work[0]; /* work vectors */ 521 G = snes->work[1]; 522 W = snes->work[2]; 523 524 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 525 snes->iter = 0; 526 snes->norm = 0.0; 527 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 528 529 ierr = SNESVIAdjustInitialGuess(X,vi->xl,vi->xu);CHKERRQ(ierr); 530 ierr = SNESComputeFunction(snes,X,vi->phi);CHKERRQ(ierr); 531 if (snes->domainerror) { 532 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 533 snes->ops->computefunction = vi->computeuserfunction; 534 PetscFunctionReturn(0); 535 } 536 /* Compute Merit function */ 537 ierr = SNESVIComputeMeritFunction(vi->phi,&vi->merit,&vi->phinorm);CHKERRQ(ierr); 538 539 ierr = VecNormBegin(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm <- ||x|| */ 540 ierr = VecNormEnd(X,NORM_2,&xnorm);CHKERRQ(ierr); 541 if PetscIsInfOrNanReal(vi->merit) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 542 543 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 544 snes->norm = vi->phinorm; 545 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 546 SNESLogConvHistory(snes,vi->phinorm,0); 547 SNESMonitor(snes,0,vi->phinorm); 548 549 /* set parameter for default relative tolerance convergence test */ 550 snes->ttol = vi->phinorm*snes->rtol; 551 /* test convergence */ 552 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,vi->phinorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 553 if (snes->reason) { 554 snes->ops->computefunction = vi->computeuserfunction; 555 PetscFunctionReturn(0); 556 } 557 558 for (i=0; i<maxits; i++) { 559 560 /* Call general purpose update function */ 561 if (snes->ops->update) { 562 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 563 } 564 565 /* Solve J Y = Phi, where J is the semismooth jacobian */ 566 /* Get the nonlinear function jacobian */ 567 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 568 /* Get the diagonal shift and row scaling vectors */ 569 ierr = SNESVIComputeBsubdifferentialVectors(snes,X,F,snes->jacobian,vi->Da,vi->Db);CHKERRQ(ierr); 570 /* Compute the semismooth jacobian */ 571 ierr = SNESVIComputeJacobian(snes->jacobian,snes->jacobian_pre,vi->Da,vi->Db);CHKERRQ(ierr); 572 /* Compute the merit function gradient */ 573 ierr = SNESVIComputeMeritFunctionGradient(snes->jacobian,vi->phi,vi->dpsi);CHKERRQ(ierr); 574 ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,flg);CHKERRQ(ierr); 575 ierr = SNES_KSPSolve(snes,snes->ksp,vi->phi,Y);CHKERRQ(ierr); 576 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 577 578 if (kspreason < 0) { 579 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 580 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 581 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 582 break; 583 } 584 } 585 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 586 snes->linear_its += lits; 587 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 588 /* 589 if (vi->precheckstep) { 590 PetscBool changed_y = PETSC_FALSE; 591 ierr = (*vi->precheckstep)(snes,X,Y,vi->precheck,&changed_y);CHKERRQ(ierr); 592 } 593 594 if (PetscLogPrintInfo){ 595 ierr = SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr); 596 } 597 */ 598 /* Compute a (scaled) negative update in the line search routine: 599 Y <- X - lambda*Y 600 and evaluate G = function(Y) (depends on the line search). 601 */ 602 ierr = VecCopy(Y,snes->vec_sol_update);CHKERRQ(ierr); 603 ynorm = 1; gnorm = vi->phinorm; 604 ierr = (*vi->LineSearch)(snes,vi->lsP,X,vi->phi,G,Y,W,vi->phinorm,xnorm,&ynorm,&gnorm,&lssucceed);CHKERRQ(ierr); 605 ierr = PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",vi->phinorm,gnorm,ynorm,(int)lssucceed);CHKERRQ(ierr); 606 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 607 if (snes->domainerror) { 608 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 609 snes->ops->computefunction = vi->computeuserfunction; 610 PetscFunctionReturn(0); 611 } 612 if (!lssucceed) { 613 if (++snes->numFailures >= snes->maxFailures) { 614 PetscBool ismin; 615 snes->reason = SNES_DIVERGED_LINE_SEARCH; 616 ierr = SNESVICheckLocalMin_Private(snes,snes->jacobian,G,W,gnorm,&ismin);CHKERRQ(ierr); 617 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 618 break; 619 } 620 } 621 /* Update function and solution vectors */ 622 vi->phinorm = gnorm; 623 vi->merit = 0.5*vi->phinorm*vi->phinorm; 624 ierr = VecCopy(G,vi->phi);CHKERRQ(ierr); 625 ierr = VecCopy(W,X);CHKERRQ(ierr); 626 /* Monitor convergence */ 627 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 628 snes->iter = i+1; 629 snes->norm = vi->phinorm; 630 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 631 SNESLogConvHistory(snes,snes->norm,lits); 632 SNESMonitor(snes,snes->iter,snes->norm); 633 /* Test for convergence, xnorm = || X || */ 634 if (snes->ops->converged != SNESSkipConverged) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 635 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,vi->phinorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 636 if (snes->reason) break; 637 } 638 if (i == maxits) { 639 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 640 if(!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 641 } 642 snes->ops->computefunction = vi->computeuserfunction; 643 PetscFunctionReturn(0); 644 } 645 646 #undef __FUNCT__ 647 #define __FUNCT__ "SNESVICreateIndexSets_AS" 648 PetscErrorCode SNESVICreateIndexSets_AS(SNES snes,Vec Db,PetscReal thresh,IS* ISact,IS* ISinact) 649 { 650 PetscErrorCode ierr; 651 PetscInt i,nlocal,ilow,ihigh,nloc_isact=0,nloc_isinact=0; 652 PetscInt *idx_act,*idx_inact,i1=0,i2=0; 653 PetscScalar *db; 654 655 PetscFunctionBegin; 656 657 ierr = VecGetLocalSize(Db,&nlocal);CHKERRQ(ierr); 658 ierr = VecGetOwnershipRange(Db,&ilow,&ihigh);CHKERRQ(ierr); 659 ierr = VecGetArray(Db,&db);CHKERRQ(ierr); 660 /* Compute the sizes of the active and inactive sets */ 661 for (i=0; i < nlocal;i++) { 662 if (PetscAbsScalar(db[i]) <= thresh) nloc_isact++; 663 else nloc_isinact++; 664 } 665 ierr = PetscMalloc(nloc_isact*sizeof(PetscInt),&idx_act);CHKERRQ(ierr); 666 ierr = PetscMalloc(nloc_isinact*sizeof(PetscInt),&idx_inact);CHKERRQ(ierr); 667 668 /* Creating the indexing arrays */ 669 for(i=0; i < nlocal; i++) { 670 if (PetscAbsScalar(db[i]) <= thresh) idx_act[i1++] = ilow+i; 671 else idx_inact[i2++] = ilow+i; 672 } 673 674 /* Create the index sets */ 675 ierr = ISCreateGeneral(((PetscObject)snes)->comm,nloc_isact,idx_act,PETSC_COPY_VALUES,ISact);CHKERRQ(ierr); 676 ierr = ISCreateGeneral(((PetscObject)snes)->comm,nloc_isinact,idx_inact,PETSC_COPY_VALUES,ISinact);CHKERRQ(ierr); 677 678 ierr = VecRestoreArray(Db,&db);CHKERRQ(ierr); 679 ierr = PetscFree(idx_act);CHKERRQ(ierr); 680 ierr = PetscFree(idx_inact);CHKERRQ(ierr); 681 PetscFunctionReturn(0); 682 } 683 684 #undef __FUNCT__ 685 #define __FUNCT__ "SNESVICreateIndexSets_RS" 686 PetscErrorCode SNESVICreateIndexSets_RS(SNES snes,Vec X,Vec Xl, Vec Xu,IS* ISact,IS* ISinact) 687 { 688 PetscErrorCode ierr; 689 PetscInt i,nlocal,ilow,ihigh,nloc_isact=0,nloc_isinact=0; 690 PetscInt *idx_act,*idx_inact,i1=0,i2=0; 691 const PetscScalar *x,*xl,*xu,*f; 692 Vec F = snes->vec_func; 693 694 PetscFunctionBegin; 695 696 ierr = VecGetLocalSize(X,&nlocal);CHKERRQ(ierr); 697 ierr = VecGetOwnershipRange(X,&ilow,&ihigh);CHKERRQ(ierr); 698 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 699 ierr = VecGetArrayRead(Xl,&xl);CHKERRQ(ierr); 700 ierr = VecGetArrayRead(Xu,&xu);CHKERRQ(ierr); 701 ierr = VecGetArrayRead(F,&f);CHKERRQ(ierr); 702 /* Compute the sizes of the active and inactive sets */ 703 for (i=0; i < nlocal;i++) { 704 if (((x[i] > xl[i] + 1.e-8 || (f[i] < 0.0)) && ((x[i] < xu[i] - 1.e-8) || f[i] > 0.0))) nloc_isinact++; 705 else nloc_isact++; 706 } 707 ierr = PetscMalloc(nloc_isact*sizeof(PetscInt),&idx_act);CHKERRQ(ierr); 708 ierr = PetscMalloc(nloc_isinact*sizeof(PetscInt),&idx_inact);CHKERRQ(ierr); 709 710 /* Creating the indexing arrays */ 711 for(i=0; i < nlocal; i++) { 712 if (((x[i] > xl[i] + 1.e-8 || (f[i] < 0.0)) && ((x[i] < xu[i] - 1.e-8) || f[i] > 0.0))) idx_inact[i2++] = ilow+i; 713 else idx_act[i1++] = ilow+i; 714 } 715 716 /* Create the index sets */ 717 ierr = ISCreateGeneral(((PetscObject)snes)->comm,nloc_isact,idx_act,PETSC_COPY_VALUES,ISact);CHKERRQ(ierr); 718 ierr = ISCreateGeneral(((PetscObject)snes)->comm,nloc_isinact,idx_inact,PETSC_COPY_VALUES,ISinact);CHKERRQ(ierr); 719 720 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 721 ierr = VecRestoreArrayRead(Xl,&xl);CHKERRQ(ierr); 722 ierr = VecRestoreArrayRead(Xu,&xu);CHKERRQ(ierr); 723 ierr = VecRestoreArrayRead(F,&f);CHKERRQ(ierr); 724 ierr = PetscFree(idx_act);CHKERRQ(ierr); 725 ierr = PetscFree(idx_inact);CHKERRQ(ierr); 726 PetscFunctionReturn(0); 727 } 728 729 /* Create active and inactive set vectors. The local size of this vector is set and petsc computes the global size */ 730 #undef __FUNCT__ 731 #define __FUNCT__ "SNESVICreateSubVectors" 732 PetscErrorCode SNESVICreateSubVectors(SNES snes,PetscInt n,Vec* newv) 733 { 734 PetscErrorCode ierr; 735 Vec v; 736 737 PetscFunctionBegin; 738 ierr = VecCreate(((PetscObject)snes)->comm,&v);CHKERRQ(ierr); 739 ierr = VecSetSizes(v,n,PETSC_DECIDE);CHKERRQ(ierr); 740 ierr = VecSetFromOptions(v);CHKERRQ(ierr); 741 *newv = v; 742 743 PetscFunctionReturn(0); 744 } 745 746 /* Resets the snes PC and KSP when the active set sizes change */ 747 #undef __FUNCT__ 748 #define __FUNCT__ "SNESVIResetPCandKSP" 749 PetscErrorCode SNESVIResetPCandKSP(SNES snes,Mat Amat,Mat Pmat) 750 { 751 PetscErrorCode ierr; 752 KSP kspnew,snesksp; 753 PC pcnew; 754 const MatSolverPackage stype; 755 756 PetscFunctionBegin; 757 /* The active and inactive set sizes have changed so need to create a new snes->ksp object */ 758 ierr = SNESGetKSP(snes,&snesksp);CHKERRQ(ierr); 759 ierr = KSPCreate(((PetscObject)snes)->comm,&kspnew);CHKERRQ(ierr); 760 /* Copy over snes->ksp info */ 761 kspnew->pc_side = snesksp->pc_side; 762 kspnew->rtol = snesksp->rtol; 763 kspnew->abstol = snesksp->abstol; 764 kspnew->max_it = snesksp->max_it; 765 ierr = KSPSetType(kspnew,((PetscObject)snesksp)->type_name);CHKERRQ(ierr); 766 ierr = KSPGetPC(kspnew,&pcnew);CHKERRQ(ierr); 767 ierr = PCSetType(kspnew->pc,((PetscObject)snesksp->pc)->type_name);CHKERRQ(ierr); 768 ierr = PCSetOperators(kspnew->pc,Amat,Pmat,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 769 ierr = PCFactorGetMatSolverPackage(snesksp->pc,&stype);CHKERRQ(ierr); 770 ierr = PCFactorSetMatSolverPackage(kspnew->pc,stype);CHKERRQ(ierr); 771 ierr = KSPDestroy(snesksp);CHKERRQ(ierr); 772 snes->ksp = kspnew; 773 ierr = PetscLogObjectParent(snes,kspnew);CHKERRQ(ierr); 774 ierr = KSPSetFromOptions(kspnew);CHKERRQ(ierr); 775 PetscFunctionReturn(0); 776 } 777 778 779 #undef __FUNCT__ 780 #define __FUNCT__ "SNESVIComputeInactiveSetFnorm" 781 PetscErrorCode SNESVIComputeInactiveSetFnorm(SNES snes,Vec F,Vec X,PetscScalar *fnorm) 782 { 783 PetscErrorCode ierr; 784 SNES_VI *vi = (SNES_VI*)snes->data; 785 const PetscScalar *x,*xl,*xu,*f; 786 PetscInt i,n; 787 PetscReal rnorm; 788 789 PetscFunctionBegin; 790 ierr = VecGetLocalSize(X,&n);CHKERRQ(ierr); 791 ierr = VecGetArrayRead(vi->xl,&xl);CHKERRQ(ierr); 792 ierr = VecGetArrayRead(vi->xu,&xu);CHKERRQ(ierr); 793 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 794 ierr = VecGetArrayRead(F,&f);CHKERRQ(ierr); 795 rnorm = 0.0; 796 for (i=0; i<n; i++) { 797 if (((x[i] > xl[i] + 1.e-8 || (f[i] < 0.0)) && ((x[i] < xu[i] - 1.e-8) || f[i] > 0.0))) rnorm += f[i]*f[i]; 798 } 799 ierr = VecRestoreArrayRead(F,&f);CHKERRQ(ierr); 800 ierr = VecRestoreArrayRead(vi->xl,&xl);CHKERRQ(ierr); 801 ierr = VecRestoreArrayRead(vi->xu,&xu);CHKERRQ(ierr); 802 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 803 ierr = MPI_Allreduce(&rnorm,fnorm,1,MPIU_REAL,MPI_SUM,((PetscObject)snes)->comm);CHKERRQ(ierr); 804 *fnorm = sqrt(*fnorm); 805 PetscFunctionReturn(0); 806 } 807 808 /* Variational Inequality solver using reduce space method. No semismooth algorithm is 809 implemented in this algorithm. It basically identifies the active variables and does 810 a linear solve on the inactive variables. */ 811 #undef __FUNCT__ 812 #define __FUNCT__ "SNESSolveVI_RS" 813 PetscErrorCode SNESSolveVI_RS(SNES snes) 814 { 815 SNES_VI *vi = (SNES_VI*)snes->data; 816 PetscErrorCode ierr; 817 PetscInt maxits,i,lits,Nis_inact; 818 PetscBool lssucceed; 819 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 820 PetscReal fnorm,gnorm,xnorm=0,ynorm; 821 Vec Y,X,F,G,W; 822 KSPConvergedReason kspreason; 823 824 PetscFunctionBegin; 825 snes->numFailures = 0; 826 snes->numLinearSolveFailures = 0; 827 snes->reason = SNES_CONVERGED_ITERATING; 828 829 maxits = snes->max_its; /* maximum number of iterations */ 830 X = snes->vec_sol; /* solution vector */ 831 F = snes->vec_func; /* residual vector */ 832 Y = snes->work[0]; /* work vectors */ 833 G = snes->work[1]; 834 W = snes->work[2]; 835 836 Nis_inact = F->map->N; 837 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 838 snes->iter = 0; 839 snes->norm = 0.0; 840 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 841 842 ierr = SNESVIAdjustInitialGuess(X,vi->xl,vi->xu);CHKERRQ(ierr); 843 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 844 if (snes->domainerror) { 845 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 846 PetscFunctionReturn(0); 847 } 848 ierr = SNESVIComputeInactiveSetFnorm(snes,F,X,&fnorm);CHKERRQ(ierr); 849 ierr = VecNormBegin(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm <- ||x|| */ 850 ierr = VecNormEnd(X,NORM_2,&xnorm);CHKERRQ(ierr); 851 if PetscIsInfOrNanReal(fnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 852 853 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 854 snes->norm = fnorm; 855 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 856 SNESLogConvHistory(snes,fnorm,0); 857 SNESMonitor(snes,0,fnorm); 858 859 /* set parameter for default relative tolerance convergence test */ 860 snes->ttol = fnorm*snes->rtol; 861 /* test convergence */ 862 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 863 if (snes->reason) PetscFunctionReturn(0); 864 865 for (i=0; i<maxits; i++) { 866 867 IS IS_act,IS_inact; /* _act -> active set _inact -> inactive set */ 868 VecScatter scat_act,scat_inact; 869 PetscInt nis_act,nis_inact,Nis_inact_prev; 870 Vec Y_act,Y_inact,F_inact; 871 Mat jac_inact_inact,prejac_inact_inact; 872 873 /* Call general purpose update function */ 874 if (snes->ops->update) { 875 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 876 } 877 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 878 /* Create active and inactive index sets */ 879 ierr = SNESVICreateIndexSets_RS(snes,X,vi->xl,vi->xu,&IS_act,&IS_inact);CHKERRQ(ierr); 880 881 Nis_inact_prev = Nis_inact; 882 /* Get sizes of active and inactive sets */ 883 ierr = ISGetLocalSize(IS_act,&nis_act);CHKERRQ(ierr); 884 ierr = ISGetLocalSize(IS_inact,&nis_inact);CHKERRQ(ierr); 885 ierr = ISGetSize(IS_inact,&Nis_inact);CHKERRQ(ierr); 886 887 ierr = PetscPrintf(PETSC_COMM_WORLD,"Size of active set = %d, size of inactive set = %d\n",nis_act,nis_inact);CHKERRQ(ierr); 888 889 /* Create active and inactive set vectors */ 890 ierr = SNESVICreateSubVectors(snes,nis_inact,&F_inact);CHKERRQ(ierr); 891 ierr = SNESVICreateSubVectors(snes,nis_act,&Y_act);CHKERRQ(ierr); 892 ierr = SNESVICreateSubVectors(snes,nis_inact,&Y_inact);CHKERRQ(ierr); 893 894 /* Create inactive set submatrices */ 895 ierr = MatGetSubMatrix(snes->jacobian,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);CHKERRQ(ierr); 896 897 /* Create scatter contexts */ 898 ierr = VecScatterCreate(Y,IS_act,Y_act,PETSC_NULL,&scat_act);CHKERRQ(ierr); 899 ierr = VecScatterCreate(Y,IS_inact,Y_inact,PETSC_NULL,&scat_inact);CHKERRQ(ierr); 900 901 /* Do a vec scatter to active and inactive set vectors */ 902 ierr = VecScatterBegin(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 903 ierr = VecScatterEnd(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 904 905 ierr = VecScatterBegin(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 906 ierr = VecScatterEnd(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 907 908 ierr = VecScatterBegin(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 909 ierr = VecScatterEnd(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 910 911 /* Active set direction = 0 */ 912 ierr = VecSet(Y_act,0);CHKERRQ(ierr); 913 if (snes->jacobian != snes->jacobian_pre) { 914 ierr = MatGetSubMatrix(snes->jacobian_pre,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&prejac_inact_inact);CHKERRQ(ierr); 915 } else prejac_inact_inact = jac_inact_inact; 916 917 if (Nis_inact != Nis_inact_prev) { 918 ierr = SNESVIResetPCandKSP(snes,jac_inact_inact,prejac_inact_inact);CHKERRQ(ierr); 919 } 920 921 ierr = KSPSetOperators(snes->ksp,jac_inact_inact,prejac_inact_inact,flg);CHKERRQ(ierr); 922 ierr = SNES_KSPSolve(snes,snes->ksp,F_inact,Y_inact);CHKERRQ(ierr); 923 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 924 if (kspreason < 0) { 925 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 926 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 927 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 928 break; 929 } 930 } 931 932 ierr = VecScatterBegin(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 933 ierr = VecScatterEnd(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 934 ierr = VecScatterBegin(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 935 ierr = VecScatterEnd(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 936 937 ierr = VecDestroy(F_inact);CHKERRQ(ierr); 938 ierr = VecDestroy(Y_act);CHKERRQ(ierr); 939 ierr = VecDestroy(Y_inact);CHKERRQ(ierr); 940 ierr = VecScatterDestroy(scat_act);CHKERRQ(ierr); 941 ierr = VecScatterDestroy(scat_inact);CHKERRQ(ierr); 942 ierr = ISDestroy(IS_act);CHKERRQ(ierr); 943 ierr = ISDestroy(IS_inact);CHKERRQ(ierr); 944 ierr = MatDestroy(jac_inact_inact);CHKERRQ(ierr); 945 if (snes->jacobian != snes->jacobian_pre) { 946 ierr = MatDestroy(prejac_inact_inact);CHKERRQ(ierr); 947 } 948 949 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 950 snes->linear_its += lits; 951 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 952 /* 953 if (vi->precheckstep) { 954 PetscBool changed_y = PETSC_FALSE; 955 ierr = (*vi->precheckstep)(snes,X,Y,vi->precheck,&changed_y);CHKERRQ(ierr); 956 } 957 958 if (PetscLogPrintInfo){ 959 ierr = SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr); 960 } 961 */ 962 /* Compute a (scaled) negative update in the line search routine: 963 Y <- X - lambda*Y 964 and evaluate G = function(Y) (depends on the line search). 965 */ 966 ierr = VecCopy(Y,snes->vec_sol_update);CHKERRQ(ierr); 967 ynorm = 1; gnorm = fnorm; 968 ierr = (*vi->LineSearch)(snes,vi->lsP,X,F,G,Y,W,fnorm,xnorm,&ynorm,&gnorm,&lssucceed);CHKERRQ(ierr); 969 ierr = PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",fnorm,gnorm,ynorm,(int)lssucceed);CHKERRQ(ierr); 970 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 971 if (snes->domainerror) { 972 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 973 PetscFunctionReturn(0); 974 } 975 if (!lssucceed) { 976 if (++snes->numFailures >= snes->maxFailures) { 977 PetscBool ismin; 978 snes->reason = SNES_DIVERGED_LINE_SEARCH; 979 ierr = SNESVICheckLocalMin_Private(snes,snes->jacobian,G,W,gnorm,&ismin);CHKERRQ(ierr); 980 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 981 break; 982 } 983 } 984 /* Update function and solution vectors */ 985 fnorm = gnorm; 986 ierr = VecCopy(G,F);CHKERRQ(ierr); 987 ierr = VecCopy(W,X);CHKERRQ(ierr); 988 /* Monitor convergence */ 989 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 990 snes->iter = i+1; 991 snes->norm = fnorm; 992 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 993 SNESLogConvHistory(snes,snes->norm,lits); 994 SNESMonitor(snes,snes->iter,snes->norm); 995 /* Test for convergence, xnorm = || X || */ 996 if (snes->ops->converged != SNESSkipConverged) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 997 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 998 if (snes->reason) break; 999 } 1000 if (i == maxits) { 1001 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 1002 if(!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 1003 } 1004 PetscFunctionReturn(0); 1005 } 1006 1007 /* -------------------------------------------------------------------------- */ 1008 /* 1009 SNESSetUp_VI - Sets up the internal data structures for the later use 1010 of the SNESVI nonlinear solver. 1011 1012 Input Parameter: 1013 . snes - the SNES context 1014 . x - the solution vector 1015 1016 Application Interface Routine: SNESSetUp() 1017 1018 Notes: 1019 For basic use of the SNES solvers, the user need not explicitly call 1020 SNESSetUp(), since these actions will automatically occur during 1021 the call to SNESSolve(). 1022 */ 1023 #undef __FUNCT__ 1024 #define __FUNCT__ "SNESSetUp_VI" 1025 PetscErrorCode SNESSetUp_VI(SNES snes) 1026 { 1027 PetscErrorCode ierr; 1028 SNES_VI *vi = (SNES_VI*) snes->data; 1029 PetscInt i_start[3],i_end[3]; 1030 1031 PetscFunctionBegin; 1032 if (!snes->vec_sol_update) { 1033 ierr = VecDuplicate(snes->vec_sol,&snes->vec_sol_update);CHKERRQ(ierr); 1034 ierr = PetscLogObjectParent(snes,snes->vec_sol_update);CHKERRQ(ierr); 1035 } 1036 if (!snes->work) { 1037 snes->nwork = 3; 1038 ierr = VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work);CHKERRQ(ierr); 1039 ierr = PetscLogObjectParents(snes,snes->nwork,snes->work);CHKERRQ(ierr); 1040 } 1041 1042 /* If the lower and upper bound on variables are not set, set it to 1043 -Inf and Inf */ 1044 if (!vi->xl && !vi->xu) { 1045 vi->usersetxbounds = PETSC_FALSE; 1046 ierr = VecDuplicate(snes->vec_sol, &vi->xl); CHKERRQ(ierr); 1047 ierr = VecSet(vi->xl,PETSC_VI_NINF);CHKERRQ(ierr); 1048 ierr = VecDuplicate(snes->vec_sol, &vi->xu); CHKERRQ(ierr); 1049 ierr = VecSet(vi->xu,PETSC_VI_INF);CHKERRQ(ierr); 1050 } else { 1051 /* Check if lower bound, upper bound and solution vector distribution across the processors is identical */ 1052 ierr = VecGetOwnershipRange(snes->vec_sol,i_start,i_end);CHKERRQ(ierr); 1053 ierr = VecGetOwnershipRange(vi->xl,i_start+1,i_end+1);CHKERRQ(ierr); 1054 ierr = VecGetOwnershipRange(vi->xu,i_start+2,i_end+2);CHKERRQ(ierr); 1055 if ((i_start[0] != i_start[1]) || (i_start[0] != i_start[2]) || (i_end[0] != i_end[1]) || (i_end[0] != i_end[2])) 1056 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Distribution of lower bound, upper bound and the solution vector should be identical across all the processors."); 1057 } 1058 1059 if (snes->ops->solve != SNESSolveVI_RS) { 1060 ierr = VecDuplicate(snes->vec_sol, &vi->dpsi);CHKERRQ(ierr); 1061 ierr = VecDuplicate(snes->vec_sol, &vi->phi); CHKERRQ(ierr); 1062 ierr = VecDuplicate(snes->vec_sol, &vi->Da); CHKERRQ(ierr); 1063 ierr = VecDuplicate(snes->vec_sol, &vi->Db); CHKERRQ(ierr); 1064 ierr = VecDuplicate(snes->vec_sol, &vi->z);CHKERRQ(ierr); 1065 ierr = VecDuplicate(snes->vec_sol, &vi->t); CHKERRQ(ierr); 1066 1067 } 1068 PetscFunctionReturn(0); 1069 } 1070 /* -------------------------------------------------------------------------- */ 1071 /* 1072 SNESDestroy_VI - Destroys the private SNES_VI context that was created 1073 with SNESCreate_VI(). 1074 1075 Input Parameter: 1076 . snes - the SNES context 1077 1078 Application Interface Routine: SNESDestroy() 1079 */ 1080 #undef __FUNCT__ 1081 #define __FUNCT__ "SNESDestroy_VI" 1082 PetscErrorCode SNESDestroy_VI(SNES snes) 1083 { 1084 SNES_VI *vi = (SNES_VI*) snes->data; 1085 PetscErrorCode ierr; 1086 1087 PetscFunctionBegin; 1088 if (snes->vec_sol_update) { 1089 ierr = VecDestroy(snes->vec_sol_update);CHKERRQ(ierr); 1090 snes->vec_sol_update = PETSC_NULL; 1091 } 1092 if (snes->nwork) { 1093 ierr = VecDestroyVecs(snes->work,snes->nwork);CHKERRQ(ierr); 1094 snes->nwork = 0; 1095 snes->work = PETSC_NULL; 1096 } 1097 1098 if (snes->ops->solve != SNESSolveVI_RS) { 1099 /* clear vectors */ 1100 ierr = VecDestroy(vi->dpsi);CHKERRQ(ierr); 1101 ierr = VecDestroy(vi->phi); CHKERRQ(ierr); 1102 ierr = VecDestroy(vi->Da); CHKERRQ(ierr); 1103 ierr = VecDestroy(vi->Db); CHKERRQ(ierr); 1104 ierr = VecDestroy(vi->z); CHKERRQ(ierr); 1105 ierr = VecDestroy(vi->t); CHKERRQ(ierr); 1106 if (!vi->usersetxbounds) { 1107 ierr = VecDestroy(vi->xl); CHKERRQ(ierr); 1108 ierr = VecDestroy(vi->xu); CHKERRQ(ierr); 1109 } 1110 } 1111 if (vi->lsmonitor) { 1112 ierr = PetscViewerASCIIMonitorDestroy(vi->lsmonitor);CHKERRQ(ierr); 1113 } 1114 ierr = PetscFree(snes->data);CHKERRQ(ierr); 1115 1116 /* clear composed functions */ 1117 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSet_C","",PETSC_NULL);CHKERRQ(ierr); 1118 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetMonitor_C","",PETSC_NULL);CHKERRQ(ierr); 1119 PetscFunctionReturn(0); 1120 } 1121 /* -------------------------------------------------------------------------- */ 1122 #undef __FUNCT__ 1123 #define __FUNCT__ "SNESLineSearchNo_VI" 1124 1125 /* 1126 This routine is a copy of SNESLineSearchNo routine in snes/impls/ls/ls.c 1127 1128 */ 1129 PetscErrorCode SNESLineSearchNo_VI(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal xnorm,PetscReal *ynorm,PetscReal *gnorm,PetscBool *flag) 1130 { 1131 PetscErrorCode ierr; 1132 SNES_VI *vi = (SNES_VI*)snes->data; 1133 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1134 1135 PetscFunctionBegin; 1136 *flag = PETSC_TRUE; 1137 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1138 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); /* ynorm = || y || */ 1139 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1140 if (vi->postcheckstep) { 1141 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1142 } 1143 if (changed_y) { 1144 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1145 } 1146 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1147 if (!snes->domainerror) { 1148 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); /* gnorm = || g || */ 1149 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1150 } 1151 if (vi->lsmonitor) { 1152 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1153 } 1154 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1155 PetscFunctionReturn(0); 1156 } 1157 1158 /* -------------------------------------------------------------------------- */ 1159 #undef __FUNCT__ 1160 #define __FUNCT__ "SNESLineSearchBoundProjectedNo_VI" 1161 1162 /* 1163 This routine does not actually do a line search but it takes a full newton 1164 step while ensuring that the new iterates remain within the constraints. 1165 1166 */ 1167 PetscErrorCode SNESLineSearchBoundProjectedNo_VI(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal xnorm,PetscReal *ynorm,PetscReal *gnorm,PetscBool *flag) 1168 { 1169 PetscErrorCode ierr; 1170 SNES_VI *vi = (SNES_VI*)snes->data; 1171 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1172 1173 PetscFunctionBegin; 1174 *flag = PETSC_TRUE; 1175 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1176 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); /* ynorm = || y || */ 1177 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1178 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1179 if (vi->postcheckstep) { 1180 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1181 } 1182 if (changed_y) { 1183 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1184 } 1185 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1186 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1187 if (!snes->domainerror) { 1188 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); /* gnorm = || g || */ 1189 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1190 } 1191 if (vi->lsmonitor) { 1192 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1193 } 1194 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1195 PetscFunctionReturn(0); 1196 } 1197 1198 /* -------------------------------------------------------------------------- */ 1199 #undef __FUNCT__ 1200 #define __FUNCT__ "SNESLineSearchNoNorms_VI" 1201 1202 /* 1203 This routine is a copy of SNESLineSearchNoNorms in snes/impls/ls/ls.c 1204 */ 1205 PetscErrorCode SNESLineSearchNoNorms_VI(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal xnorm,PetscReal *ynorm,PetscReal *gnorm,PetscBool *flag) 1206 { 1207 PetscErrorCode ierr; 1208 SNES_VI *vi = (SNES_VI*)snes->data; 1209 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1210 1211 PetscFunctionBegin; 1212 *flag = PETSC_TRUE; 1213 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1214 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1215 if (vi->postcheckstep) { 1216 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1217 } 1218 if (changed_y) { 1219 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1220 } 1221 1222 /* don't evaluate function the last time through */ 1223 if (snes->iter < snes->max_its-1) { 1224 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1225 } 1226 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1227 PetscFunctionReturn(0); 1228 } 1229 /* -------------------------------------------------------------------------- */ 1230 #undef __FUNCT__ 1231 #define __FUNCT__ "SNESLineSearchCubic_VI" 1232 /* 1233 This routine is a copy of SNESLineSearchCubic in snes/impls/ls/ls.c 1234 */ 1235 PetscErrorCode SNESLineSearchCubic_VI(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal xnorm,PetscReal *ynorm,PetscReal *gnorm,PetscBool *flag) 1236 { 1237 /* 1238 Note that for line search purposes we work with with the related 1239 minimization problem: 1240 min z(x): R^n -> R, 1241 where z(x) = .5 * fnorm*fnorm, and fnorm = || f ||_2. 1242 This function z(x) is same as the merit function for the complementarity problem. 1243 */ 1244 1245 PetscReal initslope,lambdaprev,gnormprev,a,b,d,t1,t2,rellength; 1246 PetscReal minlambda,lambda,lambdatemp; 1247 #if defined(PETSC_USE_COMPLEX) 1248 PetscScalar cinitslope; 1249 #endif 1250 PetscErrorCode ierr; 1251 PetscInt count; 1252 SNES_VI *vi = (SNES_VI*)snes->data; 1253 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1254 MPI_Comm comm; 1255 1256 PetscFunctionBegin; 1257 ierr = PetscObjectGetComm((PetscObject)snes,&comm);CHKERRQ(ierr); 1258 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1259 *flag = PETSC_TRUE; 1260 1261 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); 1262 if (*ynorm == 0.0) { 1263 if (vi->lsmonitor) { 1264 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Initial direction and size is 0\n");CHKERRQ(ierr); 1265 } 1266 *gnorm = fnorm; 1267 ierr = VecCopy(x,w);CHKERRQ(ierr); 1268 ierr = VecCopy(f,g);CHKERRQ(ierr); 1269 *flag = PETSC_FALSE; 1270 goto theend1; 1271 } 1272 if (*ynorm > vi->maxstep) { /* Step too big, so scale back */ 1273 if (vi->lsmonitor) { 1274 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Scaling step by %G old ynorm %G\n",vi->maxstep/(*ynorm),*ynorm);CHKERRQ(ierr); 1275 } 1276 ierr = VecScale(y,vi->maxstep/(*ynorm));CHKERRQ(ierr); 1277 *ynorm = vi->maxstep; 1278 } 1279 ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr); 1280 minlambda = vi->minlambda/rellength; 1281 #if defined(PETSC_USE_COMPLEX) 1282 ierr = VecDot(vi->dpsi,y,&cinitslope);CHKERRQ(ierr); 1283 initslope = PetscRealPart(cinitslope); 1284 #else 1285 ierr = VecDot(vi->dpsi,y,&initslope);CHKERRQ(ierr); 1286 #endif 1287 if (initslope > 0.0) initslope = -initslope; 1288 if (initslope == 0.0) initslope = -1.0; 1289 1290 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1291 if (snes->nfuncs >= snes->max_funcs) { 1292 ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); 1293 *flag = PETSC_FALSE; 1294 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1295 goto theend1; 1296 } 1297 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1298 if (snes->domainerror) { 1299 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1300 PetscFunctionReturn(0); 1301 } 1302 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1303 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1304 ierr = PetscInfo2(snes,"Initial fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1305 if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + vi->alpha*initslope) { /* Sufficient reduction */ 1306 if (vi->lsmonitor) { 1307 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1308 } 1309 goto theend1; 1310 } 1311 1312 /* Fit points with quadratic */ 1313 lambda = 1.0; 1314 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 1315 lambdaprev = lambda; 1316 gnormprev = *gnorm; 1317 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1318 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1319 else lambda = lambdatemp; 1320 1321 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1322 if (snes->nfuncs >= snes->max_funcs) { 1323 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);CHKERRQ(ierr); 1324 *flag = PETSC_FALSE; 1325 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1326 goto theend1; 1327 } 1328 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1329 if (snes->domainerror) { 1330 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1331 PetscFunctionReturn(0); 1332 } 1333 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1334 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1335 if (vi->lsmonitor) { 1336 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: gnorm after quadratic fit %G\n",*gnorm);CHKERRQ(ierr); 1337 } 1338 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* sufficient reduction */ 1339 if (vi->lsmonitor) { 1340 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Quadratically determined step, lambda=%18.16e\n",lambda);CHKERRQ(ierr); 1341 } 1342 goto theend1; 1343 } 1344 1345 /* Fit points with cubic */ 1346 count = 1; 1347 while (PETSC_TRUE) { 1348 if (lambda <= minlambda) { 1349 if (vi->lsmonitor) { 1350 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: unable to find good step length! After %D tries \n",count);CHKERRQ(ierr); 1351 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",fnorm,*gnorm,*ynorm,minlambda,lambda,initslope);CHKERRQ(ierr); 1352 } 1353 *flag = PETSC_FALSE; 1354 break; 1355 } 1356 t1 = .5*((*gnorm)*(*gnorm) - fnorm*fnorm) - lambda*initslope; 1357 t2 = .5*(gnormprev*gnormprev - fnorm*fnorm) - lambdaprev*initslope; 1358 a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 1359 b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 1360 d = b*b - 3*a*initslope; 1361 if (d < 0.0) d = 0.0; 1362 if (a == 0.0) { 1363 lambdatemp = -initslope/(2.0*b); 1364 } else { 1365 lambdatemp = (-b + sqrt(d))/(3.0*a); 1366 } 1367 lambdaprev = lambda; 1368 gnormprev = *gnorm; 1369 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1370 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1371 else lambda = lambdatemp; 1372 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1373 if (snes->nfuncs >= snes->max_funcs) { 1374 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 1375 ierr = PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",fnorm,*gnorm,*ynorm,lambda,initslope);CHKERRQ(ierr); 1376 *flag = PETSC_FALSE; 1377 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1378 break; 1379 } 1380 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1381 if (snes->domainerror) { 1382 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1383 PetscFunctionReturn(0); 1384 } 1385 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1386 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1387 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* is reduction enough? */ 1388 if (vi->lsmonitor) { 1389 ierr = PetscPrintf(comm," Line search: Cubically determined step, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr); 1390 } 1391 break; 1392 } else { 1393 if (vi->lsmonitor) { 1394 ierr = PetscPrintf(comm," Line search: Cubic step no good, shrinking lambda, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr); 1395 } 1396 } 1397 count++; 1398 } 1399 theend1: 1400 /* Optional user-defined check for line search step validity */ 1401 if (vi->postcheckstep && *flag) { 1402 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1403 if (changed_y) { 1404 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1405 } 1406 if (changed_y || changed_w) { /* recompute the function if the step has changed */ 1407 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1408 if (snes->domainerror) { 1409 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1410 PetscFunctionReturn(0); 1411 } 1412 ierr = VecNormBegin(g,NORM_2,gnorm);CHKERRQ(ierr); 1413 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1414 ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr); 1415 ierr = VecNormEnd(g,NORM_2,gnorm);CHKERRQ(ierr); 1416 ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr); 1417 } 1418 } 1419 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1420 PetscFunctionReturn(0); 1421 } 1422 1423 /* -------------------------------------------------------------------------- */ 1424 #undef __FUNCT__ 1425 #define __FUNCT__ "SNESLineSearchBoundProjectedCubic_VI" 1426 /* 1427 This routine implements a cubic line search while keeping the iterates within bounds 1428 */ 1429 PetscErrorCode SNESLineSearchBoundProjectedCubic_VI(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal xnorm,PetscReal *ynorm,PetscReal *gnorm,PetscBool *flag) 1430 { 1431 PetscReal initslope,lambdaprev,gnormprev,a,b,d,t1,t2,rellength; 1432 PetscReal minlambda,lambda,lambdatemp; 1433 #if defined(PETSC_USE_COMPLEX) 1434 PetscScalar cinitslope; 1435 #endif 1436 PetscErrorCode ierr; 1437 PetscInt count; 1438 SNES_VI *vi = (SNES_VI*)snes->data; 1439 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1440 MPI_Comm comm; 1441 1442 PetscFunctionBegin; 1443 ierr = PetscObjectGetComm((PetscObject)snes,&comm);CHKERRQ(ierr); 1444 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1445 *flag = PETSC_TRUE; 1446 1447 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); 1448 if (*ynorm == 0.0) { 1449 if (vi->lsmonitor) { 1450 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Initial direction and size is 0\n");CHKERRQ(ierr); 1451 } 1452 *gnorm = fnorm; 1453 ierr = VecCopy(x,w);CHKERRQ(ierr); 1454 ierr = VecCopy(f,g);CHKERRQ(ierr); 1455 *flag = PETSC_FALSE; 1456 goto theend1; 1457 } 1458 if (*ynorm > vi->maxstep) { /* Step too big, so scale back */ 1459 if (vi->lsmonitor) { 1460 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Scaling step by %G old ynorm %G\n",vi->maxstep/(*ynorm),*ynorm);CHKERRQ(ierr); 1461 } 1462 ierr = VecScale(y,vi->maxstep/(*ynorm));CHKERRQ(ierr); 1463 *ynorm = vi->maxstep; 1464 } 1465 ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr); 1466 minlambda = vi->minlambda/rellength; 1467 ierr = MatMult(snes->jacobian,y,w);CHKERRQ(ierr); 1468 #if defined(PETSC_USE_COMPLEX) 1469 ierr = VecDot(f,w,&cinitslope);CHKERRQ(ierr); 1470 initslope = PetscRealPart(cinitslope); 1471 #else 1472 ierr = VecDot(f,w,&initslope);CHKERRQ(ierr); 1473 #endif 1474 if (initslope > 0.0) initslope = -initslope; 1475 if (initslope == 0.0) initslope = -1.0; 1476 1477 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1478 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1479 if (snes->nfuncs >= snes->max_funcs) { 1480 ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); 1481 *flag = PETSC_FALSE; 1482 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1483 goto theend1; 1484 } 1485 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1486 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1487 if (snes->domainerror) { 1488 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1489 PetscFunctionReturn(0); 1490 } 1491 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1492 ierr = PetscInfo2(snes,"Initial fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1493 if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + vi->alpha*initslope) { /* Sufficient reduction */ 1494 if (vi->lsmonitor) { 1495 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1496 } 1497 goto theend1; 1498 } 1499 1500 /* Fit points with quadratic */ 1501 lambda = 1.0; 1502 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 1503 lambdaprev = lambda; 1504 gnormprev = *gnorm; 1505 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1506 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1507 else lambda = lambdatemp; 1508 1509 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1510 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1511 if (snes->nfuncs >= snes->max_funcs) { 1512 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);CHKERRQ(ierr); 1513 *flag = PETSC_FALSE; 1514 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1515 goto theend1; 1516 } 1517 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1518 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1519 if (snes->domainerror) { 1520 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1521 PetscFunctionReturn(0); 1522 } 1523 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1524 if (vi->lsmonitor) { 1525 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: gnorm after quadratic fit %G\n",*gnorm);CHKERRQ(ierr); 1526 } 1527 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* sufficient reduction */ 1528 if (vi->lsmonitor) { 1529 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Quadratically determined step, lambda=%18.16e\n",lambda);CHKERRQ(ierr); 1530 } 1531 goto theend1; 1532 } 1533 1534 /* Fit points with cubic */ 1535 count = 1; 1536 while (PETSC_TRUE) { 1537 if (lambda <= minlambda) { 1538 if (vi->lsmonitor) { 1539 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: unable to find good step length! After %D tries \n",count);CHKERRQ(ierr); 1540 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",fnorm,*gnorm,*ynorm,minlambda,lambda,initslope);CHKERRQ(ierr); 1541 } 1542 *flag = PETSC_FALSE; 1543 break; 1544 } 1545 t1 = .5*((*gnorm)*(*gnorm) - fnorm*fnorm) - lambda*initslope; 1546 t2 = .5*(gnormprev*gnormprev - fnorm*fnorm) - lambdaprev*initslope; 1547 a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 1548 b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 1549 d = b*b - 3*a*initslope; 1550 if (d < 0.0) d = 0.0; 1551 if (a == 0.0) { 1552 lambdatemp = -initslope/(2.0*b); 1553 } else { 1554 lambdatemp = (-b + sqrt(d))/(3.0*a); 1555 } 1556 lambdaprev = lambda; 1557 gnormprev = *gnorm; 1558 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1559 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1560 else lambda = lambdatemp; 1561 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1562 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1563 if (snes->nfuncs >= snes->max_funcs) { 1564 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 1565 ierr = PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",fnorm,*gnorm,*ynorm,lambda,initslope);CHKERRQ(ierr); 1566 *flag = PETSC_FALSE; 1567 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1568 break; 1569 } 1570 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1571 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1572 if (snes->domainerror) { 1573 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1574 PetscFunctionReturn(0); 1575 } 1576 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1577 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* is reduction enough? */ 1578 if (vi->lsmonitor) { 1579 ierr = PetscPrintf(comm," Line search: Cubically determined step, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr); 1580 } 1581 break; 1582 } else { 1583 if (vi->lsmonitor) { 1584 ierr = PetscPrintf(comm," Line search: Cubic step no good, shrinking lambda, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr); 1585 } 1586 } 1587 count++; 1588 } 1589 theend1: 1590 /* Optional user-defined check for line search step validity */ 1591 if (vi->postcheckstep && *flag) { 1592 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1593 if (changed_y) { 1594 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1595 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1596 } 1597 if (changed_y || changed_w) { /* recompute the function if the step has changed */ 1598 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1599 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1600 if (snes->domainerror) { 1601 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1602 PetscFunctionReturn(0); 1603 } 1604 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1605 ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr); 1606 ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr); 1607 } 1608 } 1609 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1610 PetscFunctionReturn(0); 1611 } 1612 1613 #undef __FUNCT__ 1614 #define __FUNCT__ "SNESLineSearchQuadratic_VI" 1615 /* 1616 This routine is a copy of SNESLineSearchQuadratic in snes/impls/ls/ls.c 1617 */ 1618 PetscErrorCode SNESLineSearchQuadratic_VI(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal xnorm,PetscReal *ynorm,PetscReal *gnorm,PetscBool *flag) 1619 { 1620 /* 1621 Note that for line search purposes we work with with the related 1622 minimization problem: 1623 min z(x): R^n -> R, 1624 where z(x) = .5 * fnorm*fnorm,and fnorm = || f ||_2. 1625 z(x) is the same as the merit function for the complementarity problem 1626 */ 1627 PetscReal initslope,minlambda,lambda,lambdatemp,rellength; 1628 #if defined(PETSC_USE_COMPLEX) 1629 PetscScalar cinitslope; 1630 #endif 1631 PetscErrorCode ierr; 1632 PetscInt count; 1633 SNES_VI *vi = (SNES_VI*)snes->data; 1634 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1635 1636 PetscFunctionBegin; 1637 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1638 *flag = PETSC_TRUE; 1639 1640 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); 1641 if (*ynorm == 0.0) { 1642 if (vi->lsmonitor) { 1643 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor,"Line search: Direction and size is 0\n");CHKERRQ(ierr); 1644 } 1645 *gnorm = fnorm; 1646 ierr = VecCopy(x,w);CHKERRQ(ierr); 1647 ierr = VecCopy(f,g);CHKERRQ(ierr); 1648 *flag = PETSC_FALSE; 1649 goto theend2; 1650 } 1651 if (*ynorm > vi->maxstep) { /* Step too big, so scale back */ 1652 ierr = VecScale(y,vi->maxstep/(*ynorm));CHKERRQ(ierr); 1653 *ynorm = vi->maxstep; 1654 } 1655 ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr); 1656 minlambda = vi->minlambda/rellength; 1657 #if defined(PETSC_USE_COMPLEX) 1658 ierr = VecDot(vi->dpsi,y,&cinitslope);CHKERRQ(ierr); 1659 initslope = PetscRealPart(cinitslope); 1660 #else 1661 ierr = VecDot(vi->dpsi,y,&initslope);CHKERRQ(ierr); 1662 #endif 1663 if (initslope > 0.0) initslope = -initslope; 1664 if (initslope == 0.0) initslope = -1.0; 1665 ierr = PetscInfo1(snes,"Initslope %G \n",initslope);CHKERRQ(ierr); 1666 1667 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1668 if (snes->nfuncs >= snes->max_funcs) { 1669 ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); 1670 *flag = PETSC_FALSE; 1671 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1672 goto theend2; 1673 } 1674 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1675 if (snes->domainerror) { 1676 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1677 PetscFunctionReturn(0); 1678 } 1679 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1680 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1681 if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + vi->alpha*initslope) { /* Sufficient reduction */ 1682 if (vi->lsmonitor) { 1683 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1684 } 1685 goto theend2; 1686 } 1687 1688 /* Fit points with quadratic */ 1689 lambda = 1.0; 1690 count = 1; 1691 while (PETSC_TRUE) { 1692 if (lambda <= minlambda) { /* bad luck; use full step */ 1693 if (vi->lsmonitor) { 1694 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor,"Line search: Unable to find good step length! %D \n",count);CHKERRQ(ierr); 1695 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor,"Line search: fnorm=%G, gnorm=%G, ynorm=%G, lambda=%G, initial slope=%G\n",fnorm,*gnorm,*ynorm,lambda,initslope);CHKERRQ(ierr); 1696 } 1697 ierr = VecCopy(x,w);CHKERRQ(ierr); 1698 *flag = PETSC_FALSE; 1699 break; 1700 } 1701 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 1702 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1703 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1704 else lambda = lambdatemp; 1705 1706 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1707 if (snes->nfuncs >= snes->max_funcs) { 1708 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 1709 ierr = PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",fnorm,*gnorm,*ynorm,lambda,initslope);CHKERRQ(ierr); 1710 *flag = PETSC_FALSE; 1711 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1712 break; 1713 } 1714 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1715 if (snes->domainerror) { 1716 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1717 PetscFunctionReturn(0); 1718 } 1719 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1720 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1721 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* sufficient reduction */ 1722 if (vi->lsmonitor) { 1723 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line Search: Quadratically determined step, lambda=%G\n",lambda);CHKERRQ(ierr); 1724 } 1725 break; 1726 } 1727 count++; 1728 } 1729 theend2: 1730 /* Optional user-defined check for line search step validity */ 1731 if (vi->postcheckstep) { 1732 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1733 if (changed_y) { 1734 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1735 } 1736 if (changed_y || changed_w) { /* recompute the function if the step has changed */ 1737 ierr = SNESComputeFunction(snes,w,g); 1738 if (snes->domainerror) { 1739 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1740 PetscFunctionReturn(0); 1741 } 1742 ierr = VecNormBegin(g,NORM_2,gnorm);CHKERRQ(ierr); 1743 ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr); 1744 ierr = VecNormEnd(g,NORM_2,gnorm);CHKERRQ(ierr); 1745 ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr); 1746 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1747 } 1748 } 1749 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1750 PetscFunctionReturn(0); 1751 } 1752 1753 typedef PetscErrorCode (*FCN2)(SNES,void*,Vec,Vec,Vec,Vec,Vec,PetscReal,PetscReal,PetscReal*,PetscReal*,PetscBool*); /* force argument to next function to not be extern C*/ 1754 /* -------------------------------------------------------------------------- */ 1755 EXTERN_C_BEGIN 1756 #undef __FUNCT__ 1757 #define __FUNCT__ "SNESLineSearchSet_VI" 1758 PetscErrorCode SNESLineSearchSet_VI(SNES snes,FCN2 func,void *lsctx) 1759 { 1760 PetscFunctionBegin; 1761 ((SNES_VI *)(snes->data))->LineSearch = func; 1762 ((SNES_VI *)(snes->data))->lsP = lsctx; 1763 PetscFunctionReturn(0); 1764 } 1765 EXTERN_C_END 1766 1767 /* -------------------------------------------------------------------------- */ 1768 EXTERN_C_BEGIN 1769 #undef __FUNCT__ 1770 #define __FUNCT__ "SNESLineSearchSetMonitor_VI" 1771 PetscErrorCode SNESLineSearchSetMonitor_VI(SNES snes,PetscBool flg) 1772 { 1773 SNES_VI *vi = (SNES_VI*)snes->data; 1774 PetscErrorCode ierr; 1775 1776 PetscFunctionBegin; 1777 if (flg && !vi->lsmonitor) { 1778 ierr = PetscViewerASCIIMonitorCreate(((PetscObject)snes)->comm,"stdout",((PetscObject)snes)->tablevel,&vi->lsmonitor);CHKERRQ(ierr); 1779 } else if (!flg && vi->lsmonitor) { 1780 ierr = PetscViewerASCIIMonitorDestroy(vi->lsmonitor);CHKERRQ(ierr); 1781 } 1782 PetscFunctionReturn(0); 1783 } 1784 EXTERN_C_END 1785 1786 /* 1787 SNESView_VI - Prints info from the SNESVI data structure. 1788 1789 Input Parameters: 1790 . SNES - the SNES context 1791 . viewer - visualization context 1792 1793 Application Interface Routine: SNESView() 1794 */ 1795 #undef __FUNCT__ 1796 #define __FUNCT__ "SNESView_VI" 1797 static PetscErrorCode SNESView_VI(SNES snes,PetscViewer viewer) 1798 { 1799 SNES_VI *vi = (SNES_VI *)snes->data; 1800 const char *cstr,*tstr; 1801 PetscErrorCode ierr; 1802 PetscBool iascii; 1803 1804 PetscFunctionBegin; 1805 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1806 if (iascii) { 1807 if (vi->LineSearch == SNESLineSearchNo_VI) cstr = "SNESLineSearchNo"; 1808 else if (vi->LineSearch == SNESLineSearchQuadratic_VI) cstr = "SNESLineSearchQuadratic"; 1809 else if (vi->LineSearch == SNESLineSearchCubic_VI) cstr = "SNESLineSearchCubic"; 1810 else if (vi->LineSearch == SNESLineSearchBoundProjectedCubic_VI) cstr = "SNESLineSearchBoundProjectedCubic"; 1811 else if(vi->LineSearch == SNESLineSearchBoundProjectedNo_VI) cstr = "SNESLineSearchBoundProjectedBasic"; 1812 else cstr = "unknown"; 1813 if (snes->ops->solve == SNESSolveVI_SS) tstr = "Semismooth"; 1814 else if (snes->ops->solve == SNESSolveVI_RS) tstr = "Reduced Space"; 1815 else tstr = "unknown"; 1816 ierr = PetscViewerASCIIPrintf(viewer," VI algorithm: %s\n",tstr);CHKERRQ(ierr); 1817 ierr = PetscViewerASCIIPrintf(viewer," line search variant: %s\n",cstr);CHKERRQ(ierr); 1818 ierr = PetscViewerASCIIPrintf(viewer," alpha=%G, maxstep=%G, minlambda=%G\n",vi->alpha,vi->maxstep,vi->minlambda);CHKERRQ(ierr); 1819 } else { 1820 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for SNES EQ VI",((PetscObject)viewer)->type_name); 1821 } 1822 PetscFunctionReturn(0); 1823 } 1824 1825 1826 #undef __FUNCT__ 1827 #define __FUNCT__ "SNESVISetVariableBounds" 1828 /*@ 1829 SNESVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu. 1830 1831 Input Parameters: 1832 . snes - the SNES context. 1833 . xl - lower bound. 1834 . xu - upper bound. 1835 1836 Notes: 1837 If this routine is not called then the lower and upper bounds are set to 1838 PETSC_VI_INF and PETSC_VI_NINF respectively during SNESSetUp(). 1839 1840 @*/ 1841 PetscErrorCode SNESVISetVariableBounds(SNES snes, Vec xl, Vec xu) 1842 { 1843 SNES_VI *vi = (SNES_VI*)snes->data; 1844 1845 PetscFunctionBegin; 1846 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1847 PetscValidHeaderSpecific(xl,VEC_CLASSID,2); 1848 PetscValidHeaderSpecific(xu,VEC_CLASSID,3); 1849 if (!snes->vec_func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call SNESSetFunction() first"); 1850 if (xl->map->N != snes->vec_func->map->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector lengths lower bound = %D solution vector = %D",xl->map->N,snes->vec_func->map->N); 1851 if (xu->map->N != snes->vec_func->map->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector lengths: upper bound = %D solution vector = %D",xu->map->N,snes->vec_func->map->N); 1852 1853 vi->usersetxbounds = PETSC_TRUE; 1854 vi->xl = xl; 1855 vi->xu = xu; 1856 PetscFunctionReturn(0); 1857 } 1858 /* -------------------------------------------------------------------------- */ 1859 /* 1860 SNESSetFromOptions_VI - Sets various parameters for the SNESVI method. 1861 1862 Input Parameter: 1863 . snes - the SNES context 1864 1865 Application Interface Routine: SNESSetFromOptions() 1866 */ 1867 #undef __FUNCT__ 1868 #define __FUNCT__ "SNESSetFromOptions_VI" 1869 static PetscErrorCode SNESSetFromOptions_VI(SNES snes) 1870 { 1871 SNES_VI *vi = (SNES_VI *)snes->data; 1872 const char *lses[] = {"basic","basicnonorms","quadratic","cubic","boundproj_basic"}; 1873 const char *vies[] = {"ss","rs"}; 1874 PetscErrorCode ierr; 1875 PetscInt indx; 1876 PetscBool flg,set,flg2; 1877 1878 PetscFunctionBegin; 1879 ierr = PetscOptionsHead("SNES semismooth method options");CHKERRQ(ierr); 1880 ierr = PetscOptionsBool("-snes_vi_monitor","Monitor all non-active variables","None",PETSC_FALSE,&flg,0);CHKERRQ(ierr); 1881 if (flg) { 1882 ierr = SNESMonitorSet(snes,SNESMonitorVI,0,0);CHKERRQ(ierr); 1883 } 1884 ierr = PetscOptionsReal("-snes_vi_alpha","Function norm must decrease by","None",vi->alpha,&vi->alpha,0);CHKERRQ(ierr); 1885 ierr = PetscOptionsReal("-snes_vi_maxstep","Step must be less than","None",vi->maxstep,&vi->maxstep,0);CHKERRQ(ierr); 1886 ierr = PetscOptionsReal("-snes_vi_minlambda","Minimum lambda allowed","None",vi->minlambda,&vi->minlambda,0);CHKERRQ(ierr); 1887 ierr = PetscOptionsReal("-snes_vi_const_tol","constraint tolerance","None",vi->const_tol,&vi->const_tol,0);CHKERRQ(ierr); 1888 ierr = PetscOptionsBool("-snes_vi_lsmonitor","Print progress of line searches","SNESLineSearchSetMonitor",vi->lsmonitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);CHKERRQ(ierr); 1889 if (set) {ierr = SNESLineSearchSetMonitor(snes,flg);CHKERRQ(ierr);} 1890 ierr = PetscOptionsEList("-snes_vi_type","Semismooth algorithm used","",vies,2,"ss",&indx,&flg2);CHKERRQ(ierr); 1891 if (flg2) { 1892 switch (indx) { 1893 case 0: 1894 snes->ops->solve = SNESSolveVI_SS; 1895 break; 1896 case 1: 1897 snes->ops->solve = SNESSolveVI_RS; 1898 break; 1899 } 1900 } 1901 if (snes->ops->solve == SNESSolveVI_RS) { 1902 ierr = SNESLineSearchSet(snes,SNESLineSearchBoundProjectedCubic_VI,PETSC_NULL);CHKERRQ(ierr); 1903 } 1904 else { 1905 ierr = PetscOptionsEList("-snes_vi_ls","Line search used","SNESLineSearchSet",lses,5,"boundproj_basic",&indx,&flg);CHKERRQ(ierr); 1906 if (flg) { 1907 switch (indx) { 1908 case 0: 1909 ierr = SNESLineSearchSet(snes,SNESLineSearchNo_VI,PETSC_NULL);CHKERRQ(ierr); 1910 break; 1911 case 1: 1912 ierr = SNESLineSearchSet(snes,SNESLineSearchNoNorms_VI,PETSC_NULL);CHKERRQ(ierr); 1913 break; 1914 case 2: 1915 ierr = SNESLineSearchSet(snes,SNESLineSearchQuadratic_VI,PETSC_NULL);CHKERRQ(ierr); 1916 break; 1917 case 3: 1918 ierr = SNESLineSearchSet(snes,SNESLineSearchCubic_VI,PETSC_NULL);CHKERRQ(ierr); 1919 break; 1920 case 4: 1921 ierr = SNESLineSearchSet(snes,SNESLineSearchBoundProjectedNo_VI,PETSC_NULL);CHKERRQ(ierr); 1922 break; 1923 } 1924 } 1925 } 1926 ierr = PetscOptionsTail();CHKERRQ(ierr); 1927 PetscFunctionReturn(0); 1928 } 1929 /* -------------------------------------------------------------------------- */ 1930 /*MC 1931 SNESVI - Semismooth newton method based nonlinear solver that uses a line search 1932 1933 Options Database: 1934 + -snes_vi [cubic,quadratic,basic,basicnonorms] - Selects line search 1935 . -snes_vi_alpha <alpha> - Sets alpha 1936 . -snes_vi_maxstep <maxstep> - Sets the maximum stepsize the line search will use (if the 2-norm(y) > maxstep then scale y to be y = (maxstep/2-norm(y)) *y) 1937 . -snes_vi_minlambda <minlambda> - Sets the minimum lambda the line search will use minlambda / max_i ( y[i]/x[i] ) 1938 - -snes_vi_monitor - print information about progress of line searches 1939 1940 1941 Level: beginner 1942 1943 .seealso: SNESCreate(), SNES, SNESSetType(), SNESTR, SNESLineSearchSet(), 1944 SNESLineSearchSetPostCheck(), SNESLineSearchNo(), SNESLineSearchCubic(), SNESLineSearchQuadratic(), 1945 SNESLineSearchSet(), SNESLineSearchNoNorms(), SNESLineSearchSetPreCheck(), SNESLineSearchSetParams(), SNESLineSearchGetParams() 1946 1947 M*/ 1948 EXTERN_C_BEGIN 1949 #undef __FUNCT__ 1950 #define __FUNCT__ "SNESCreate_VI" 1951 PetscErrorCode SNESCreate_VI(SNES snes) 1952 { 1953 PetscErrorCode ierr; 1954 SNES_VI *vi; 1955 1956 PetscFunctionBegin; 1957 snes->ops->setup = SNESSetUp_VI; 1958 snes->ops->solve = SNESSolveVI_SS; 1959 snes->ops->destroy = SNESDestroy_VI; 1960 snes->ops->setfromoptions = SNESSetFromOptions_VI; 1961 snes->ops->view = SNESView_VI; 1962 snes->ops->converged = SNESDefaultConverged_VI; 1963 1964 ierr = PetscNewLog(snes,SNES_VI,&vi);CHKERRQ(ierr); 1965 snes->data = (void*)vi; 1966 vi->alpha = 1.e-4; 1967 vi->maxstep = 1.e8; 1968 vi->minlambda = 1.e-12; 1969 vi->LineSearch = SNESLineSearchBoundProjectedNo_VI; 1970 vi->lsP = PETSC_NULL; 1971 vi->postcheckstep = PETSC_NULL; 1972 vi->postcheck = PETSC_NULL; 1973 vi->precheckstep = PETSC_NULL; 1974 vi->precheck = PETSC_NULL; 1975 vi->const_tol = 2.2204460492503131e-16; 1976 1977 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetMonitor_C","SNESLineSearchSetMonitor_VI",SNESLineSearchSetMonitor_VI);CHKERRQ(ierr); 1978 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSet_C","SNESLineSearchSet_VI",SNESLineSearchSet_VI);CHKERRQ(ierr); 1979 1980 PetscFunctionReturn(0); 1981 } 1982 EXTERN_C_END 1983