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