1 2 #include <../src/snes/impls/vi/viimpl.h> /*I "petscsnes.h" I*/ 3 #include <../include/private/kspimpl.h> 4 #include <../include/private/matimpl.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,MPIU_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 fnorm,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 = fnorm*snes->rtol; 137 } 138 if (fnorm != fnorm) { 139 ierr = PetscInfo(snes,"Failed to converged, function norm is NaN\n");CHKERRQ(ierr); 140 *reason = SNES_DIVERGED_FNORM_NAN; 141 } else if (fnorm < snes->abstol) { 142 ierr = PetscInfo2(snes,"Converged due to function norm %G < %G\n",fnorm,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 (fnorm < snes->ttol) { 151 ierr = PetscInfo2(snes,"Converged due to function norm %G < %G (relative tolerance)\n",fnorm,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);CHKERRQ(ierr); 180 ierr = VecNormEnd(phi,NORM_2,phinorm);CHKERRQ(ierr); 181 182 *merit = 0.5*(*phinorm)*(*phinorm); 183 PetscFunctionReturn(0); 184 } 185 186 PETSC_STATIC_INLINE PetscScalar Phi(PetscScalar a,PetscScalar b) 187 { 188 return a + b - sqrt(a*a + b*b); 189 } 190 191 PETSC_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);CHKERRQ(ierr); 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 This file implements a semismooth truncated Newton method with a line search, 405 for solving a system of nonlinear equations in complementarity form, using the KSP, Vec, 406 and Mat interfaces for linear solvers, vectors, and matrices, 407 respectively. 408 409 The following basic routines are required for each nonlinear solver: 410 SNESCreate_XXX() - Creates a nonlinear solver context 411 SNESSetFromOptions_XXX() - Sets runtime options 412 SNESSolve_XXX() - Solves the nonlinear system 413 SNESDestroy_XXX() - Destroys the nonlinear solver context 414 The suffix "_XXX" denotes a particular implementation, in this case 415 we use _VI (e.g., SNESCreate_VI, SNESSolve_VI) for solving 416 systems of nonlinear equations with a line search (LS) method. 417 These routines are actually called via the common user interface 418 routines SNESCreate(), SNESSetFromOptions(), SNESSolve(), and 419 SNESDestroy(), so the application code interface remains identical 420 for all nonlinear solvers. 421 422 Another key routine is: 423 SNESSetUp_XXX() - Prepares for the use of a nonlinear solver 424 by setting data structures and options. The interface routine SNESSetUp() 425 is not usually called directly by the user, but instead is called by 426 SNESSolve() if necessary. 427 428 Additional basic routines are: 429 SNESView_XXX() - Prints details of runtime options that 430 have actually been used. 431 These are called by application codes via the interface routines 432 SNESView(). 433 434 The various types of solvers (preconditioners, Krylov subspace methods, 435 nonlinear solvers, timesteppers) are all organized similarly, so the 436 above description applies to these categories also. 437 438 -------------------------------------------------------------------- */ 439 /* 440 SNESSolveVI_SS - Solves the complementarity problem with a semismooth Newton 441 method using a line search. 442 443 Input Parameters: 444 . snes - the SNES context 445 446 Output Parameter: 447 . outits - number of iterations until termination 448 449 Application Interface Routine: SNESSolve() 450 451 Notes: 452 This implements essentially a semismooth Newton method with a 453 line search. The default line search does not do any line seach 454 but rather takes a full newton step. 455 */ 456 #undef __FUNCT__ 457 #define __FUNCT__ "SNESSolveVI_SS" 458 PetscErrorCode SNESSolveVI_SS(SNES snes) 459 { 460 SNES_VI *vi = (SNES_VI*)snes->data; 461 PetscErrorCode ierr; 462 PetscInt maxits,i,lits; 463 PetscBool lssucceed; 464 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 465 PetscReal gnorm,xnorm=0,ynorm; 466 Vec Y,X,F,G,W; 467 KSPConvergedReason kspreason; 468 469 PetscFunctionBegin; 470 vi->computeuserfunction = snes->ops->computefunction; 471 snes->ops->computefunction = SNESVIComputeFunction; 472 473 snes->numFailures = 0; 474 snes->numLinearSolveFailures = 0; 475 snes->reason = SNES_CONVERGED_ITERATING; 476 477 maxits = snes->max_its; /* maximum number of iterations */ 478 X = snes->vec_sol; /* solution vector */ 479 F = snes->vec_func; /* residual vector */ 480 Y = snes->work[0]; /* work vectors */ 481 G = snes->work[1]; 482 W = snes->work[2]; 483 484 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 485 snes->iter = 0; 486 snes->norm = 0.0; 487 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 488 489 ierr = SNESVIProjectOntoBounds(snes,X);CHKERRQ(ierr); 490 ierr = SNESComputeFunction(snes,X,vi->phi);CHKERRQ(ierr); 491 if (snes->domainerror) { 492 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 493 snes->ops->computefunction = vi->computeuserfunction; 494 PetscFunctionReturn(0); 495 } 496 /* Compute Merit function */ 497 ierr = SNESVIComputeMeritFunction(vi->phi,&vi->merit,&vi->phinorm);CHKERRQ(ierr); 498 499 ierr = VecNormBegin(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm <- ||x|| */ 500 ierr = VecNormEnd(X,NORM_2,&xnorm);CHKERRQ(ierr); 501 if (PetscIsInfOrNanReal(vi->merit)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 502 503 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 504 snes->norm = vi->phinorm; 505 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 506 SNESLogConvHistory(snes,vi->phinorm,0); 507 SNESMonitor(snes,0,vi->phinorm); 508 509 /* set parameter for default relative tolerance convergence test */ 510 snes->ttol = vi->phinorm*snes->rtol; 511 /* test convergence */ 512 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,vi->phinorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 513 if (snes->reason) { 514 snes->ops->computefunction = vi->computeuserfunction; 515 PetscFunctionReturn(0); 516 } 517 518 for (i=0; i<maxits; i++) { 519 520 /* Call general purpose update function */ 521 if (snes->ops->update) { 522 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 523 } 524 525 /* Solve J Y = Phi, where J is the semismooth jacobian */ 526 /* Get the nonlinear function jacobian */ 527 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 528 /* Get the diagonal shift and row scaling vectors */ 529 ierr = SNESVIComputeBsubdifferentialVectors(snes,X,F,snes->jacobian,vi->Da,vi->Db);CHKERRQ(ierr); 530 /* Compute the semismooth jacobian */ 531 ierr = SNESVIComputeJacobian(snes->jacobian,snes->jacobian_pre,vi->Da,vi->Db);CHKERRQ(ierr); 532 /* Compute the merit function gradient */ 533 ierr = SNESVIComputeMeritFunctionGradient(snes->jacobian,vi->phi,vi->dpsi);CHKERRQ(ierr); 534 ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,flg);CHKERRQ(ierr); 535 ierr = SNES_KSPSolve(snes,snes->ksp,vi->phi,Y);CHKERRQ(ierr); 536 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 537 538 if (kspreason < 0) { 539 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 540 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 541 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 542 break; 543 } 544 } 545 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 546 snes->linear_its += lits; 547 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 548 /* 549 if (vi->precheckstep) { 550 PetscBool changed_y = PETSC_FALSE; 551 ierr = (*vi->precheckstep)(snes,X,Y,vi->precheck,&changed_y);CHKERRQ(ierr); 552 } 553 554 if (PetscLogPrintInfo){ 555 ierr = SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr); 556 } 557 */ 558 /* Compute a (scaled) negative update in the line search routine: 559 Y <- X - lambda*Y 560 and evaluate G = function(Y) (depends on the line search). 561 */ 562 ierr = VecCopy(Y,snes->vec_sol_update);CHKERRQ(ierr); 563 ynorm = 1; gnorm = vi->phinorm; 564 ierr = (*vi->LineSearch)(snes,vi->lsP,X,vi->phi,G,Y,W,vi->phinorm,xnorm,&ynorm,&gnorm,&lssucceed);CHKERRQ(ierr); 565 ierr = PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",vi->phinorm,gnorm,ynorm,(int)lssucceed);CHKERRQ(ierr); 566 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 567 if (snes->domainerror) { 568 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 569 snes->ops->computefunction = vi->computeuserfunction; 570 PetscFunctionReturn(0); 571 } 572 if (!lssucceed) { 573 if (++snes->numFailures >= snes->maxFailures) { 574 PetscBool ismin; 575 snes->reason = SNES_DIVERGED_LINE_SEARCH; 576 ierr = SNESVICheckLocalMin_Private(snes,snes->jacobian,G,W,gnorm,&ismin);CHKERRQ(ierr); 577 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 578 break; 579 } 580 } 581 /* Update function and solution vectors */ 582 vi->phinorm = gnorm; 583 vi->merit = 0.5*vi->phinorm*vi->phinorm; 584 ierr = VecCopy(G,vi->phi);CHKERRQ(ierr); 585 ierr = VecCopy(W,X);CHKERRQ(ierr); 586 /* Monitor convergence */ 587 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 588 snes->iter = i+1; 589 snes->norm = vi->phinorm; 590 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 591 SNESLogConvHistory(snes,snes->norm,lits); 592 SNESMonitor(snes,snes->iter,snes->norm); 593 /* Test for convergence, xnorm = || X || */ 594 if (snes->ops->converged != SNESSkipConverged) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 595 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,vi->phinorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 596 if (snes->reason) break; 597 } 598 if (i == maxits) { 599 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 600 if(!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 601 } 602 snes->ops->computefunction = vi->computeuserfunction; 603 PetscFunctionReturn(0); 604 } 605 606 #undef __FUNCT__ 607 #define __FUNCT__ "SNESVIGetActiveSetIS" 608 /* 609 SNESVIGetActiveSetIndices - Gets the global indices for the active set variables 610 611 Input parameter 612 . snes - the SNES context 613 . X - the snes solution vector 614 . F - the nonlinear function vector 615 616 Output parameter 617 . ISact - active set index set 618 */ 619 PetscErrorCode SNESVIGetActiveSetIS(SNES snes,Vec X,Vec F,IS* ISact) 620 { 621 PetscErrorCode ierr; 622 SNES_VI *vi = (SNES_VI*)snes->data; 623 Vec Xl=vi->xl,Xu=vi->xu; 624 const PetscScalar *x,*f,*xl,*xu; 625 PetscInt *idx_act,i,nlocal,nloc_isact=0,ilow,ihigh,i1=0; 626 627 PetscFunctionBegin; 628 ierr = VecGetLocalSize(X,&nlocal);CHKERRQ(ierr); 629 ierr = VecGetOwnershipRange(X,&ilow,&ihigh);CHKERRQ(ierr); 630 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 631 ierr = VecGetArrayRead(Xl,&xl);CHKERRQ(ierr); 632 ierr = VecGetArrayRead(Xu,&xu);CHKERRQ(ierr); 633 ierr = VecGetArrayRead(F,&f);CHKERRQ(ierr); 634 /* Compute active set size */ 635 for (i=0; i < nlocal;i++) { 636 if (!((x[i] > xl[i] + 1.e-8 || (f[i] < 0.0)) && ((x[i] < xu[i] - 1.e-8) || f[i] > 0.0))) nloc_isact++; 637 } 638 639 ierr = PetscMalloc(nloc_isact*sizeof(PetscInt),&idx_act);CHKERRQ(ierr); 640 641 /* Set active set indices */ 642 for(i=0; i < nlocal; i++) { 643 if (!((x[i] > xl[i] + 1.e-8 || (f[i] < 0.0)) && ((x[i] < xu[i] - 1.e-8) || f[i] > 0.0))) idx_act[i1++] = ilow+i; 644 } 645 646 /* Create active set IS */ 647 ierr = ISCreateGeneral(((PetscObject)snes)->comm,nloc_isact,idx_act,PETSC_OWN_POINTER,ISact);CHKERRQ(ierr); 648 649 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 650 ierr = VecRestoreArrayRead(Xl,&xl);CHKERRQ(ierr); 651 ierr = VecRestoreArrayRead(Xu,&xu);CHKERRQ(ierr); 652 ierr = VecRestoreArrayRead(F,&f);CHKERRQ(ierr); 653 PetscFunctionReturn(0); 654 } 655 656 #undef __FUNCT__ 657 #define __FUNCT__ "SNESVICreateIndexSets_RS" 658 PetscErrorCode SNESVICreateIndexSets_RS(SNES snes,Vec X,Vec F,IS* ISact,IS* ISinact) 659 { 660 PetscErrorCode ierr; 661 662 PetscFunctionBegin; 663 ierr = SNESVIGetActiveSetIS(snes,X,F,ISact);CHKERRQ(ierr); 664 ierr = ISComplement(*ISact,X->map->rstart,X->map->rend,ISinact);CHKERRQ(ierr); 665 PetscFunctionReturn(0); 666 } 667 668 /* Create active and inactive set vectors. The local size of this vector is set and petsc computes the global size */ 669 #undef __FUNCT__ 670 #define __FUNCT__ "SNESVICreateSubVectors" 671 PetscErrorCode SNESVICreateSubVectors(SNES snes,PetscInt n,Vec* newv) 672 { 673 PetscErrorCode ierr; 674 Vec v; 675 676 PetscFunctionBegin; 677 ierr = VecCreate(((PetscObject)snes)->comm,&v);CHKERRQ(ierr); 678 ierr = VecSetSizes(v,n,PETSC_DECIDE);CHKERRQ(ierr); 679 ierr = VecSetFromOptions(v);CHKERRQ(ierr); 680 *newv = v; 681 682 PetscFunctionReturn(0); 683 } 684 685 /* Resets the snes PC and KSP when the active set sizes change */ 686 #undef __FUNCT__ 687 #define __FUNCT__ "SNESVIResetPCandKSP" 688 PetscErrorCode SNESVIResetPCandKSP(SNES snes,Mat Amat,Mat Pmat) 689 { 690 PetscErrorCode ierr; 691 KSP kspnew,snesksp; 692 PC pcnew; 693 const MatSolverPackage stype; 694 695 PetscFunctionBegin; 696 /* The active and inactive set sizes have changed so need to create a new snes->ksp object */ 697 ierr = SNESGetKSP(snes,&snesksp);CHKERRQ(ierr); 698 ierr = KSPCreate(((PetscObject)snes)->comm,&kspnew);CHKERRQ(ierr); 699 /* Copy over snes->ksp info */ 700 kspnew->pc_side = snesksp->pc_side; 701 kspnew->rtol = snesksp->rtol; 702 kspnew->abstol = snesksp->abstol; 703 kspnew->max_it = snesksp->max_it; 704 ierr = KSPSetType(kspnew,((PetscObject)snesksp)->type_name);CHKERRQ(ierr); 705 ierr = KSPGetPC(kspnew,&pcnew);CHKERRQ(ierr); 706 ierr = PCSetType(kspnew->pc,((PetscObject)snesksp->pc)->type_name);CHKERRQ(ierr); 707 ierr = PCSetOperators(kspnew->pc,Amat,Pmat,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 708 ierr = PCFactorGetMatSolverPackage(snesksp->pc,&stype);CHKERRQ(ierr); 709 ierr = PCFactorSetMatSolverPackage(kspnew->pc,stype);CHKERRQ(ierr); 710 ierr = KSPDestroy(snesksp);CHKERRQ(ierr); 711 snes->ksp = kspnew; 712 ierr = PetscLogObjectParent(snes,kspnew);CHKERRQ(ierr); 713 ierr = KSPSetFromOptions(kspnew);CHKERRQ(ierr); 714 PetscFunctionReturn(0); 715 } 716 717 718 #undef __FUNCT__ 719 #define __FUNCT__ "SNESVIComputeInactiveSetFnorm" 720 PetscErrorCode SNESVIComputeInactiveSetFnorm(SNES snes,Vec F,Vec X,PetscScalar *fnorm) 721 { 722 PetscErrorCode ierr; 723 SNES_VI *vi = (SNES_VI*)snes->data; 724 const PetscScalar *x,*xl,*xu,*f; 725 PetscInt i,n; 726 PetscReal rnorm; 727 728 PetscFunctionBegin; 729 ierr = VecGetLocalSize(X,&n);CHKERRQ(ierr); 730 ierr = VecGetArrayRead(vi->xl,&xl);CHKERRQ(ierr); 731 ierr = VecGetArrayRead(vi->xu,&xu);CHKERRQ(ierr); 732 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 733 ierr = VecGetArrayRead(F,&f);CHKERRQ(ierr); 734 rnorm = 0.0; 735 for (i=0; i<n; i++) { 736 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]; 737 } 738 ierr = VecRestoreArrayRead(F,&f);CHKERRQ(ierr); 739 ierr = VecRestoreArrayRead(vi->xl,&xl);CHKERRQ(ierr); 740 ierr = VecRestoreArrayRead(vi->xu,&xu);CHKERRQ(ierr); 741 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 742 ierr = MPI_Allreduce(&rnorm,fnorm,1,MPIU_REAL,MPIU_SUM,((PetscObject)snes)->comm);CHKERRQ(ierr); 743 *fnorm = sqrt(*fnorm); 744 PetscFunctionReturn(0); 745 } 746 747 /* Variational Inequality solver using reduce space method. No semismooth algorithm is 748 implemented in this algorithm. It basically identifies the active variables and does 749 a linear solve on the inactive variables. */ 750 #undef __FUNCT__ 751 #define __FUNCT__ "SNESSolveVI_RS" 752 PetscErrorCode SNESSolveVI_RS(SNES snes) 753 { 754 SNES_VI *vi = (SNES_VI*)snes->data; 755 PetscErrorCode ierr; 756 PetscInt maxits,i,lits; 757 PetscBool lssucceed; 758 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 759 PetscReal fnorm,gnorm,xnorm=0,ynorm; 760 Vec Y,X,F,G,W; 761 KSPConvergedReason kspreason; 762 763 PetscFunctionBegin; 764 snes->numFailures = 0; 765 snes->numLinearSolveFailures = 0; 766 snes->reason = SNES_CONVERGED_ITERATING; 767 768 maxits = snes->max_its; /* maximum number of iterations */ 769 X = snes->vec_sol; /* solution vector */ 770 F = snes->vec_func; /* residual vector */ 771 Y = snes->work[0]; /* work vectors */ 772 G = snes->work[1]; 773 W = snes->work[2]; 774 775 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 776 snes->iter = 0; 777 snes->norm = 0.0; 778 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 779 780 ierr = SNESVIProjectOntoBounds(snes,X);CHKERRQ(ierr); 781 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 782 if (snes->domainerror) { 783 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 784 PetscFunctionReturn(0); 785 } 786 ierr = SNESVIComputeInactiveSetFnorm(snes,F,X,&fnorm);CHKERRQ(ierr); 787 ierr = VecNormBegin(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm <- ||x|| */ 788 ierr = VecNormEnd(X,NORM_2,&xnorm);CHKERRQ(ierr); 789 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 790 791 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 792 snes->norm = fnorm; 793 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 794 SNESLogConvHistory(snes,fnorm,0); 795 SNESMonitor(snes,0,fnorm); 796 797 /* set parameter for default relative tolerance convergence test */ 798 snes->ttol = fnorm*snes->rtol; 799 /* test convergence */ 800 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 801 if (snes->reason) PetscFunctionReturn(0); 802 803 for (i=0; i<maxits; i++) { 804 805 IS IS_act,IS_inact; /* _act -> active set _inact -> inactive set */ 806 VecScatter scat_act,scat_inact; 807 PetscInt nis_act,nis_inact; 808 Vec Y_act,Y_inact,F_inact; 809 Mat jac_inact_inact,prejac_inact_inact; 810 IS keptrows; 811 PetscBool isequal; 812 813 /* Call general purpose update function */ 814 if (snes->ops->update) { 815 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 816 } 817 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 818 819 /* Create active and inactive index sets */ 820 ierr = SNESVICreateIndexSets_RS(snes,X,F,&IS_act,&IS_inact);CHKERRQ(ierr); 821 822 /* Create inactive set submatrix */ 823 ierr = MatGetSubMatrix(snes->jacobian,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);CHKERRQ(ierr); 824 ierr = MatFindNonzeroRows(jac_inact_inact,&keptrows);CHKERRQ(ierr); 825 if (keptrows) { 826 PetscInt cnt,*nrows,k; 827 const PetscInt *krows,*inact; 828 PetscInt rstart=jac_inact_inact->rmap->rstart; 829 830 ierr = MatDestroy(jac_inact_inact);CHKERRQ(ierr); 831 ierr = ISDestroy(IS_act);CHKERRQ(ierr); 832 833 ierr = ISGetLocalSize(keptrows,&cnt);CHKERRQ(ierr); 834 ierr = ISGetIndices(keptrows,&krows);CHKERRQ(ierr); 835 ierr = ISGetIndices(IS_inact,&inact);CHKERRQ(ierr); 836 ierr = PetscMalloc(cnt*sizeof(PetscInt),&nrows);CHKERRQ(ierr); 837 for (k=0; k<cnt; k++) { 838 nrows[k] = inact[krows[k]-rstart]; 839 } 840 ierr = ISRestoreIndices(keptrows,&krows);CHKERRQ(ierr); 841 ierr = ISRestoreIndices(IS_inact,&inact);CHKERRQ(ierr); 842 ierr = ISDestroy(keptrows);CHKERRQ(ierr); 843 ierr = ISDestroy(IS_inact);CHKERRQ(ierr); 844 845 ierr = ISCreateGeneral(PETSC_COMM_WORLD,cnt,nrows,PETSC_OWN_POINTER,&IS_inact);CHKERRQ(ierr); 846 ierr = ISComplement(IS_inact,F->map->rstart,F->map->rend,&IS_act);CHKERRQ(ierr); 847 ierr = MatGetSubMatrix(snes->jacobian,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);CHKERRQ(ierr); 848 } 849 850 /* Get sizes of active and inactive sets */ 851 ierr = ISGetLocalSize(IS_act,&nis_act);CHKERRQ(ierr); 852 ierr = ISGetLocalSize(IS_inact,&nis_inact);CHKERRQ(ierr); 853 854 /* Create active and inactive set vectors */ 855 ierr = SNESVICreateSubVectors(snes,nis_inact,&F_inact);CHKERRQ(ierr); 856 ierr = SNESVICreateSubVectors(snes,nis_act,&Y_act);CHKERRQ(ierr); 857 ierr = SNESVICreateSubVectors(snes,nis_inact,&Y_inact);CHKERRQ(ierr); 858 859 /* Create scatter contexts */ 860 ierr = VecScatterCreate(Y,IS_act,Y_act,PETSC_NULL,&scat_act);CHKERRQ(ierr); 861 ierr = VecScatterCreate(Y,IS_inact,Y_inact,PETSC_NULL,&scat_inact);CHKERRQ(ierr); 862 863 /* Do a vec scatter to active and inactive set vectors */ 864 ierr = VecScatterBegin(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 865 ierr = VecScatterEnd(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 866 867 ierr = VecScatterBegin(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 868 ierr = VecScatterEnd(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 869 870 ierr = VecScatterBegin(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 871 ierr = VecScatterEnd(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 872 873 /* Active set direction = 0 */ 874 ierr = VecSet(Y_act,0);CHKERRQ(ierr); 875 if (snes->jacobian != snes->jacobian_pre) { 876 ierr = MatGetSubMatrix(snes->jacobian_pre,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&prejac_inact_inact);CHKERRQ(ierr); 877 } else prejac_inact_inact = jac_inact_inact; 878 879 ierr = ISEqual(vi->IS_inact_prev,IS_inact,&isequal);CHKERRQ(ierr); 880 if (!isequal) { 881 ierr = SNESVIResetPCandKSP(snes,jac_inact_inact,prejac_inact_inact);CHKERRQ(ierr); 882 } 883 884 /* ierr = ISView(IS_inact,0);CHKERRQ(ierr); */ 885 /* ierr = ISView(IS_act,0);CHKERRQ(ierr);*/ 886 /* ierr = MatView(snes->jacobian_pre,0); */ 887 888 ierr = KSPSetOperators(snes->ksp,jac_inact_inact,prejac_inact_inact,flg);CHKERRQ(ierr); 889 ierr = KSPSetUp(snes->ksp);CHKERRQ(ierr); 890 { 891 PC pc; 892 PetscBool flg; 893 ierr = KSPGetPC(snes->ksp,&pc);CHKERRQ(ierr); 894 ierr = PetscTypeCompare((PetscObject)pc,PCFIELDSPLIT,&flg);CHKERRQ(ierr); 895 if (flg) { 896 KSP *subksps; 897 ierr = PCFieldSplitGetSubKSP(pc,PETSC_NULL,&subksps);CHKERRQ(ierr); 898 ierr = KSPGetPC(subksps[0],&pc);CHKERRQ(ierr); 899 ierr = PetscFree(subksps);CHKERRQ(ierr); 900 ierr = PetscTypeCompare((PetscObject)pc,PCBJACOBI,&flg);CHKERRQ(ierr); 901 if (flg) { 902 PetscInt n,N = 101*101,j,cnts[3] = {0,0,0}; 903 const PetscInt *ii; 904 905 ierr = ISGetSize(IS_inact,&n);CHKERRQ(ierr); 906 ierr = ISGetIndices(IS_inact,&ii);CHKERRQ(ierr); 907 for (j=0; j<n; j++) { 908 if (ii[j] < N) cnts[0]++; 909 else if (ii[j] < 2*N) cnts[1]++; 910 else if (ii[j] < 3*N) cnts[2]++; 911 } 912 ierr = ISRestoreIndices(IS_inact,&ii);CHKERRQ(ierr); 913 914 ierr = PCBJacobiSetTotalBlocks(pc,3,cnts);CHKERRQ(ierr); 915 } 916 } 917 } 918 919 ierr = SNES_KSPSolve(snes,snes->ksp,F_inact,Y_inact);CHKERRQ(ierr); 920 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 921 if (kspreason < 0) { 922 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 923 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 924 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 925 break; 926 } 927 } 928 929 ierr = VecScatterBegin(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 930 ierr = VecScatterEnd(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 931 ierr = VecScatterBegin(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 932 ierr = VecScatterEnd(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 933 934 ierr = VecDestroy(F_inact);CHKERRQ(ierr); 935 ierr = VecDestroy(Y_act);CHKERRQ(ierr); 936 ierr = VecDestroy(Y_inact);CHKERRQ(ierr); 937 ierr = VecScatterDestroy(scat_act);CHKERRQ(ierr); 938 ierr = VecScatterDestroy(scat_inact);CHKERRQ(ierr); 939 ierr = ISDestroy(IS_act);CHKERRQ(ierr); 940 if (!isequal) { 941 ierr = ISDestroy(vi->IS_inact_prev);CHKERRQ(ierr); 942 ierr = ISDuplicate(IS_inact,&vi->IS_inact_prev);CHKERRQ(ierr); 943 } 944 ierr = ISDestroy(IS_inact);CHKERRQ(ierr); 945 ierr = MatDestroy(jac_inact_inact);CHKERRQ(ierr); 946 if (snes->jacobian != snes->jacobian_pre) { 947 ierr = MatDestroy(prejac_inact_inact);CHKERRQ(ierr); 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 #undef __FUNCT__ 1008 #define __FUNCT__ "SNESVISetRedundancyCheck" 1009 PetscErrorCode SNESVISetRedundancyCheck(SNES snes,PetscErrorCode (*func)(SNES,IS,IS*,void*)) 1010 { 1011 SNES_VI *vi = (SNES_VI*)snes->data; 1012 1013 PetscFunctionBegin; 1014 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1015 vi->checkredundancy = func; 1016 PetscFunctionReturn(0); 1017 } 1018 1019 #if defined(PETSC_HAVE_MATLAB_ENGINE) 1020 #include <engine.h> 1021 #include <mex.h> 1022 typedef struct {char *funcname; mxArray *ctx;} SNESVIMatlabContext; 1023 1024 #undef __FUNCT__ 1025 #define __FUNCT__ "SNESVIRedundancyCheck_Matlab" 1026 PetscErrorCode SNESVIRedundancyCheck_Matlab(SNES snes,IS is_act,IS* is_redact,void* ctx) 1027 { 1028 PetscErrorCode ierr; 1029 SNESVIMatlabContext *sctx = (SNESVIMatlabContext*)ctx; 1030 int nlhs = 2, nrhs = 4; 1031 mxArray *plhs[2], *prhs[4]; 1032 long long int l1 = 0,ls = 0; 1033 1034 PetscFunctionBegin; 1035 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1036 PetscValidHeaderSpecific(is_act,IS_CLASSID,2); 1037 PetscValidPointer(is_redact,3); 1038 PetscCheckSameComm(snes,1,is_act,2); 1039 1040 /* call Matlab function in ctx with arguments is_act */ 1041 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 1042 ierr = PetscMemcpy(&l1,&is_act,sizeof(is_act));CHKERRQ(ierr); 1043 prhs[0] = mxCreateDoubleScalar((double)ls); 1044 prhs[1] = mxCreateDoubleScalar((double)l1); 1045 prhs[2] = mxCreateString(sctx->funcname); 1046 prhs[3] = sctx->ctx; 1047 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESVIRedundancyCheckInternal");CHKERRQ(ierr); 1048 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 1049 *is_redact = (IS)mxGetPr(plhs[1]); 1050 mxDestroyArray(prhs[0]); 1051 mxDestroyArray(prhs[1]); 1052 mxDestroyArray(prhs[2]); 1053 mxDestroyArray(prhs[3]); 1054 mxDestroyArray(plhs[0]); 1055 PetscFunctionReturn(0); 1056 } 1057 1058 #undef __FUNCT__ 1059 #define __FUNCT__ "SNESVISetRedundancyCheckMatlab" 1060 PetscErrorCode SNESVISetRedundancyCheckMatlab(SNES snes,const char* func,mxArray* ctx) 1061 { 1062 PetscErrorCode ierr; 1063 SNESVIMatlabContext *sctx; 1064 SNES_VI *vi = (SNES_VI*)snes->data; 1065 1066 PetscFunctionBegin; 1067 /* currently sctx is memory bleed */ 1068 ierr = PetscMalloc(sizeof(SNESVIMatlabContext),&sctx);CHKERRQ(ierr); 1069 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 1070 1071 sctx->ctx = mxDuplicateArray(ctx); 1072 vi->checkredundancy = SNESVIRedundancyCheck_Matlab; 1073 PetscFunctionReturn(0); 1074 } 1075 1076 #endif 1077 1078 1079 /* Variational Inequality solver using augmented space method. It does the opposite of the 1080 reduced space method i.e. it identifies the active set variables and instead of discarding 1081 them it augments the original system by introducing additional equality 1082 constraint equations for active set variables. The user can optionally provide an IS for 1083 a subset of 'redundant' active set variables which will treated as free variables. 1084 Specific implementation for Allen-Cahn problem 1085 */ 1086 #undef __FUNCT__ 1087 #define __FUNCT__ "SNESSolveVI_RS2" 1088 PetscErrorCode SNESSolveVI_RS2(SNES snes) 1089 { 1090 SNES_VI *vi = (SNES_VI*)snes->data; 1091 PetscErrorCode ierr; 1092 PetscInt maxits,i,lits; 1093 PetscBool lssucceed; 1094 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 1095 PetscReal fnorm,gnorm,xnorm=0,ynorm; 1096 Vec Y,X,F,G,W; 1097 KSPConvergedReason kspreason; 1098 1099 PetscFunctionBegin; 1100 snes->numFailures = 0; 1101 snes->numLinearSolveFailures = 0; 1102 snes->reason = SNES_CONVERGED_ITERATING; 1103 1104 maxits = snes->max_its; /* maximum number of iterations */ 1105 X = snes->vec_sol; /* solution vector */ 1106 F = snes->vec_func; /* residual vector */ 1107 Y = snes->work[0]; /* work vectors */ 1108 G = snes->work[1]; 1109 W = snes->work[2]; 1110 1111 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1112 snes->iter = 0; 1113 snes->norm = 0.0; 1114 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1115 1116 ierr = SNESVIProjectOntoBounds(snes,X);CHKERRQ(ierr); 1117 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 1118 if (snes->domainerror) { 1119 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 1120 PetscFunctionReturn(0); 1121 } 1122 ierr = SNESVIComputeInactiveSetFnorm(snes,F,X,&fnorm);CHKERRQ(ierr); 1123 ierr = VecNormBegin(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm <- ||x|| */ 1124 ierr = VecNormEnd(X,NORM_2,&xnorm);CHKERRQ(ierr); 1125 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1126 1127 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1128 snes->norm = fnorm; 1129 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1130 SNESLogConvHistory(snes,fnorm,0); 1131 SNESMonitor(snes,0,fnorm); 1132 1133 /* set parameter for default relative tolerance convergence test */ 1134 snes->ttol = fnorm*snes->rtol; 1135 /* test convergence */ 1136 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 1137 if (snes->reason) PetscFunctionReturn(0); 1138 1139 for (i=0; i<maxits; i++) { 1140 IS IS_act,IS_inact; /* _act -> active set _inact -> inactive set */ 1141 IS IS_redact=PETSC_NULL; /* redundant active set */ 1142 Mat J_aug,Jpre_aug; 1143 Vec F_aug,Y_aug; 1144 PetscInt nis_redact,nis_act; 1145 const PetscInt *idx_redact,*idx_act; 1146 PetscInt k; 1147 PetscInt *idx_actkept=PETSC_NULL,nkept=0; /* list of kept active set */ 1148 PetscScalar *f,*f2; 1149 PetscBool isequal; 1150 PetscInt i1=0,j1=0; 1151 1152 /* Call general purpose update function */ 1153 if (snes->ops->update) { 1154 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 1155 } 1156 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 1157 1158 /* Create active and inactive index sets */ 1159 ierr = SNESVICreateIndexSets_RS(snes,X,F,&IS_act,&IS_inact);CHKERRQ(ierr); 1160 1161 /* Get local active set size */ 1162 ierr = ISGetLocalSize(IS_act,&nis_act);CHKERRQ(ierr); 1163 ierr = ISGetIndices(IS_act,&idx_act);CHKERRQ(ierr); 1164 if(nis_act) { 1165 if(vi->checkredundancy) { 1166 (*vi->checkredundancy)(snes,IS_act,&IS_redact,snes->funP); 1167 } 1168 1169 if(!IS_redact) { 1170 /* User called checkredundancy function but didn't create IS_redact because 1171 there were no redundant active set variables */ 1172 /* Copy over all active set indices to the list */ 1173 ierr = PetscMalloc(nis_act*sizeof(PetscInt),&idx_actkept);CHKERRQ(ierr); 1174 for(k=0;k < nis_act;k++) idx_actkept[k] = idx_act[k]; 1175 nkept = nis_act; 1176 } else { 1177 ierr = ISGetLocalSize(IS_redact,&nis_redact);CHKERRQ(ierr); 1178 ierr = PetscMalloc((nis_act-nis_redact)*sizeof(PetscInt),&idx_actkept);CHKERRQ(ierr); 1179 1180 /* Create reduced active set list */ 1181 ierr = ISGetIndices(IS_act,&idx_act);CHKERRQ(ierr); 1182 ierr = ISGetIndices(IS_redact,&idx_redact);CHKERRQ(ierr); 1183 j1 = 0; 1184 for(k=0;k<nis_act;k++) { 1185 if(j1 < nis_redact && idx_act[k] == idx_redact[j1]) j1++; 1186 else idx_actkept[nkept++] = idx_act[k]; 1187 } 1188 ierr = ISRestoreIndices(IS_act,&idx_act);CHKERRQ(ierr); 1189 ierr = ISRestoreIndices(IS_redact,&idx_redact);CHKERRQ(ierr); 1190 } 1191 1192 /* Create augmented F and Y */ 1193 ierr = VecCreate(((PetscObject)snes)->comm,&F_aug);CHKERRQ(ierr); 1194 ierr = VecSetSizes(F_aug,F->map->n+nkept,PETSC_DECIDE);CHKERRQ(ierr); 1195 ierr = VecSetFromOptions(F_aug);CHKERRQ(ierr); 1196 ierr = VecDuplicate(F_aug,&Y_aug);CHKERRQ(ierr); 1197 1198 /* Copy over F to F_aug in the first n locations */ 1199 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 1200 ierr = VecGetArray(F_aug,&f2);CHKERRQ(ierr); 1201 ierr = PetscMemcpy(f2,f,F->map->n*sizeof(PetscScalar));CHKERRQ(ierr); 1202 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 1203 ierr = VecRestoreArray(F_aug,&f2);CHKERRQ(ierr); 1204 1205 /* Create the augmented jacobian matrix */ 1206 ierr = MatCreate(((PetscObject)snes)->comm,&J_aug);CHKERRQ(ierr); 1207 ierr = MatSetSizes(J_aug,X->map->n+nkept,X->map->n+nkept,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 1208 ierr = MatSetFromOptions(J_aug);CHKERRQ(ierr); 1209 1210 /* set preallocation info..will add this later */ 1211 1212 /* Set values in the augmented matrix...Doing only sequential case first*/ 1213 PetscInt ncols; 1214 const PetscInt *cols; 1215 const PetscScalar *vals; 1216 PetscScalar value=1.0; 1217 PetscInt row,col; 1218 1219 /* Set the original jacobian matrix in J_aug */ 1220 for(row=0;row<snes->jacobian->rmap->n;row++) { 1221 ierr = MatGetRow(snes->jacobian,row,&ncols,&cols,&vals);CHKERRQ(ierr); 1222 ierr = MatSetValues(J_aug,1,&row,ncols,cols,vals,INSERT_VALUES);CHKERRQ(ierr); 1223 ierr = MatRestoreRow(snes->jacobian,row,&ncols,&cols,&vals);CHKERRQ(ierr); 1224 } 1225 /* Add the augmented part */ 1226 for(k=0;k<nkept;k++) { 1227 row = idx_actkept[k]; 1228 col = snes->jacobian->rmap->n+k; 1229 ierr = MatSetValues(J_aug,1,&row,1,&col,&value,INSERT_VALUES);CHKERRQ(ierr); 1230 ierr = MatSetValues(J_aug,1,&col,1,&row,&value,INSERT_VALUES);CHKERRQ(ierr); 1231 } 1232 ierr = MatAssemblyBegin(J_aug,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1233 ierr = MatAssemblyEnd(J_aug,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1234 /* Only considering prejac = jac for now */ 1235 Jpre_aug = J_aug; 1236 1237 } else { 1238 F_aug = F; J_aug = snes->jacobian; Y_aug = Y; Jpre_aug = snes->jacobian_pre; 1239 } 1240 1241 ierr = ISEqual(vi->IS_inact_prev,IS_inact,&isequal);CHKERRQ(ierr); 1242 if (!isequal) { 1243 ierr = SNESVIResetPCandKSP(snes,J_aug,Jpre_aug);CHKERRQ(ierr); 1244 } 1245 ierr = KSPSetOperators(snes->ksp,J_aug,Jpre_aug,flg);CHKERRQ(ierr); 1246 ierr = KSPSetUp(snes->ksp);CHKERRQ(ierr); 1247 /* { 1248 PC pc; 1249 PetscBool flg; 1250 ierr = KSPGetPC(snes->ksp,&pc);CHKERRQ(ierr); 1251 ierr = PetscTypeCompare((PetscObject)pc,PCFIELDSPLIT,&flg);CHKERRQ(ierr); 1252 if (flg) { 1253 KSP *subksps; 1254 ierr = PCFieldSplitGetSubKSP(pc,PETSC_NULL,&subksps);CHKERRQ(ierr); 1255 ierr = KSPGetPC(subksps[0],&pc);CHKERRQ(ierr); 1256 ierr = PetscFree(subksps);CHKERRQ(ierr); 1257 ierr = PetscTypeCompare((PetscObject)pc,PCBJACOBI,&flg);CHKERRQ(ierr); 1258 if (flg) { 1259 PetscInt n,N = 101*101,j,cnts[3] = {0,0,0}; 1260 const PetscInt *ii; 1261 1262 ierr = ISGetSize(IS_inact,&n);CHKERRQ(ierr); 1263 ierr = ISGetIndices(IS_inact,&ii);CHKERRQ(ierr); 1264 for (j=0; j<n; j++) { 1265 if (ii[j] < N) cnts[0]++; 1266 else if (ii[j] < 2*N) cnts[1]++; 1267 else if (ii[j] < 3*N) cnts[2]++; 1268 } 1269 ierr = ISRestoreIndices(IS_inact,&ii);CHKERRQ(ierr); 1270 1271 ierr = PCBJacobiSetTotalBlocks(pc,3,cnts);CHKERRQ(ierr); 1272 } 1273 } 1274 } 1275 */ 1276 ierr = SNES_KSPSolve(snes,snes->ksp,F_aug,Y_aug);CHKERRQ(ierr); 1277 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 1278 if (kspreason < 0) { 1279 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 1280 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 1281 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 1282 break; 1283 } 1284 } 1285 1286 if(nis_act) { 1287 PetscScalar *y1,*y2; 1288 ierr = VecGetArray(Y,&y1);CHKERRQ(ierr); 1289 ierr = VecGetArray(Y_aug,&y2);CHKERRQ(ierr); 1290 /* Copy over inactive Y_aug to Y */ 1291 j1 = 0; 1292 for(i1=Y->map->rstart;i1 < Y->map->rend;i1++) { 1293 if(j1 < nkept && idx_actkept[j1] == i1) j1++; 1294 else y1[i1-Y->map->rstart] = y2[i1-Y->map->rstart]; 1295 } 1296 ierr = VecRestoreArray(Y,&y1);CHKERRQ(ierr); 1297 ierr = VecRestoreArray(Y_aug,&y2);CHKERRQ(ierr); 1298 1299 ierr = VecDestroy(F_aug);CHKERRQ(ierr); 1300 ierr = VecDestroy(Y_aug);CHKERRQ(ierr); 1301 ierr = MatDestroy(J_aug);CHKERRQ(ierr); 1302 ierr = PetscFree(idx_actkept);CHKERRQ(ierr); 1303 } 1304 1305 if (!isequal) { 1306 ierr = ISDestroy(vi->IS_inact_prev);CHKERRQ(ierr); 1307 ierr = ISDuplicate(IS_inact,&vi->IS_inact_prev);CHKERRQ(ierr); 1308 } 1309 ierr = ISDestroy(IS_inact);CHKERRQ(ierr); 1310 1311 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 1312 snes->linear_its += lits; 1313 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 1314 /* 1315 if (vi->precheckstep) { 1316 PetscBool changed_y = PETSC_FALSE; 1317 ierr = (*vi->precheckstep)(snes,X,Y,vi->precheck,&changed_y);CHKERRQ(ierr); 1318 } 1319 1320 if (PetscLogPrintInfo){ 1321 ierr = SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr); 1322 } 1323 */ 1324 /* Compute a (scaled) negative update in the line search routine: 1325 Y <- X - lambda*Y 1326 and evaluate G = function(Y) (depends on the line search). 1327 */ 1328 ierr = VecCopy(Y,snes->vec_sol_update);CHKERRQ(ierr); 1329 ynorm = 1; gnorm = fnorm; 1330 ierr = (*vi->LineSearch)(snes,vi->lsP,X,F,G,Y,W,fnorm,xnorm,&ynorm,&gnorm,&lssucceed);CHKERRQ(ierr); 1331 ierr = PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",fnorm,gnorm,ynorm,(int)lssucceed);CHKERRQ(ierr); 1332 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 1333 if (snes->domainerror) { 1334 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 1335 PetscFunctionReturn(0); 1336 } 1337 if (!lssucceed) { 1338 if (++snes->numFailures >= snes->maxFailures) { 1339 PetscBool ismin; 1340 snes->reason = SNES_DIVERGED_LINE_SEARCH; 1341 ierr = SNESVICheckLocalMin_Private(snes,snes->jacobian,G,W,gnorm,&ismin);CHKERRQ(ierr); 1342 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 1343 break; 1344 } 1345 } 1346 /* Update function and solution vectors */ 1347 fnorm = gnorm; 1348 ierr = VecCopy(G,F);CHKERRQ(ierr); 1349 ierr = VecCopy(W,X);CHKERRQ(ierr); 1350 /* Monitor convergence */ 1351 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1352 snes->iter = i+1; 1353 snes->norm = fnorm; 1354 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1355 SNESLogConvHistory(snes,snes->norm,lits); 1356 SNESMonitor(snes,snes->iter,snes->norm); 1357 /* Test for convergence, xnorm = || X || */ 1358 if (snes->ops->converged != SNESSkipConverged) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 1359 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 1360 if (snes->reason) break; 1361 } 1362 if (i == maxits) { 1363 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 1364 if(!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 1365 } 1366 PetscFunctionReturn(0); 1367 } 1368 1369 /* -------------------------------------------------------------------------- */ 1370 /* 1371 SNESSetUp_VI - Sets up the internal data structures for the later use 1372 of the SNESVI nonlinear solver. 1373 1374 Input Parameter: 1375 . snes - the SNES context 1376 . x - the solution vector 1377 1378 Application Interface Routine: SNESSetUp() 1379 1380 Notes: 1381 For basic use of the SNES solvers, the user need not explicitly call 1382 SNESSetUp(), since these actions will automatically occur during 1383 the call to SNESSolve(). 1384 */ 1385 #undef __FUNCT__ 1386 #define __FUNCT__ "SNESSetUp_VI" 1387 PetscErrorCode SNESSetUp_VI(SNES snes) 1388 { 1389 PetscErrorCode ierr; 1390 SNES_VI *vi = (SNES_VI*) snes->data; 1391 PetscInt i_start[3],i_end[3]; 1392 1393 PetscFunctionBegin; 1394 if (!snes->vec_sol_update) { 1395 ierr = VecDuplicate(snes->vec_sol,&snes->vec_sol_update);CHKERRQ(ierr); 1396 ierr = PetscLogObjectParent(snes,snes->vec_sol_update);CHKERRQ(ierr); 1397 } 1398 if (!snes->work) { 1399 snes->nwork = 3; 1400 ierr = VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work);CHKERRQ(ierr); 1401 ierr = PetscLogObjectParents(snes,snes->nwork,snes->work);CHKERRQ(ierr); 1402 } 1403 1404 /* If the lower and upper bound on variables are not set, set it to 1405 -Inf and Inf */ 1406 if (!vi->xl && !vi->xu) { 1407 vi->usersetxbounds = PETSC_FALSE; 1408 ierr = VecDuplicate(snes->vec_sol, &vi->xl); CHKERRQ(ierr); 1409 ierr = VecSet(vi->xl,PETSC_VI_NINF);CHKERRQ(ierr); 1410 ierr = VecDuplicate(snes->vec_sol, &vi->xu); CHKERRQ(ierr); 1411 ierr = VecSet(vi->xu,PETSC_VI_INF);CHKERRQ(ierr); 1412 } else { 1413 /* Check if lower bound, upper bound and solution vector distribution across the processors is identical */ 1414 ierr = VecGetOwnershipRange(snes->vec_sol,i_start,i_end);CHKERRQ(ierr); 1415 ierr = VecGetOwnershipRange(vi->xl,i_start+1,i_end+1);CHKERRQ(ierr); 1416 ierr = VecGetOwnershipRange(vi->xu,i_start+2,i_end+2);CHKERRQ(ierr); 1417 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])) 1418 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."); 1419 } 1420 if (snes->ops->solve != SNESSolveVI_SS) { 1421 /* Set up previous active index set for the first snes solve 1422 vi->IS_inact_prev = 0,1,2,....N */ 1423 PetscInt *indices; 1424 PetscInt i,n,rstart,rend; 1425 1426 ierr = VecGetOwnershipRange(snes->vec_sol,&rstart,&rend);CHKERRQ(ierr); 1427 ierr = VecGetLocalSize(snes->vec_sol,&n);CHKERRQ(ierr); 1428 ierr = PetscMalloc(n*sizeof(PetscInt),&indices);CHKERRQ(ierr); 1429 for(i=0;i < n; i++) indices[i] = rstart + i; 1430 ierr = ISCreateGeneral(((PetscObject)snes)->comm,n,indices,PETSC_OWN_POINTER,&vi->IS_inact_prev); 1431 } 1432 1433 if (snes->ops->solve == SNESSolveVI_SS) { 1434 ierr = VecDuplicate(snes->vec_sol, &vi->dpsi);CHKERRQ(ierr); 1435 ierr = VecDuplicate(snes->vec_sol, &vi->phi); CHKERRQ(ierr); 1436 ierr = VecDuplicate(snes->vec_sol, &vi->Da); CHKERRQ(ierr); 1437 ierr = VecDuplicate(snes->vec_sol, &vi->Db); CHKERRQ(ierr); 1438 ierr = VecDuplicate(snes->vec_sol, &vi->z);CHKERRQ(ierr); 1439 ierr = VecDuplicate(snes->vec_sol, &vi->t); CHKERRQ(ierr); 1440 } 1441 PetscFunctionReturn(0); 1442 } 1443 /* -------------------------------------------------------------------------- */ 1444 /* 1445 SNESDestroy_VI - Destroys the private SNES_VI context that was created 1446 with SNESCreate_VI(). 1447 1448 Input Parameter: 1449 . snes - the SNES context 1450 1451 Application Interface Routine: SNESDestroy() 1452 */ 1453 #undef __FUNCT__ 1454 #define __FUNCT__ "SNESDestroy_VI" 1455 PetscErrorCode SNESDestroy_VI(SNES snes) 1456 { 1457 SNES_VI *vi = (SNES_VI*) snes->data; 1458 PetscErrorCode ierr; 1459 1460 PetscFunctionBegin; 1461 if (snes->vec_sol_update) { 1462 ierr = VecDestroy(snes->vec_sol_update);CHKERRQ(ierr); 1463 snes->vec_sol_update = PETSC_NULL; 1464 } 1465 if (snes->work) { 1466 ierr = VecDestroyVecs(snes->nwork,&snes->work);CHKERRQ(ierr); 1467 } 1468 if (snes->ops->solve != SNESSolveVI_SS) { 1469 ierr = ISDestroy(vi->IS_inact_prev); 1470 } 1471 1472 if (snes->ops->solve == SNESSolveVI_SS) { 1473 /* clear vectors */ 1474 ierr = VecDestroy(vi->dpsi);CHKERRQ(ierr); 1475 ierr = VecDestroy(vi->phi); CHKERRQ(ierr); 1476 ierr = VecDestroy(vi->Da); CHKERRQ(ierr); 1477 ierr = VecDestroy(vi->Db); CHKERRQ(ierr); 1478 ierr = VecDestroy(vi->z); CHKERRQ(ierr); 1479 ierr = VecDestroy(vi->t); CHKERRQ(ierr); 1480 } 1481 1482 if (!vi->usersetxbounds) { 1483 ierr = VecDestroy(vi->xl); CHKERRQ(ierr); 1484 ierr = VecDestroy(vi->xu); CHKERRQ(ierr); 1485 } 1486 1487 if (vi->lsmonitor) { 1488 ierr = PetscViewerASCIIMonitorDestroy(vi->lsmonitor);CHKERRQ(ierr); 1489 } 1490 ierr = PetscFree(snes->data);CHKERRQ(ierr); 1491 1492 /* clear composed functions */ 1493 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSet_C","",PETSC_NULL);CHKERRQ(ierr); 1494 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetMonitor_C","",PETSC_NULL);CHKERRQ(ierr); 1495 PetscFunctionReturn(0); 1496 } 1497 1498 /* -------------------------------------------------------------------------- */ 1499 #undef __FUNCT__ 1500 #define __FUNCT__ "SNESLineSearchNo_VI" 1501 1502 /* 1503 This routine does not actually do a line search but it takes a full newton 1504 step while ensuring that the new iterates remain within the constraints. 1505 1506 */ 1507 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) 1508 { 1509 PetscErrorCode ierr; 1510 SNES_VI *vi = (SNES_VI*)snes->data; 1511 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1512 1513 PetscFunctionBegin; 1514 *flag = PETSC_TRUE; 1515 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1516 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); /* ynorm = || y || */ 1517 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1518 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1519 if (vi->postcheckstep) { 1520 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1521 } 1522 if (changed_y) { 1523 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1524 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1525 } 1526 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1527 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1528 if (!snes->domainerror) { 1529 if (snes->ops->solve != SNESSolveVI_SS) { 1530 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1531 } else { 1532 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); /* gnorm = || g || */ 1533 } 1534 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1535 } 1536 if (vi->lsmonitor) { 1537 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1538 } 1539 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1540 PetscFunctionReturn(0); 1541 } 1542 1543 /* -------------------------------------------------------------------------- */ 1544 #undef __FUNCT__ 1545 #define __FUNCT__ "SNESLineSearchNoNorms_VI" 1546 1547 /* 1548 This routine is a copy of SNESLineSearchNoNorms in snes/impls/ls/ls.c 1549 */ 1550 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) 1551 { 1552 PetscErrorCode ierr; 1553 SNES_VI *vi = (SNES_VI*)snes->data; 1554 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1555 1556 PetscFunctionBegin; 1557 *flag = PETSC_TRUE; 1558 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1559 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1560 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1561 if (vi->postcheckstep) { 1562 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1563 } 1564 if (changed_y) { 1565 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1566 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1567 } 1568 1569 /* don't evaluate function the last time through */ 1570 if (snes->iter < snes->max_its-1) { 1571 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1572 } 1573 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1574 PetscFunctionReturn(0); 1575 } 1576 1577 /* -------------------------------------------------------------------------- */ 1578 #undef __FUNCT__ 1579 #define __FUNCT__ "SNESLineSearchCubic_VI" 1580 /* 1581 This routine implements a cubic line search while doing a projection on the variable bounds 1582 */ 1583 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) 1584 { 1585 PetscReal initslope,lambdaprev,gnormprev,a,b,d,t1,t2,rellength; 1586 PetscReal minlambda,lambda,lambdatemp; 1587 #if defined(PETSC_USE_COMPLEX) 1588 PetscScalar cinitslope; 1589 #endif 1590 PetscErrorCode ierr; 1591 PetscInt count; 1592 SNES_VI *vi = (SNES_VI*)snes->data; 1593 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1594 MPI_Comm comm; 1595 1596 PetscFunctionBegin; 1597 ierr = PetscObjectGetComm((PetscObject)snes,&comm);CHKERRQ(ierr); 1598 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1599 *flag = PETSC_TRUE; 1600 1601 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); 1602 if (*ynorm == 0.0) { 1603 if (vi->lsmonitor) { 1604 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Initial direction and size is 0\n");CHKERRQ(ierr); 1605 } 1606 *gnorm = fnorm; 1607 ierr = VecCopy(x,w);CHKERRQ(ierr); 1608 ierr = VecCopy(f,g);CHKERRQ(ierr); 1609 *flag = PETSC_FALSE; 1610 goto theend1; 1611 } 1612 if (*ynorm > vi->maxstep) { /* Step too big, so scale back */ 1613 if (vi->lsmonitor) { 1614 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Scaling step by %G old ynorm %G\n",vi->maxstep/(*ynorm),*ynorm);CHKERRQ(ierr); 1615 } 1616 ierr = VecScale(y,vi->maxstep/(*ynorm));CHKERRQ(ierr); 1617 *ynorm = vi->maxstep; 1618 } 1619 ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr); 1620 minlambda = vi->minlambda/rellength; 1621 ierr = MatMult(snes->jacobian,y,w);CHKERRQ(ierr); 1622 #if defined(PETSC_USE_COMPLEX) 1623 ierr = VecDot(f,w,&cinitslope);CHKERRQ(ierr); 1624 initslope = PetscRealPart(cinitslope); 1625 #else 1626 ierr = VecDot(f,w,&initslope);CHKERRQ(ierr); 1627 #endif 1628 if (initslope > 0.0) initslope = -initslope; 1629 if (initslope == 0.0) initslope = -1.0; 1630 1631 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1632 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1633 if (snes->nfuncs >= snes->max_funcs) { 1634 ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); 1635 *flag = PETSC_FALSE; 1636 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1637 goto theend1; 1638 } 1639 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1640 if (snes->ops->solve != SNESSolveVI_SS) { 1641 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1642 } else { 1643 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1644 } 1645 if (snes->domainerror) { 1646 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1647 PetscFunctionReturn(0); 1648 } 1649 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1650 ierr = PetscInfo2(snes,"Initial fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1651 if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + vi->alpha*initslope) { /* Sufficient reduction */ 1652 if (vi->lsmonitor) { 1653 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1654 } 1655 goto theend1; 1656 } 1657 1658 /* Fit points with quadratic */ 1659 lambda = 1.0; 1660 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 1661 lambdaprev = lambda; 1662 gnormprev = *gnorm; 1663 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1664 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1665 else lambda = lambdatemp; 1666 1667 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1668 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1669 if (snes->nfuncs >= snes->max_funcs) { 1670 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);CHKERRQ(ierr); 1671 *flag = PETSC_FALSE; 1672 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1673 goto theend1; 1674 } 1675 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1676 if (snes->ops->solve != SNESSolveVI_SS) { 1677 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1678 } else { 1679 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1680 } 1681 if (snes->domainerror) { 1682 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1683 PetscFunctionReturn(0); 1684 } 1685 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1686 if (vi->lsmonitor) { 1687 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: gnorm after quadratic fit %G\n",*gnorm);CHKERRQ(ierr); 1688 } 1689 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* sufficient reduction */ 1690 if (vi->lsmonitor) { 1691 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Quadratically determined step, lambda=%18.16e\n",lambda);CHKERRQ(ierr); 1692 } 1693 goto theend1; 1694 } 1695 1696 /* Fit points with cubic */ 1697 count = 1; 1698 while (PETSC_TRUE) { 1699 if (lambda <= minlambda) { 1700 if (vi->lsmonitor) { 1701 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: unable to find good step length! After %D tries \n",count);CHKERRQ(ierr); 1702 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); 1703 } 1704 *flag = PETSC_FALSE; 1705 break; 1706 } 1707 t1 = .5*((*gnorm)*(*gnorm) - fnorm*fnorm) - lambda*initslope; 1708 t2 = .5*(gnormprev*gnormprev - fnorm*fnorm) - lambdaprev*initslope; 1709 a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 1710 b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 1711 d = b*b - 3*a*initslope; 1712 if (d < 0.0) d = 0.0; 1713 if (a == 0.0) { 1714 lambdatemp = -initslope/(2.0*b); 1715 } else { 1716 lambdatemp = (-b + sqrt(d))/(3.0*a); 1717 } 1718 lambdaprev = lambda; 1719 gnormprev = *gnorm; 1720 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1721 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1722 else lambda = lambdatemp; 1723 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1724 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1725 if (snes->nfuncs >= snes->max_funcs) { 1726 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 1727 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); 1728 *flag = PETSC_FALSE; 1729 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1730 break; 1731 } 1732 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1733 if (snes->ops->solve != SNESSolveVI_SS) { 1734 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1735 } else { 1736 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1737 } 1738 if (snes->domainerror) { 1739 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1740 PetscFunctionReturn(0); 1741 } 1742 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1743 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* is reduction enough? */ 1744 if (vi->lsmonitor) { 1745 ierr = PetscPrintf(comm," Line search: Cubically determined step, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr); 1746 } 1747 break; 1748 } else { 1749 if (vi->lsmonitor) { 1750 ierr = PetscPrintf(comm," Line search: Cubic step no good, shrinking lambda, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr); 1751 } 1752 } 1753 count++; 1754 } 1755 theend1: 1756 /* Optional user-defined check for line search step validity */ 1757 if (vi->postcheckstep && *flag) { 1758 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1759 if (changed_y) { 1760 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1761 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1762 } 1763 if (changed_y || changed_w) { /* recompute the function if the step has changed */ 1764 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1765 if (snes->ops->solve != SNESSolveVI_SS) { 1766 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1767 } else { 1768 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1769 } 1770 if (snes->domainerror) { 1771 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1772 PetscFunctionReturn(0); 1773 } 1774 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1775 ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr); 1776 ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr); 1777 } 1778 } 1779 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1780 PetscFunctionReturn(0); 1781 } 1782 1783 #undef __FUNCT__ 1784 #define __FUNCT__ "SNESLineSearchQuadratic_VI" 1785 /* 1786 This routine does a quadratic line search while keeping the iterates within the variable bounds 1787 */ 1788 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) 1789 { 1790 /* 1791 Note that for line search purposes we work with with the related 1792 minimization problem: 1793 min z(x): R^n -> R, 1794 where z(x) = .5 * fnorm*fnorm,and fnorm = || f ||_2. 1795 */ 1796 PetscReal initslope,minlambda,lambda,lambdatemp,rellength; 1797 #if defined(PETSC_USE_COMPLEX) 1798 PetscScalar cinitslope; 1799 #endif 1800 PetscErrorCode ierr; 1801 PetscInt count; 1802 SNES_VI *vi = (SNES_VI*)snes->data; 1803 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1804 1805 PetscFunctionBegin; 1806 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1807 *flag = PETSC_TRUE; 1808 1809 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); 1810 if (*ynorm == 0.0) { 1811 if (vi->lsmonitor) { 1812 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor,"Line search: Direction and size is 0\n");CHKERRQ(ierr); 1813 } 1814 *gnorm = fnorm; 1815 ierr = VecCopy(x,w);CHKERRQ(ierr); 1816 ierr = VecCopy(f,g);CHKERRQ(ierr); 1817 *flag = PETSC_FALSE; 1818 goto theend2; 1819 } 1820 if (*ynorm > vi->maxstep) { /* Step too big, so scale back */ 1821 ierr = VecScale(y,vi->maxstep/(*ynorm));CHKERRQ(ierr); 1822 *ynorm = vi->maxstep; 1823 } 1824 ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr); 1825 minlambda = vi->minlambda/rellength; 1826 ierr = MatMult(snes->jacobian,y,w);CHKERRQ(ierr); 1827 #if defined(PETSC_USE_COMPLEX) 1828 ierr = VecDot(f,w,&cinitslope);CHKERRQ(ierr); 1829 initslope = PetscRealPart(cinitslope); 1830 #else 1831 ierr = VecDot(f,w,&initslope);CHKERRQ(ierr); 1832 #endif 1833 if (initslope > 0.0) initslope = -initslope; 1834 if (initslope == 0.0) initslope = -1.0; 1835 ierr = PetscInfo1(snes,"Initslope %G \n",initslope);CHKERRQ(ierr); 1836 1837 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1838 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1839 if (snes->nfuncs >= snes->max_funcs) { 1840 ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); 1841 *flag = PETSC_FALSE; 1842 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1843 goto theend2; 1844 } 1845 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1846 if (snes->ops->solve != SNESSolveVI_SS) { 1847 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1848 } else { 1849 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1850 } 1851 if (snes->domainerror) { 1852 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1853 PetscFunctionReturn(0); 1854 } 1855 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1856 if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + vi->alpha*initslope) { /* Sufficient reduction */ 1857 if (vi->lsmonitor) { 1858 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1859 } 1860 goto theend2; 1861 } 1862 1863 /* Fit points with quadratic */ 1864 lambda = 1.0; 1865 count = 1; 1866 while (PETSC_TRUE) { 1867 if (lambda <= minlambda) { /* bad luck; use full step */ 1868 if (vi->lsmonitor) { 1869 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor,"Line search: Unable to find good step length! %D \n",count);CHKERRQ(ierr); 1870 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); 1871 } 1872 ierr = VecCopy(x,w);CHKERRQ(ierr); 1873 *flag = PETSC_FALSE; 1874 break; 1875 } 1876 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 1877 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1878 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1879 else lambda = lambdatemp; 1880 1881 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1882 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1883 if (snes->nfuncs >= snes->max_funcs) { 1884 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 1885 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); 1886 *flag = PETSC_FALSE; 1887 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1888 break; 1889 } 1890 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1891 if (snes->domainerror) { 1892 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1893 PetscFunctionReturn(0); 1894 } 1895 if (snes->ops->solve != SNESSolveVI_SS) { 1896 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1897 } else { 1898 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1899 } 1900 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1901 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* sufficient reduction */ 1902 if (vi->lsmonitor) { 1903 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line Search: Quadratically determined step, lambda=%G\n",lambda);CHKERRQ(ierr); 1904 } 1905 break; 1906 } 1907 count++; 1908 } 1909 theend2: 1910 /* Optional user-defined check for line search step validity */ 1911 if (vi->postcheckstep) { 1912 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1913 if (changed_y) { 1914 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1915 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1916 } 1917 if (changed_y || changed_w) { /* recompute the function if the step has changed */ 1918 ierr = SNESComputeFunction(snes,w,g); 1919 if (snes->domainerror) { 1920 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1921 PetscFunctionReturn(0); 1922 } 1923 if (snes->ops->solve != SNESSolveVI_SS) { 1924 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1925 } else { 1926 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1927 } 1928 1929 ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr); 1930 ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr); 1931 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1932 } 1933 } 1934 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1935 PetscFunctionReturn(0); 1936 } 1937 1938 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*/ 1939 /* -------------------------------------------------------------------------- */ 1940 EXTERN_C_BEGIN 1941 #undef __FUNCT__ 1942 #define __FUNCT__ "SNESLineSearchSet_VI" 1943 PetscErrorCode SNESLineSearchSet_VI(SNES snes,FCN2 func,void *lsctx) 1944 { 1945 PetscFunctionBegin; 1946 ((SNES_VI *)(snes->data))->LineSearch = func; 1947 ((SNES_VI *)(snes->data))->lsP = lsctx; 1948 PetscFunctionReturn(0); 1949 } 1950 EXTERN_C_END 1951 1952 /* -------------------------------------------------------------------------- */ 1953 EXTERN_C_BEGIN 1954 #undef __FUNCT__ 1955 #define __FUNCT__ "SNESLineSearchSetMonitor_VI" 1956 PetscErrorCode SNESLineSearchSetMonitor_VI(SNES snes,PetscBool flg) 1957 { 1958 SNES_VI *vi = (SNES_VI*)snes->data; 1959 PetscErrorCode ierr; 1960 1961 PetscFunctionBegin; 1962 if (flg && !vi->lsmonitor) { 1963 ierr = PetscViewerASCIIMonitorCreate(((PetscObject)snes)->comm,"stdout",((PetscObject)snes)->tablevel,&vi->lsmonitor);CHKERRQ(ierr); 1964 } else if (!flg && vi->lsmonitor) { 1965 ierr = PetscViewerASCIIMonitorDestroy(vi->lsmonitor);CHKERRQ(ierr); 1966 } 1967 PetscFunctionReturn(0); 1968 } 1969 EXTERN_C_END 1970 1971 /* 1972 SNESView_VI - Prints info from the SNESVI data structure. 1973 1974 Input Parameters: 1975 . SNES - the SNES context 1976 . viewer - visualization context 1977 1978 Application Interface Routine: SNESView() 1979 */ 1980 #undef __FUNCT__ 1981 #define __FUNCT__ "SNESView_VI" 1982 static PetscErrorCode SNESView_VI(SNES snes,PetscViewer viewer) 1983 { 1984 SNES_VI *vi = (SNES_VI *)snes->data; 1985 const char *cstr,*tstr; 1986 PetscErrorCode ierr; 1987 PetscBool iascii; 1988 1989 PetscFunctionBegin; 1990 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1991 if (iascii) { 1992 if (vi->LineSearch == SNESLineSearchNo_VI) cstr = "SNESLineSearchNo"; 1993 else if (vi->LineSearch == SNESLineSearchQuadratic_VI) cstr = "SNESLineSearchQuadratic"; 1994 else if (vi->LineSearch == SNESLineSearchCubic_VI) cstr = "SNESLineSearchCubic"; 1995 else cstr = "unknown"; 1996 if (snes->ops->solve == SNESSolveVI_SS) tstr = "Semismooth"; 1997 else if (snes->ops->solve == SNESSolveVI_RS) tstr = "Reduced Space"; 1998 else if (snes->ops->solve == SNESSolveVI_RS2) tstr = "Augmented Space"; 1999 else tstr = "unknown"; 2000 ierr = PetscViewerASCIIPrintf(viewer," VI algorithm: %s\n",tstr);CHKERRQ(ierr); 2001 ierr = PetscViewerASCIIPrintf(viewer," line search variant: %s\n",cstr);CHKERRQ(ierr); 2002 ierr = PetscViewerASCIIPrintf(viewer," alpha=%G, maxstep=%G, minlambda=%G\n",vi->alpha,vi->maxstep,vi->minlambda);CHKERRQ(ierr); 2003 } else { 2004 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for SNES EQ VI",((PetscObject)viewer)->type_name); 2005 } 2006 PetscFunctionReturn(0); 2007 } 2008 2009 #undef __FUNCT__ 2010 #define __FUNCT__ "SNESVISetVariableBounds" 2011 /*@ 2012 SNESVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu. 2013 2014 Input Parameters: 2015 . snes - the SNES context. 2016 . xl - lower bound. 2017 . xu - upper bound. 2018 2019 Notes: 2020 If this routine is not called then the lower and upper bounds are set to 2021 PETSC_VI_INF and PETSC_VI_NINF respectively during SNESSetUp(). 2022 2023 @*/ 2024 PetscErrorCode SNESVISetVariableBounds(SNES snes, Vec xl, Vec xu) 2025 { 2026 SNES_VI *vi = (SNES_VI*)snes->data; 2027 2028 PetscFunctionBegin; 2029 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2030 PetscValidHeaderSpecific(xl,VEC_CLASSID,2); 2031 PetscValidHeaderSpecific(xu,VEC_CLASSID,3); 2032 if (!snes->vec_func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call SNESSetFunction() first"); 2033 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); 2034 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); 2035 2036 vi->usersetxbounds = PETSC_TRUE; 2037 vi->xl = xl; 2038 vi->xu = xu; 2039 PetscFunctionReturn(0); 2040 } 2041 2042 /* -------------------------------------------------------------------------- */ 2043 /* 2044 SNESSetFromOptions_VI - Sets various parameters for the SNESVI method. 2045 2046 Input Parameter: 2047 . snes - the SNES context 2048 2049 Application Interface Routine: SNESSetFromOptions() 2050 */ 2051 #undef __FUNCT__ 2052 #define __FUNCT__ "SNESSetFromOptions_VI" 2053 static PetscErrorCode SNESSetFromOptions_VI(SNES snes) 2054 { 2055 SNES_VI *vi = (SNES_VI *)snes->data; 2056 const char *lses[] = {"basic","basicnonorms","quadratic","cubic"}; 2057 const char *vies[] = {"ss","rs","rs2"}; 2058 PetscErrorCode ierr; 2059 PetscInt indx; 2060 PetscBool flg,set,flg2; 2061 2062 PetscFunctionBegin; 2063 ierr = PetscOptionsHead("SNES semismooth method options");CHKERRQ(ierr); 2064 ierr = PetscOptionsBool("-snes_vi_monitor","Monitor all non-active variables","None",PETSC_FALSE,&flg,0);CHKERRQ(ierr); 2065 if (flg) { 2066 ierr = SNESMonitorSet(snes,SNESMonitorVI,0,0);CHKERRQ(ierr); 2067 } 2068 ierr = PetscOptionsReal("-snes_ls_alpha","Function norm must decrease by","None",vi->alpha,&vi->alpha,0);CHKERRQ(ierr); 2069 ierr = PetscOptionsReal("-snes_ls_maxstep","Step must be less than","None",vi->maxstep,&vi->maxstep,0);CHKERRQ(ierr); 2070 ierr = PetscOptionsReal("-snes_ls_minlambda","Minimum lambda allowed","None",vi->minlambda,&vi->minlambda,0);CHKERRQ(ierr); 2071 ierr = PetscOptionsReal("-snes_vi_const_tol","constraint tolerance","None",vi->const_tol,&vi->const_tol,0);CHKERRQ(ierr); 2072 ierr = PetscOptionsBool("-snes_ls_monitor","Print progress of line searches","SNESLineSearchSetMonitor",vi->lsmonitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);CHKERRQ(ierr); 2073 if (set) {ierr = SNESLineSearchSetMonitor(snes,flg);CHKERRQ(ierr);} 2074 ierr = PetscOptionsEList("-snes_vi_type","Semismooth algorithm used","",vies,3,"ss",&indx,&flg2);CHKERRQ(ierr); 2075 if (flg2) { 2076 switch (indx) { 2077 case 0: 2078 snes->ops->solve = SNESSolveVI_SS; 2079 break; 2080 case 1: 2081 snes->ops->solve = SNESSolveVI_RS; 2082 break; 2083 case 2: 2084 snes->ops->solve = SNESSolveVI_RS2; 2085 } 2086 } 2087 ierr = PetscOptionsEList("-snes_ls","Line search used","SNESLineSearchSet",lses,4,"basic",&indx,&flg);CHKERRQ(ierr); 2088 if (flg) { 2089 switch (indx) { 2090 case 0: 2091 ierr = SNESLineSearchSet(snes,SNESLineSearchNo_VI,PETSC_NULL);CHKERRQ(ierr); 2092 break; 2093 case 1: 2094 ierr = SNESLineSearchSet(snes,SNESLineSearchNoNorms_VI,PETSC_NULL);CHKERRQ(ierr); 2095 break; 2096 case 2: 2097 ierr = SNESLineSearchSet(snes,SNESLineSearchQuadratic_VI,PETSC_NULL);CHKERRQ(ierr); 2098 break; 2099 case 3: 2100 ierr = SNESLineSearchSet(snes,SNESLineSearchCubic_VI,PETSC_NULL);CHKERRQ(ierr); 2101 break; 2102 } 2103 } 2104 ierr = PetscOptionsTail();CHKERRQ(ierr); 2105 PetscFunctionReturn(0); 2106 } 2107 /* -------------------------------------------------------------------------- */ 2108 /*MC 2109 SNESVI - Semismooth newton method based nonlinear solver that uses a line search 2110 2111 Options Database: 2112 + -snes_ls [cubic,quadratic,basic,basicnonorms] - Selects line search 2113 . -snes_ls_alpha <alpha> - Sets alpha 2114 . -snes_ls_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) 2115 . -snes_ls_minlambda <minlambda> - Sets the minimum lambda the line search will use minlambda / max_i ( y[i]/x[i] ) 2116 - -snes_ls_monitor - print information about progress of line searches 2117 2118 2119 Level: beginner 2120 2121 .seealso: SNESCreate(), SNES, SNESSetType(), SNESTR, SNESLineSearchSet(), 2122 SNESLineSearchSetPostCheck(), SNESLineSearchNo(), SNESLineSearchCubic(), SNESLineSearchQuadratic(), 2123 SNESLineSearchSet(), SNESLineSearchNoNorms(), SNESLineSearchSetPreCheck(), SNESLineSearchSetParams(), SNESLineSearchGetParams() 2124 2125 M*/ 2126 EXTERN_C_BEGIN 2127 #undef __FUNCT__ 2128 #define __FUNCT__ "SNESCreate_VI" 2129 PetscErrorCode SNESCreate_VI(SNES snes) 2130 { 2131 PetscErrorCode ierr; 2132 SNES_VI *vi; 2133 2134 PetscFunctionBegin; 2135 snes->ops->setup = SNESSetUp_VI; 2136 snes->ops->solve = SNESSolveVI_SS; 2137 snes->ops->destroy = SNESDestroy_VI; 2138 snes->ops->setfromoptions = SNESSetFromOptions_VI; 2139 snes->ops->view = SNESView_VI; 2140 snes->ops->converged = SNESDefaultConverged_VI; 2141 2142 ierr = PetscNewLog(snes,SNES_VI,&vi);CHKERRQ(ierr); 2143 snes->data = (void*)vi; 2144 vi->alpha = 1.e-4; 2145 vi->maxstep = 1.e8; 2146 vi->minlambda = 1.e-12; 2147 vi->LineSearch = SNESLineSearchNo_VI; 2148 vi->lsP = PETSC_NULL; 2149 vi->postcheckstep = PETSC_NULL; 2150 vi->postcheck = PETSC_NULL; 2151 vi->precheckstep = PETSC_NULL; 2152 vi->precheck = PETSC_NULL; 2153 vi->const_tol = 2.2204460492503131e-16; 2154 vi->checkredundancy = PETSC_NULL; 2155 2156 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetMonitor_C","SNESLineSearchSetMonitor_VI",SNESLineSearchSetMonitor_VI);CHKERRQ(ierr); 2157 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSet_C","SNESLineSearchSet_VI",SNESLineSearchSet_VI);CHKERRQ(ierr); 2158 2159 PetscFunctionReturn(0); 2160 } 2161 EXTERN_C_END 2162