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