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