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