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