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