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 /* Resets the snes PC and KSP when the active set sizes change */ 737 #undef __FUNCT__ 738 #define __FUNCT__ "SNESVIResetPCandKSP" 739 PetscErrorCode SNESVIResetPCandKSP(SNES snes,Mat Amat,Mat Pmat) 740 { 741 PetscErrorCode ierr; 742 KSP kspnew,snesksp; 743 PC pcnew; 744 const MatSolverPackage stype; 745 746 PetscFunctionBegin; 747 /* The active and inactive set sizes have changed so need to create a new snes->ksp object */ 748 ierr = SNESGetKSP(snes,&snesksp);CHKERRQ(ierr); 749 ierr = KSPCreate(((PetscObject)snes)->comm,&kspnew);CHKERRQ(ierr); 750 /* Copy over snes->ksp info */ 751 kspnew->pc_side = snesksp->pc_side; 752 kspnew->rtol = snesksp->rtol; 753 kspnew->abstol = snesksp->abstol; 754 kspnew->max_it = snesksp->max_it; 755 ierr = KSPSetType(kspnew,((PetscObject)snesksp)->type_name);CHKERRQ(ierr); 756 ierr = KSPGetPC(kspnew,&pcnew);CHKERRQ(ierr); 757 ierr = PCSetType(kspnew->pc,((PetscObject)snesksp->pc)->type_name);CHKERRQ(ierr); 758 ierr = PCSetOperators(kspnew->pc,Amat,Pmat,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 759 ierr = PCFactorGetMatSolverPackage(snesksp->pc,&stype);CHKERRQ(ierr); 760 ierr = PCFactorSetMatSolverPackage(kspnew->pc,stype);CHKERRQ(ierr); 761 ierr = KSPDestroy(snesksp);CHKERRQ(ierr); 762 snes->ksp = kspnew; 763 ierr = PetscLogObjectParent(snes,kspnew);CHKERRQ(ierr); 764 ierr = KSPSetFromOptions(kspnew);CHKERRQ(ierr); 765 PetscFunctionReturn(0); 766 } 767 /* Variational Inequality solver using active set method */ 768 #undef __FUNCT__ 769 #define __FUNCT__ "SNESSolveVI_AS" 770 PetscErrorCode SNESSolveVI_AS(SNES snes) 771 { 772 SNES_VI *vi = (SNES_VI*)snes->data; 773 PetscErrorCode ierr; 774 PetscInt maxits,i,lits,Nis_act=0; 775 PetscBool lssucceed; 776 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 777 PetscReal gnorm,xnorm=0,ynorm; 778 Vec Y,X,F,G,W; 779 KSPConvergedReason kspreason; 780 781 PetscFunctionBegin; 782 snes->numFailures = 0; 783 snes->numLinearSolveFailures = 0; 784 snes->reason = SNES_CONVERGED_ITERATING; 785 786 maxits = snes->max_its; /* maximum number of iterations */ 787 X = snes->vec_sol; /* solution vector */ 788 F = snes->vec_func; /* residual vector */ 789 Y = snes->work[0]; /* work vectors */ 790 G = snes->work[1]; 791 W = snes->work[2]; 792 793 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 794 snes->iter = 0; 795 snes->norm = 0.0; 796 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 797 798 ierr = SNESVIAdjustInitialGuess(X,vi->xl,vi->xu);CHKERRQ(ierr); 799 ierr = SNESComputeFunction(snes,X,vi->phi);CHKERRQ(ierr); 800 if (snes->domainerror) { 801 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 802 PetscFunctionReturn(0); 803 } 804 /* Compute Merit function */ 805 ierr = SNESVIComputeMeritFunction(vi->phi,&vi->merit,&vi->phinorm);CHKERRQ(ierr); 806 807 ierr = VecNormBegin(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm <- ||x|| */ 808 ierr = VecNormEnd(X,NORM_2,&xnorm);CHKERRQ(ierr); 809 if PetscIsInfOrNanReal(vi->merit) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 810 811 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 812 snes->norm = vi->phinorm; 813 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 814 SNESLogConvHistory(snes,vi->phinorm,0); 815 SNESMonitor(snes,0,vi->phinorm); 816 817 /* set parameter for default relative tolerance convergence test */ 818 snes->ttol = vi->phinorm*snes->rtol; 819 /* test convergence */ 820 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,vi->phinorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 821 if (snes->reason) PetscFunctionReturn(0); 822 823 for (i=0; i<maxits; i++) { 824 825 IS IS_act,IS_inact; /* _act -> active set _inact -> inactive set */ 826 PetscReal thresh,J_norm1; 827 VecScatter scat_act,scat_inact; 828 PetscInt nis_act,nis_inact,Nis_act_prev; 829 Vec Da_act,Da_inact,Db_inact; 830 Vec Y_act,Y_inact,phi_act,phi_inact; 831 Mat jac_inact_inact,jac_inact_act,prejac_inact_inact; 832 833 /* Call general purpose update function */ 834 if (snes->ops->update) { 835 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 836 } 837 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 838 /* Compute the threshold value for creating active and inactive sets */ 839 ierr = MatNorm(snes->jacobian,NORM_1,&J_norm1);CHKERRQ(ierr); 840 thresh = PetscMin(vi->merit,1e-2)/(1+J_norm1); 841 842 /* Compute B-subdifferential vectors Da and Db */ 843 ierr = SNESVIComputeBsubdifferentialVectors(snes,X,F,snes->jacobian,vi->Da,vi->Db);CHKERRQ(ierr); 844 845 /* Create active and inactive index sets */ 846 ierr = SNESVICreateIndexSets_AS(snes,vi->Db,thresh,&IS_act,&IS_inact);CHKERRQ(ierr); 847 848 Nis_act_prev = Nis_act; 849 /* Get sizes of active and inactive sets */ 850 ierr = ISGetLocalSize(IS_act,&nis_act);CHKERRQ(ierr); 851 ierr = ISGetLocalSize(IS_inact,&nis_inact);CHKERRQ(ierr); 852 ierr = ISGetSize(IS_act,&Nis_act);CHKERRQ(ierr); 853 854 ierr = PetscPrintf(PETSC_COMM_WORLD,"Size of active set = %d, size of inactive set = %d\n",nis_act,nis_inact);CHKERRQ(ierr); 855 856 /* Create active and inactive set vectors */ 857 ierr = SNESVICreateVectors_AS(snes,nis_act,&Da_act);CHKERRQ(ierr); 858 ierr = SNESVICreateVectors_AS(snes,nis_inact,&Da_inact);CHKERRQ(ierr); 859 ierr = SNESVICreateVectors_AS(snes,nis_inact,&Db_inact);CHKERRQ(ierr); 860 ierr = SNESVICreateVectors_AS(snes,nis_act,&phi_act);CHKERRQ(ierr); 861 ierr = SNESVICreateVectors_AS(snes,nis_inact,&phi_inact);CHKERRQ(ierr); 862 ierr = SNESVICreateVectors_AS(snes,nis_act,&Y_act);CHKERRQ(ierr); 863 ierr = SNESVICreateVectors_AS(snes,nis_inact,&Y_inact);CHKERRQ(ierr); 864 865 /* Create inactive set submatrices */ 866 ierr = MatGetSubMatrix(snes->jacobian,IS_inact,IS_act,MAT_INITIAL_MATRIX,&jac_inact_act);CHKERRQ(ierr); 867 ierr = MatGetSubMatrix(snes->jacobian,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);CHKERRQ(ierr); 868 869 /* Create scatter contexts */ 870 ierr = VecScatterCreate(vi->Da,IS_act,Da_act,PETSC_NULL,&scat_act);CHKERRQ(ierr); 871 ierr = VecScatterCreate(vi->Da,IS_inact,Da_inact,PETSC_NULL,&scat_inact);CHKERRQ(ierr); 872 873 /* Do a vec scatter to active and inactive set vectors */ 874 ierr = VecScatterBegin(scat_act,vi->Da,Da_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 875 ierr = VecScatterEnd(scat_act,vi->Da,Da_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 876 877 ierr = VecScatterBegin(scat_inact,vi->Da,Da_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 878 ierr = VecScatterEnd(scat_inact,vi->Da,Da_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 879 880 ierr = VecScatterBegin(scat_inact,vi->Db,Db_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 881 ierr = VecScatterEnd(scat_inact,vi->Db,Db_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 882 883 ierr = VecScatterBegin(scat_act,vi->phi,phi_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 884 ierr = VecScatterEnd(scat_act,vi->phi,phi_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 885 886 ierr = VecScatterBegin(scat_inact,vi->phi,phi_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 887 ierr = VecScatterEnd(scat_inact,vi->phi,phi_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 888 889 ierr = VecScatterBegin(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 890 ierr = VecScatterEnd(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 891 892 ierr = VecScatterBegin(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 893 ierr = VecScatterEnd(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 894 895 /* Active set direction */ 896 ierr = VecPointwiseDivide(Y_act,phi_act,Da_act);CHKERRQ(ierr); 897 /* inactive set jacobian and preconditioner */ 898 ierr = VecPointwiseDivide(Da_inact,Da_inact,Db_inact);CHKERRQ(ierr); 899 ierr = MatDiagonalSet(jac_inact_inact,Da_inact,ADD_VALUES);CHKERRQ(ierr); 900 if (snes->jacobian != snes->jacobian_pre) { 901 ierr = MatGetSubMatrix(snes->jacobian_pre,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&prejac_inact_inact);CHKERRQ(ierr); 902 ierr = MatDiagonalSet(prejac_inact_inact,Da_inact,ADD_VALUES);CHKERRQ(ierr); 903 } else prejac_inact_inact = jac_inact_inact; 904 905 /* right hand side */ 906 ierr = VecPointwiseDivide(phi_inact,phi_inact,Db_inact);CHKERRQ(ierr); 907 ierr = MatMult(jac_inact_act,Y_act,Db_inact);CHKERRQ(ierr); 908 ierr = VecAXPY(phi_inact,-1.0,Db_inact);CHKERRQ(ierr); 909 910 if ((i != 0) && (Nis_act != Nis_act_prev)) { 911 ierr = SNESVIResetPCandKSP(snes,jac_inact_inact,prejac_inact_inact);CHKERRQ(ierr); 912 } 913 914 ierr = KSPSetOperators(snes->ksp,jac_inact_inact,prejac_inact_inact,flg);CHKERRQ(ierr); 915 ierr = SNES_KSPSolve(snes,snes->ksp,phi_inact,Y_inact);CHKERRQ(ierr); 916 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 917 if (kspreason < 0) { 918 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 919 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 920 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 921 break; 922 } 923 } 924 925 ierr = VecScatterBegin(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 926 ierr = VecScatterEnd(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 927 ierr = VecScatterBegin(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 928 ierr = VecScatterEnd(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 929 930 ierr = VecDestroy(Da_act);CHKERRQ(ierr); 931 ierr = VecDestroy(Da_inact);CHKERRQ(ierr); 932 ierr = VecDestroy(Db_inact);CHKERRQ(ierr); 933 ierr = VecDestroy(phi_act);CHKERRQ(ierr); 934 ierr = VecDestroy(phi_inact);CHKERRQ(ierr); 935 ierr = VecDestroy(Y_act);CHKERRQ(ierr); 936 ierr = VecDestroy(Y_inact);CHKERRQ(ierr); 937 ierr = VecScatterDestroy(scat_act);CHKERRQ(ierr); 938 ierr = VecScatterDestroy(scat_inact);CHKERRQ(ierr); 939 ierr = ISDestroy(IS_act);CHKERRQ(ierr); 940 ierr = ISDestroy(IS_inact);CHKERRQ(ierr); 941 ierr = MatDestroy(jac_inact_act);CHKERRQ(ierr); 942 ierr = MatDestroy(jac_inact_inact);CHKERRQ(ierr); 943 if (snes->jacobian != snes->jacobian_pre) { 944 ierr = MatDestroy(prejac_inact_inact);CHKERRQ(ierr); 945 } 946 947 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 948 snes->linear_its += lits; 949 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 950 /* 951 if (vi->precheckstep) { 952 PetscBool changed_y = PETSC_FALSE; 953 ierr = (*vi->precheckstep)(snes,X,Y,vi->precheck,&changed_y);CHKERRQ(ierr); 954 } 955 956 if (PetscLogPrintInfo){ 957 ierr = SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr); 958 } 959 */ 960 /* Compute a (scaled) negative update in the line search routine: 961 Y <- X - lambda*Y 962 and evaluate G = function(Y) (depends on the line search). 963 */ 964 ierr = VecCopy(Y,snes->vec_sol_update);CHKERRQ(ierr); 965 ynorm = 1; gnorm = vi->phinorm; 966 ierr = SNESVIComputeJacobian(snes->jacobian,snes->jacobian_pre,vi->Da,vi->Db);CHKERRQ(ierr); 967 ierr = (*vi->LineSearch)(snes,vi->lsP,X,vi->phi,G,Y,W,vi->phinorm,xnorm,&ynorm,&gnorm,&lssucceed);CHKERRQ(ierr); 968 ierr = PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",vi->phinorm,gnorm,ynorm,(int)lssucceed);CHKERRQ(ierr); 969 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 970 if (snes->domainerror) { 971 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 972 PetscFunctionReturn(0); 973 } 974 if (!lssucceed) { 975 if (++snes->numFailures >= snes->maxFailures) { 976 PetscBool ismin; 977 snes->reason = SNES_DIVERGED_LINE_SEARCH; 978 ierr = SNESVICheckLocalMin_Private(snes,snes->jacobian,G,W,gnorm,&ismin);CHKERRQ(ierr); 979 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 980 break; 981 } 982 } 983 /* Update function and solution vectors */ 984 vi->phinorm = gnorm; 985 vi->merit = 0.5*vi->phinorm*vi->phinorm; 986 ierr = VecCopy(G,vi->phi);CHKERRQ(ierr); 987 ierr = VecCopy(W,X);CHKERRQ(ierr); 988 /* Monitor convergence */ 989 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 990 snes->iter = i+1; 991 snes->norm = vi->phinorm; 992 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 993 SNESLogConvHistory(snes,snes->norm,lits); 994 SNESMonitor(snes,snes->iter,snes->norm); 995 /* Test for convergence, xnorm = || X || */ 996 if (snes->ops->converged != SNESSkipConverged) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 997 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,vi->phinorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 998 if (snes->reason) break; 999 } 1000 if (i == maxits) { 1001 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 1002 if(!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 1003 } 1004 PetscFunctionReturn(0); 1005 } 1006 1007 /* Variational Inequality solver using active set method */ 1008 #undef __FUNCT__ 1009 #define __FUNCT__ "SNESSolveVI_RS" 1010 PetscErrorCode SNESSolveVI_RS(SNES snes) 1011 { 1012 SNES_VI *vi = (SNES_VI*)snes->data; 1013 PetscErrorCode ierr; 1014 PetscInt maxits,i,lits,Nis_act=0; 1015 PetscBool lssucceed; 1016 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 1017 PetscReal gnorm,xnorm=0,ynorm; 1018 Vec Y,X,F,G,W; 1019 KSPConvergedReason kspreason; 1020 1021 PetscFunctionBegin; 1022 snes->numFailures = 0; 1023 snes->numLinearSolveFailures = 0; 1024 snes->reason = SNES_CONVERGED_ITERATING; 1025 1026 maxits = snes->max_its; /* maximum number of iterations */ 1027 X = snes->vec_sol; /* solution vector */ 1028 F = snes->vec_func; /* residual vector */ 1029 Y = snes->work[0]; /* work vectors */ 1030 G = snes->work[1]; 1031 W = snes->work[2]; 1032 1033 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1034 snes->iter = 0; 1035 snes->norm = 0.0; 1036 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1037 1038 ierr = SNESVIAdjustInitialGuess(X,vi->xl,vi->xu);CHKERRQ(ierr); 1039 ierr = SNESComputeFunction(snes,X,vi->phi);CHKERRQ(ierr); 1040 if (snes->domainerror) { 1041 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 1042 PetscFunctionReturn(0); 1043 } 1044 /* Compute Merit function */ 1045 ierr = SNESVIComputeMeritFunction(vi->phi,&vi->merit,&vi->phinorm);CHKERRQ(ierr); 1046 1047 ierr = VecNormBegin(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm <- ||x|| */ 1048 ierr = VecNormEnd(X,NORM_2,&xnorm);CHKERRQ(ierr); 1049 if PetscIsInfOrNanReal(vi->merit) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1050 1051 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1052 snes->norm = vi->phinorm; 1053 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1054 SNESLogConvHistory(snes,vi->phinorm,0); 1055 SNESMonitor(snes,0,vi->phinorm); 1056 1057 /* set parameter for default relative tolerance convergence test */ 1058 snes->ttol = vi->phinorm*snes->rtol; 1059 /* test convergence */ 1060 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,vi->phinorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 1061 if (snes->reason) PetscFunctionReturn(0); 1062 1063 for (i=0; i<maxits; i++) { 1064 1065 IS IS_act,IS_inact; /* _act -> active set _inact -> inactive set */ 1066 VecScatter scat_act,scat_inact; 1067 PetscInt nis_act,nis_inact,Nis_act_prev; 1068 Vec Y_act,Y_inact,phi_inact; 1069 Mat jac_inact_inact,prejac_inact_inact; 1070 1071 /* Call general purpose update function */ 1072 if (snes->ops->update) { 1073 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 1074 } 1075 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 1076 ierr = SNESVIComputeBsubdifferentialVectors(snes,X,F,snes->jacobian,vi->Da,vi->Db);CHKERRQ(ierr); 1077 ierr = SNESVIComputeJacobian(snes->jacobian,snes->jacobian_pre,vi->Da,vi->Db);CHKERRQ(ierr); 1078 /* Create active and inactive index sets */ 1079 ierr = SNESVICreateIndexSets_RS(snes,X,vi->xl,vi->xu,&IS_act,&IS_inact);CHKERRQ(ierr); 1080 1081 Nis_act_prev = Nis_act; 1082 /* Get sizes of active and inactive sets */ 1083 ierr = ISGetLocalSize(IS_act,&nis_act);CHKERRQ(ierr); 1084 ierr = ISGetLocalSize(IS_inact,&nis_inact);CHKERRQ(ierr); 1085 ierr = ISGetSize(IS_act,&Nis_act);CHKERRQ(ierr); 1086 1087 ierr = PetscPrintf(PETSC_COMM_WORLD,"Size of active set = %d, size of inactive set = %d\n",nis_act,nis_inact);CHKERRQ(ierr); 1088 1089 /* Create active and inactive set vectors */ 1090 ierr = SNESVICreateVectors_AS(snes,nis_inact,&phi_inact);CHKERRQ(ierr); 1091 ierr = SNESVICreateVectors_AS(snes,nis_act,&Y_act);CHKERRQ(ierr); 1092 ierr = SNESVICreateVectors_AS(snes,nis_inact,&Y_inact);CHKERRQ(ierr); 1093 1094 /* Create inactive set submatrices */ 1095 ierr = MatGetSubMatrix(snes->jacobian,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);CHKERRQ(ierr); 1096 1097 /* Create scatter contexts */ 1098 ierr = VecScatterCreate(Y,IS_act,Y_act,PETSC_NULL,&scat_act);CHKERRQ(ierr); 1099 ierr = VecScatterCreate(Y,IS_inact,Y_inact,PETSC_NULL,&scat_inact);CHKERRQ(ierr); 1100 1101 /* Do a vec scatter to active and inactive set vectors */ 1102 ierr = VecScatterBegin(scat_inact,vi->phi,phi_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1103 ierr = VecScatterEnd(scat_inact,vi->phi,phi_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1104 1105 ierr = VecScatterBegin(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1106 ierr = VecScatterEnd(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1107 1108 ierr = VecScatterBegin(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1109 ierr = VecScatterEnd(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1110 1111 /* Active set direction = 0*/ 1112 ierr = VecSet(Y_act,0);CHKERRQ(ierr); 1113 if (snes->jacobian != snes->jacobian_pre) { 1114 ierr = MatGetSubMatrix(snes->jacobian_pre,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&prejac_inact_inact);CHKERRQ(ierr); 1115 } else prejac_inact_inact = jac_inact_inact; 1116 1117 if ((i != 0) && (Nis_act != Nis_act_prev)) { 1118 ierr = SNESVIResetPCandKSP(snes,jac_inact_inact,prejac_inact_inact);CHKERRQ(ierr); 1119 } 1120 1121 ierr = KSPSetOperators(snes->ksp,jac_inact_inact,prejac_inact_inact,flg);CHKERRQ(ierr); 1122 ierr = SNES_KSPSolve(snes,snes->ksp,phi_inact,Y_inact);CHKERRQ(ierr); 1123 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 1124 if (kspreason < 0) { 1125 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 1126 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 1127 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 1128 break; 1129 } 1130 } 1131 1132 ierr = VecScatterBegin(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1133 ierr = VecScatterEnd(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1134 ierr = VecScatterBegin(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1135 ierr = VecScatterEnd(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1136 1137 ierr = VecDestroy(phi_inact);CHKERRQ(ierr); 1138 ierr = VecDestroy(Y_act);CHKERRQ(ierr); 1139 ierr = VecDestroy(Y_inact);CHKERRQ(ierr); 1140 ierr = VecScatterDestroy(scat_act);CHKERRQ(ierr); 1141 ierr = VecScatterDestroy(scat_inact);CHKERRQ(ierr); 1142 ierr = ISDestroy(IS_act);CHKERRQ(ierr); 1143 ierr = ISDestroy(IS_inact);CHKERRQ(ierr); 1144 ierr = MatDestroy(jac_inact_inact);CHKERRQ(ierr); 1145 if (snes->jacobian != snes->jacobian_pre) { 1146 ierr = MatDestroy(prejac_inact_inact);CHKERRQ(ierr); 1147 } 1148 1149 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 1150 snes->linear_its += lits; 1151 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 1152 /* 1153 if (vi->precheckstep) { 1154 PetscBool changed_y = PETSC_FALSE; 1155 ierr = (*vi->precheckstep)(snes,X,Y,vi->precheck,&changed_y);CHKERRQ(ierr); 1156 } 1157 1158 if (PetscLogPrintInfo){ 1159 ierr = SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr); 1160 } 1161 */ 1162 /* Compute a (scaled) negative update in the line search routine: 1163 Y <- X - lambda*Y 1164 and evaluate G = function(Y) (depends on the line search). 1165 */ 1166 ierr = VecCopy(Y,snes->vec_sol_update);CHKERRQ(ierr); 1167 ynorm = 1; gnorm = vi->phinorm; 1168 ierr = (*vi->LineSearch)(snes,vi->lsP,X,vi->phi,G,Y,W,vi->phinorm,xnorm,&ynorm,&gnorm,&lssucceed);CHKERRQ(ierr); 1169 ierr = PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",vi->phinorm,gnorm,ynorm,(int)lssucceed);CHKERRQ(ierr); 1170 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 1171 if (snes->domainerror) { 1172 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 1173 PetscFunctionReturn(0); 1174 } 1175 if (!lssucceed) { 1176 if (++snes->numFailures >= snes->maxFailures) { 1177 PetscBool ismin; 1178 snes->reason = SNES_DIVERGED_LINE_SEARCH; 1179 ierr = SNESVICheckLocalMin_Private(snes,snes->jacobian,G,W,gnorm,&ismin);CHKERRQ(ierr); 1180 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 1181 break; 1182 } 1183 } 1184 /* Update function and solution vectors */ 1185 vi->phinorm = gnorm; 1186 vi->merit = 0.5*vi->phinorm*vi->phinorm; 1187 ierr = VecCopy(G,vi->phi);CHKERRQ(ierr); 1188 ierr = VecCopy(W,X);CHKERRQ(ierr); 1189 /* Monitor convergence */ 1190 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1191 snes->iter = i+1; 1192 snes->norm = vi->phinorm; 1193 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1194 SNESLogConvHistory(snes,snes->norm,lits); 1195 SNESMonitor(snes,snes->iter,snes->norm); 1196 /* Test for convergence, xnorm = || X || */ 1197 if (snes->ops->converged != SNESSkipConverged) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 1198 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,vi->phinorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 1199 if (snes->reason) break; 1200 } 1201 if (i == maxits) { 1202 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 1203 if(!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 1204 } 1205 PetscFunctionReturn(0); 1206 } 1207 1208 /* -------------------------------------------------------------------------- */ 1209 /* 1210 SNESSetUp_VI - Sets up the internal data structures for the later use 1211 of the SNESVI nonlinear solver. 1212 1213 Input Parameter: 1214 . snes - the SNES context 1215 . x - the solution vector 1216 1217 Application Interface Routine: SNESSetUp() 1218 1219 Notes: 1220 For basic use of the SNES solvers, the user need not explicitly call 1221 SNESSetUp(), since these actions will automatically occur during 1222 the call to SNESSolve(). 1223 */ 1224 #undef __FUNCT__ 1225 #define __FUNCT__ "SNESSetUp_VI" 1226 PetscErrorCode SNESSetUp_VI(SNES snes) 1227 { 1228 PetscErrorCode ierr; 1229 SNES_VI *vi = (SNES_VI*) snes->data; 1230 PetscInt i_start[3],i_end[3]; 1231 1232 PetscFunctionBegin; 1233 if (!snes->vec_sol_update) { 1234 ierr = VecDuplicate(snes->vec_sol,&snes->vec_sol_update);CHKERRQ(ierr); 1235 ierr = PetscLogObjectParent(snes,snes->vec_sol_update);CHKERRQ(ierr); 1236 } 1237 if (!snes->work) { 1238 snes->nwork = 3; 1239 ierr = VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work);CHKERRQ(ierr); 1240 ierr = PetscLogObjectParents(snes,snes->nwork,snes->work);CHKERRQ(ierr); 1241 } 1242 1243 ierr = VecDuplicate(snes->vec_sol, &vi->phi); CHKERRQ(ierr); 1244 ierr = VecDuplicate(snes->vec_sol, &vi->Da); CHKERRQ(ierr); 1245 ierr = VecDuplicate(snes->vec_sol, &vi->Db); CHKERRQ(ierr); 1246 ierr = VecDuplicate(snes->vec_sol, &vi->z);CHKERRQ(ierr); 1247 ierr = VecDuplicate(snes->vec_sol, &vi->t); CHKERRQ(ierr); 1248 1249 /* If the lower and upper bound on variables are not set, set it to 1250 -Inf and Inf */ 1251 if (!vi->xl && !vi->xu) { 1252 vi->usersetxbounds = PETSC_FALSE; 1253 ierr = VecDuplicate(snes->vec_sol, &vi->xl); CHKERRQ(ierr); 1254 ierr = VecSet(vi->xl,PETSC_VI_NINF);CHKERRQ(ierr); 1255 ierr = VecDuplicate(snes->vec_sol, &vi->xu); CHKERRQ(ierr); 1256 ierr = VecSet(vi->xu,PETSC_VI_INF);CHKERRQ(ierr); 1257 } else { 1258 /* Check if lower bound, upper bound and solution vector distribution across the processors is identical */ 1259 ierr = VecGetOwnershipRange(snes->vec_sol,i_start,i_end);CHKERRQ(ierr); 1260 ierr = VecGetOwnershipRange(vi->xl,i_start+1,i_end+1);CHKERRQ(ierr); 1261 ierr = VecGetOwnershipRange(vi->xu,i_start+2,i_end+2);CHKERRQ(ierr); 1262 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])) 1263 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."); 1264 } 1265 1266 vi->computeuserfunction = snes->ops->computefunction; 1267 snes->ops->computefunction = SNESVIComputeFunction; 1268 1269 PetscFunctionReturn(0); 1270 } 1271 /* -------------------------------------------------------------------------- */ 1272 /* 1273 SNESDestroy_VI - Destroys the private SNES_VI context that was created 1274 with SNESCreate_VI(). 1275 1276 Input Parameter: 1277 . snes - the SNES context 1278 1279 Application Interface Routine: SNESDestroy() 1280 */ 1281 #undef __FUNCT__ 1282 #define __FUNCT__ "SNESDestroy_VI" 1283 PetscErrorCode SNESDestroy_VI(SNES snes) 1284 { 1285 SNES_VI *vi = (SNES_VI*) snes->data; 1286 PetscErrorCode ierr; 1287 1288 PetscFunctionBegin; 1289 if (snes->vec_sol_update) { 1290 ierr = VecDestroy(snes->vec_sol_update);CHKERRQ(ierr); 1291 snes->vec_sol_update = PETSC_NULL; 1292 } 1293 if (snes->nwork) { 1294 ierr = VecDestroyVecs(snes->work,snes->nwork);CHKERRQ(ierr); 1295 snes->nwork = 0; 1296 snes->work = PETSC_NULL; 1297 } 1298 1299 /* clear vectors */ 1300 ierr = VecDestroy(vi->phi); CHKERRQ(ierr); 1301 ierr = VecDestroy(vi->Da); CHKERRQ(ierr); 1302 ierr = VecDestroy(vi->Db); CHKERRQ(ierr); 1303 ierr = VecDestroy(vi->z); CHKERRQ(ierr); 1304 ierr = VecDestroy(vi->t); CHKERRQ(ierr); 1305 if (!vi->usersetxbounds) { 1306 ierr = VecDestroy(vi->xl); CHKERRQ(ierr); 1307 ierr = VecDestroy(vi->xu); CHKERRQ(ierr); 1308 } 1309 if (vi->lsmonitor) { 1310 ierr = PetscViewerASCIIMonitorDestroy(vi->lsmonitor);CHKERRQ(ierr); 1311 } 1312 ierr = PetscFree(snes->data);CHKERRQ(ierr); 1313 1314 /* clear composed functions */ 1315 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSet_C","",PETSC_NULL);CHKERRQ(ierr); 1316 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetMonitor_C","",PETSC_NULL);CHKERRQ(ierr); 1317 PetscFunctionReturn(0); 1318 } 1319 /* -------------------------------------------------------------------------- */ 1320 #undef __FUNCT__ 1321 #define __FUNCT__ "SNESLineSearchNo_VI" 1322 1323 /* 1324 This routine is a copy of SNESLineSearchNo routine in snes/impls/ls/ls.c 1325 1326 */ 1327 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) 1328 { 1329 PetscErrorCode ierr; 1330 SNES_VI *vi = (SNES_VI*)snes->data; 1331 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1332 1333 PetscFunctionBegin; 1334 *flag = PETSC_TRUE; 1335 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1336 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); /* ynorm = || y || */ 1337 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1338 if (vi->postcheckstep) { 1339 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1340 } 1341 if (changed_y) { 1342 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1343 } 1344 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1345 if (!snes->domainerror) { 1346 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); /* gnorm = || g || */ 1347 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1348 } 1349 if (vi->lsmonitor) { 1350 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1351 } 1352 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1353 PetscFunctionReturn(0); 1354 } 1355 1356 /* -------------------------------------------------------------------------- */ 1357 #undef __FUNCT__ 1358 #define __FUNCT__ "SNESLineSearchNoNorms_VI" 1359 1360 /* 1361 This routine is a copy of SNESLineSearchNoNorms in snes/impls/ls/ls.c 1362 */ 1363 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) 1364 { 1365 PetscErrorCode ierr; 1366 SNES_VI *vi = (SNES_VI*)snes->data; 1367 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1368 1369 PetscFunctionBegin; 1370 *flag = PETSC_TRUE; 1371 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1372 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1373 if (vi->postcheckstep) { 1374 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1375 } 1376 if (changed_y) { 1377 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1378 } 1379 1380 /* don't evaluate function the last time through */ 1381 if (snes->iter < snes->max_its-1) { 1382 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1383 } 1384 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1385 PetscFunctionReturn(0); 1386 } 1387 /* -------------------------------------------------------------------------- */ 1388 #undef __FUNCT__ 1389 #define __FUNCT__ "SNESLineSearchCubic_VI" 1390 /* 1391 This routine is a copy of SNESLineSearchCubic in snes/impls/ls/ls.c 1392 */ 1393 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) 1394 { 1395 /* 1396 Note that for line search purposes we work with with the related 1397 minimization problem: 1398 min z(x): R^n -> R, 1399 where z(x) = .5 * fnorm*fnorm, and fnorm = || f ||_2. 1400 This function z(x) is same as the merit function for the complementarity problem. 1401 */ 1402 1403 PetscReal initslope,lambdaprev,gnormprev,a,b,d,t1,t2,rellength; 1404 PetscReal minlambda,lambda,lambdatemp; 1405 #if defined(PETSC_USE_COMPLEX) 1406 PetscScalar cinitslope; 1407 #endif 1408 PetscErrorCode ierr; 1409 PetscInt count; 1410 SNES_VI *vi = (SNES_VI*)snes->data; 1411 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1412 MPI_Comm comm; 1413 1414 PetscFunctionBegin; 1415 ierr = PetscObjectGetComm((PetscObject)snes,&comm);CHKERRQ(ierr); 1416 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1417 *flag = PETSC_TRUE; 1418 1419 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); 1420 if (*ynorm == 0.0) { 1421 if (vi->lsmonitor) { 1422 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Initial direction and size is 0\n");CHKERRQ(ierr); 1423 } 1424 *gnorm = fnorm; 1425 ierr = VecCopy(x,w);CHKERRQ(ierr); 1426 ierr = VecCopy(f,g);CHKERRQ(ierr); 1427 *flag = PETSC_FALSE; 1428 goto theend1; 1429 } 1430 if (*ynorm > vi->maxstep) { /* Step too big, so scale back */ 1431 if (vi->lsmonitor) { 1432 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Scaling step by %G old ynorm %G\n",vi->maxstep/(*ynorm),*ynorm);CHKERRQ(ierr); 1433 } 1434 ierr = VecScale(y,vi->maxstep/(*ynorm));CHKERRQ(ierr); 1435 *ynorm = vi->maxstep; 1436 } 1437 ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr); 1438 minlambda = vi->minlambda/rellength; 1439 ierr = MatMult(snes->jacobian,y,w);CHKERRQ(ierr); 1440 #if defined(PETSC_USE_COMPLEX) 1441 ierr = VecDot(f,w,&cinitslope);CHKERRQ(ierr); 1442 initslope = PetscRealPart(cinitslope); 1443 #else 1444 ierr = VecDot(f,w,&initslope);CHKERRQ(ierr); 1445 #endif 1446 if (initslope > 0.0) initslope = -initslope; 1447 if (initslope == 0.0) initslope = -1.0; 1448 1449 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1450 if (snes->nfuncs >= snes->max_funcs) { 1451 ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); 1452 *flag = PETSC_FALSE; 1453 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1454 goto theend1; 1455 } 1456 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1457 if (snes->domainerror) { 1458 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1459 PetscFunctionReturn(0); 1460 } 1461 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1462 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1463 ierr = PetscInfo2(snes,"Initial fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1464 if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + vi->alpha*initslope) { /* Sufficient reduction */ 1465 if (vi->lsmonitor) { 1466 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1467 } 1468 goto theend1; 1469 } 1470 1471 /* Fit points with quadratic */ 1472 lambda = 1.0; 1473 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 1474 lambdaprev = lambda; 1475 gnormprev = *gnorm; 1476 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1477 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1478 else lambda = lambdatemp; 1479 1480 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1481 if (snes->nfuncs >= snes->max_funcs) { 1482 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);CHKERRQ(ierr); 1483 *flag = PETSC_FALSE; 1484 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1485 goto theend1; 1486 } 1487 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1488 if (snes->domainerror) { 1489 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1490 PetscFunctionReturn(0); 1491 } 1492 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1493 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1494 if (vi->lsmonitor) { 1495 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: gnorm after quadratic fit %G\n",*gnorm);CHKERRQ(ierr); 1496 } 1497 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* sufficient reduction */ 1498 if (vi->lsmonitor) { 1499 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Quadratically determined step, lambda=%18.16e\n",lambda);CHKERRQ(ierr); 1500 } 1501 goto theend1; 1502 } 1503 1504 /* Fit points with cubic */ 1505 count = 1; 1506 while (PETSC_TRUE) { 1507 if (lambda <= minlambda) { 1508 if (vi->lsmonitor) { 1509 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: unable to find good step length! After %D tries \n",count);CHKERRQ(ierr); 1510 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); 1511 } 1512 *flag = PETSC_FALSE; 1513 break; 1514 } 1515 t1 = .5*((*gnorm)*(*gnorm) - fnorm*fnorm) - lambda*initslope; 1516 t2 = .5*(gnormprev*gnormprev - fnorm*fnorm) - lambdaprev*initslope; 1517 a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 1518 b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 1519 d = b*b - 3*a*initslope; 1520 if (d < 0.0) d = 0.0; 1521 if (a == 0.0) { 1522 lambdatemp = -initslope/(2.0*b); 1523 } else { 1524 lambdatemp = (-b + sqrt(d))/(3.0*a); 1525 } 1526 lambdaprev = lambda; 1527 gnormprev = *gnorm; 1528 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1529 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1530 else lambda = lambdatemp; 1531 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1532 if (snes->nfuncs >= snes->max_funcs) { 1533 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 1534 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); 1535 *flag = PETSC_FALSE; 1536 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1537 break; 1538 } 1539 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1540 if (snes->domainerror) { 1541 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1542 PetscFunctionReturn(0); 1543 } 1544 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1545 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1546 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* is reduction enough? */ 1547 if (vi->lsmonitor) { 1548 ierr = PetscPrintf(comm," Line search: Cubically determined step, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr); 1549 } 1550 break; 1551 } else { 1552 if (vi->lsmonitor) { 1553 ierr = PetscPrintf(comm," Line search: Cubic step no good, shrinking lambda, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr); 1554 } 1555 } 1556 count++; 1557 } 1558 theend1: 1559 /* Optional user-defined check for line search step validity */ 1560 if (vi->postcheckstep && *flag) { 1561 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1562 if (changed_y) { 1563 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1564 } 1565 if (changed_y || changed_w) { /* recompute the function if the step has changed */ 1566 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1567 if (snes->domainerror) { 1568 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1569 PetscFunctionReturn(0); 1570 } 1571 ierr = VecNormBegin(g,NORM_2,gnorm);CHKERRQ(ierr); 1572 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1573 ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr); 1574 ierr = VecNormEnd(g,NORM_2,gnorm);CHKERRQ(ierr); 1575 ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr); 1576 } 1577 } 1578 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1579 PetscFunctionReturn(0); 1580 } 1581 /* -------------------------------------------------------------------------- */ 1582 #undef __FUNCT__ 1583 #define __FUNCT__ "SNESLineSearchQuadratic_VI" 1584 /* 1585 This routine is a copy of SNESLineSearchQuadratic in snes/impls/ls/ls.c 1586 */ 1587 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) 1588 { 1589 /* 1590 Note that for line search purposes we work with with the related 1591 minimization problem: 1592 min z(x): R^n -> R, 1593 where z(x) = .5 * fnorm*fnorm,and fnorm = || f ||_2. 1594 z(x) is the same as the merit function for the complementarity problem 1595 */ 1596 PetscReal initslope,minlambda,lambda,lambdatemp,rellength; 1597 #if defined(PETSC_USE_COMPLEX) 1598 PetscScalar cinitslope; 1599 #endif 1600 PetscErrorCode ierr; 1601 PetscInt count; 1602 SNES_VI *vi = (SNES_VI*)snes->data; 1603 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1604 1605 PetscFunctionBegin; 1606 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1607 *flag = PETSC_TRUE; 1608 1609 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); 1610 if (*ynorm == 0.0) { 1611 if (vi->lsmonitor) { 1612 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor,"Line search: Direction and size is 0\n");CHKERRQ(ierr); 1613 } 1614 *gnorm = fnorm; 1615 ierr = VecCopy(x,w);CHKERRQ(ierr); 1616 ierr = VecCopy(f,g);CHKERRQ(ierr); 1617 *flag = PETSC_FALSE; 1618 goto theend2; 1619 } 1620 if (*ynorm > vi->maxstep) { /* Step too big, so scale back */ 1621 ierr = VecScale(y,vi->maxstep/(*ynorm));CHKERRQ(ierr); 1622 *ynorm = vi->maxstep; 1623 } 1624 ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr); 1625 minlambda = vi->minlambda/rellength; 1626 ierr = MatMult(snes->jacobian,y,w);CHKERRQ(ierr); 1627 #if defined(PETSC_USE_COMPLEX) 1628 ierr = VecDot(f,w,&cinitslope);CHKERRQ(ierr); 1629 initslope = PetscRealPart(cinitslope); 1630 #else 1631 ierr = VecDot(f,w,&initslope);CHKERRQ(ierr); 1632 #endif 1633 if (initslope > 0.0) initslope = -initslope; 1634 if (initslope == 0.0) initslope = -1.0; 1635 ierr = PetscInfo1(snes,"Initslope %G \n",initslope);CHKERRQ(ierr); 1636 1637 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1638 if (snes->nfuncs >= snes->max_funcs) { 1639 ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); 1640 *flag = PETSC_FALSE; 1641 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1642 goto theend2; 1643 } 1644 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1645 if (snes->domainerror) { 1646 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1647 PetscFunctionReturn(0); 1648 } 1649 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1650 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1651 if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + vi->alpha*initslope) { /* Sufficient reduction */ 1652 if (vi->lsmonitor) { 1653 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1654 } 1655 goto theend2; 1656 } 1657 1658 /* Fit points with quadratic */ 1659 lambda = 1.0; 1660 count = 1; 1661 while (PETSC_TRUE) { 1662 if (lambda <= minlambda) { /* bad luck; use full step */ 1663 if (vi->lsmonitor) { 1664 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor,"Line search: Unable to find good step length! %D \n",count);CHKERRQ(ierr); 1665 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); 1666 } 1667 ierr = VecCopy(x,w);CHKERRQ(ierr); 1668 *flag = PETSC_FALSE; 1669 break; 1670 } 1671 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 1672 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1673 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1674 else lambda = lambdatemp; 1675 1676 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1677 if (snes->nfuncs >= snes->max_funcs) { 1678 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 1679 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); 1680 *flag = PETSC_FALSE; 1681 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1682 break; 1683 } 1684 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1685 if (snes->domainerror) { 1686 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1687 PetscFunctionReturn(0); 1688 } 1689 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1690 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1691 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* sufficient reduction */ 1692 if (vi->lsmonitor) { 1693 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line Search: Quadratically determined step, lambda=%G\n",lambda);CHKERRQ(ierr); 1694 } 1695 break; 1696 } 1697 count++; 1698 } 1699 theend2: 1700 /* Optional user-defined check for line search step validity */ 1701 if (vi->postcheckstep) { 1702 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1703 if (changed_y) { 1704 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1705 } 1706 if (changed_y || changed_w) { /* recompute the function if the step has changed */ 1707 ierr = SNESComputeFunction(snes,w,g); 1708 if (snes->domainerror) { 1709 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1710 PetscFunctionReturn(0); 1711 } 1712 ierr = VecNormBegin(g,NORM_2,gnorm);CHKERRQ(ierr); 1713 ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr); 1714 ierr = VecNormEnd(g,NORM_2,gnorm);CHKERRQ(ierr); 1715 ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr); 1716 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1717 } 1718 } 1719 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1720 PetscFunctionReturn(0); 1721 } 1722 1723 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*/ 1724 /* -------------------------------------------------------------------------- */ 1725 EXTERN_C_BEGIN 1726 #undef __FUNCT__ 1727 #define __FUNCT__ "SNESLineSearchSet_VI" 1728 PetscErrorCode PETSCSNES_DLLEXPORT SNESLineSearchSet_VI(SNES snes,FCN2 func,void *lsctx) 1729 { 1730 PetscFunctionBegin; 1731 ((SNES_VI *)(snes->data))->LineSearch = func; 1732 ((SNES_VI *)(snes->data))->lsP = lsctx; 1733 PetscFunctionReturn(0); 1734 } 1735 EXTERN_C_END 1736 1737 /* -------------------------------------------------------------------------- */ 1738 EXTERN_C_BEGIN 1739 #undef __FUNCT__ 1740 #define __FUNCT__ "SNESLineSearchSetMonitor_VI" 1741 PetscErrorCode PETSCSNES_DLLEXPORT SNESLineSearchSetMonitor_VI(SNES snes,PetscBool flg) 1742 { 1743 SNES_VI *vi = (SNES_VI*)snes->data; 1744 PetscErrorCode ierr; 1745 1746 PetscFunctionBegin; 1747 if (flg && !vi->lsmonitor) { 1748 ierr = PetscViewerASCIIMonitorCreate(((PetscObject)snes)->comm,"stdout",((PetscObject)snes)->tablevel,&vi->lsmonitor);CHKERRQ(ierr); 1749 } else if (!flg && vi->lsmonitor) { 1750 ierr = PetscViewerASCIIMonitorDestroy(vi->lsmonitor);CHKERRQ(ierr); 1751 } 1752 PetscFunctionReturn(0); 1753 } 1754 EXTERN_C_END 1755 1756 /* 1757 SNESView_VI - Prints info from the SNESVI data structure. 1758 1759 Input Parameters: 1760 . SNES - the SNES context 1761 . viewer - visualization context 1762 1763 Application Interface Routine: SNESView() 1764 */ 1765 #undef __FUNCT__ 1766 #define __FUNCT__ "SNESView_VI" 1767 static PetscErrorCode SNESView_VI(SNES snes,PetscViewer viewer) 1768 { 1769 SNES_VI *vi = (SNES_VI *)snes->data; 1770 const char *cstr,*tstr; 1771 PetscErrorCode ierr; 1772 PetscBool iascii; 1773 1774 PetscFunctionBegin; 1775 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1776 if (iascii) { 1777 if (vi->LineSearch == SNESLineSearchNo_VI) cstr = "SNESLineSearchNo"; 1778 else if (vi->LineSearch == SNESLineSearchQuadratic_VI) cstr = "SNESLineSearchQuadratic"; 1779 else if (vi->LineSearch == SNESLineSearchCubic_VI) cstr = "SNESLineSearchCubic"; 1780 else cstr = "unknown"; 1781 if (snes->ops->solve == SNESSolveVI_SS) tstr = "Semismooth"; 1782 else if (snes->ops->solve == SNESSolveVI_AS) tstr = "Active Set"; 1783 else if (snes->ops->solve == SNESSolveVI_RS) tstr = "Reduced Space"; 1784 else tstr = "unknown"; 1785 ierr = PetscViewerASCIIPrintf(viewer," VI algorithm: %s\n",tstr);CHKERRQ(ierr); 1786 ierr = PetscViewerASCIIPrintf(viewer," line search variant: %s\n",cstr);CHKERRQ(ierr); 1787 ierr = PetscViewerASCIIPrintf(viewer," alpha=%G, maxstep=%G, minlambda=%G\n",vi->alpha,vi->maxstep,vi->minlambda);CHKERRQ(ierr); 1788 } else { 1789 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for SNES EQ VI",((PetscObject)viewer)->type_name); 1790 } 1791 PetscFunctionReturn(0); 1792 } 1793 1794 /* 1795 SNESVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu. 1796 1797 Input Parameters: 1798 . snes - the SNES context. 1799 . xl - lower bound. 1800 . xu - upper bound. 1801 1802 Notes: 1803 If this routine is not called then the lower and upper bounds are set to 1804 PETSC_VI_INF and PETSC_VI_NINF respectively during SNESSetUp(). 1805 1806 */ 1807 1808 #undef __FUNCT__ 1809 #define __FUNCT__ "SNESVISetVariableBounds" 1810 PetscErrorCode SNESVISetVariableBounds(SNES snes, Vec xl, Vec xu) 1811 { 1812 SNES_VI *vi = (SNES_VI*)snes->data; 1813 1814 PetscFunctionBegin; 1815 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1816 PetscValidHeaderSpecific(xl,VEC_CLASSID,2); 1817 PetscValidHeaderSpecific(xu,VEC_CLASSID,3); 1818 1819 /* Check if SNESSetFunction is called */ 1820 if(!snes->vec_func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call SNESSetFunction() first"); 1821 1822 /* Check if the vector sizes are compatible for lower and upper bounds */ 1823 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); 1824 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); 1825 vi->usersetxbounds = PETSC_TRUE; 1826 vi->xl = xl; 1827 vi->xu = xu; 1828 1829 PetscFunctionReturn(0); 1830 } 1831 /* -------------------------------------------------------------------------- */ 1832 /* 1833 SNESSetFromOptions_VI - Sets various parameters for the SNESVI method. 1834 1835 Input Parameter: 1836 . snes - the SNES context 1837 1838 Application Interface Routine: SNESSetFromOptions() 1839 */ 1840 #undef __FUNCT__ 1841 #define __FUNCT__ "SNESSetFromOptions_VI" 1842 static PetscErrorCode SNESSetFromOptions_VI(SNES snes) 1843 { 1844 SNES_VI *vi = (SNES_VI *)snes->data; 1845 const char *lses[] = {"basic","basicnonorms","quadratic","cubic"}; 1846 const char *vies[] = {"ss","as","rs"}; 1847 PetscErrorCode ierr; 1848 PetscInt indx; 1849 PetscBool flg,set,flg2; 1850 1851 PetscFunctionBegin; 1852 ierr = PetscOptionsHead("SNES semismooth method options");CHKERRQ(ierr); 1853 ierr = PetscOptionsReal("-snes_vi_alpha","Function norm must decrease by","None",vi->alpha,&vi->alpha,0);CHKERRQ(ierr); 1854 ierr = PetscOptionsReal("-snes_vi_maxstep","Step must be less than","None",vi->maxstep,&vi->maxstep,0);CHKERRQ(ierr); 1855 ierr = PetscOptionsReal("-snes_vi_minlambda","Minimum lambda allowed","None",vi->minlambda,&vi->minlambda,0);CHKERRQ(ierr); 1856 ierr = PetscOptionsReal("-snes_vi_const_tol","constraint tolerance","None",vi->const_tol,&vi->const_tol,0);CHKERRQ(ierr); 1857 ierr = PetscOptionsBool("-snes_vi_lsmonitor","Print progress of line searches","SNESLineSearchSetMonitor",vi->lsmonitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);CHKERRQ(ierr); 1858 if (set) {ierr = SNESLineSearchSetMonitor(snes,flg);CHKERRQ(ierr);} 1859 ierr = PetscOptionsEList("-snes_vi_type","Semismooth algorithm used","",vies,3,"ss",&indx,&flg2);CHKERRQ(ierr); 1860 if (flg2) { 1861 switch (indx) { 1862 case 0: 1863 snes->ops->solve = SNESSolveVI_SS; 1864 break; 1865 case 1: 1866 snes->ops->solve = SNESSolveVI_AS; 1867 break; 1868 case 2: 1869 snes->ops->solve = SNESSolveVI_RS; 1870 break; 1871 } 1872 } 1873 ierr = PetscOptionsEList("-snes_vi_ls","Line search used","SNESLineSearchSet",lses,4,"cubic",&indx,&flg);CHKERRQ(ierr); 1874 if (flg) { 1875 switch (indx) { 1876 case 0: 1877 ierr = SNESLineSearchSet(snes,SNESLineSearchNo_VI,PETSC_NULL);CHKERRQ(ierr); 1878 break; 1879 case 1: 1880 ierr = SNESLineSearchSet(snes,SNESLineSearchNoNorms_VI,PETSC_NULL);CHKERRQ(ierr); 1881 break; 1882 case 2: 1883 ierr = SNESLineSearchSet(snes,SNESLineSearchQuadratic_VI,PETSC_NULL);CHKERRQ(ierr); 1884 break; 1885 case 3: 1886 ierr = SNESLineSearchSet(snes,SNESLineSearchCubic_VI,PETSC_NULL);CHKERRQ(ierr); 1887 break; 1888 } 1889 } 1890 ierr = PetscOptionsTail();CHKERRQ(ierr); 1891 PetscFunctionReturn(0); 1892 } 1893 /* -------------------------------------------------------------------------- */ 1894 /*MC 1895 SNESVI - Semismooth newton method based nonlinear solver that uses a line search 1896 1897 Options Database: 1898 + -snes_vi [cubic,quadratic,basic,basicnonorms] - Selects line search 1899 . -snes_vi_alpha <alpha> - Sets alpha 1900 . -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) 1901 . -snes_vi_minlambda <minlambda> - Sets the minimum lambda the line search will use minlambda / max_i ( y[i]/x[i] ) 1902 - -snes_vi_monitor - print information about progress of line searches 1903 1904 1905 Level: beginner 1906 1907 .seealso: SNESCreate(), SNES, SNESSetType(), SNESTR, SNESLineSearchSet(), 1908 SNESLineSearchSetPostCheck(), SNESLineSearchNo(), SNESLineSearchCubic(), SNESLineSearchQuadratic(), 1909 SNESLineSearchSet(), SNESLineSearchNoNorms(), SNESLineSearchSetPreCheck(), SNESLineSearchSetParams(), SNESLineSearchGetParams() 1910 1911 M*/ 1912 EXTERN_C_BEGIN 1913 #undef __FUNCT__ 1914 #define __FUNCT__ "SNESCreate_VI" 1915 PetscErrorCode PETSCSNES_DLLEXPORT SNESCreate_VI(SNES snes) 1916 { 1917 PetscErrorCode ierr; 1918 SNES_VI *vi; 1919 1920 PetscFunctionBegin; 1921 snes->ops->setup = SNESSetUp_VI; 1922 snes->ops->solve = SNESSolveVI_SS; 1923 snes->ops->destroy = SNESDestroy_VI; 1924 snes->ops->setfromoptions = SNESSetFromOptions_VI; 1925 snes->ops->view = SNESView_VI; 1926 snes->ops->converged = SNESDefaultConverged_VI; 1927 1928 ierr = PetscNewLog(snes,SNES_VI,&vi);CHKERRQ(ierr); 1929 snes->data = (void*)vi; 1930 vi->alpha = 1.e-4; 1931 vi->maxstep = 1.e8; 1932 vi->minlambda = 1.e-12; 1933 vi->LineSearch = SNESLineSearchCubic_VI; 1934 vi->lsP = PETSC_NULL; 1935 vi->postcheckstep = PETSC_NULL; 1936 vi->postcheck = PETSC_NULL; 1937 vi->precheckstep = PETSC_NULL; 1938 vi->precheck = PETSC_NULL; 1939 vi->const_tol = 2.2204460492503131e-16; 1940 vi->computessfunction = ComputeFischerFunction; 1941 1942 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetMonitor_C","SNESLineSearchSetMonitor_VI",SNESLineSearchSetMonitor_VI);CHKERRQ(ierr); 1943 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSet_C","SNESLineSearchSet_VI",SNESLineSearchSet_VI);CHKERRQ(ierr); 1944 1945 PetscFunctionReturn(0); 1946 } 1947 EXTERN_C_END 1948