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