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 (((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + 1.e-8 || (PetscRealPart(f[i]) < 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - 1.e-8) || PetscRealPart(f[i]) > 0.0))) rnorm += PetscRealPart(PetscConj(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 ((PetscRealPart(l[i]) <= SNES_VI_NINF) && (PetscRealPart(u[i]) >= SNES_VI_INF)) { /* no constraints on variable */ 230 phi_arr[i] = f_arr[i]; 231 } else if (PetscRealPart(l[i]) <= SNES_VI_NINF) { /* upper bound on variable only */ 232 phi_arr[i] = -Phi(u[i] - x_arr[i],-f_arr[i]); 233 } else if (PetscRealPart(u[i]) >= SNES_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 ((PetscRealPart(l[i]) <= SNES_VI_NINF) && (PetscRealPart(u[i]) >= SNES_VI_INF)) {/* no constraints on variable */ 275 da[i] = 0; 276 db[i] = 1; 277 } else if (PetscRealPart(l[i]) <= SNES_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 (PetscRealPart(u[i]) >= SNES_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 (PetscRealPart(x[i]) < PetscRealPart(xl[i])) x[i] = xl[i]; 394 else if (PetscRealPart(x[i]) > PetscRealPart(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 ierr = SNESMonitor(snes,0,vi->phinorm);CHKERRQ(ierr); 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 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 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 (!((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + 1.e-8 || (PetscRealPart(f[i]) < 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - 1.e-8) || PetscRealPart(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 (!((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + 1.e-8 || (PetscRealPart(f[i]) < 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - 1.e-8) || PetscRealPart(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 ierr = SNESGetKSP(snes,&snesksp);CHKERRQ(ierr); 697 698 /* 699 ierr = KSPReset(snesksp);CHKERRQ(ierr); 700 */ 701 702 ierr = KSPCreate(((PetscObject)snes)->comm,&kspnew);CHKERRQ(ierr); 703 kspnew->pc_side = snesksp->pc_side; 704 kspnew->rtol = snesksp->rtol; 705 kspnew->abstol = snesksp->abstol; 706 kspnew->max_it = snesksp->max_it; 707 ierr = KSPSetType(kspnew,((PetscObject)snesksp)->type_name);CHKERRQ(ierr); 708 ierr = KSPGetPC(kspnew,&pcnew);CHKERRQ(ierr); 709 ierr = PCSetType(kspnew->pc,((PetscObject)snesksp->pc)->type_name);CHKERRQ(ierr); 710 ierr = PCSetOperators(kspnew->pc,Amat,Pmat,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 711 ierr = PCFactorGetMatSolverPackage(snesksp->pc,&stype);CHKERRQ(ierr); 712 ierr = PCFactorSetMatSolverPackage(kspnew->pc,stype);CHKERRQ(ierr); 713 ierr = KSPDestroy(&snesksp);CHKERRQ(ierr); 714 snes->ksp = kspnew; 715 ierr = PetscLogObjectParent(snes,kspnew);CHKERRQ(ierr); 716 ierr = KSPSetFromOptions(kspnew);CHKERRQ(ierr); 717 PetscFunctionReturn(0); 718 } 719 720 721 #undef __FUNCT__ 722 #define __FUNCT__ "SNESVIComputeInactiveSetFnorm" 723 PetscErrorCode SNESVIComputeInactiveSetFnorm(SNES snes,Vec F,Vec X,PetscReal *fnorm) 724 { 725 PetscErrorCode ierr; 726 SNES_VI *vi = (SNES_VI*)snes->data; 727 const PetscScalar *x,*xl,*xu,*f; 728 PetscInt i,n; 729 PetscReal rnorm; 730 731 PetscFunctionBegin; 732 ierr = VecGetLocalSize(X,&n);CHKERRQ(ierr); 733 ierr = VecGetArrayRead(vi->xl,&xl);CHKERRQ(ierr); 734 ierr = VecGetArrayRead(vi->xu,&xu);CHKERRQ(ierr); 735 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 736 ierr = VecGetArrayRead(F,&f);CHKERRQ(ierr); 737 rnorm = 0.0; 738 for (i=0; i<n; i++) { 739 if (((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + 1.e-8 || (PetscRealPart(f[i]) < 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - 1.e-8) || PetscRealPart(f[i]) > 0.0))) rnorm += PetscRealPart(PetscConj(f[i])*f[i]); 740 } 741 ierr = VecRestoreArrayRead(F,&f);CHKERRQ(ierr); 742 ierr = VecRestoreArrayRead(vi->xl,&xl);CHKERRQ(ierr); 743 ierr = VecRestoreArrayRead(vi->xu,&xu);CHKERRQ(ierr); 744 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 745 ierr = MPI_Allreduce(&rnorm,fnorm,1,MPIU_REAL,MPIU_SUM,((PetscObject)snes)->comm);CHKERRQ(ierr); 746 *fnorm = sqrt(*fnorm); 747 PetscFunctionReturn(0); 748 } 749 750 /* Variational Inequality solver using reduce space method. No semismooth algorithm is 751 implemented in this algorithm. It basically identifies the active constraints and does 752 a linear solve on the other variables (those not associated with the active constraints). */ 753 #undef __FUNCT__ 754 #define __FUNCT__ "SNESSolveVI_RS" 755 PetscErrorCode SNESSolveVI_RS(SNES snes) 756 { 757 SNES_VI *vi = (SNES_VI*)snes->data; 758 PetscErrorCode ierr; 759 PetscInt maxits,i,lits; 760 PetscBool lssucceed; 761 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 762 PetscReal fnorm,gnorm,xnorm=0,ynorm; 763 Vec Y,X,F,G,W; 764 KSPConvergedReason kspreason; 765 766 PetscFunctionBegin; 767 snes->numFailures = 0; 768 snes->numLinearSolveFailures = 0; 769 snes->reason = SNES_CONVERGED_ITERATING; 770 771 maxits = snes->max_its; /* maximum number of iterations */ 772 X = snes->vec_sol; /* solution vector */ 773 F = snes->vec_func; /* residual vector */ 774 Y = snes->work[0]; /* work vectors */ 775 G = snes->work[1]; 776 W = snes->work[2]; 777 778 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 779 snes->iter = 0; 780 snes->norm = 0.0; 781 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 782 783 ierr = SNESVIProjectOntoBounds(snes,X);CHKERRQ(ierr); 784 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 785 if (snes->domainerror) { 786 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 787 PetscFunctionReturn(0); 788 } 789 ierr = SNESVIComputeInactiveSetFnorm(snes,F,X,&fnorm);CHKERRQ(ierr); 790 ierr = VecNormBegin(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm <- ||x|| */ 791 ierr = VecNormEnd(X,NORM_2,&xnorm);CHKERRQ(ierr); 792 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 793 794 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 795 snes->norm = fnorm; 796 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 797 SNESLogConvHistory(snes,fnorm,0); 798 ierr = SNESMonitor(snes,0,fnorm);CHKERRQ(ierr); 799 800 /* set parameter for default relative tolerance convergence test */ 801 snes->ttol = fnorm*snes->rtol; 802 /* test convergence */ 803 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 804 if (snes->reason) PetscFunctionReturn(0); 805 806 for (i=0; i<maxits; i++) { 807 808 IS IS_act,IS_inact; /* _act -> active set _inact -> inactive set */ 809 VecScatter scat_act,scat_inact; 810 PetscInt nis_act,nis_inact; 811 Vec Y_act,Y_inact,F_inact; 812 Mat jac_inact_inact,prejac_inact_inact; 813 IS keptrows; 814 PetscBool isequal; 815 816 /* Call general purpose update function */ 817 if (snes->ops->update) { 818 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 819 } 820 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 821 822 /* Create active and inactive index sets */ 823 ierr = SNESVICreateIndexSets_RS(snes,X,F,&IS_act,&IS_inact);CHKERRQ(ierr); 824 825 /* Create inactive set submatrix */ 826 ierr = MatGetSubMatrix(snes->jacobian,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);CHKERRQ(ierr); 827 ierr = MatFindNonzeroRows(jac_inact_inact,&keptrows);CHKERRQ(ierr); 828 if (keptrows) { 829 PetscInt cnt,*nrows,k; 830 const PetscInt *krows,*inact; 831 PetscInt rstart=jac_inact_inact->rmap->rstart; 832 833 ierr = MatDestroy(&jac_inact_inact);CHKERRQ(ierr); 834 ierr = ISDestroy(&IS_act);CHKERRQ(ierr); 835 836 ierr = ISGetLocalSize(keptrows,&cnt);CHKERRQ(ierr); 837 ierr = ISGetIndices(keptrows,&krows);CHKERRQ(ierr); 838 ierr = ISGetIndices(IS_inact,&inact);CHKERRQ(ierr); 839 ierr = PetscMalloc(cnt*sizeof(PetscInt),&nrows);CHKERRQ(ierr); 840 for (k=0; k<cnt; k++) { 841 nrows[k] = inact[krows[k]-rstart]; 842 } 843 ierr = ISRestoreIndices(keptrows,&krows);CHKERRQ(ierr); 844 ierr = ISRestoreIndices(IS_inact,&inact);CHKERRQ(ierr); 845 ierr = ISDestroy(&keptrows);CHKERRQ(ierr); 846 ierr = ISDestroy(&IS_inact);CHKERRQ(ierr); 847 848 ierr = ISCreateGeneral(PETSC_COMM_WORLD,cnt,nrows,PETSC_OWN_POINTER,&IS_inact);CHKERRQ(ierr); 849 ierr = ISComplement(IS_inact,F->map->rstart,F->map->rend,&IS_act);CHKERRQ(ierr); 850 ierr = MatGetSubMatrix(snes->jacobian,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);CHKERRQ(ierr); 851 } 852 853 /* Get sizes of active and inactive sets */ 854 ierr = ISGetLocalSize(IS_act,&nis_act);CHKERRQ(ierr); 855 ierr = ISGetLocalSize(IS_inact,&nis_inact);CHKERRQ(ierr); 856 857 /* Create active and inactive set vectors */ 858 ierr = SNESVICreateSubVectors(snes,nis_inact,&F_inact);CHKERRQ(ierr); 859 ierr = SNESVICreateSubVectors(snes,nis_act,&Y_act);CHKERRQ(ierr); 860 ierr = SNESVICreateSubVectors(snes,nis_inact,&Y_inact);CHKERRQ(ierr); 861 862 /* Create scatter contexts */ 863 ierr = VecScatterCreate(Y,IS_act,Y_act,PETSC_NULL,&scat_act);CHKERRQ(ierr); 864 ierr = VecScatterCreate(Y,IS_inact,Y_inact,PETSC_NULL,&scat_inact);CHKERRQ(ierr); 865 866 /* Do a vec scatter to active and inactive set vectors */ 867 ierr = VecScatterBegin(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 868 ierr = VecScatterEnd(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 869 870 ierr = VecScatterBegin(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 871 ierr = VecScatterEnd(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 872 873 ierr = VecScatterBegin(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 874 ierr = VecScatterEnd(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 875 876 /* Active set direction = 0 */ 877 ierr = VecSet(Y_act,0);CHKERRQ(ierr); 878 if (snes->jacobian != snes->jacobian_pre) { 879 ierr = MatGetSubMatrix(snes->jacobian_pre,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&prejac_inact_inact);CHKERRQ(ierr); 880 } else prejac_inact_inact = jac_inact_inact; 881 882 ierr = ISEqual(vi->IS_inact_prev,IS_inact,&isequal);CHKERRQ(ierr); 883 if (!isequal) { 884 ierr = SNESVIResetPCandKSP(snes,jac_inact_inact,prejac_inact_inact);CHKERRQ(ierr); 885 } 886 887 /* ierr = ISView(IS_inact,0);CHKERRQ(ierr); */ 888 /* ierr = ISView(IS_act,0);CHKERRQ(ierr);*/ 889 /* ierr = MatView(snes->jacobian_pre,0); */ 890 891 ierr = KSPSetOperators(snes->ksp,jac_inact_inact,prejac_inact_inact,flg);CHKERRQ(ierr); 892 ierr = KSPSetUp(snes->ksp);CHKERRQ(ierr); 893 { 894 PC pc; 895 PetscBool flg; 896 ierr = KSPGetPC(snes->ksp,&pc);CHKERRQ(ierr); 897 ierr = PetscTypeCompare((PetscObject)pc,PCFIELDSPLIT,&flg);CHKERRQ(ierr); 898 if (flg) { 899 KSP *subksps; 900 ierr = PCFieldSplitGetSubKSP(pc,PETSC_NULL,&subksps);CHKERRQ(ierr); 901 ierr = KSPGetPC(subksps[0],&pc);CHKERRQ(ierr); 902 ierr = PetscFree(subksps);CHKERRQ(ierr); 903 ierr = PetscTypeCompare((PetscObject)pc,PCBJACOBI,&flg);CHKERRQ(ierr); 904 if (flg) { 905 PetscInt n,N = 101*101,j,cnts[3] = {0,0,0}; 906 const PetscInt *ii; 907 908 ierr = ISGetSize(IS_inact,&n);CHKERRQ(ierr); 909 ierr = ISGetIndices(IS_inact,&ii);CHKERRQ(ierr); 910 for (j=0; j<n; j++) { 911 if (ii[j] < N) cnts[0]++; 912 else if (ii[j] < 2*N) cnts[1]++; 913 else if (ii[j] < 3*N) cnts[2]++; 914 } 915 ierr = ISRestoreIndices(IS_inact,&ii);CHKERRQ(ierr); 916 917 ierr = PCBJacobiSetTotalBlocks(pc,3,cnts);CHKERRQ(ierr); 918 } 919 } 920 } 921 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 if (!isequal) { 944 ierr = ISDestroy(&vi->IS_inact_prev);CHKERRQ(ierr); 945 ierr = ISDuplicate(IS_inact,&vi->IS_inact_prev);CHKERRQ(ierr); 946 } 947 ierr = ISDestroy(&IS_inact);CHKERRQ(ierr); 948 ierr = MatDestroy(&jac_inact_inact);CHKERRQ(ierr); 949 if (snes->jacobian != snes->jacobian_pre) { 950 ierr = MatDestroy(&prejac_inact_inact);CHKERRQ(ierr); 951 } 952 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 953 snes->linear_its += lits; 954 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 955 /* 956 if (vi->precheckstep) { 957 PetscBool changed_y = PETSC_FALSE; 958 ierr = (*vi->precheckstep)(snes,X,Y,vi->precheck,&changed_y);CHKERRQ(ierr); 959 } 960 961 if (PetscLogPrintInfo){ 962 ierr = SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr); 963 } 964 */ 965 /* Compute a (scaled) negative update in the line search routine: 966 Y <- X - lambda*Y 967 and evaluate G = function(Y) (depends on the line search). 968 */ 969 ierr = VecCopy(Y,snes->vec_sol_update);CHKERRQ(ierr); 970 ynorm = 1; gnorm = fnorm; 971 ierr = (*vi->LineSearch)(snes,vi->lsP,X,F,G,Y,W,fnorm,xnorm,&ynorm,&gnorm,&lssucceed);CHKERRQ(ierr); 972 ierr = PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",fnorm,gnorm,ynorm,(int)lssucceed);CHKERRQ(ierr); 973 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 974 if (snes->domainerror) { 975 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 976 PetscFunctionReturn(0); 977 } 978 if (!lssucceed) { 979 if (++snes->numFailures >= snes->maxFailures) { 980 PetscBool ismin; 981 snes->reason = SNES_DIVERGED_LINE_SEARCH; 982 ierr = SNESVICheckLocalMin_Private(snes,snes->jacobian,G,W,gnorm,&ismin);CHKERRQ(ierr); 983 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 984 break; 985 } 986 } 987 /* Update function and solution vectors */ 988 fnorm = gnorm; 989 ierr = VecCopy(G,F);CHKERRQ(ierr); 990 ierr = VecCopy(W,X);CHKERRQ(ierr); 991 /* Monitor convergence */ 992 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 993 snes->iter = i+1; 994 snes->norm = fnorm; 995 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 996 SNESLogConvHistory(snes,snes->norm,lits); 997 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 998 /* Test for convergence, xnorm = || X || */ 999 if (snes->ops->converged != SNESSkipConverged) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 1000 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 1001 if (snes->reason) break; 1002 } 1003 if (i == maxits) { 1004 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 1005 if(!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 1006 } 1007 PetscFunctionReturn(0); 1008 } 1009 1010 #undef __FUNCT__ 1011 #define __FUNCT__ "SNESVISetRedundancyCheck" 1012 PetscErrorCode SNESVISetRedundancyCheck(SNES snes,PetscErrorCode (*func)(SNES,IS,IS*,void*),void *ctx) 1013 { 1014 SNES_VI *vi = (SNES_VI*)snes->data; 1015 1016 PetscFunctionBegin; 1017 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1018 vi->checkredundancy = func; 1019 vi->ctxP = ctx; 1020 PetscFunctionReturn(0); 1021 } 1022 1023 #if defined(PETSC_HAVE_MATLAB_ENGINE) 1024 #include <engine.h> 1025 #include <mex.h> 1026 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 1027 1028 #undef __FUNCT__ 1029 #define __FUNCT__ "SNESVIRedundancyCheck_Matlab" 1030 PetscErrorCode SNESVIRedundancyCheck_Matlab(SNES snes,IS is_act,IS* is_redact,void* ctx) 1031 { 1032 PetscErrorCode ierr; 1033 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 1034 int nlhs = 1, nrhs = 5; 1035 mxArray *plhs[1], *prhs[5]; 1036 long long int l1 = 0, l2 = 0, ls = 0; 1037 PetscInt *indices=PETSC_NULL; 1038 1039 PetscFunctionBegin; 1040 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1041 PetscValidHeaderSpecific(is_act,IS_CLASSID,2); 1042 PetscValidPointer(is_redact,3); 1043 PetscCheckSameComm(snes,1,is_act,2); 1044 1045 /* Create IS for reduced active set of size 0, its size and indices will 1046 bet set by the Matlab function */ 1047 ierr = ISCreateGeneral(((PetscObject)snes)->comm,0,indices,PETSC_OWN_POINTER,is_redact);CHKERRQ(ierr); 1048 /* call Matlab function in ctx */ 1049 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 1050 ierr = PetscMemcpy(&l1,&is_act,sizeof(is_act));CHKERRQ(ierr); 1051 ierr = PetscMemcpy(&l2,is_redact,sizeof(is_act));CHKERRQ(ierr); 1052 prhs[0] = mxCreateDoubleScalar((double)ls); 1053 prhs[1] = mxCreateDoubleScalar((double)l1); 1054 prhs[2] = mxCreateDoubleScalar((double)l2); 1055 prhs[3] = mxCreateString(sctx->funcname); 1056 prhs[4] = sctx->ctx; 1057 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESVIRedundancyCheckInternal");CHKERRQ(ierr); 1058 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 1059 mxDestroyArray(prhs[0]); 1060 mxDestroyArray(prhs[1]); 1061 mxDestroyArray(prhs[2]); 1062 mxDestroyArray(prhs[3]); 1063 mxDestroyArray(plhs[0]); 1064 PetscFunctionReturn(0); 1065 } 1066 1067 #undef __FUNCT__ 1068 #define __FUNCT__ "SNESVISetRedundancyCheckMatlab" 1069 PetscErrorCode SNESVISetRedundancyCheckMatlab(SNES snes,const char* func,mxArray* ctx) 1070 { 1071 PetscErrorCode ierr; 1072 SNESMatlabContext *sctx; 1073 1074 PetscFunctionBegin; 1075 /* currently sctx is memory bleed */ 1076 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 1077 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 1078 sctx->ctx = mxDuplicateArray(ctx); 1079 ierr = SNESVISetRedundancyCheck(snes,SNESVIRedundancyCheck_Matlab,sctx);CHKERRQ(ierr); 1080 PetscFunctionReturn(0); 1081 } 1082 1083 #endif 1084 1085 1086 /* Variational Inequality solver using augmented space method. It does the opposite of the 1087 reduced space method i.e. it identifies the active set variables and instead of discarding 1088 them it augments the original system by introducing additional equality 1089 constraint equations for active set variables. The user can optionally provide an IS for 1090 a subset of 'redundant' active set variables which will treated as free variables. 1091 Specific implementation for Allen-Cahn problem 1092 */ 1093 #undef __FUNCT__ 1094 #define __FUNCT__ "SNESSolveVI_RSAUG" 1095 PetscErrorCode SNESSolveVI_RSAUG(SNES snes) 1096 { 1097 SNES_VI *vi = (SNES_VI*)snes->data; 1098 PetscErrorCode ierr; 1099 PetscInt maxits,i,lits; 1100 PetscBool lssucceed; 1101 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 1102 PetscReal fnorm,gnorm,xnorm=0,ynorm; 1103 Vec Y,X,F,G,W; 1104 KSPConvergedReason kspreason; 1105 1106 PetscFunctionBegin; 1107 snes->numFailures = 0; 1108 snes->numLinearSolveFailures = 0; 1109 snes->reason = SNES_CONVERGED_ITERATING; 1110 1111 maxits = snes->max_its; /* maximum number of iterations */ 1112 X = snes->vec_sol; /* solution vector */ 1113 F = snes->vec_func; /* residual vector */ 1114 Y = snes->work[0]; /* work vectors */ 1115 G = snes->work[1]; 1116 W = snes->work[2]; 1117 1118 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1119 snes->iter = 0; 1120 snes->norm = 0.0; 1121 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1122 1123 ierr = SNESVIProjectOntoBounds(snes,X);CHKERRQ(ierr); 1124 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 1125 if (snes->domainerror) { 1126 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 1127 PetscFunctionReturn(0); 1128 } 1129 ierr = SNESVIComputeInactiveSetFnorm(snes,F,X,&fnorm);CHKERRQ(ierr); 1130 ierr = VecNormBegin(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm <- ||x|| */ 1131 ierr = VecNormEnd(X,NORM_2,&xnorm);CHKERRQ(ierr); 1132 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1133 1134 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1135 snes->norm = fnorm; 1136 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1137 SNESLogConvHistory(snes,fnorm,0); 1138 ierr = SNESMonitor(snes,0,fnorm);CHKERRQ(ierr); 1139 1140 /* set parameter for default relative tolerance convergence test */ 1141 snes->ttol = fnorm*snes->rtol; 1142 /* test convergence */ 1143 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 1144 if (snes->reason) PetscFunctionReturn(0); 1145 1146 for (i=0; i<maxits; i++) { 1147 IS IS_act,IS_inact; /* _act -> active set _inact -> inactive set */ 1148 IS IS_redact; /* redundant active set */ 1149 Mat J_aug,Jpre_aug; 1150 Vec F_aug,Y_aug; 1151 PetscInt nis_redact,nis_act; 1152 const PetscInt *idx_redact,*idx_act; 1153 PetscInt k; 1154 PetscInt *idx_actkept=PETSC_NULL,nkept=0; /* list of kept active set */ 1155 PetscScalar *f,*f2; 1156 PetscBool isequal; 1157 PetscInt i1=0,j1=0; 1158 1159 /* Call general purpose update function */ 1160 if (snes->ops->update) { 1161 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 1162 } 1163 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 1164 1165 /* Create active and inactive index sets */ 1166 ierr = SNESVICreateIndexSets_RS(snes,X,F,&IS_act,&IS_inact);CHKERRQ(ierr); 1167 1168 /* Get local active set size */ 1169 ierr = ISGetLocalSize(IS_act,&nis_act);CHKERRQ(ierr); 1170 if (nis_act) { 1171 ierr = ISGetIndices(IS_act,&idx_act);CHKERRQ(ierr); 1172 IS_redact = PETSC_NULL; 1173 if(vi->checkredundancy) { 1174 (*vi->checkredundancy)(snes,IS_act,&IS_redact,vi->ctxP); 1175 } 1176 1177 if(!IS_redact) { 1178 /* User called checkredundancy function but didn't create IS_redact because 1179 there were no redundant active set variables */ 1180 /* Copy over all active set indices to the list */ 1181 ierr = PetscMalloc(nis_act*sizeof(PetscInt),&idx_actkept);CHKERRQ(ierr); 1182 for(k=0;k < nis_act;k++) idx_actkept[k] = idx_act[k]; 1183 nkept = nis_act; 1184 } else { 1185 ierr = ISGetLocalSize(IS_redact,&nis_redact);CHKERRQ(ierr); 1186 ierr = PetscMalloc((nis_act-nis_redact)*sizeof(PetscInt),&idx_actkept);CHKERRQ(ierr); 1187 1188 /* Create reduced active set list */ 1189 ierr = ISGetIndices(IS_act,&idx_act);CHKERRQ(ierr); 1190 ierr = ISGetIndices(IS_redact,&idx_redact);CHKERRQ(ierr); 1191 j1 = 0;nkept = 0; 1192 for(k=0;k<nis_act;k++) { 1193 if(j1 < nis_redact && idx_act[k] == idx_redact[j1]) j1++; 1194 else idx_actkept[nkept++] = idx_act[k]; 1195 } 1196 ierr = ISRestoreIndices(IS_act,&idx_act);CHKERRQ(ierr); 1197 ierr = ISRestoreIndices(IS_redact,&idx_redact);CHKERRQ(ierr); 1198 1199 ierr = ISDestroy(&IS_redact);CHKERRQ(ierr); 1200 } 1201 1202 /* Create augmented F and Y */ 1203 ierr = VecCreate(((PetscObject)snes)->comm,&F_aug);CHKERRQ(ierr); 1204 ierr = VecSetSizes(F_aug,F->map->n+nkept,PETSC_DECIDE);CHKERRQ(ierr); 1205 ierr = VecSetFromOptions(F_aug);CHKERRQ(ierr); 1206 ierr = VecDuplicate(F_aug,&Y_aug);CHKERRQ(ierr); 1207 1208 /* Copy over F to F_aug in the first n locations */ 1209 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 1210 ierr = VecGetArray(F_aug,&f2);CHKERRQ(ierr); 1211 ierr = PetscMemcpy(f2,f,F->map->n*sizeof(PetscScalar));CHKERRQ(ierr); 1212 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 1213 ierr = VecRestoreArray(F_aug,&f2);CHKERRQ(ierr); 1214 1215 /* Create the augmented jacobian matrix */ 1216 ierr = MatCreate(((PetscObject)snes)->comm,&J_aug);CHKERRQ(ierr); 1217 ierr = MatSetSizes(J_aug,X->map->n+nkept,X->map->n+nkept,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 1218 ierr = MatSetFromOptions(J_aug);CHKERRQ(ierr); 1219 1220 1221 { /* local vars */ 1222 /* Preallocate augmented matrix and set values in it...Doing only sequential case first*/ 1223 PetscInt ncols; 1224 const PetscInt *cols; 1225 const PetscScalar *vals; 1226 PetscScalar value[2]; 1227 PetscInt row,col[2]; 1228 PetscInt *d_nnz; 1229 value[0] = 1.0; value[1] = 0.0; 1230 ierr = PetscMalloc((X->map->n+nkept)*sizeof(PetscInt),&d_nnz);CHKERRQ(ierr); 1231 ierr = PetscMemzero(d_nnz,(X->map->n+nkept)*sizeof(PetscInt));CHKERRQ(ierr); 1232 for(row=0;row<snes->jacobian->rmap->n;row++) { 1233 ierr = MatGetRow(snes->jacobian,row,&ncols,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 1234 d_nnz[row] += ncols; 1235 ierr = MatRestoreRow(snes->jacobian,row,&ncols,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 1236 } 1237 1238 for(k=0;k<nkept;k++) { 1239 d_nnz[idx_actkept[k]] += 1; 1240 d_nnz[snes->jacobian->rmap->n+k] += 2; 1241 } 1242 ierr = MatSeqAIJSetPreallocation(J_aug,PETSC_NULL,d_nnz);CHKERRQ(ierr); 1243 1244 ierr = PetscFree(d_nnz);CHKERRQ(ierr); 1245 1246 /* Set the original jacobian matrix in J_aug */ 1247 for(row=0;row<snes->jacobian->rmap->n;row++) { 1248 ierr = MatGetRow(snes->jacobian,row,&ncols,&cols,&vals);CHKERRQ(ierr); 1249 ierr = MatSetValues(J_aug,1,&row,ncols,cols,vals,INSERT_VALUES);CHKERRQ(ierr); 1250 ierr = MatRestoreRow(snes->jacobian,row,&ncols,&cols,&vals);CHKERRQ(ierr); 1251 } 1252 /* Add the augmented part */ 1253 for(k=0;k<nkept;k++) { 1254 row = snes->jacobian->rmap->n+k; 1255 col[0] = idx_actkept[k]; col[1] = snes->jacobian->rmap->n+k; 1256 ierr = MatSetValues(J_aug,1,&row,1,col,value,INSERT_VALUES);CHKERRQ(ierr); 1257 ierr = MatSetValues(J_aug,1,&col[0],1,&row,&value[0],INSERT_VALUES);CHKERRQ(ierr); 1258 } 1259 ierr = MatAssemblyBegin(J_aug,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1260 ierr = MatAssemblyEnd(J_aug,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1261 /* Only considering prejac = jac for now */ 1262 Jpre_aug = J_aug; 1263 } /* local vars*/ 1264 ierr = ISRestoreIndices(IS_act,&idx_act);CHKERRQ(ierr); 1265 ierr = ISDestroy(&IS_act);CHKERRQ(ierr); 1266 } else { 1267 F_aug = F; J_aug = snes->jacobian; Y_aug = Y; Jpre_aug = snes->jacobian_pre; 1268 } 1269 1270 ierr = ISEqual(vi->IS_inact_prev,IS_inact,&isequal);CHKERRQ(ierr); 1271 if (!isequal) { 1272 ierr = SNESVIResetPCandKSP(snes,J_aug,Jpre_aug);CHKERRQ(ierr); 1273 } 1274 ierr = KSPSetOperators(snes->ksp,J_aug,Jpre_aug,flg);CHKERRQ(ierr); 1275 ierr = KSPSetUp(snes->ksp);CHKERRQ(ierr); 1276 /* { 1277 PC pc; 1278 PetscBool flg; 1279 ierr = KSPGetPC(snes->ksp,&pc);CHKERRQ(ierr); 1280 ierr = PetscTypeCompare((PetscObject)pc,PCFIELDSPLIT,&flg);CHKERRQ(ierr); 1281 if (flg) { 1282 KSP *subksps; 1283 ierr = PCFieldSplitGetSubKSP(pc,PETSC_NULL,&subksps);CHKERRQ(ierr); 1284 ierr = KSPGetPC(subksps[0],&pc);CHKERRQ(ierr); 1285 ierr = PetscFree(subksps);CHKERRQ(ierr); 1286 ierr = PetscTypeCompare((PetscObject)pc,PCBJACOBI,&flg);CHKERRQ(ierr); 1287 if (flg) { 1288 PetscInt n,N = 101*101,j,cnts[3] = {0,0,0}; 1289 const PetscInt *ii; 1290 1291 ierr = ISGetSize(IS_inact,&n);CHKERRQ(ierr); 1292 ierr = ISGetIndices(IS_inact,&ii);CHKERRQ(ierr); 1293 for (j=0; j<n; j++) { 1294 if (ii[j] < N) cnts[0]++; 1295 else if (ii[j] < 2*N) cnts[1]++; 1296 else if (ii[j] < 3*N) cnts[2]++; 1297 } 1298 ierr = ISRestoreIndices(IS_inact,&ii);CHKERRQ(ierr); 1299 1300 ierr = PCBJacobiSetTotalBlocks(pc,3,cnts);CHKERRQ(ierr); 1301 } 1302 } 1303 } 1304 */ 1305 ierr = SNES_KSPSolve(snes,snes->ksp,F_aug,Y_aug);CHKERRQ(ierr); 1306 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 1307 if (kspreason < 0) { 1308 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 1309 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 1310 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 1311 break; 1312 } 1313 } 1314 1315 if(nis_act) { 1316 PetscScalar *y1,*y2; 1317 ierr = VecGetArray(Y,&y1);CHKERRQ(ierr); 1318 ierr = VecGetArray(Y_aug,&y2);CHKERRQ(ierr); 1319 /* Copy over inactive Y_aug to Y */ 1320 j1 = 0; 1321 for(i1=Y->map->rstart;i1 < Y->map->rend;i1++) { 1322 if(j1 < nkept && idx_actkept[j1] == i1) j1++; 1323 else y1[i1-Y->map->rstart] = y2[i1-Y->map->rstart]; 1324 } 1325 ierr = VecRestoreArray(Y,&y1);CHKERRQ(ierr); 1326 ierr = VecRestoreArray(Y_aug,&y2);CHKERRQ(ierr); 1327 1328 ierr = VecDestroy(&F_aug);CHKERRQ(ierr); 1329 ierr = VecDestroy(&Y_aug);CHKERRQ(ierr); 1330 ierr = MatDestroy(&J_aug);CHKERRQ(ierr); 1331 ierr = PetscFree(idx_actkept);CHKERRQ(ierr); 1332 } 1333 1334 if (!isequal) { 1335 ierr = ISDestroy(&vi->IS_inact_prev);CHKERRQ(ierr); 1336 ierr = ISDuplicate(IS_inact,&vi->IS_inact_prev);CHKERRQ(ierr); 1337 } 1338 ierr = ISDestroy(&IS_inact);CHKERRQ(ierr); 1339 1340 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 1341 snes->linear_its += lits; 1342 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 1343 /* 1344 if (vi->precheckstep) { 1345 PetscBool changed_y = PETSC_FALSE; 1346 ierr = (*vi->precheckstep)(snes,X,Y,vi->precheck,&changed_y);CHKERRQ(ierr); 1347 } 1348 1349 if (PetscLogPrintInfo){ 1350 ierr = SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr); 1351 } 1352 */ 1353 /* Compute a (scaled) negative update in the line search routine: 1354 Y <- X - lambda*Y 1355 and evaluate G = function(Y) (depends on the line search). 1356 */ 1357 ierr = VecCopy(Y,snes->vec_sol_update);CHKERRQ(ierr); 1358 ynorm = 1; gnorm = fnorm; 1359 ierr = (*vi->LineSearch)(snes,vi->lsP,X,F,G,Y,W,fnorm,xnorm,&ynorm,&gnorm,&lssucceed);CHKERRQ(ierr); 1360 ierr = PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",fnorm,gnorm,ynorm,(int)lssucceed);CHKERRQ(ierr); 1361 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 1362 if (snes->domainerror) { 1363 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 1364 PetscFunctionReturn(0); 1365 } 1366 if (!lssucceed) { 1367 if (++snes->numFailures >= snes->maxFailures) { 1368 PetscBool ismin; 1369 snes->reason = SNES_DIVERGED_LINE_SEARCH; 1370 ierr = SNESVICheckLocalMin_Private(snes,snes->jacobian,G,W,gnorm,&ismin);CHKERRQ(ierr); 1371 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 1372 break; 1373 } 1374 } 1375 /* Update function and solution vectors */ 1376 fnorm = gnorm; 1377 ierr = VecCopy(G,F);CHKERRQ(ierr); 1378 ierr = VecCopy(W,X);CHKERRQ(ierr); 1379 /* Monitor convergence */ 1380 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1381 snes->iter = i+1; 1382 snes->norm = fnorm; 1383 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1384 SNESLogConvHistory(snes,snes->norm,lits); 1385 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 1386 /* Test for convergence, xnorm = || X || */ 1387 if (snes->ops->converged != SNESSkipConverged) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 1388 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 1389 if (snes->reason) break; 1390 } 1391 if (i == maxits) { 1392 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 1393 if(!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 1394 } 1395 PetscFunctionReturn(0); 1396 } 1397 1398 /* -------------------------------------------------------------------------- */ 1399 /* 1400 SNESSetUp_VI - Sets up the internal data structures for the later use 1401 of the SNESVI nonlinear solver. 1402 1403 Input Parameter: 1404 . snes - the SNES context 1405 . x - the solution vector 1406 1407 Application Interface Routine: SNESSetUp() 1408 1409 Notes: 1410 For basic use of the SNES solvers, the user need not explicitly call 1411 SNESSetUp(), since these actions will automatically occur during 1412 the call to SNESSolve(). 1413 */ 1414 #undef __FUNCT__ 1415 #define __FUNCT__ "SNESSetUp_VI" 1416 PetscErrorCode SNESSetUp_VI(SNES snes) 1417 { 1418 PetscErrorCode ierr; 1419 SNES_VI *vi = (SNES_VI*) snes->data; 1420 PetscInt i_start[3],i_end[3]; 1421 1422 PetscFunctionBegin; 1423 1424 ierr = SNESDefaultGetWork(snes,3);CHKERRQ(ierr); 1425 1426 /* If the lower and upper bound on variables are not set, set it to 1427 -Inf and Inf */ 1428 if (!vi->xl && !vi->xu) { 1429 vi->usersetxbounds = PETSC_FALSE; 1430 ierr = VecDuplicate(snes->vec_sol, &vi->xl); CHKERRQ(ierr); 1431 ierr = VecSet(vi->xl,SNES_VI_NINF);CHKERRQ(ierr); 1432 ierr = VecDuplicate(snes->vec_sol, &vi->xu); CHKERRQ(ierr); 1433 ierr = VecSet(vi->xu,SNES_VI_INF);CHKERRQ(ierr); 1434 } else { 1435 /* Check if lower bound, upper bound and solution vector distribution across the processors is identical */ 1436 ierr = VecGetOwnershipRange(snes->vec_sol,i_start,i_end);CHKERRQ(ierr); 1437 ierr = VecGetOwnershipRange(vi->xl,i_start+1,i_end+1);CHKERRQ(ierr); 1438 ierr = VecGetOwnershipRange(vi->xu,i_start+2,i_end+2);CHKERRQ(ierr); 1439 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])) 1440 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."); 1441 } 1442 if (snes->ops->solve != SNESSolveVI_SS) { 1443 /* Set up previous active index set for the first snes solve 1444 vi->IS_inact_prev = 0,1,2,....N */ 1445 PetscInt *indices; 1446 PetscInt i,n,rstart,rend; 1447 1448 ierr = VecGetOwnershipRange(snes->vec_sol,&rstart,&rend);CHKERRQ(ierr); 1449 ierr = VecGetLocalSize(snes->vec_sol,&n);CHKERRQ(ierr); 1450 ierr = PetscMalloc(n*sizeof(PetscInt),&indices);CHKERRQ(ierr); 1451 for(i=0;i < n; i++) indices[i] = rstart + i; 1452 ierr = ISCreateGeneral(((PetscObject)snes)->comm,n,indices,PETSC_OWN_POINTER,&vi->IS_inact_prev); 1453 } 1454 1455 if (snes->ops->solve == SNESSolveVI_SS) { 1456 ierr = VecDuplicate(snes->vec_sol, &vi->dpsi);CHKERRQ(ierr); 1457 ierr = VecDuplicate(snes->vec_sol, &vi->phi); CHKERRQ(ierr); 1458 ierr = VecDuplicate(snes->vec_sol, &vi->Da); CHKERRQ(ierr); 1459 ierr = VecDuplicate(snes->vec_sol, &vi->Db); CHKERRQ(ierr); 1460 ierr = VecDuplicate(snes->vec_sol, &vi->z);CHKERRQ(ierr); 1461 ierr = VecDuplicate(snes->vec_sol, &vi->t); CHKERRQ(ierr); 1462 } 1463 PetscFunctionReturn(0); 1464 } 1465 /* -------------------------------------------------------------------------- */ 1466 /* 1467 SNESDestroy_VI - Destroys the private SNES_VI context that was created 1468 with SNESCreate_VI(). 1469 1470 Input Parameter: 1471 . snes - the SNES context 1472 1473 Application Interface Routine: SNESDestroy() 1474 */ 1475 #undef __FUNCT__ 1476 #define __FUNCT__ "SNESDestroy_VI" 1477 PetscErrorCode SNESDestroy_VI(SNES snes) 1478 { 1479 SNES_VI *vi = (SNES_VI*) snes->data; 1480 PetscErrorCode ierr; 1481 1482 PetscFunctionBegin; 1483 if (snes->ops->solve != SNESSolveVI_SS) { 1484 ierr = ISDestroy(&vi->IS_inact_prev); 1485 } 1486 1487 if (snes->ops->solve == SNESSolveVI_SS) { 1488 /* clear vectors */ 1489 ierr = VecDestroy(&vi->dpsi);CHKERRQ(ierr); 1490 ierr = VecDestroy(&vi->phi); CHKERRQ(ierr); 1491 ierr = VecDestroy(&vi->Da); CHKERRQ(ierr); 1492 ierr = VecDestroy(&vi->Db); CHKERRQ(ierr); 1493 ierr = VecDestroy(&vi->z); CHKERRQ(ierr); 1494 ierr = VecDestroy(&vi->t); CHKERRQ(ierr); 1495 } 1496 1497 if (!vi->usersetxbounds) { 1498 ierr = VecDestroy(&vi->xl); CHKERRQ(ierr); 1499 ierr = VecDestroy(&vi->xu); CHKERRQ(ierr); 1500 } 1501 1502 ierr = PetscViewerASCIIMonitorDestroy(&vi->lsmonitor);CHKERRQ(ierr); 1503 ierr = PetscFree(snes->data);CHKERRQ(ierr); 1504 1505 /* clear composed functions */ 1506 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSet_C","",PETSC_NULL);CHKERRQ(ierr); 1507 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetMonitor_C","",PETSC_NULL);CHKERRQ(ierr); 1508 PetscFunctionReturn(0); 1509 } 1510 1511 /* -------------------------------------------------------------------------- */ 1512 #undef __FUNCT__ 1513 #define __FUNCT__ "SNESLineSearchNo_VI" 1514 1515 /* 1516 This routine does not actually do a line search but it takes a full newton 1517 step while ensuring that the new iterates remain within the constraints. 1518 1519 */ 1520 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) 1521 { 1522 PetscErrorCode ierr; 1523 SNES_VI *vi = (SNES_VI*)snes->data; 1524 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1525 1526 PetscFunctionBegin; 1527 *flag = PETSC_TRUE; 1528 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1529 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); /* ynorm = || y || */ 1530 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1531 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1532 if (vi->postcheckstep) { 1533 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1534 } 1535 if (changed_y) { 1536 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1537 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1538 } 1539 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1540 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1541 if (!snes->domainerror) { 1542 if (snes->ops->solve != SNESSolveVI_SS) { 1543 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1544 } else { 1545 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); /* gnorm = || g || */ 1546 } 1547 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1548 } 1549 if (vi->lsmonitor) { 1550 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1551 } 1552 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1553 PetscFunctionReturn(0); 1554 } 1555 1556 /* -------------------------------------------------------------------------- */ 1557 #undef __FUNCT__ 1558 #define __FUNCT__ "SNESLineSearchNoNorms_VI" 1559 1560 /* 1561 This routine is a copy of SNESLineSearchNoNorms in snes/impls/ls/ls.c 1562 */ 1563 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) 1564 { 1565 PetscErrorCode ierr; 1566 SNES_VI *vi = (SNES_VI*)snes->data; 1567 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1568 1569 PetscFunctionBegin; 1570 *flag = PETSC_TRUE; 1571 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1572 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1573 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1574 if (vi->postcheckstep) { 1575 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1576 } 1577 if (changed_y) { 1578 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1579 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1580 } 1581 1582 /* don't evaluate function the last time through */ 1583 if (snes->iter < snes->max_its-1) { 1584 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1585 } 1586 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1587 PetscFunctionReturn(0); 1588 } 1589 1590 /* -------------------------------------------------------------------------- */ 1591 #undef __FUNCT__ 1592 #define __FUNCT__ "SNESLineSearchCubic_VI" 1593 /* 1594 This routine implements a cubic line search while doing a projection on the variable bounds 1595 */ 1596 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) 1597 { 1598 PetscReal initslope,lambdaprev,gnormprev,a,b,d,t1,t2,rellength; 1599 PetscReal minlambda,lambda,lambdatemp; 1600 #if defined(PETSC_USE_COMPLEX) 1601 PetscScalar cinitslope; 1602 #endif 1603 PetscErrorCode ierr; 1604 PetscInt count; 1605 SNES_VI *vi = (SNES_VI*)snes->data; 1606 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1607 MPI_Comm comm; 1608 1609 PetscFunctionBegin; 1610 ierr = PetscObjectGetComm((PetscObject)snes,&comm);CHKERRQ(ierr); 1611 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1612 *flag = PETSC_TRUE; 1613 1614 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); 1615 if (*ynorm == 0.0) { 1616 if (vi->lsmonitor) { 1617 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Initial direction and size is 0\n");CHKERRQ(ierr); 1618 } 1619 *gnorm = fnorm; 1620 ierr = VecCopy(x,w);CHKERRQ(ierr); 1621 ierr = VecCopy(f,g);CHKERRQ(ierr); 1622 *flag = PETSC_FALSE; 1623 goto theend1; 1624 } 1625 if (*ynorm > vi->maxstep) { /* Step too big, so scale back */ 1626 if (vi->lsmonitor) { 1627 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Scaling step by %G old ynorm %G\n",vi->maxstep/(*ynorm),*ynorm);CHKERRQ(ierr); 1628 } 1629 ierr = VecScale(y,vi->maxstep/(*ynorm));CHKERRQ(ierr); 1630 *ynorm = vi->maxstep; 1631 } 1632 ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr); 1633 minlambda = vi->minlambda/rellength; 1634 ierr = MatMult(snes->jacobian,y,w);CHKERRQ(ierr); 1635 #if defined(PETSC_USE_COMPLEX) 1636 ierr = VecDot(f,w,&cinitslope);CHKERRQ(ierr); 1637 initslope = PetscRealPart(cinitslope); 1638 #else 1639 ierr = VecDot(f,w,&initslope);CHKERRQ(ierr); 1640 #endif 1641 if (initslope > 0.0) initslope = -initslope; 1642 if (initslope == 0.0) initslope = -1.0; 1643 1644 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1645 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1646 if (snes->nfuncs >= snes->max_funcs) { 1647 ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); 1648 *flag = PETSC_FALSE; 1649 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1650 goto theend1; 1651 } 1652 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1653 if (snes->ops->solve != SNESSolveVI_SS) { 1654 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1655 } else { 1656 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1657 } 1658 if (snes->domainerror) { 1659 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1660 PetscFunctionReturn(0); 1661 } 1662 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1663 ierr = PetscInfo2(snes,"Initial fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1664 if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + vi->alpha*initslope) { /* Sufficient reduction */ 1665 if (vi->lsmonitor) { 1666 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1667 } 1668 goto theend1; 1669 } 1670 1671 /* Fit points with quadratic */ 1672 lambda = 1.0; 1673 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 1674 lambdaprev = lambda; 1675 gnormprev = *gnorm; 1676 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1677 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1678 else lambda = lambdatemp; 1679 1680 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1681 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1682 if (snes->nfuncs >= snes->max_funcs) { 1683 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);CHKERRQ(ierr); 1684 *flag = PETSC_FALSE; 1685 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1686 goto theend1; 1687 } 1688 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1689 if (snes->ops->solve != SNESSolveVI_SS) { 1690 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1691 } else { 1692 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1693 } 1694 if (snes->domainerror) { 1695 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1696 PetscFunctionReturn(0); 1697 } 1698 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1699 if (vi->lsmonitor) { 1700 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: gnorm after quadratic fit %G\n",*gnorm);CHKERRQ(ierr); 1701 } 1702 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* sufficient reduction */ 1703 if (vi->lsmonitor) { 1704 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Quadratically determined step, lambda=%18.16e\n",lambda);CHKERRQ(ierr); 1705 } 1706 goto theend1; 1707 } 1708 1709 /* Fit points with cubic */ 1710 count = 1; 1711 while (PETSC_TRUE) { 1712 if (lambda <= minlambda) { 1713 if (vi->lsmonitor) { 1714 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: unable to find good step length! After %D tries \n",count);CHKERRQ(ierr); 1715 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); 1716 } 1717 *flag = PETSC_FALSE; 1718 break; 1719 } 1720 t1 = .5*((*gnorm)*(*gnorm) - fnorm*fnorm) - lambda*initslope; 1721 t2 = .5*(gnormprev*gnormprev - fnorm*fnorm) - lambdaprev*initslope; 1722 a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 1723 b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 1724 d = b*b - 3*a*initslope; 1725 if (d < 0.0) d = 0.0; 1726 if (a == 0.0) { 1727 lambdatemp = -initslope/(2.0*b); 1728 } else { 1729 lambdatemp = (-b + sqrt(d))/(3.0*a); 1730 } 1731 lambdaprev = lambda; 1732 gnormprev = *gnorm; 1733 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1734 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1735 else lambda = lambdatemp; 1736 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1737 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1738 if (snes->nfuncs >= snes->max_funcs) { 1739 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 1740 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); 1741 *flag = PETSC_FALSE; 1742 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1743 break; 1744 } 1745 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1746 if (snes->ops->solve != SNESSolveVI_SS) { 1747 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1748 } else { 1749 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1750 } 1751 if (snes->domainerror) { 1752 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1753 PetscFunctionReturn(0); 1754 } 1755 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1756 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* is reduction enough? */ 1757 if (vi->lsmonitor) { 1758 ierr = PetscPrintf(comm," Line search: Cubically determined step, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr); 1759 } 1760 break; 1761 } else { 1762 if (vi->lsmonitor) { 1763 ierr = PetscPrintf(comm," Line search: Cubic step no good, shrinking lambda, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr); 1764 } 1765 } 1766 count++; 1767 } 1768 theend1: 1769 /* Optional user-defined check for line search step validity */ 1770 if (vi->postcheckstep && *flag) { 1771 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1772 if (changed_y) { 1773 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1774 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1775 } 1776 if (changed_y || changed_w) { /* recompute the function if the step has changed */ 1777 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1778 if (snes->ops->solve != SNESSolveVI_SS) { 1779 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1780 } else { 1781 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1782 } 1783 if (snes->domainerror) { 1784 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1785 PetscFunctionReturn(0); 1786 } 1787 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1788 ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr); 1789 ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr); 1790 } 1791 } 1792 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1793 PetscFunctionReturn(0); 1794 } 1795 1796 #undef __FUNCT__ 1797 #define __FUNCT__ "SNESLineSearchQuadratic_VI" 1798 /* 1799 This routine does a quadratic line search while keeping the iterates within the variable bounds 1800 */ 1801 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) 1802 { 1803 /* 1804 Note that for line search purposes we work with with the related 1805 minimization problem: 1806 min z(x): R^n -> R, 1807 where z(x) = .5 * fnorm*fnorm,and fnorm = || f ||_2. 1808 */ 1809 PetscReal initslope,minlambda,lambda,lambdatemp,rellength; 1810 #if defined(PETSC_USE_COMPLEX) 1811 PetscScalar cinitslope; 1812 #endif 1813 PetscErrorCode ierr; 1814 PetscInt count; 1815 SNES_VI *vi = (SNES_VI*)snes->data; 1816 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1817 1818 PetscFunctionBegin; 1819 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1820 *flag = PETSC_TRUE; 1821 1822 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); 1823 if (*ynorm == 0.0) { 1824 if (vi->lsmonitor) { 1825 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor,"Line search: Direction and size is 0\n");CHKERRQ(ierr); 1826 } 1827 *gnorm = fnorm; 1828 ierr = VecCopy(x,w);CHKERRQ(ierr); 1829 ierr = VecCopy(f,g);CHKERRQ(ierr); 1830 *flag = PETSC_FALSE; 1831 goto theend2; 1832 } 1833 if (*ynorm > vi->maxstep) { /* Step too big, so scale back */ 1834 ierr = VecScale(y,vi->maxstep/(*ynorm));CHKERRQ(ierr); 1835 *ynorm = vi->maxstep; 1836 } 1837 ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr); 1838 minlambda = vi->minlambda/rellength; 1839 ierr = MatMult(snes->jacobian,y,w);CHKERRQ(ierr); 1840 #if defined(PETSC_USE_COMPLEX) 1841 ierr = VecDot(f,w,&cinitslope);CHKERRQ(ierr); 1842 initslope = PetscRealPart(cinitslope); 1843 #else 1844 ierr = VecDot(f,w,&initslope);CHKERRQ(ierr); 1845 #endif 1846 if (initslope > 0.0) initslope = -initslope; 1847 if (initslope == 0.0) initslope = -1.0; 1848 ierr = PetscInfo1(snes,"Initslope %G \n",initslope);CHKERRQ(ierr); 1849 1850 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1851 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1852 if (snes->nfuncs >= snes->max_funcs) { 1853 ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); 1854 *flag = PETSC_FALSE; 1855 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1856 goto theend2; 1857 } 1858 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1859 if (snes->ops->solve != SNESSolveVI_SS) { 1860 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1861 } else { 1862 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1863 } 1864 if (snes->domainerror) { 1865 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1866 PetscFunctionReturn(0); 1867 } 1868 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1869 if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + vi->alpha*initslope) { /* Sufficient reduction */ 1870 if (vi->lsmonitor) { 1871 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1872 } 1873 goto theend2; 1874 } 1875 1876 /* Fit points with quadratic */ 1877 lambda = 1.0; 1878 count = 1; 1879 while (PETSC_TRUE) { 1880 if (lambda <= minlambda) { /* bad luck; use full step */ 1881 if (vi->lsmonitor) { 1882 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor,"Line search: Unable to find good step length! %D \n",count);CHKERRQ(ierr); 1883 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); 1884 } 1885 ierr = VecCopy(x,w);CHKERRQ(ierr); 1886 *flag = PETSC_FALSE; 1887 break; 1888 } 1889 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 1890 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1891 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1892 else lambda = lambdatemp; 1893 1894 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1895 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1896 if (snes->nfuncs >= snes->max_funcs) { 1897 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 1898 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); 1899 *flag = PETSC_FALSE; 1900 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1901 break; 1902 } 1903 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1904 if (snes->domainerror) { 1905 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1906 PetscFunctionReturn(0); 1907 } 1908 if (snes->ops->solve != SNESSolveVI_SS) { 1909 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1910 } else { 1911 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1912 } 1913 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1914 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* sufficient reduction */ 1915 if (vi->lsmonitor) { 1916 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line Search: Quadratically determined step, lambda=%G\n",lambda);CHKERRQ(ierr); 1917 } 1918 break; 1919 } 1920 count++; 1921 } 1922 theend2: 1923 /* Optional user-defined check for line search step validity */ 1924 if (vi->postcheckstep) { 1925 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1926 if (changed_y) { 1927 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1928 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1929 } 1930 if (changed_y || changed_w) { /* recompute the function if the step has changed */ 1931 ierr = SNESComputeFunction(snes,w,g); 1932 if (snes->domainerror) { 1933 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1934 PetscFunctionReturn(0); 1935 } 1936 if (snes->ops->solve != SNESSolveVI_SS) { 1937 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1938 } else { 1939 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1940 } 1941 1942 ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr); 1943 ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr); 1944 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1945 } 1946 } 1947 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1948 PetscFunctionReturn(0); 1949 } 1950 1951 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*/ 1952 /* -------------------------------------------------------------------------- */ 1953 EXTERN_C_BEGIN 1954 #undef __FUNCT__ 1955 #define __FUNCT__ "SNESLineSearchSet_VI" 1956 PetscErrorCode SNESLineSearchSet_VI(SNES snes,FCN2 func,void *lsctx) 1957 { 1958 PetscFunctionBegin; 1959 ((SNES_VI *)(snes->data))->LineSearch = func; 1960 ((SNES_VI *)(snes->data))->lsP = lsctx; 1961 PetscFunctionReturn(0); 1962 } 1963 EXTERN_C_END 1964 1965 /* -------------------------------------------------------------------------- */ 1966 EXTERN_C_BEGIN 1967 #undef __FUNCT__ 1968 #define __FUNCT__ "SNESLineSearchSetMonitor_VI" 1969 PetscErrorCode SNESLineSearchSetMonitor_VI(SNES snes,PetscBool flg) 1970 { 1971 SNES_VI *vi = (SNES_VI*)snes->data; 1972 PetscErrorCode ierr; 1973 1974 PetscFunctionBegin; 1975 if (flg && !vi->lsmonitor) { 1976 ierr = PetscViewerASCIIMonitorCreate(((PetscObject)snes)->comm,"stdout",((PetscObject)snes)->tablevel,&vi->lsmonitor);CHKERRQ(ierr); 1977 } else if (!flg && vi->lsmonitor) { 1978 ierr = PetscViewerASCIIMonitorDestroy(&vi->lsmonitor);CHKERRQ(ierr); 1979 } 1980 PetscFunctionReturn(0); 1981 } 1982 EXTERN_C_END 1983 1984 /* 1985 SNESView_VI - Prints info from the SNESVI data structure. 1986 1987 Input Parameters: 1988 . SNES - the SNES context 1989 . viewer - visualization context 1990 1991 Application Interface Routine: SNESView() 1992 */ 1993 #undef __FUNCT__ 1994 #define __FUNCT__ "SNESView_VI" 1995 static PetscErrorCode SNESView_VI(SNES snes,PetscViewer viewer) 1996 { 1997 SNES_VI *vi = (SNES_VI *)snes->data; 1998 const char *cstr,*tstr; 1999 PetscErrorCode ierr; 2000 PetscBool iascii; 2001 2002 PetscFunctionBegin; 2003 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 2004 if (iascii) { 2005 if (vi->LineSearch == SNESLineSearchNo_VI) cstr = "SNESLineSearchNo"; 2006 else if (vi->LineSearch == SNESLineSearchQuadratic_VI) cstr = "SNESLineSearchQuadratic"; 2007 else if (vi->LineSearch == SNESLineSearchCubic_VI) cstr = "SNESLineSearchCubic"; 2008 else cstr = "unknown"; 2009 if (snes->ops->solve == SNESSolveVI_SS) tstr = "Semismooth"; 2010 else if (snes->ops->solve == SNESSolveVI_RS) tstr = "Reduced Space"; 2011 else if (snes->ops->solve == SNESSolveVI_RSAUG) tstr = "Reduced space with augmented variables"; 2012 else tstr = "unknown"; 2013 ierr = PetscViewerASCIIPrintf(viewer," VI algorithm: %s\n",tstr);CHKERRQ(ierr); 2014 ierr = PetscViewerASCIIPrintf(viewer," line search variant: %s\n",cstr);CHKERRQ(ierr); 2015 ierr = PetscViewerASCIIPrintf(viewer," alpha=%G, maxstep=%G, minlambda=%G\n",vi->alpha,vi->maxstep,vi->minlambda);CHKERRQ(ierr); 2016 } else { 2017 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for SNES EQ VI",((PetscObject)viewer)->type_name); 2018 } 2019 PetscFunctionReturn(0); 2020 } 2021 2022 #undef __FUNCT__ 2023 #define __FUNCT__ "SNESVISetVariableBounds" 2024 /*@ 2025 SNESVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu. 2026 2027 Input Parameters: 2028 . snes - the SNES context. 2029 . xl - lower bound. 2030 . xu - upper bound. 2031 2032 Notes: 2033 If this routine is not called then the lower and upper bounds are set to 2034 SNES_VI_INF and SNES_VI_NINF respectively during SNESSetUp(). 2035 2036 @*/ 2037 PetscErrorCode SNESVISetVariableBounds(SNES snes, Vec xl, Vec xu) 2038 { 2039 SNES_VI *vi; 2040 PetscErrorCode ierr; 2041 2042 PetscFunctionBegin; 2043 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2044 PetscValidHeaderSpecific(xl,VEC_CLASSID,2); 2045 PetscValidHeaderSpecific(xu,VEC_CLASSID,3); 2046 if (!snes->vec_func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call SNESSetFunction() first"); 2047 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); 2048 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); 2049 ierr = SNESSetType(snes,SNESVI);CHKERRQ(ierr); 2050 vi = (SNES_VI*)snes->data; 2051 vi->usersetxbounds = PETSC_TRUE; 2052 vi->xl = xl; 2053 vi->xu = xu; 2054 PetscFunctionReturn(0); 2055 } 2056 2057 /* -------------------------------------------------------------------------- */ 2058 /* 2059 SNESSetFromOptions_VI - Sets various parameters for the SNESVI method. 2060 2061 Input Parameter: 2062 . snes - the SNES context 2063 2064 Application Interface Routine: SNESSetFromOptions() 2065 */ 2066 #undef __FUNCT__ 2067 #define __FUNCT__ "SNESSetFromOptions_VI" 2068 static PetscErrorCode SNESSetFromOptions_VI(SNES snes) 2069 { 2070 SNES_VI *vi = (SNES_VI *)snes->data; 2071 const char *lses[] = {"basic","basicnonorms","quadratic","cubic"}; 2072 const char *vies[] = {"ss","rs","rsaug"}; 2073 PetscErrorCode ierr; 2074 PetscInt indx; 2075 PetscBool flg,set,flg2; 2076 2077 PetscFunctionBegin; 2078 ierr = PetscOptionsHead("SNES semismooth method options");CHKERRQ(ierr); 2079 ierr = PetscOptionsBool("-snes_vi_monitor","Monitor all non-active variables","None",PETSC_FALSE,&flg,0);CHKERRQ(ierr); 2080 if (flg) { 2081 ierr = SNESMonitorSet(snes,SNESMonitorVI,0,0);CHKERRQ(ierr); 2082 } 2083 ierr = PetscOptionsReal("-snes_ls_alpha","Function norm must decrease by","None",vi->alpha,&vi->alpha,0);CHKERRQ(ierr); 2084 ierr = PetscOptionsReal("-snes_ls_maxstep","Step must be less than","None",vi->maxstep,&vi->maxstep,0);CHKERRQ(ierr); 2085 ierr = PetscOptionsReal("-snes_ls_minlambda","Minimum lambda allowed","None",vi->minlambda,&vi->minlambda,0);CHKERRQ(ierr); 2086 ierr = PetscOptionsReal("-snes_vi_const_tol","constraint tolerance","None",vi->const_tol,&vi->const_tol,0);CHKERRQ(ierr); 2087 ierr = PetscOptionsBool("-snes_ls_monitor","Print progress of line searches","SNESLineSearchSetMonitor",vi->lsmonitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);CHKERRQ(ierr); 2088 if (set) {ierr = SNESLineSearchSetMonitor(snes,flg);CHKERRQ(ierr);} 2089 ierr = PetscOptionsEList("-snes_vi_type","Semismooth algorithm used","",vies,3,"ss",&indx,&flg2);CHKERRQ(ierr); 2090 if (flg2) { 2091 switch (indx) { 2092 case 0: 2093 snes->ops->solve = SNESSolveVI_SS; 2094 break; 2095 case 1: 2096 snes->ops->solve = SNESSolveVI_RS; 2097 break; 2098 case 2: 2099 snes->ops->solve = SNESSolveVI_RSAUG; 2100 } 2101 } 2102 ierr = PetscOptionsEList("-snes_ls","Line search used","SNESLineSearchSet",lses,4,"basic",&indx,&flg);CHKERRQ(ierr); 2103 if (flg) { 2104 switch (indx) { 2105 case 0: 2106 ierr = SNESLineSearchSet(snes,SNESLineSearchNo_VI,PETSC_NULL);CHKERRQ(ierr); 2107 break; 2108 case 1: 2109 ierr = SNESLineSearchSet(snes,SNESLineSearchNoNorms_VI,PETSC_NULL);CHKERRQ(ierr); 2110 break; 2111 case 2: 2112 ierr = SNESLineSearchSet(snes,SNESLineSearchQuadratic_VI,PETSC_NULL);CHKERRQ(ierr); 2113 break; 2114 case 3: 2115 ierr = SNESLineSearchSet(snes,SNESLineSearchCubic_VI,PETSC_NULL);CHKERRQ(ierr); 2116 break; 2117 } 2118 } 2119 ierr = PetscOptionsTail();CHKERRQ(ierr); 2120 PetscFunctionReturn(0); 2121 } 2122 /* -------------------------------------------------------------------------- */ 2123 /*MC 2124 SNESVI - Various solvers for variational inequalities based on Newton's method 2125 2126 Options Database: 2127 + -snes_vi_type <ss,rs,rsaug> a semi-smooth solver, a reduced space active set method and a reduced space active set method that does not eliminate the active constraints from the Jacobian instead augments the Jacobian with 2128 additional variables that enforce the constraints 2129 - -snes_vi_monitor - prints the number of active constraints at each iteration. 2130 2131 2132 Level: beginner 2133 2134 .seealso: SNESCreate(), SNES, SNESSetType(), SNESTR, SNESLineSearchSet(), 2135 SNESLineSearchSetPostCheck(), SNESLineSearchNo(), SNESLineSearchCubic(), SNESLineSearchQuadratic(), 2136 SNESLineSearchSet(), SNESLineSearchNoNorms(), SNESLineSearchSetPreCheck(), SNESLineSearchSetParams(), SNESLineSearchGetParams() 2137 2138 M*/ 2139 EXTERN_C_BEGIN 2140 #undef __FUNCT__ 2141 #define __FUNCT__ "SNESCreate_VI" 2142 PetscErrorCode SNESCreate_VI(SNES snes) 2143 { 2144 PetscErrorCode ierr; 2145 SNES_VI *vi; 2146 2147 PetscFunctionBegin; 2148 snes->ops->setup = SNESSetUp_VI; 2149 snes->ops->solve = SNESSolveVI_RS; 2150 snes->ops->destroy = SNESDestroy_VI; 2151 snes->ops->setfromoptions = SNESSetFromOptions_VI; 2152 snes->ops->view = SNESView_VI; 2153 snes->ops->reset = 0; /* XXX Implement!!! */ 2154 snes->ops->converged = SNESDefaultConverged_VI; 2155 2156 ierr = PetscNewLog(snes,SNES_VI,&vi);CHKERRQ(ierr); 2157 snes->data = (void*)vi; 2158 vi->alpha = 1.e-4; 2159 vi->maxstep = 1.e8; 2160 vi->minlambda = 1.e-12; 2161 vi->LineSearch = SNESLineSearchNo_VI; 2162 vi->lsP = PETSC_NULL; 2163 vi->postcheckstep = PETSC_NULL; 2164 vi->postcheck = PETSC_NULL; 2165 vi->precheckstep = PETSC_NULL; 2166 vi->precheck = PETSC_NULL; 2167 vi->const_tol = 2.2204460492503131e-16; 2168 vi->checkredundancy = PETSC_NULL; 2169 2170 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetMonitor_C","SNESLineSearchSetMonitor_VI",SNESLineSearchSetMonitor_VI);CHKERRQ(ierr); 2171 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSet_C","SNESLineSearchSet_VI",SNESLineSearchSet_VI);CHKERRQ(ierr); 2172 2173 PetscFunctionReturn(0); 2174 } 2175 EXTERN_C_END 2176