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