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 active set variables. The user can optionally provide an IS for 1024 a subset of 'redundant' active set variables which will treated as free 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 IS IS_act,IS_inact; /* _act -> active set _inact -> inactive set */ 1082 IS IS_redact=PETSC_NULL; /* redundant active set */ 1083 Mat J_aug,Jpre_aug; 1084 Vec F_aug,Y_aug; 1085 PetscInt nis_redact,nis_act; 1086 const PetscInt *idx_redact,*idx_act; 1087 PetscInt k; 1088 PetscInt *idx_actkept=PETSC_NULL,nkept=0; /* list of kept active set */ 1089 PetscScalar *f,*f2; 1090 PetscBool isequal; 1091 PetscInt i1=0,j1=0; 1092 1093 /* Call general purpose update function */ 1094 if (snes->ops->update) { 1095 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 1096 } 1097 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 1098 1099 /* Create active and inactive index sets */ 1100 ierr = SNESVICreateIndexSets_RS(snes,X,F,&IS_act,&IS_inact);CHKERRQ(ierr); 1101 1102 /* Get local active set size */ 1103 ierr = ISGetLocalSize(IS_act,&nis_act);CHKERRQ(ierr); 1104 ierr = ISGetIndices(IS_act,&idx_act);CHKERRQ(ierr); 1105 if(nis_act) { 1106 if(vi->checkredundancy) { 1107 (*vi->checkredundancy)(snes,IS_act,&IS_redact,snes->funP); 1108 } 1109 1110 if(!IS_redact) { 1111 /* User called checkredundancy function but didn't create IS_redact because 1112 there were no redundant active set variables */ 1113 /* Copy over all active set indices to the list */ 1114 ierr = PetscMalloc(nis_act*sizeof(PetscInt),&idx_actkept);CHKERRQ(ierr); 1115 for(k=0;k < nis_act;k++) idx_actkept[k] = idx_act[k]; 1116 nkept = nis_act; 1117 } else { 1118 ierr = ISGetLocalSize(IS_redact,&nis_redact);CHKERRQ(ierr); 1119 ierr = PetscMalloc((nis_act-nis_redact)*sizeof(PetscInt),&idx_actkept);CHKERRQ(ierr); 1120 1121 /* Create reduced active set list */ 1122 ierr = ISGetIndices(IS_act,&idx_act);CHKERRQ(ierr); 1123 ierr = ISGetIndices(IS_redact,&idx_redact);CHKERRQ(ierr); 1124 j1 = 0; 1125 for(k=0;k<nis_act;k++) { 1126 if(j1 < nis_redact && idx_act[k] == idx_redact[j1]) j1++; 1127 else idx_actkept[nkept++] = idx_act[k]; 1128 } 1129 ierr = ISRestoreIndices(IS_act,&idx_act);CHKERRQ(ierr); 1130 ierr = ISRestoreIndices(IS_redact,&idx_redact);CHKERRQ(ierr); 1131 } 1132 1133 /* Create augmented F and Y */ 1134 ierr = VecCreate(((PetscObject)snes)->comm,&F_aug);CHKERRQ(ierr); 1135 ierr = VecSetSizes(F_aug,F->map->n+nkept,PETSC_DECIDE);CHKERRQ(ierr); 1136 ierr = VecSetFromOptions(F_aug);CHKERRQ(ierr); 1137 ierr = VecDuplicate(F_aug,&Y_aug);CHKERRQ(ierr); 1138 1139 /* Copy over F to F_aug in the first n locations */ 1140 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 1141 ierr = VecGetArray(F_aug,&f2);CHKERRQ(ierr); 1142 ierr = PetscMemcpy(f2,f,F->map->n*sizeof(PetscScalar));CHKERRQ(ierr); 1143 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 1144 ierr = VecRestoreArray(F_aug,&f2);CHKERRQ(ierr); 1145 1146 /* Create the augmented jacobian matrix */ 1147 ierr = MatCreate(((PetscObject)snes)->comm,&J_aug);CHKERRQ(ierr); 1148 ierr = MatSetSizes(J_aug,X->map->n+nkept,X->map->n+nkept,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 1149 ierr = MatSetFromOptions(J_aug);CHKERRQ(ierr); 1150 1151 /* set preallocation info..will add this later */ 1152 1153 /* Set values in the augmented matrix...Doing only sequential case first*/ 1154 PetscInt ncols; 1155 const PetscInt *cols; 1156 const PetscScalar *vals; 1157 PetscScalar value=1.0; 1158 PetscInt row,col; 1159 1160 /* Set the original jacobian matrix in J_aug */ 1161 for(row=0;row<snes->jacobian->rmap->n;row++) { 1162 ierr = MatGetRow(snes->jacobian,row,&ncols,&cols,&vals);CHKERRQ(ierr); 1163 ierr = MatSetValues(J_aug,1,&row,ncols,cols,vals,INSERT_VALUES);CHKERRQ(ierr); 1164 ierr = MatRestoreRow(snes->jacobian,row,&ncols,&cols,&vals);CHKERRQ(ierr); 1165 } 1166 /* Add the augmented part */ 1167 for(k=0;k<nkept;k++) { 1168 row = idx_actkept[k]; 1169 col = snes->jacobian->rmap->n+k; 1170 ierr = MatSetValues(J_aug,1,&row,1,&col,&value,INSERT_VALUES);CHKERRQ(ierr); 1171 ierr = MatSetValues(J_aug,1,&col,1,&row,&value,INSERT_VALUES);CHKERRQ(ierr); 1172 } 1173 ierr = MatAssemblyBegin(J_aug,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1174 ierr = MatAssemblyEnd(J_aug,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1175 /* Only considering prejac = jac for now */ 1176 Jpre_aug = J_aug; 1177 1178 } else { 1179 F_aug = F; J_aug = snes->jacobian; Y_aug = Y; Jpre_aug = snes->jacobian_pre; 1180 } 1181 1182 ierr = ISEqual(vi->IS_inact_prev,IS_inact,&isequal);CHKERRQ(ierr); 1183 if (!isequal) { 1184 ierr = SNESVIResetPCandKSP(snes,J_aug,Jpre_aug);CHKERRQ(ierr); 1185 } 1186 ierr = KSPSetOperators(snes->ksp,J_aug,Jpre_aug,flg);CHKERRQ(ierr); 1187 ierr = KSPSetUp(snes->ksp);CHKERRQ(ierr); 1188 /* { 1189 PC pc; 1190 PetscBool flg; 1191 ierr = KSPGetPC(snes->ksp,&pc);CHKERRQ(ierr); 1192 ierr = PetscTypeCompare((PetscObject)pc,PCFIELDSPLIT,&flg);CHKERRQ(ierr); 1193 if (flg) { 1194 KSP *subksps; 1195 ierr = PCFieldSplitGetSubKSP(pc,PETSC_NULL,&subksps);CHKERRQ(ierr); 1196 ierr = KSPGetPC(subksps[0],&pc);CHKERRQ(ierr); 1197 ierr = PetscFree(subksps);CHKERRQ(ierr); 1198 ierr = PetscTypeCompare((PetscObject)pc,PCBJACOBI,&flg);CHKERRQ(ierr); 1199 if (flg) { 1200 PetscInt n,N = 101*101,j,cnts[3] = {0,0,0}; 1201 const PetscInt *ii; 1202 1203 ierr = ISGetSize(IS_inact,&n);CHKERRQ(ierr); 1204 ierr = ISGetIndices(IS_inact,&ii);CHKERRQ(ierr); 1205 for (j=0; j<n; j++) { 1206 if (ii[j] < N) cnts[0]++; 1207 else if (ii[j] < 2*N) cnts[1]++; 1208 else if (ii[j] < 3*N) cnts[2]++; 1209 } 1210 ierr = ISRestoreIndices(IS_inact,&ii);CHKERRQ(ierr); 1211 1212 ierr = PCBJacobiSetTotalBlocks(pc,3,cnts);CHKERRQ(ierr); 1213 } 1214 } 1215 } 1216 */ 1217 ierr = SNES_KSPSolve(snes,snes->ksp,F_aug,Y_aug);CHKERRQ(ierr); 1218 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 1219 if (kspreason < 0) { 1220 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 1221 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 1222 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 1223 break; 1224 } 1225 } 1226 1227 if(nis_act) { 1228 PetscScalar *y1,*y2; 1229 ierr = VecGetArray(Y,&y1);CHKERRQ(ierr); 1230 ierr = VecGetArray(Y_aug,&y2);CHKERRQ(ierr); 1231 /* Copy over inactive Y_aug to Y */ 1232 j1 = 0; 1233 for(i1=Y->map->rstart;i1 < Y->map->rend;i1++) { 1234 if(j1 < nkept && idx_actkept[j1] == i1) j1++; 1235 else y1[i1-Y->map->rstart] = y2[i1-Y->map->rstart]; 1236 } 1237 ierr = VecRestoreArray(Y,&y1);CHKERRQ(ierr); 1238 ierr = VecRestoreArray(Y_aug,&y2);CHKERRQ(ierr); 1239 1240 ierr = VecDestroy(F_aug);CHKERRQ(ierr); 1241 ierr = VecDestroy(Y_aug);CHKERRQ(ierr); 1242 ierr = MatDestroy(J_aug);CHKERRQ(ierr); 1243 ierr = PetscFree(idx_actkept);CHKERRQ(ierr); 1244 } 1245 1246 if (!isequal) { 1247 ierr = ISDestroy(vi->IS_inact_prev);CHKERRQ(ierr); 1248 ierr = ISDuplicate(IS_inact,&vi->IS_inact_prev);CHKERRQ(ierr); 1249 } 1250 ierr = ISDestroy(IS_inact);CHKERRQ(ierr); 1251 1252 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 1253 snes->linear_its += lits; 1254 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 1255 /* 1256 if (vi->precheckstep) { 1257 PetscBool changed_y = PETSC_FALSE; 1258 ierr = (*vi->precheckstep)(snes,X,Y,vi->precheck,&changed_y);CHKERRQ(ierr); 1259 } 1260 1261 if (PetscLogPrintInfo){ 1262 ierr = SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr); 1263 } 1264 */ 1265 /* Compute a (scaled) negative update in the line search routine: 1266 Y <- X - lambda*Y 1267 and evaluate G = function(Y) (depends on the line search). 1268 */ 1269 ierr = VecCopy(Y,snes->vec_sol_update);CHKERRQ(ierr); 1270 ynorm = 1; gnorm = fnorm; 1271 ierr = (*vi->LineSearch)(snes,vi->lsP,X,F,G,Y,W,fnorm,xnorm,&ynorm,&gnorm,&lssucceed);CHKERRQ(ierr); 1272 ierr = PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",fnorm,gnorm,ynorm,(int)lssucceed);CHKERRQ(ierr); 1273 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 1274 if (snes->domainerror) { 1275 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 1276 PetscFunctionReturn(0); 1277 } 1278 if (!lssucceed) { 1279 if (++snes->numFailures >= snes->maxFailures) { 1280 PetscBool ismin; 1281 snes->reason = SNES_DIVERGED_LINE_SEARCH; 1282 ierr = SNESVICheckLocalMin_Private(snes,snes->jacobian,G,W,gnorm,&ismin);CHKERRQ(ierr); 1283 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 1284 break; 1285 } 1286 } 1287 /* Update function and solution vectors */ 1288 fnorm = gnorm; 1289 ierr = VecCopy(G,F);CHKERRQ(ierr); 1290 ierr = VecCopy(W,X);CHKERRQ(ierr); 1291 /* Monitor convergence */ 1292 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1293 snes->iter = i+1; 1294 snes->norm = fnorm; 1295 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1296 SNESLogConvHistory(snes,snes->norm,lits); 1297 SNESMonitor(snes,snes->iter,snes->norm); 1298 /* Test for convergence, xnorm = || X || */ 1299 if (snes->ops->converged != SNESSkipConverged) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 1300 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 1301 if (snes->reason) break; 1302 } 1303 if (i == maxits) { 1304 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 1305 if(!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 1306 } 1307 PetscFunctionReturn(0); 1308 } 1309 1310 /* -------------------------------------------------------------------------- */ 1311 /* 1312 SNESSetUp_VI - Sets up the internal data structures for the later use 1313 of the SNESVI nonlinear solver. 1314 1315 Input Parameter: 1316 . snes - the SNES context 1317 . x - the solution vector 1318 1319 Application Interface Routine: SNESSetUp() 1320 1321 Notes: 1322 For basic use of the SNES solvers, the user need not explicitly call 1323 SNESSetUp(), since these actions will automatically occur during 1324 the call to SNESSolve(). 1325 */ 1326 #undef __FUNCT__ 1327 #define __FUNCT__ "SNESSetUp_VI" 1328 PetscErrorCode SNESSetUp_VI(SNES snes) 1329 { 1330 PetscErrorCode ierr; 1331 SNES_VI *vi = (SNES_VI*) snes->data; 1332 PetscInt i_start[3],i_end[3]; 1333 1334 PetscFunctionBegin; 1335 if (!snes->vec_sol_update) { 1336 ierr = VecDuplicate(snes->vec_sol,&snes->vec_sol_update);CHKERRQ(ierr); 1337 ierr = PetscLogObjectParent(snes,snes->vec_sol_update);CHKERRQ(ierr); 1338 } 1339 if (!snes->work) { 1340 snes->nwork = 3; 1341 ierr = VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work);CHKERRQ(ierr); 1342 ierr = PetscLogObjectParents(snes,snes->nwork,snes->work);CHKERRQ(ierr); 1343 } 1344 1345 /* If the lower and upper bound on variables are not set, set it to 1346 -Inf and Inf */ 1347 if (!vi->xl && !vi->xu) { 1348 vi->usersetxbounds = PETSC_FALSE; 1349 ierr = VecDuplicate(snes->vec_sol, &vi->xl); CHKERRQ(ierr); 1350 ierr = VecSet(vi->xl,PETSC_VI_NINF);CHKERRQ(ierr); 1351 ierr = VecDuplicate(snes->vec_sol, &vi->xu); CHKERRQ(ierr); 1352 ierr = VecSet(vi->xu,PETSC_VI_INF);CHKERRQ(ierr); 1353 } else { 1354 /* Check if lower bound, upper bound and solution vector distribution across the processors is identical */ 1355 ierr = VecGetOwnershipRange(snes->vec_sol,i_start,i_end);CHKERRQ(ierr); 1356 ierr = VecGetOwnershipRange(vi->xl,i_start+1,i_end+1);CHKERRQ(ierr); 1357 ierr = VecGetOwnershipRange(vi->xu,i_start+2,i_end+2);CHKERRQ(ierr); 1358 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])) 1359 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."); 1360 } 1361 if (snes->ops->solve != SNESSolveVI_SS) { 1362 /* Set up previous active index set for the first snes solve 1363 vi->IS_inact_prev = 0,1,2,....N */ 1364 PetscInt *indices; 1365 PetscInt i,n,rstart,rend; 1366 1367 ierr = VecGetOwnershipRange(snes->vec_sol,&rstart,&rend);CHKERRQ(ierr); 1368 ierr = VecGetLocalSize(snes->vec_sol,&n);CHKERRQ(ierr); 1369 ierr = PetscMalloc(n*sizeof(PetscInt),&indices);CHKERRQ(ierr); 1370 for(i=0;i < n; i++) indices[i] = rstart + i; 1371 ierr = ISCreateGeneral(((PetscObject)snes)->comm,n,indices,PETSC_OWN_POINTER,&vi->IS_inact_prev); 1372 } 1373 1374 if (snes->ops->solve == SNESSolveVI_SS) { 1375 ierr = VecDuplicate(snes->vec_sol, &vi->dpsi);CHKERRQ(ierr); 1376 ierr = VecDuplicate(snes->vec_sol, &vi->phi); CHKERRQ(ierr); 1377 ierr = VecDuplicate(snes->vec_sol, &vi->Da); CHKERRQ(ierr); 1378 ierr = VecDuplicate(snes->vec_sol, &vi->Db); CHKERRQ(ierr); 1379 ierr = VecDuplicate(snes->vec_sol, &vi->z);CHKERRQ(ierr); 1380 ierr = VecDuplicate(snes->vec_sol, &vi->t); CHKERRQ(ierr); 1381 } 1382 PetscFunctionReturn(0); 1383 } 1384 /* -------------------------------------------------------------------------- */ 1385 /* 1386 SNESDestroy_VI - Destroys the private SNES_VI context that was created 1387 with SNESCreate_VI(). 1388 1389 Input Parameter: 1390 . snes - the SNES context 1391 1392 Application Interface Routine: SNESDestroy() 1393 */ 1394 #undef __FUNCT__ 1395 #define __FUNCT__ "SNESDestroy_VI" 1396 PetscErrorCode SNESDestroy_VI(SNES snes) 1397 { 1398 SNES_VI *vi = (SNES_VI*) snes->data; 1399 PetscErrorCode ierr; 1400 1401 PetscFunctionBegin; 1402 if (snes->vec_sol_update) { 1403 ierr = VecDestroy(snes->vec_sol_update);CHKERRQ(ierr); 1404 snes->vec_sol_update = PETSC_NULL; 1405 } 1406 if (snes->work) { 1407 ierr = VecDestroyVecs(snes->nwork,&snes->work);CHKERRQ(ierr); 1408 } 1409 if (snes->ops->solve != SNESSolveVI_SS) { 1410 ierr = ISDestroy(vi->IS_inact_prev); 1411 } 1412 1413 if (snes->ops->solve == SNESSolveVI_SS) { 1414 /* clear vectors */ 1415 ierr = VecDestroy(vi->dpsi);CHKERRQ(ierr); 1416 ierr = VecDestroy(vi->phi); CHKERRQ(ierr); 1417 ierr = VecDestroy(vi->Da); CHKERRQ(ierr); 1418 ierr = VecDestroy(vi->Db); CHKERRQ(ierr); 1419 ierr = VecDestroy(vi->z); CHKERRQ(ierr); 1420 ierr = VecDestroy(vi->t); CHKERRQ(ierr); 1421 } 1422 1423 if (!vi->usersetxbounds) { 1424 ierr = VecDestroy(vi->xl); CHKERRQ(ierr); 1425 ierr = VecDestroy(vi->xu); CHKERRQ(ierr); 1426 } 1427 1428 if (vi->lsmonitor) { 1429 ierr = PetscViewerASCIIMonitorDestroy(vi->lsmonitor);CHKERRQ(ierr); 1430 } 1431 ierr = PetscFree(snes->data);CHKERRQ(ierr); 1432 1433 /* clear composed functions */ 1434 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSet_C","",PETSC_NULL);CHKERRQ(ierr); 1435 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetMonitor_C","",PETSC_NULL);CHKERRQ(ierr); 1436 PetscFunctionReturn(0); 1437 } 1438 1439 /* -------------------------------------------------------------------------- */ 1440 #undef __FUNCT__ 1441 #define __FUNCT__ "SNESLineSearchNo_VI" 1442 1443 /* 1444 This routine does not actually do a line search but it takes a full newton 1445 step while ensuring that the new iterates remain within the constraints. 1446 1447 */ 1448 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) 1449 { 1450 PetscErrorCode ierr; 1451 SNES_VI *vi = (SNES_VI*)snes->data; 1452 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1453 1454 PetscFunctionBegin; 1455 *flag = PETSC_TRUE; 1456 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1457 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); /* ynorm = || y || */ 1458 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1459 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1460 if (vi->postcheckstep) { 1461 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1462 } 1463 if (changed_y) { 1464 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1465 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1466 } 1467 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1468 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1469 if (!snes->domainerror) { 1470 if (snes->ops->solve != SNESSolveVI_SS) { 1471 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1472 } else { 1473 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); /* gnorm = || g || */ 1474 } 1475 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1476 } 1477 if (vi->lsmonitor) { 1478 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1479 } 1480 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1481 PetscFunctionReturn(0); 1482 } 1483 1484 /* -------------------------------------------------------------------------- */ 1485 #undef __FUNCT__ 1486 #define __FUNCT__ "SNESLineSearchNoNorms_VI" 1487 1488 /* 1489 This routine is a copy of SNESLineSearchNoNorms in snes/impls/ls/ls.c 1490 */ 1491 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) 1492 { 1493 PetscErrorCode ierr; 1494 SNES_VI *vi = (SNES_VI*)snes->data; 1495 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1496 1497 PetscFunctionBegin; 1498 *flag = PETSC_TRUE; 1499 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1500 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1501 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1502 if (vi->postcheckstep) { 1503 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1504 } 1505 if (changed_y) { 1506 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1507 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1508 } 1509 1510 /* don't evaluate function the last time through */ 1511 if (snes->iter < snes->max_its-1) { 1512 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1513 } 1514 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1515 PetscFunctionReturn(0); 1516 } 1517 1518 /* -------------------------------------------------------------------------- */ 1519 #undef __FUNCT__ 1520 #define __FUNCT__ "SNESLineSearchCubic_VI" 1521 /* 1522 This routine implements a cubic line search while doing a projection on the variable bounds 1523 */ 1524 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) 1525 { 1526 PetscReal initslope,lambdaprev,gnormprev,a,b,d,t1,t2,rellength; 1527 PetscReal minlambda,lambda,lambdatemp; 1528 #if defined(PETSC_USE_COMPLEX) 1529 PetscScalar cinitslope; 1530 #endif 1531 PetscErrorCode ierr; 1532 PetscInt count; 1533 SNES_VI *vi = (SNES_VI*)snes->data; 1534 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1535 MPI_Comm comm; 1536 1537 PetscFunctionBegin; 1538 ierr = PetscObjectGetComm((PetscObject)snes,&comm);CHKERRQ(ierr); 1539 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1540 *flag = PETSC_TRUE; 1541 1542 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); 1543 if (*ynorm == 0.0) { 1544 if (vi->lsmonitor) { 1545 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Initial direction and size is 0\n");CHKERRQ(ierr); 1546 } 1547 *gnorm = fnorm; 1548 ierr = VecCopy(x,w);CHKERRQ(ierr); 1549 ierr = VecCopy(f,g);CHKERRQ(ierr); 1550 *flag = PETSC_FALSE; 1551 goto theend1; 1552 } 1553 if (*ynorm > vi->maxstep) { /* Step too big, so scale back */ 1554 if (vi->lsmonitor) { 1555 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Scaling step by %G old ynorm %G\n",vi->maxstep/(*ynorm),*ynorm);CHKERRQ(ierr); 1556 } 1557 ierr = VecScale(y,vi->maxstep/(*ynorm));CHKERRQ(ierr); 1558 *ynorm = vi->maxstep; 1559 } 1560 ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr); 1561 minlambda = vi->minlambda/rellength; 1562 ierr = MatMult(snes->jacobian,y,w);CHKERRQ(ierr); 1563 #if defined(PETSC_USE_COMPLEX) 1564 ierr = VecDot(f,w,&cinitslope);CHKERRQ(ierr); 1565 initslope = PetscRealPart(cinitslope); 1566 #else 1567 ierr = VecDot(f,w,&initslope);CHKERRQ(ierr); 1568 #endif 1569 if (initslope > 0.0) initslope = -initslope; 1570 if (initslope == 0.0) initslope = -1.0; 1571 1572 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1573 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1574 if (snes->nfuncs >= snes->max_funcs) { 1575 ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); 1576 *flag = PETSC_FALSE; 1577 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1578 goto theend1; 1579 } 1580 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1581 if (snes->ops->solve != SNESSolveVI_SS) { 1582 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1583 } else { 1584 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1585 } 1586 if (snes->domainerror) { 1587 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1588 PetscFunctionReturn(0); 1589 } 1590 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1591 ierr = PetscInfo2(snes,"Initial fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1592 if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + vi->alpha*initslope) { /* Sufficient reduction */ 1593 if (vi->lsmonitor) { 1594 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1595 } 1596 goto theend1; 1597 } 1598 1599 /* Fit points with quadratic */ 1600 lambda = 1.0; 1601 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 1602 lambdaprev = lambda; 1603 gnormprev = *gnorm; 1604 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1605 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1606 else lambda = lambdatemp; 1607 1608 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1609 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1610 if (snes->nfuncs >= snes->max_funcs) { 1611 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);CHKERRQ(ierr); 1612 *flag = PETSC_FALSE; 1613 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1614 goto theend1; 1615 } 1616 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1617 if (snes->ops->solve != SNESSolveVI_SS) { 1618 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1619 } else { 1620 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1621 } 1622 if (snes->domainerror) { 1623 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1624 PetscFunctionReturn(0); 1625 } 1626 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1627 if (vi->lsmonitor) { 1628 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: gnorm after quadratic fit %G\n",*gnorm);CHKERRQ(ierr); 1629 } 1630 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* sufficient reduction */ 1631 if (vi->lsmonitor) { 1632 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Quadratically determined step, lambda=%18.16e\n",lambda);CHKERRQ(ierr); 1633 } 1634 goto theend1; 1635 } 1636 1637 /* Fit points with cubic */ 1638 count = 1; 1639 while (PETSC_TRUE) { 1640 if (lambda <= minlambda) { 1641 if (vi->lsmonitor) { 1642 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: unable to find good step length! After %D tries \n",count);CHKERRQ(ierr); 1643 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); 1644 } 1645 *flag = PETSC_FALSE; 1646 break; 1647 } 1648 t1 = .5*((*gnorm)*(*gnorm) - fnorm*fnorm) - lambda*initslope; 1649 t2 = .5*(gnormprev*gnormprev - fnorm*fnorm) - lambdaprev*initslope; 1650 a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 1651 b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 1652 d = b*b - 3*a*initslope; 1653 if (d < 0.0) d = 0.0; 1654 if (a == 0.0) { 1655 lambdatemp = -initslope/(2.0*b); 1656 } else { 1657 lambdatemp = (-b + sqrt(d))/(3.0*a); 1658 } 1659 lambdaprev = lambda; 1660 gnormprev = *gnorm; 1661 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1662 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1663 else lambda = lambdatemp; 1664 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1665 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1666 if (snes->nfuncs >= snes->max_funcs) { 1667 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 1668 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); 1669 *flag = PETSC_FALSE; 1670 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1671 break; 1672 } 1673 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1674 if (snes->ops->solve != SNESSolveVI_SS) { 1675 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1676 } else { 1677 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1678 } 1679 if (snes->domainerror) { 1680 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1681 PetscFunctionReturn(0); 1682 } 1683 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1684 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* is reduction enough? */ 1685 if (vi->lsmonitor) { 1686 ierr = PetscPrintf(comm," Line search: Cubically determined step, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr); 1687 } 1688 break; 1689 } else { 1690 if (vi->lsmonitor) { 1691 ierr = PetscPrintf(comm," Line search: Cubic step no good, shrinking lambda, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr); 1692 } 1693 } 1694 count++; 1695 } 1696 theend1: 1697 /* Optional user-defined check for line search step validity */ 1698 if (vi->postcheckstep && *flag) { 1699 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1700 if (changed_y) { 1701 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1702 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1703 } 1704 if (changed_y || changed_w) { /* recompute the function if the step has changed */ 1705 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1706 if (snes->ops->solve != SNESSolveVI_SS) { 1707 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1708 } else { 1709 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1710 } 1711 if (snes->domainerror) { 1712 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1713 PetscFunctionReturn(0); 1714 } 1715 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1716 ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr); 1717 ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr); 1718 } 1719 } 1720 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1721 PetscFunctionReturn(0); 1722 } 1723 1724 #undef __FUNCT__ 1725 #define __FUNCT__ "SNESLineSearchQuadratic_VI" 1726 /* 1727 This routine does a quadratic line search while keeping the iterates within the variable bounds 1728 */ 1729 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) 1730 { 1731 /* 1732 Note that for line search purposes we work with with the related 1733 minimization problem: 1734 min z(x): R^n -> R, 1735 where z(x) = .5 * fnorm*fnorm,and fnorm = || f ||_2. 1736 */ 1737 PetscReal initslope,minlambda,lambda,lambdatemp,rellength; 1738 #if defined(PETSC_USE_COMPLEX) 1739 PetscScalar cinitslope; 1740 #endif 1741 PetscErrorCode ierr; 1742 PetscInt count; 1743 SNES_VI *vi = (SNES_VI*)snes->data; 1744 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1745 1746 PetscFunctionBegin; 1747 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1748 *flag = PETSC_TRUE; 1749 1750 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); 1751 if (*ynorm == 0.0) { 1752 if (vi->lsmonitor) { 1753 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor,"Line search: Direction and size is 0\n");CHKERRQ(ierr); 1754 } 1755 *gnorm = fnorm; 1756 ierr = VecCopy(x,w);CHKERRQ(ierr); 1757 ierr = VecCopy(f,g);CHKERRQ(ierr); 1758 *flag = PETSC_FALSE; 1759 goto theend2; 1760 } 1761 if (*ynorm > vi->maxstep) { /* Step too big, so scale back */ 1762 ierr = VecScale(y,vi->maxstep/(*ynorm));CHKERRQ(ierr); 1763 *ynorm = vi->maxstep; 1764 } 1765 ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr); 1766 minlambda = vi->minlambda/rellength; 1767 ierr = MatMult(snes->jacobian,y,w);CHKERRQ(ierr); 1768 #if defined(PETSC_USE_COMPLEX) 1769 ierr = VecDot(f,w,&cinitslope);CHKERRQ(ierr); 1770 initslope = PetscRealPart(cinitslope); 1771 #else 1772 ierr = VecDot(f,w,&initslope);CHKERRQ(ierr); 1773 #endif 1774 if (initslope > 0.0) initslope = -initslope; 1775 if (initslope == 0.0) initslope = -1.0; 1776 ierr = PetscInfo1(snes,"Initslope %G \n",initslope);CHKERRQ(ierr); 1777 1778 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1779 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1780 if (snes->nfuncs >= snes->max_funcs) { 1781 ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); 1782 *flag = PETSC_FALSE; 1783 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1784 goto theend2; 1785 } 1786 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1787 if (snes->ops->solve != SNESSolveVI_SS) { 1788 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1789 } else { 1790 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1791 } 1792 if (snes->domainerror) { 1793 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1794 PetscFunctionReturn(0); 1795 } 1796 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1797 if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + vi->alpha*initslope) { /* Sufficient reduction */ 1798 if (vi->lsmonitor) { 1799 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1800 } 1801 goto theend2; 1802 } 1803 1804 /* Fit points with quadratic */ 1805 lambda = 1.0; 1806 count = 1; 1807 while (PETSC_TRUE) { 1808 if (lambda <= minlambda) { /* bad luck; use full step */ 1809 if (vi->lsmonitor) { 1810 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor,"Line search: Unable to find good step length! %D \n",count);CHKERRQ(ierr); 1811 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); 1812 } 1813 ierr = VecCopy(x,w);CHKERRQ(ierr); 1814 *flag = PETSC_FALSE; 1815 break; 1816 } 1817 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 1818 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1819 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1820 else lambda = lambdatemp; 1821 1822 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1823 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1824 if (snes->nfuncs >= snes->max_funcs) { 1825 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 1826 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); 1827 *flag = PETSC_FALSE; 1828 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1829 break; 1830 } 1831 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1832 if (snes->domainerror) { 1833 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1834 PetscFunctionReturn(0); 1835 } 1836 if (snes->ops->solve != SNESSolveVI_SS) { 1837 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1838 } else { 1839 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1840 } 1841 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1842 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* sufficient reduction */ 1843 if (vi->lsmonitor) { 1844 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line Search: Quadratically determined step, lambda=%G\n",lambda);CHKERRQ(ierr); 1845 } 1846 break; 1847 } 1848 count++; 1849 } 1850 theend2: 1851 /* Optional user-defined check for line search step validity */ 1852 if (vi->postcheckstep) { 1853 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1854 if (changed_y) { 1855 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1856 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1857 } 1858 if (changed_y || changed_w) { /* recompute the function if the step has changed */ 1859 ierr = SNESComputeFunction(snes,w,g); 1860 if (snes->domainerror) { 1861 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1862 PetscFunctionReturn(0); 1863 } 1864 if (snes->ops->solve != SNESSolveVI_SS) { 1865 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1866 } else { 1867 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1868 } 1869 1870 ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr); 1871 ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr); 1872 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1873 } 1874 } 1875 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1876 PetscFunctionReturn(0); 1877 } 1878 1879 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*/ 1880 /* -------------------------------------------------------------------------- */ 1881 EXTERN_C_BEGIN 1882 #undef __FUNCT__ 1883 #define __FUNCT__ "SNESLineSearchSet_VI" 1884 PetscErrorCode SNESLineSearchSet_VI(SNES snes,FCN2 func,void *lsctx) 1885 { 1886 PetscFunctionBegin; 1887 ((SNES_VI *)(snes->data))->LineSearch = func; 1888 ((SNES_VI *)(snes->data))->lsP = lsctx; 1889 PetscFunctionReturn(0); 1890 } 1891 EXTERN_C_END 1892 1893 /* -------------------------------------------------------------------------- */ 1894 EXTERN_C_BEGIN 1895 #undef __FUNCT__ 1896 #define __FUNCT__ "SNESLineSearchSetMonitor_VI" 1897 PetscErrorCode SNESLineSearchSetMonitor_VI(SNES snes,PetscBool flg) 1898 { 1899 SNES_VI *vi = (SNES_VI*)snes->data; 1900 PetscErrorCode ierr; 1901 1902 PetscFunctionBegin; 1903 if (flg && !vi->lsmonitor) { 1904 ierr = PetscViewerASCIIMonitorCreate(((PetscObject)snes)->comm,"stdout",((PetscObject)snes)->tablevel,&vi->lsmonitor);CHKERRQ(ierr); 1905 } else if (!flg && vi->lsmonitor) { 1906 ierr = PetscViewerASCIIMonitorDestroy(vi->lsmonitor);CHKERRQ(ierr); 1907 } 1908 PetscFunctionReturn(0); 1909 } 1910 EXTERN_C_END 1911 1912 /* 1913 SNESView_VI - Prints info from the SNESVI data structure. 1914 1915 Input Parameters: 1916 . SNES - the SNES context 1917 . viewer - visualization context 1918 1919 Application Interface Routine: SNESView() 1920 */ 1921 #undef __FUNCT__ 1922 #define __FUNCT__ "SNESView_VI" 1923 static PetscErrorCode SNESView_VI(SNES snes,PetscViewer viewer) 1924 { 1925 SNES_VI *vi = (SNES_VI *)snes->data; 1926 const char *cstr,*tstr; 1927 PetscErrorCode ierr; 1928 PetscBool iascii; 1929 1930 PetscFunctionBegin; 1931 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1932 if (iascii) { 1933 if (vi->LineSearch == SNESLineSearchNo_VI) cstr = "SNESLineSearchNo"; 1934 else if (vi->LineSearch == SNESLineSearchQuadratic_VI) cstr = "SNESLineSearchQuadratic"; 1935 else if (vi->LineSearch == SNESLineSearchCubic_VI) cstr = "SNESLineSearchCubic"; 1936 else cstr = "unknown"; 1937 if (snes->ops->solve == SNESSolveVI_SS) tstr = "Semismooth"; 1938 else if (snes->ops->solve == SNESSolveVI_RS) tstr = "Reduced Space"; 1939 else if (snes->ops->solve == SNESSolveVI_RS2) tstr = "Augmented Space"; 1940 else tstr = "unknown"; 1941 ierr = PetscViewerASCIIPrintf(viewer," VI algorithm: %s\n",tstr);CHKERRQ(ierr); 1942 ierr = PetscViewerASCIIPrintf(viewer," line search variant: %s\n",cstr);CHKERRQ(ierr); 1943 ierr = PetscViewerASCIIPrintf(viewer," alpha=%G, maxstep=%G, minlambda=%G\n",vi->alpha,vi->maxstep,vi->minlambda);CHKERRQ(ierr); 1944 } else { 1945 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for SNES EQ VI",((PetscObject)viewer)->type_name); 1946 } 1947 PetscFunctionReturn(0); 1948 } 1949 1950 #undef __FUNCT__ 1951 #define __FUNCT__ "SNESVISetVariableBounds" 1952 /*@ 1953 SNESVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu. 1954 1955 Input Parameters: 1956 . snes - the SNES context. 1957 . xl - lower bound. 1958 . xu - upper bound. 1959 1960 Notes: 1961 If this routine is not called then the lower and upper bounds are set to 1962 PETSC_VI_INF and PETSC_VI_NINF respectively during SNESSetUp(). 1963 1964 @*/ 1965 PetscErrorCode SNESVISetVariableBounds(SNES snes, Vec xl, Vec xu) 1966 { 1967 SNES_VI *vi = (SNES_VI*)snes->data; 1968 1969 PetscFunctionBegin; 1970 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1971 PetscValidHeaderSpecific(xl,VEC_CLASSID,2); 1972 PetscValidHeaderSpecific(xu,VEC_CLASSID,3); 1973 if (!snes->vec_func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call SNESSetFunction() first"); 1974 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); 1975 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); 1976 1977 vi->usersetxbounds = PETSC_TRUE; 1978 vi->xl = xl; 1979 vi->xu = xu; 1980 PetscFunctionReturn(0); 1981 } 1982 1983 /* -------------------------------------------------------------------------- */ 1984 /* 1985 SNESSetFromOptions_VI - Sets various parameters for the SNESVI method. 1986 1987 Input Parameter: 1988 . snes - the SNES context 1989 1990 Application Interface Routine: SNESSetFromOptions() 1991 */ 1992 #undef __FUNCT__ 1993 #define __FUNCT__ "SNESSetFromOptions_VI" 1994 static PetscErrorCode SNESSetFromOptions_VI(SNES snes) 1995 { 1996 SNES_VI *vi = (SNES_VI *)snes->data; 1997 const char *lses[] = {"basic","basicnonorms","quadratic","cubic"}; 1998 const char *vies[] = {"ss","rs","rs2"}; 1999 PetscErrorCode ierr; 2000 PetscInt indx; 2001 PetscBool flg,set,flg2; 2002 2003 PetscFunctionBegin; 2004 ierr = PetscOptionsHead("SNES semismooth method options");CHKERRQ(ierr); 2005 ierr = PetscOptionsBool("-snes_vi_monitor","Monitor all non-active variables","None",PETSC_FALSE,&flg,0);CHKERRQ(ierr); 2006 if (flg) { 2007 ierr = SNESMonitorSet(snes,SNESMonitorVI,0,0);CHKERRQ(ierr); 2008 } 2009 ierr = PetscOptionsReal("-snes_ls_alpha","Function norm must decrease by","None",vi->alpha,&vi->alpha,0);CHKERRQ(ierr); 2010 ierr = PetscOptionsReal("-snes_ls_maxstep","Step must be less than","None",vi->maxstep,&vi->maxstep,0);CHKERRQ(ierr); 2011 ierr = PetscOptionsReal("-snes_ls_minlambda","Minimum lambda allowed","None",vi->minlambda,&vi->minlambda,0);CHKERRQ(ierr); 2012 ierr = PetscOptionsReal("-snes_vi_const_tol","constraint tolerance","None",vi->const_tol,&vi->const_tol,0);CHKERRQ(ierr); 2013 ierr = PetscOptionsBool("-snes_ls_monitor","Print progress of line searches","SNESLineSearchSetMonitor",vi->lsmonitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);CHKERRQ(ierr); 2014 if (set) {ierr = SNESLineSearchSetMonitor(snes,flg);CHKERRQ(ierr);} 2015 ierr = PetscOptionsEList("-snes_vi_type","Semismooth algorithm used","",vies,3,"ss",&indx,&flg2);CHKERRQ(ierr); 2016 if (flg2) { 2017 switch (indx) { 2018 case 0: 2019 snes->ops->solve = SNESSolveVI_SS; 2020 break; 2021 case 1: 2022 snes->ops->solve = SNESSolveVI_RS; 2023 break; 2024 case 2: 2025 snes->ops->solve = SNESSolveVI_RS2; 2026 } 2027 } 2028 ierr = PetscOptionsEList("-snes_ls","Line search used","SNESLineSearchSet",lses,4,"basic",&indx,&flg);CHKERRQ(ierr); 2029 if (flg) { 2030 switch (indx) { 2031 case 0: 2032 ierr = SNESLineSearchSet(snes,SNESLineSearchNo_VI,PETSC_NULL);CHKERRQ(ierr); 2033 break; 2034 case 1: 2035 ierr = SNESLineSearchSet(snes,SNESLineSearchNoNorms_VI,PETSC_NULL);CHKERRQ(ierr); 2036 break; 2037 case 2: 2038 ierr = SNESLineSearchSet(snes,SNESLineSearchQuadratic_VI,PETSC_NULL);CHKERRQ(ierr); 2039 break; 2040 case 3: 2041 ierr = SNESLineSearchSet(snes,SNESLineSearchCubic_VI,PETSC_NULL);CHKERRQ(ierr); 2042 break; 2043 } 2044 } 2045 ierr = PetscOptionsTail();CHKERRQ(ierr); 2046 PetscFunctionReturn(0); 2047 } 2048 /* -------------------------------------------------------------------------- */ 2049 /*MC 2050 SNESVI - Semismooth newton method based nonlinear solver that uses a line search 2051 2052 Options Database: 2053 + -snes_ls [cubic,quadratic,basic,basicnonorms] - Selects line search 2054 . -snes_ls_alpha <alpha> - Sets alpha 2055 . -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) 2056 . -snes_ls_minlambda <minlambda> - Sets the minimum lambda the line search will use minlambda / max_i ( y[i]/x[i] ) 2057 - -snes_ls_monitor - print information about progress of line searches 2058 2059 2060 Level: beginner 2061 2062 .seealso: SNESCreate(), SNES, SNESSetType(), SNESTR, SNESLineSearchSet(), 2063 SNESLineSearchSetPostCheck(), SNESLineSearchNo(), SNESLineSearchCubic(), SNESLineSearchQuadratic(), 2064 SNESLineSearchSet(), SNESLineSearchNoNorms(), SNESLineSearchSetPreCheck(), SNESLineSearchSetParams(), SNESLineSearchGetParams() 2065 2066 M*/ 2067 EXTERN_C_BEGIN 2068 #undef __FUNCT__ 2069 #define __FUNCT__ "SNESCreate_VI" 2070 PetscErrorCode SNESCreate_VI(SNES snes) 2071 { 2072 PetscErrorCode ierr; 2073 SNES_VI *vi; 2074 2075 PetscFunctionBegin; 2076 snes->ops->setup = SNESSetUp_VI; 2077 snes->ops->solve = SNESSolveVI_SS; 2078 snes->ops->destroy = SNESDestroy_VI; 2079 snes->ops->setfromoptions = SNESSetFromOptions_VI; 2080 snes->ops->view = SNESView_VI; 2081 snes->ops->converged = SNESDefaultConverged_VI; 2082 2083 ierr = PetscNewLog(snes,SNES_VI,&vi);CHKERRQ(ierr); 2084 snes->data = (void*)vi; 2085 vi->alpha = 1.e-4; 2086 vi->maxstep = 1.e8; 2087 vi->minlambda = 1.e-12; 2088 vi->LineSearch = SNESLineSearchNo_VI; 2089 vi->lsP = PETSC_NULL; 2090 vi->postcheckstep = PETSC_NULL; 2091 vi->postcheck = PETSC_NULL; 2092 vi->precheckstep = PETSC_NULL; 2093 vi->precheck = PETSC_NULL; 2094 vi->const_tol = 2.2204460492503131e-16; 2095 vi->checkredundancy = PETSC_NULL; 2096 2097 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetMonitor_C","SNESLineSearchSetMonitor_VI",SNESLineSearchSetMonitor_VI);CHKERRQ(ierr); 2098 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSet_C","SNESLineSearchSet_VI",SNESLineSearchSet_VI);CHKERRQ(ierr); 2099 2100 PetscFunctionReturn(0); 2101 } 2102 EXTERN_C_END 2103