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; 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 /* Create active and inactive set vectors */ 844 ierr = SNESVICreateVectors_AS(snes,nis_act,&Da_act);CHKERRQ(ierr); 845 ierr = SNESVICreateVectors_AS(snes,nis_inact,&Da_inact);CHKERRQ(ierr); 846 ierr = SNESVICreateVectors_AS(snes,nis_inact,&Db_inact);CHKERRQ(ierr); 847 ierr = SNESVICreateVectors_AS(snes,nis_act,&phi_act);CHKERRQ(ierr); 848 ierr = SNESVICreateVectors_AS(snes,nis_inact,&phi_inact);CHKERRQ(ierr); 849 ierr = SNESVICreateVectors_AS(snes,nis_act,&Y_act);CHKERRQ(ierr); 850 ierr = SNESVICreateVectors_AS(snes,nis_inact,&Y_inact);CHKERRQ(ierr); 851 852 /* Create inactive set submatrices */ 853 ierr = MatGetSubMatrix(snes->jacobian,IS_inact,IS_act,MAT_INITIAL_MATRIX,&jac_inact_act);CHKERRQ(ierr); 854 ierr = MatGetSubMatrix(snes->jacobian,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);CHKERRQ(ierr); 855 856 /* Create scatter contexts */ 857 ierr = VecScatterCreate(vi->Da,IS_act,Da_act,PETSC_NULL,&scat_act);CHKERRQ(ierr); 858 ierr = VecScatterCreate(vi->Da,IS_inact,Da_inact,PETSC_NULL,&scat_inact);CHKERRQ(ierr); 859 860 /* Do a vec scatter to active and inactive set vectors */ 861 ierr = VecScatterBegin(scat_act,vi->Da,Da_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 862 ierr = VecScatterEnd(scat_act,vi->Da,Da_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 863 864 ierr = VecScatterBegin(scat_inact,vi->Da,Da_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 865 ierr = VecScatterEnd(scat_inact,vi->Da,Da_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 866 867 ierr = VecScatterBegin(scat_inact,vi->Db,Db_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 868 ierr = VecScatterEnd(scat_inact,vi->Db,Db_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 869 870 ierr = VecScatterBegin(scat_act,vi->phi,phi_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 871 ierr = VecScatterEnd(scat_act,vi->phi,phi_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 872 873 ierr = VecScatterBegin(scat_inact,vi->phi,phi_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 874 ierr = VecScatterEnd(scat_inact,vi->phi,phi_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 875 876 ierr = VecScatterBegin(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 877 ierr = VecScatterEnd(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 878 879 ierr = VecScatterBegin(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 880 ierr = VecScatterEnd(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 881 882 /* Active set direction */ 883 ierr = VecPointwiseDivide(Y_act,phi_act,Da_act);CHKERRQ(ierr); 884 /* inactive set jacobian */ 885 ierr = VecPointwiseDivide(Da_inact,Da_inact,Db_inact);CHKERRQ(ierr); 886 ierr = MatDiagonalSet(jac_inact_inact,Da_inact,ADD_VALUES);CHKERRQ(ierr); 887 /* right hand side */ 888 ierr = VecPointwiseDivide(phi_inact,phi_inact,Db_inact);CHKERRQ(ierr); 889 ierr = MatMult(jac_inact_act,Y_act,Db_inact);CHKERRQ(ierr); 890 ierr = VecAXPY(phi_inact,-1.0,Db_inact);CHKERRQ(ierr); 891 892 /* USING THE SAME MATRIX AS PRECONDITIONER....NEED TO CHANGE THIS */ 893 ierr = KSPSetOperators(snes->ksp,jac_inact_inact,jac_inact_inact,flg);CHKERRQ(ierr); 894 ierr = SNES_KSPSolve(snes,snes->ksp,phi_inact,Y_inact);CHKERRQ(ierr); 895 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 896 /* Compute the jacobian of the semismooth function which is needed for calculating the merit function 897 gradient */ 898 ierr = SNESVIComputeJacobian(snes->jacobian,snes->jacobian_pre,vi->Da,vi->Db);CHKERRQ(ierr); 899 ierr = SNESVIComputeMeritFunctionGradient(snes->jacobian,vi->phi,vi->dpsi);CHKERRQ(ierr); 900 901 ierr = VecScatterBegin(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 902 ierr = VecScatterEnd(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 903 ierr = VecScatterBegin(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 904 ierr = VecScatterEnd(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 905 906 ierr = VecDestroy(Da_act);CHKERRQ(ierr); 907 ierr = VecDestroy(Da_inact);CHKERRQ(ierr); 908 ierr = VecDestroy(Db_inact);CHKERRQ(ierr); 909 ierr = VecDestroy(phi_act);CHKERRQ(ierr); 910 ierr = VecDestroy(phi_inact);CHKERRQ(ierr); 911 ierr = VecDestroy(Y_act);CHKERRQ(ierr); 912 ierr = VecDestroy(Y_inact);CHKERRQ(ierr); 913 ierr = VecScatterDestroy(scat_act);CHKERRQ(ierr); 914 ierr = VecScatterDestroy(scat_inact);CHKERRQ(ierr); 915 ierr = ISDestroy(IS_act);CHKERRQ(ierr); 916 ierr = ISDestroy(IS_inact);CHKERRQ(ierr); 917 918 ierr = SNESVICheckDescentDirection(snes,vi->dpsi,Y,&changedir);CHKERRQ(ierr); 919 if (kspreason < 0 || changedir) { 920 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 921 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 922 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 923 break; 924 } 925 ierr = VecCopy(vi->dpsi,Y);CHKERRQ(ierr); 926 } 927 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 928 snes->linear_its += lits; 929 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 930 /* 931 if (vi->precheckstep) { 932 PetscBool changed_y = PETSC_FALSE; 933 ierr = (*vi->precheckstep)(snes,X,Y,vi->precheck,&changed_y);CHKERRQ(ierr); 934 } 935 936 if (PetscLogPrintInfo){ 937 ierr = SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr); 938 } 939 */ 940 /* Compute a (scaled) negative update in the line search routine: 941 Y <- X - lambda*Y 942 and evaluate G = function(Y) (depends on the line search). 943 */ 944 ierr = VecCopy(Y,snes->vec_sol_update);CHKERRQ(ierr); 945 ynorm = 1; gnorm = vi->phinorm; 946 ierr = (*vi->LineSearch)(snes,vi->lsP,X,vi->phi,G,Y,W,vi->phinorm,xnorm,&ynorm,&gnorm,&lssucceed);CHKERRQ(ierr); 947 ierr = PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",vi->phinorm,gnorm,ynorm,(int)lssucceed);CHKERRQ(ierr); 948 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 949 if (snes->domainerror) { 950 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 951 PetscFunctionReturn(0); 952 } 953 if (!lssucceed) { 954 if (++snes->numFailures >= snes->maxFailures) { 955 PetscBool ismin; 956 snes->reason = SNES_DIVERGED_LINE_SEARCH; 957 ierr = SNESVICheckLocalMin_Private(snes,snes->jacobian,G,W,gnorm,&ismin);CHKERRQ(ierr); 958 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 959 break; 960 } 961 } 962 /* Update function and solution vectors */ 963 vi->phinorm = gnorm; 964 vi->merit = 0.5*vi->phinorm*vi->phinorm; 965 ierr = VecCopy(G,vi->phi);CHKERRQ(ierr); 966 ierr = VecCopy(W,X);CHKERRQ(ierr); 967 /* Monitor convergence */ 968 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 969 snes->iter = i+1; 970 snes->norm = vi->phinorm; 971 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 972 SNESLogConvHistory(snes,snes->norm,lits); 973 SNESMonitor(snes,snes->iter,snes->norm); 974 /* Test for convergence, xnorm = || X || */ 975 if (snes->ops->converged != SNESSkipConverged) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 976 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,vi->phinorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 977 if (snes->reason) break; 978 } 979 if (i == maxits) { 980 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 981 if(!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 982 } 983 PetscFunctionReturn(0); 984 } 985 986 /* -------------------------------------------------------------------------- */ 987 /* 988 SNESSetUp_VI - Sets up the internal data structures for the later use 989 of the SNESVI nonlinear solver. 990 991 Input Parameter: 992 . snes - the SNES context 993 . x - the solution vector 994 995 Application Interface Routine: SNESSetUp() 996 997 Notes: 998 For basic use of the SNES solvers, the user need not explicitly call 999 SNESSetUp(), since these actions will automatically occur during 1000 the call to SNESSolve(). 1001 */ 1002 #undef __FUNCT__ 1003 #define __FUNCT__ "SNESSetUp_VI" 1004 PetscErrorCode SNESSetUp_VI(SNES snes) 1005 { 1006 PetscErrorCode ierr; 1007 SNES_VI *vi = (SNES_VI*) snes->data; 1008 PetscInt i_start[3],i_end[3]; 1009 1010 PetscFunctionBegin; 1011 if (!snes->vec_sol_update) { 1012 ierr = VecDuplicate(snes->vec_sol,&snes->vec_sol_update);CHKERRQ(ierr); 1013 ierr = PetscLogObjectParent(snes,snes->vec_sol_update);CHKERRQ(ierr); 1014 } 1015 if (!snes->work) { 1016 snes->nwork = 3; 1017 ierr = VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work);CHKERRQ(ierr); 1018 ierr = PetscLogObjectParents(snes,snes->nwork,snes->work);CHKERRQ(ierr); 1019 } 1020 1021 ierr = VecDuplicate(snes->vec_sol, &vi->phi); CHKERRQ(ierr); 1022 ierr = VecDuplicate(snes->vec_sol, &vi->dpsi); CHKERRQ(ierr); 1023 ierr = VecDuplicate(snes->vec_sol, &vi->Da); CHKERRQ(ierr); 1024 ierr = VecDuplicate(snes->vec_sol, &vi->Db); CHKERRQ(ierr); 1025 ierr = VecDuplicate(snes->vec_sol, &vi->z);CHKERRQ(ierr); 1026 ierr = VecDuplicate(snes->vec_sol, &vi->t); CHKERRQ(ierr); 1027 1028 /* If the lower and upper bound on variables are not set, set it to 1029 -Inf and Inf */ 1030 if (!vi->xl && !vi->xu) { 1031 vi->usersetxbounds = PETSC_FALSE; 1032 ierr = VecDuplicate(snes->vec_sol, &vi->xl); CHKERRQ(ierr); 1033 ierr = VecSet(vi->xl,PETSC_VI_NINF);CHKERRQ(ierr); 1034 ierr = VecDuplicate(snes->vec_sol, &vi->xu); CHKERRQ(ierr); 1035 ierr = VecSet(vi->xu,PETSC_VI_INF);CHKERRQ(ierr); 1036 } else { 1037 /* Check if lower bound, upper bound and solution vector distribution across the processors is identical */ 1038 ierr = VecGetOwnershipRange(snes->vec_sol,i_start,i_end);CHKERRQ(ierr); 1039 ierr = VecGetOwnershipRange(vi->xl,i_start+1,i_end+1);CHKERRQ(ierr); 1040 ierr = VecGetOwnershipRange(vi->xu,i_start+2,i_end+2);CHKERRQ(ierr); 1041 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])) 1042 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."); 1043 } 1044 1045 vi->computeuserfunction = snes->ops->computefunction; 1046 snes->ops->computefunction = SNESVIComputeFunction; 1047 1048 PetscFunctionReturn(0); 1049 } 1050 /* -------------------------------------------------------------------------- */ 1051 /* 1052 SNESDestroy_VI - Destroys the private SNES_VI context that was created 1053 with SNESCreate_VI(). 1054 1055 Input Parameter: 1056 . snes - the SNES context 1057 1058 Application Interface Routine: SNESDestroy() 1059 */ 1060 #undef __FUNCT__ 1061 #define __FUNCT__ "SNESDestroy_VI" 1062 PetscErrorCode SNESDestroy_VI(SNES snes) 1063 { 1064 SNES_VI *vi = (SNES_VI*) snes->data; 1065 PetscErrorCode ierr; 1066 1067 PetscFunctionBegin; 1068 if (snes->vec_sol_update) { 1069 ierr = VecDestroy(snes->vec_sol_update);CHKERRQ(ierr); 1070 snes->vec_sol_update = PETSC_NULL; 1071 } 1072 if (snes->nwork) { 1073 ierr = VecDestroyVecs(snes->work,snes->nwork);CHKERRQ(ierr); 1074 snes->nwork = 0; 1075 snes->work = PETSC_NULL; 1076 } 1077 1078 /* clear vectors */ 1079 ierr = VecDestroy(vi->phi); CHKERRQ(ierr); 1080 ierr = VecDestroy(vi->dpsi); CHKERRQ(ierr); 1081 ierr = VecDestroy(vi->Da); CHKERRQ(ierr); 1082 ierr = VecDestroy(vi->Db); CHKERRQ(ierr); 1083 ierr = VecDestroy(vi->z); CHKERRQ(ierr); 1084 ierr = VecDestroy(vi->t); CHKERRQ(ierr); 1085 if (!vi->usersetxbounds) { 1086 ierr = VecDestroy(vi->xl); CHKERRQ(ierr); 1087 ierr = VecDestroy(vi->xu); CHKERRQ(ierr); 1088 } 1089 if (vi->lsmonitor) { 1090 ierr = PetscViewerASCIIMonitorDestroy(vi->lsmonitor);CHKERRQ(ierr); 1091 } 1092 ierr = PetscFree(snes->data);CHKERRQ(ierr); 1093 1094 /* clear composed functions */ 1095 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSet_C","",PETSC_NULL);CHKERRQ(ierr); 1096 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetMonitor_C","",PETSC_NULL);CHKERRQ(ierr); 1097 PetscFunctionReturn(0); 1098 } 1099 /* -------------------------------------------------------------------------- */ 1100 #undef __FUNCT__ 1101 #define __FUNCT__ "SNESLineSearchNo_VI" 1102 1103 /* 1104 This routine is a copy of SNESLineSearchNo routine in snes/impls/ls/ls.c 1105 1106 */ 1107 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) 1108 { 1109 PetscErrorCode ierr; 1110 SNES_VI *vi = (SNES_VI*)snes->data; 1111 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1112 1113 PetscFunctionBegin; 1114 *flag = PETSC_TRUE; 1115 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1116 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); /* ynorm = || y || */ 1117 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1118 if (vi->postcheckstep) { 1119 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1120 } 1121 if (changed_y) { 1122 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1123 } 1124 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1125 if (!snes->domainerror) { 1126 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); /* gnorm = || g || */ 1127 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1128 } 1129 if (vi->lsmonitor) { 1130 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1131 } 1132 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1133 PetscFunctionReturn(0); 1134 } 1135 1136 /* -------------------------------------------------------------------------- */ 1137 #undef __FUNCT__ 1138 #define __FUNCT__ "SNESLineSearchNoNorms_VI" 1139 1140 /* 1141 This routine is a copy of SNESLineSearchNoNorms in snes/impls/ls/ls.c 1142 */ 1143 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) 1144 { 1145 PetscErrorCode ierr; 1146 SNES_VI *vi = (SNES_VI*)snes->data; 1147 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1148 1149 PetscFunctionBegin; 1150 *flag = PETSC_TRUE; 1151 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1152 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1153 if (vi->postcheckstep) { 1154 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1155 } 1156 if (changed_y) { 1157 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1158 } 1159 1160 /* don't evaluate function the last time through */ 1161 if (snes->iter < snes->max_its-1) { 1162 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1163 } 1164 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1165 PetscFunctionReturn(0); 1166 } 1167 /* -------------------------------------------------------------------------- */ 1168 #undef __FUNCT__ 1169 #define __FUNCT__ "SNESLineSearchCubic_VI" 1170 /* 1171 This routine is a copy of SNESLineSearchCubic in snes/impls/ls/ls.c 1172 */ 1173 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) 1174 { 1175 /* 1176 Note that for line search purposes we work with with the related 1177 minimization problem: 1178 min z(x): R^n -> R, 1179 where z(x) = .5 * fnorm*fnorm, and fnorm = || f ||_2. 1180 This function z(x) is same as the merit function for the complementarity problem. 1181 */ 1182 1183 PetscReal initslope,lambdaprev,gnormprev,a,b,d,t1,t2,rellength; 1184 PetscReal minlambda,lambda,lambdatemp; 1185 #if defined(PETSC_USE_COMPLEX) 1186 PetscScalar cinitslope; 1187 #endif 1188 PetscErrorCode ierr; 1189 PetscInt count; 1190 SNES_VI *vi = (SNES_VI*)snes->data; 1191 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1192 MPI_Comm comm; 1193 1194 PetscFunctionBegin; 1195 ierr = PetscObjectGetComm((PetscObject)snes,&comm);CHKERRQ(ierr); 1196 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1197 *flag = PETSC_TRUE; 1198 1199 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); 1200 if (*ynorm == 0.0) { 1201 if (vi->lsmonitor) { 1202 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Initial direction and size is 0\n");CHKERRQ(ierr); 1203 } 1204 *gnorm = fnorm; 1205 ierr = VecCopy(x,w);CHKERRQ(ierr); 1206 ierr = VecCopy(f,g);CHKERRQ(ierr); 1207 *flag = PETSC_FALSE; 1208 goto theend1; 1209 } 1210 if (*ynorm > vi->maxstep) { /* Step too big, so scale back */ 1211 if (vi->lsmonitor) { 1212 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Scaling step by %G old ynorm %G\n",vi->maxstep/(*ynorm),*ynorm);CHKERRQ(ierr); 1213 } 1214 ierr = VecScale(y,vi->maxstep/(*ynorm));CHKERRQ(ierr); 1215 *ynorm = vi->maxstep; 1216 } 1217 ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr); 1218 minlambda = vi->minlambda/rellength; 1219 /* ierr = MatMult(snes->jacobian,y,w);CHKERRQ(ierr); */ 1220 #if defined(PETSC_USE_COMPLEX) 1221 ierr = VecDot(vi->dpsi,y,&cinitslope);CHKERRQ(ierr); 1222 initslope = PetscRealPart(cinitslope); 1223 #else 1224 ierr = VecDot(vi->dpsi,y,&initslope);CHKERRQ(ierr); 1225 #endif 1226 if (initslope > 0.0) initslope = -initslope; 1227 if (initslope == 0.0) initslope = -1.0; 1228 1229 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1230 if (snes->nfuncs >= snes->max_funcs) { 1231 ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); 1232 *flag = PETSC_FALSE; 1233 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1234 goto theend1; 1235 } 1236 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1237 if (snes->domainerror) { 1238 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1239 PetscFunctionReturn(0); 1240 } 1241 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1242 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1243 ierr = PetscInfo2(snes,"Initial fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1244 if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + vi->alpha*initslope) { /* Sufficient reduction */ 1245 if (vi->lsmonitor) { 1246 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1247 } 1248 goto theend1; 1249 } 1250 1251 /* Fit points with quadratic */ 1252 lambda = 1.0; 1253 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 1254 lambdaprev = lambda; 1255 gnormprev = *gnorm; 1256 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1257 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1258 else lambda = lambdatemp; 1259 1260 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1261 if (snes->nfuncs >= snes->max_funcs) { 1262 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);CHKERRQ(ierr); 1263 *flag = PETSC_FALSE; 1264 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1265 goto theend1; 1266 } 1267 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1268 if (snes->domainerror) { 1269 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1270 PetscFunctionReturn(0); 1271 } 1272 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1273 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1274 if (vi->lsmonitor) { 1275 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: gnorm after quadratic fit %G\n",*gnorm);CHKERRQ(ierr); 1276 } 1277 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* sufficient reduction */ 1278 if (vi->lsmonitor) { 1279 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Quadratically determined step, lambda=%18.16e\n",lambda);CHKERRQ(ierr); 1280 } 1281 goto theend1; 1282 } 1283 1284 /* Fit points with cubic */ 1285 count = 1; 1286 while (PETSC_TRUE) { 1287 if (lambda <= minlambda) { 1288 if (vi->lsmonitor) { 1289 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: unable to find good step length! After %D tries \n",count);CHKERRQ(ierr); 1290 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); 1291 } 1292 *flag = PETSC_FALSE; 1293 break; 1294 } 1295 t1 = .5*((*gnorm)*(*gnorm) - fnorm*fnorm) - lambda*initslope; 1296 t2 = .5*(gnormprev*gnormprev - fnorm*fnorm) - lambdaprev*initslope; 1297 a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 1298 b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 1299 d = b*b - 3*a*initslope; 1300 if (d < 0.0) d = 0.0; 1301 if (a == 0.0) { 1302 lambdatemp = -initslope/(2.0*b); 1303 } else { 1304 lambdatemp = (-b + sqrt(d))/(3.0*a); 1305 } 1306 lambdaprev = lambda; 1307 gnormprev = *gnorm; 1308 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1309 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1310 else lambda = lambdatemp; 1311 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1312 if (snes->nfuncs >= snes->max_funcs) { 1313 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 1314 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); 1315 *flag = PETSC_FALSE; 1316 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1317 break; 1318 } 1319 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1320 if (snes->domainerror) { 1321 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1322 PetscFunctionReturn(0); 1323 } 1324 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1325 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1326 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* is reduction enough? */ 1327 if (vi->lsmonitor) { 1328 ierr = PetscPrintf(comm," Line search: Cubically determined step, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr); 1329 } 1330 break; 1331 } else { 1332 if (vi->lsmonitor) { 1333 ierr = PetscPrintf(comm," Line search: Cubic step no good, shrinking lambda, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr); 1334 } 1335 } 1336 count++; 1337 } 1338 theend1: 1339 /* Optional user-defined check for line search step validity */ 1340 if (vi->postcheckstep && *flag) { 1341 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1342 if (changed_y) { 1343 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1344 } 1345 if (changed_y || changed_w) { /* recompute the function if the step has changed */ 1346 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1347 if (snes->domainerror) { 1348 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1349 PetscFunctionReturn(0); 1350 } 1351 ierr = VecNormBegin(g,NORM_2,gnorm);CHKERRQ(ierr); 1352 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1353 ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr); 1354 ierr = VecNormEnd(g,NORM_2,gnorm);CHKERRQ(ierr); 1355 ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr); 1356 } 1357 } 1358 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1359 PetscFunctionReturn(0); 1360 } 1361 /* -------------------------------------------------------------------------- */ 1362 #undef __FUNCT__ 1363 #define __FUNCT__ "SNESLineSearchQuadratic_VI" 1364 /* 1365 This routine is a copy of SNESLineSearchQuadratic in snes/impls/ls/ls.c 1366 */ 1367 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) 1368 { 1369 /* 1370 Note that for line search purposes we work with with the related 1371 minimization problem: 1372 min z(x): R^n -> R, 1373 where z(x) = .5 * fnorm*fnorm,and fnorm = || f ||_2. 1374 z(x) is the same as the merit function for the complementarity problem 1375 */ 1376 PetscReal initslope,minlambda,lambda,lambdatemp,rellength; 1377 #if defined(PETSC_USE_COMPLEX) 1378 PetscScalar cinitslope; 1379 #endif 1380 PetscErrorCode ierr; 1381 PetscInt count; 1382 SNES_VI *vi = (SNES_VI*)snes->data; 1383 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1384 1385 PetscFunctionBegin; 1386 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1387 *flag = PETSC_TRUE; 1388 1389 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); 1390 if (*ynorm == 0.0) { 1391 if (vi->lsmonitor) { 1392 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor,"Line search: Direction and size is 0\n");CHKERRQ(ierr); 1393 } 1394 *gnorm = fnorm; 1395 ierr = VecCopy(x,w);CHKERRQ(ierr); 1396 ierr = VecCopy(f,g);CHKERRQ(ierr); 1397 *flag = PETSC_FALSE; 1398 goto theend2; 1399 } 1400 if (*ynorm > vi->maxstep) { /* Step too big, so scale back */ 1401 ierr = VecScale(y,vi->maxstep/(*ynorm));CHKERRQ(ierr); 1402 *ynorm = vi->maxstep; 1403 } 1404 ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr); 1405 minlambda = vi->minlambda/rellength; 1406 /* ierr = MatMult(snes->jacobian,y,w);CHKERRQ(ierr); */ 1407 #if defined(PETSC_USE_COMPLEX) 1408 ierr = VecDot(f,w,&cinitslope);CHKERRQ(ierr); 1409 initslope = PetscRealPart(cinitslope); 1410 #else 1411 ierr = VecDot(vi->dpsi,y,&initslope);CHKERRQ(ierr); 1412 #endif 1413 if (initslope > 0.0) initslope = -initslope; 1414 if (initslope == 0.0) initslope = -1.0; 1415 ierr = PetscInfo1(snes,"Initslope %G \n",initslope);CHKERRQ(ierr); 1416 1417 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1418 if (snes->nfuncs >= snes->max_funcs) { 1419 ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); 1420 *flag = PETSC_FALSE; 1421 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1422 goto theend2; 1423 } 1424 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1425 if (snes->domainerror) { 1426 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1427 PetscFunctionReturn(0); 1428 } 1429 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1430 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1431 if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + vi->alpha*initslope) { /* Sufficient reduction */ 1432 if (vi->lsmonitor) { 1433 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1434 } 1435 goto theend2; 1436 } 1437 1438 /* Fit points with quadratic */ 1439 lambda = 1.0; 1440 count = 1; 1441 while (PETSC_TRUE) { 1442 if (lambda <= minlambda) { /* bad luck; use full step */ 1443 if (vi->lsmonitor) { 1444 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor,"Line search: Unable to find good step length! %D \n",count);CHKERRQ(ierr); 1445 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); 1446 } 1447 ierr = VecCopy(x,w);CHKERRQ(ierr); 1448 *flag = PETSC_FALSE; 1449 break; 1450 } 1451 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 1452 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1453 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1454 else lambda = lambdatemp; 1455 1456 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1457 if (snes->nfuncs >= snes->max_funcs) { 1458 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 1459 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); 1460 *flag = PETSC_FALSE; 1461 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1462 break; 1463 } 1464 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1465 if (snes->domainerror) { 1466 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1467 PetscFunctionReturn(0); 1468 } 1469 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1470 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1471 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* sufficient reduction */ 1472 if (vi->lsmonitor) { 1473 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line Search: Quadratically determined step, lambda=%G\n",lambda);CHKERRQ(ierr); 1474 } 1475 break; 1476 } 1477 count++; 1478 } 1479 theend2: 1480 /* Optional user-defined check for line search step validity */ 1481 if (vi->postcheckstep) { 1482 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1483 if (changed_y) { 1484 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1485 } 1486 if (changed_y || changed_w) { /* recompute the function if the step has changed */ 1487 ierr = SNESComputeFunction(snes,w,g); 1488 if (snes->domainerror) { 1489 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1490 PetscFunctionReturn(0); 1491 } 1492 ierr = VecNormBegin(g,NORM_2,gnorm);CHKERRQ(ierr); 1493 ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr); 1494 ierr = VecNormEnd(g,NORM_2,gnorm);CHKERRQ(ierr); 1495 ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr); 1496 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1497 } 1498 } 1499 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1500 PetscFunctionReturn(0); 1501 } 1502 1503 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*/ 1504 /* -------------------------------------------------------------------------- */ 1505 EXTERN_C_BEGIN 1506 #undef __FUNCT__ 1507 #define __FUNCT__ "SNESLineSearchSet_VI" 1508 PetscErrorCode PETSCSNES_DLLEXPORT SNESLineSearchSet_VI(SNES snes,FCN2 func,void *lsctx) 1509 { 1510 PetscFunctionBegin; 1511 ((SNES_VI *)(snes->data))->LineSearch = func; 1512 ((SNES_VI *)(snes->data))->lsP = lsctx; 1513 PetscFunctionReturn(0); 1514 } 1515 EXTERN_C_END 1516 1517 /* -------------------------------------------------------------------------- */ 1518 EXTERN_C_BEGIN 1519 #undef __FUNCT__ 1520 #define __FUNCT__ "SNESLineSearchSetMonitor_VI" 1521 PetscErrorCode PETSCSNES_DLLEXPORT SNESLineSearchSetMonitor_VI(SNES snes,PetscBool flg) 1522 { 1523 SNES_VI *vi = (SNES_VI*)snes->data; 1524 PetscErrorCode ierr; 1525 1526 PetscFunctionBegin; 1527 if (flg && !vi->lsmonitor) { 1528 ierr = PetscViewerASCIIMonitorCreate(((PetscObject)snes)->comm,"stdout",((PetscObject)snes)->tablevel,&vi->lsmonitor);CHKERRQ(ierr); 1529 } else if (!flg && vi->lsmonitor) { 1530 ierr = PetscViewerASCIIMonitorDestroy(vi->lsmonitor);CHKERRQ(ierr); 1531 } 1532 PetscFunctionReturn(0); 1533 } 1534 EXTERN_C_END 1535 1536 /* 1537 SNESView_VI - Prints info from the SNESVI data structure. 1538 1539 Input Parameters: 1540 . SNES - the SNES context 1541 . viewer - visualization context 1542 1543 Application Interface Routine: SNESView() 1544 */ 1545 #undef __FUNCT__ 1546 #define __FUNCT__ "SNESView_VI" 1547 static PetscErrorCode SNESView_VI(SNES snes,PetscViewer viewer) 1548 { 1549 SNES_VI *vi = (SNES_VI *)snes->data; 1550 const char *cstr; 1551 PetscErrorCode ierr; 1552 PetscBool iascii; 1553 1554 PetscFunctionBegin; 1555 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1556 if (iascii) { 1557 if (vi->LineSearch == SNESLineSearchNo_VI) cstr = "SNESLineSearchNo"; 1558 else if (vi->LineSearch == SNESLineSearchQuadratic_VI) cstr = "SNESLineSearchQuadratic"; 1559 else if (vi->LineSearch == SNESLineSearchCubic_VI) cstr = "SNESLineSearchCubic"; 1560 else cstr = "unknown"; 1561 ierr = PetscViewerASCIIPrintf(viewer," line search variant: %s\n",cstr);CHKERRQ(ierr); 1562 ierr = PetscViewerASCIIPrintf(viewer," alpha=%G, maxstep=%G, minlambda=%G\n",vi->alpha,vi->maxstep,vi->minlambda);CHKERRQ(ierr); 1563 } else { 1564 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for SNES EQ VI",((PetscObject)viewer)->type_name); 1565 } 1566 PetscFunctionReturn(0); 1567 } 1568 1569 /* 1570 SNESVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu. 1571 1572 Input Parameters: 1573 . snes - the SNES context. 1574 . xl - lower bound. 1575 . xu - upper bound. 1576 1577 Notes: 1578 If this routine is not called then the lower and upper bounds are set to 1579 -Infinity and Infinity respectively during SNESSetUp. 1580 */ 1581 1582 #undef __FUNCT__ 1583 #define __FUNCT__ "SNESVISetVariableBounds" 1584 PetscErrorCode SNESVISetVariableBounds(SNES snes, Vec xl, Vec xu) 1585 { 1586 SNES_VI *vi = (SNES_VI*)snes->data; 1587 1588 PetscFunctionBegin; 1589 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1590 PetscValidHeaderSpecific(xl,VEC_CLASSID,2); 1591 PetscValidHeaderSpecific(xu,VEC_CLASSID,3); 1592 1593 /* Check if SNESSetFunction is called */ 1594 if(!snes->vec_func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call SNESSetFunction() first"); 1595 1596 /* Check if the vector sizes are compatible for lower and upper bounds */ 1597 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); 1598 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); 1599 vi->usersetxbounds = PETSC_TRUE; 1600 vi->xl = xl; 1601 vi->xu = xu; 1602 1603 PetscFunctionReturn(0); 1604 } 1605 /* -------------------------------------------------------------------------- */ 1606 /* 1607 SNESSetFromOptions_VI - Sets various parameters for the SNESVI method. 1608 1609 Input Parameter: 1610 . snes - the SNES context 1611 1612 Application Interface Routine: SNESSetFromOptions() 1613 */ 1614 #undef __FUNCT__ 1615 #define __FUNCT__ "SNESSetFromOptions_VI" 1616 static PetscErrorCode SNESSetFromOptions_VI(SNES snes) 1617 { 1618 SNES_VI *vi = (SNES_VI *)snes->data; 1619 const char *lses[] = {"basic","basicnonorms","quadratic","cubic"}; 1620 const char *vies[] = {"ss","as"}; 1621 PetscErrorCode ierr; 1622 PetscInt indx; 1623 PetscBool flg,set,flg2; 1624 1625 PetscFunctionBegin; 1626 ierr = PetscOptionsHead("SNES semismooth method options");CHKERRQ(ierr); 1627 ierr = PetscOptionsReal("-snes_vi_alpha","Function norm must decrease by","None",vi->alpha,&vi->alpha,0);CHKERRQ(ierr); 1628 ierr = PetscOptionsReal("-snes_vi_maxstep","Step must be less than","None",vi->maxstep,&vi->maxstep,0);CHKERRQ(ierr); 1629 ierr = PetscOptionsReal("-snes_vi_minlambda","Minimum lambda allowed","None",vi->minlambda,&vi->minlambda,0);CHKERRQ(ierr); 1630 ierr = PetscOptionsReal("-snes_vi_delta","descent test fraction","None",vi->delta,&vi->delta,0);CHKERRQ(ierr); 1631 ierr = PetscOptionsReal("-snes_vi_rho","descent test power","None",vi->rho,&vi->rho,0);CHKERRQ(ierr); 1632 ierr = PetscOptionsReal("-snes_vi_const_tol","constraint tolerance","None",vi->const_tol,&vi->const_tol,0);CHKERRQ(ierr); 1633 ierr = PetscOptionsBool("-snes_vi_lsmonitor","Print progress of line searches","SNESLineSearchSetMonitor",vi->lsmonitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);CHKERRQ(ierr); 1634 if (set) {ierr = SNESLineSearchSetMonitor(snes,flg);CHKERRQ(ierr);} 1635 ierr = PetscOptionsEList("-snes_vi_type","Semismooth algorithm used","",vies,2,"ss",&indx,&flg2);CHKERRQ(ierr); 1636 if (flg2) { 1637 switch (indx) { 1638 case 0: 1639 snes->ops->solve = SNESSolveVI_SS; 1640 break; 1641 case 1: 1642 snes->ops->solve = SNESSolveVI_AS; 1643 break; 1644 } 1645 } 1646 ierr = PetscOptionsEList("-snes_vi_ls","Line search used","SNESLineSearchSet",lses,4,"cubic",&indx,&flg);CHKERRQ(ierr); 1647 if (flg) { 1648 switch (indx) { 1649 case 0: 1650 ierr = SNESLineSearchSet(snes,SNESLineSearchNo_VI,PETSC_NULL);CHKERRQ(ierr); 1651 break; 1652 case 1: 1653 ierr = SNESLineSearchSet(snes,SNESLineSearchNoNorms_VI,PETSC_NULL);CHKERRQ(ierr); 1654 break; 1655 case 2: 1656 ierr = SNESLineSearchSet(snes,SNESLineSearchQuadratic_VI,PETSC_NULL);CHKERRQ(ierr); 1657 break; 1658 case 3: 1659 ierr = SNESLineSearchSet(snes,SNESLineSearchCubic_VI,PETSC_NULL);CHKERRQ(ierr); 1660 break; 1661 } 1662 } 1663 ierr = PetscOptionsTail();CHKERRQ(ierr); 1664 PetscFunctionReturn(0); 1665 } 1666 /* -------------------------------------------------------------------------- */ 1667 /*MC 1668 SNESVI - Semismooth newton method based nonlinear solver that uses a line search 1669 1670 Options Database: 1671 + -snes_vi [cubic,quadratic,basic,basicnonorms] - Selects line search 1672 . -snes_vi_alpha <alpha> - Sets alpha 1673 . -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) 1674 . -snes_vi_minlambda <minlambda> - Sets the minimum lambda the line search will use minlambda / max_i ( y[i]/x[i] ) 1675 -snes_vi_delta <delta> - Sets the fraction used in the descent test. 1676 -snes_vi_rho <rho> - Sets the power used in the descent test. 1677 For a descent direction to be accepted it has to satisfy the condition dpsi^T*d <= -delta*||d||^rho 1678 - -snes_vi_monitor - print information about progress of line searches 1679 1680 1681 Level: beginner 1682 1683 .seealso: SNESCreate(), SNES, SNESSetType(), SNESTR, SNESLineSearchSet(), 1684 SNESLineSearchSetPostCheck(), SNESLineSearchNo(), SNESLineSearchCubic(), SNESLineSearchQuadratic(), 1685 SNESLineSearchSet(), SNESLineSearchNoNorms(), SNESLineSearchSetPreCheck(), SNESLineSearchSetParams(), SNESLineSearchGetParams() 1686 1687 M*/ 1688 EXTERN_C_BEGIN 1689 #undef __FUNCT__ 1690 #define __FUNCT__ "SNESCreate_VI" 1691 PetscErrorCode PETSCSNES_DLLEXPORT SNESCreate_VI(SNES snes) 1692 { 1693 PetscErrorCode ierr; 1694 SNES_VI *vi; 1695 1696 PetscFunctionBegin; 1697 snes->ops->setup = SNESSetUp_VI; 1698 snes->ops->solve = SNESSolveVI_SS; 1699 snes->ops->destroy = SNESDestroy_VI; 1700 snes->ops->setfromoptions = SNESSetFromOptions_VI; 1701 snes->ops->view = SNESView_VI; 1702 snes->ops->converged = SNESDefaultConverged_VI; 1703 1704 ierr = PetscNewLog(snes,SNES_VI,&vi);CHKERRQ(ierr); 1705 snes->data = (void*)vi; 1706 vi->alpha = 1.e-4; 1707 vi->maxstep = 1.e8; 1708 vi->minlambda = 1.e-12; 1709 vi->LineSearch = SNESLineSearchCubic_VI; 1710 vi->lsP = PETSC_NULL; 1711 vi->postcheckstep = PETSC_NULL; 1712 vi->postcheck = PETSC_NULL; 1713 vi->precheckstep = PETSC_NULL; 1714 vi->precheck = PETSC_NULL; 1715 vi->rho = 2.1; 1716 vi->delta = 1e-10; 1717 vi->const_tol = 2.2204460492503131e-16; 1718 vi->computessfunction = ComputeFischerFunction; 1719 1720 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetMonitor_C","SNESLineSearchSetMonitor_VI",SNESLineSearchSetMonitor_VI);CHKERRQ(ierr); 1721 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSet_C","SNESLineSearchSet_VI",SNESLineSearchSet_VI);CHKERRQ(ierr); 1722 1723 PetscFunctionReturn(0); 1724 } 1725 EXTERN_C_END 1726