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 #if defined(PETSC_USE_COMPLEX) 1469 ierr = VecDot(vi->dpsi,y,&cinitslope);CHKERRQ(ierr); 1470 initslope = PetscRealPart(cinitslope); 1471 #else 1472 ierr = VecDot(vi->dpsi,y,&initslope);CHKERRQ(ierr); 1473 #endif 1474 if (initslope > 0.0) initslope = -initslope; 1475 if (initslope == 0.0) initslope = -1.0; 1476 1477 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1478 if (snes->nfuncs >= snes->max_funcs) { 1479 ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); 1480 *flag = PETSC_FALSE; 1481 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1482 goto theend1; 1483 } 1484 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1485 if (snes->domainerror) { 1486 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1487 PetscFunctionReturn(0); 1488 } 1489 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1490 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1491 ierr = PetscInfo2(snes,"Initial fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1492 if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + vi->alpha*initslope) { /* Sufficient reduction */ 1493 if (vi->lsmonitor) { 1494 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1495 } 1496 goto theend1; 1497 } 1498 1499 /* Fit points with quadratic */ 1500 lambda = 1.0; 1501 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 1502 lambdaprev = lambda; 1503 gnormprev = *gnorm; 1504 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1505 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1506 else lambda = lambdatemp; 1507 1508 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1509 if (snes->nfuncs >= snes->max_funcs) { 1510 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);CHKERRQ(ierr); 1511 *flag = PETSC_FALSE; 1512 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1513 goto theend1; 1514 } 1515 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1516 if (snes->domainerror) { 1517 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1518 PetscFunctionReturn(0); 1519 } 1520 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1521 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1522 if (vi->lsmonitor) { 1523 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: gnorm after quadratic fit %G\n",*gnorm);CHKERRQ(ierr); 1524 } 1525 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* sufficient reduction */ 1526 if (vi->lsmonitor) { 1527 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Quadratically determined step, lambda=%18.16e\n",lambda);CHKERRQ(ierr); 1528 } 1529 goto theend1; 1530 } 1531 1532 /* Fit points with cubic */ 1533 count = 1; 1534 while (PETSC_TRUE) { 1535 if (lambda <= minlambda) { 1536 if (vi->lsmonitor) { 1537 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: unable to find good step length! After %D tries \n",count);CHKERRQ(ierr); 1538 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); 1539 } 1540 *flag = PETSC_FALSE; 1541 break; 1542 } 1543 t1 = .5*((*gnorm)*(*gnorm) - fnorm*fnorm) - lambda*initslope; 1544 t2 = .5*(gnormprev*gnormprev - fnorm*fnorm) - lambdaprev*initslope; 1545 a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 1546 b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 1547 d = b*b - 3*a*initslope; 1548 if (d < 0.0) d = 0.0; 1549 if (a == 0.0) { 1550 lambdatemp = -initslope/(2.0*b); 1551 } else { 1552 lambdatemp = (-b + sqrt(d))/(3.0*a); 1553 } 1554 lambdaprev = lambda; 1555 gnormprev = *gnorm; 1556 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1557 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1558 else lambda = lambdatemp; 1559 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1560 if (snes->nfuncs >= snes->max_funcs) { 1561 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 1562 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); 1563 *flag = PETSC_FALSE; 1564 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1565 break; 1566 } 1567 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1568 if (snes->domainerror) { 1569 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1570 PetscFunctionReturn(0); 1571 } 1572 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1573 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1574 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* is reduction enough? */ 1575 if (vi->lsmonitor) { 1576 ierr = PetscPrintf(comm," Line search: Cubically determined step, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr); 1577 } 1578 break; 1579 } else { 1580 if (vi->lsmonitor) { 1581 ierr = PetscPrintf(comm," Line search: Cubic step no good, shrinking lambda, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr); 1582 } 1583 } 1584 count++; 1585 } 1586 theend1: 1587 /* Optional user-defined check for line search step validity */ 1588 if (vi->postcheckstep && *flag) { 1589 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1590 if (changed_y) { 1591 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1592 } 1593 if (changed_y || changed_w) { /* recompute the function if the step has changed */ 1594 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1595 if (snes->domainerror) { 1596 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1597 PetscFunctionReturn(0); 1598 } 1599 ierr = VecNormBegin(g,NORM_2,gnorm);CHKERRQ(ierr); 1600 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1601 ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr); 1602 ierr = VecNormEnd(g,NORM_2,gnorm);CHKERRQ(ierr); 1603 ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr); 1604 } 1605 } 1606 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1607 PetscFunctionReturn(0); 1608 } 1609 /* -------------------------------------------------------------------------- */ 1610 #undef __FUNCT__ 1611 #define __FUNCT__ "SNESLineSearchQuadratic_VI" 1612 /* 1613 This routine is a copy of SNESLineSearchQuadratic in snes/impls/ls/ls.c 1614 */ 1615 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) 1616 { 1617 /* 1618 Note that for line search purposes we work with with the related 1619 minimization problem: 1620 min z(x): R^n -> R, 1621 where z(x) = .5 * fnorm*fnorm,and fnorm = || f ||_2. 1622 z(x) is the same as the merit function for the complementarity problem 1623 */ 1624 PetscReal initslope,minlambda,lambda,lambdatemp,rellength; 1625 #if defined(PETSC_USE_COMPLEX) 1626 PetscScalar cinitslope; 1627 #endif 1628 PetscErrorCode ierr; 1629 PetscInt count; 1630 SNES_VI *vi = (SNES_VI*)snes->data; 1631 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1632 1633 PetscFunctionBegin; 1634 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1635 *flag = PETSC_TRUE; 1636 1637 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); 1638 if (*ynorm == 0.0) { 1639 if (vi->lsmonitor) { 1640 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor,"Line search: Direction and size is 0\n");CHKERRQ(ierr); 1641 } 1642 *gnorm = fnorm; 1643 ierr = VecCopy(x,w);CHKERRQ(ierr); 1644 ierr = VecCopy(f,g);CHKERRQ(ierr); 1645 *flag = PETSC_FALSE; 1646 goto theend2; 1647 } 1648 if (*ynorm > vi->maxstep) { /* Step too big, so scale back */ 1649 ierr = VecScale(y,vi->maxstep/(*ynorm));CHKERRQ(ierr); 1650 *ynorm = vi->maxstep; 1651 } 1652 ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr); 1653 minlambda = vi->minlambda/rellength; 1654 #if defined(PETSC_USE_COMPLEX) 1655 ierr = VecDot(vi->dpsi,y,&cinitslope);CHKERRQ(ierr); 1656 initslope = PetscRealPart(cinitslope); 1657 #else 1658 ierr = VecDot(vi->dpsi,y,&initslope);CHKERRQ(ierr); 1659 #endif 1660 if (initslope > 0.0) initslope = -initslope; 1661 if (initslope == 0.0) initslope = -1.0; 1662 ierr = PetscInfo1(snes,"Initslope %G \n",initslope);CHKERRQ(ierr); 1663 1664 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1665 if (snes->nfuncs >= snes->max_funcs) { 1666 ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); 1667 *flag = PETSC_FALSE; 1668 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1669 goto theend2; 1670 } 1671 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1672 if (snes->domainerror) { 1673 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1674 PetscFunctionReturn(0); 1675 } 1676 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1677 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1678 if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + vi->alpha*initslope) { /* Sufficient reduction */ 1679 if (vi->lsmonitor) { 1680 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1681 } 1682 goto theend2; 1683 } 1684 1685 /* Fit points with quadratic */ 1686 lambda = 1.0; 1687 count = 1; 1688 while (PETSC_TRUE) { 1689 if (lambda <= minlambda) { /* bad luck; use full step */ 1690 if (vi->lsmonitor) { 1691 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor,"Line search: Unable to find good step length! %D \n",count);CHKERRQ(ierr); 1692 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); 1693 } 1694 ierr = VecCopy(x,w);CHKERRQ(ierr); 1695 *flag = PETSC_FALSE; 1696 break; 1697 } 1698 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 1699 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1700 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1701 else lambda = lambdatemp; 1702 1703 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1704 if (snes->nfuncs >= snes->max_funcs) { 1705 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 1706 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); 1707 *flag = PETSC_FALSE; 1708 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1709 break; 1710 } 1711 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1712 if (snes->domainerror) { 1713 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1714 PetscFunctionReturn(0); 1715 } 1716 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1717 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1718 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* sufficient reduction */ 1719 if (vi->lsmonitor) { 1720 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line Search: Quadratically determined step, lambda=%G\n",lambda);CHKERRQ(ierr); 1721 } 1722 break; 1723 } 1724 count++; 1725 } 1726 theend2: 1727 /* Optional user-defined check for line search step validity */ 1728 if (vi->postcheckstep) { 1729 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1730 if (changed_y) { 1731 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1732 } 1733 if (changed_y || changed_w) { /* recompute the function if the step has changed */ 1734 ierr = SNESComputeFunction(snes,w,g); 1735 if (snes->domainerror) { 1736 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1737 PetscFunctionReturn(0); 1738 } 1739 ierr = VecNormBegin(g,NORM_2,gnorm);CHKERRQ(ierr); 1740 ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr); 1741 ierr = VecNormEnd(g,NORM_2,gnorm);CHKERRQ(ierr); 1742 ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr); 1743 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1744 } 1745 } 1746 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1747 PetscFunctionReturn(0); 1748 } 1749 1750 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*/ 1751 /* -------------------------------------------------------------------------- */ 1752 EXTERN_C_BEGIN 1753 #undef __FUNCT__ 1754 #define __FUNCT__ "SNESLineSearchSet_VI" 1755 PetscErrorCode PETSCSNES_DLLEXPORT SNESLineSearchSet_VI(SNES snes,FCN2 func,void *lsctx) 1756 { 1757 PetscFunctionBegin; 1758 ((SNES_VI *)(snes->data))->LineSearch = func; 1759 ((SNES_VI *)(snes->data))->lsP = lsctx; 1760 PetscFunctionReturn(0); 1761 } 1762 EXTERN_C_END 1763 1764 /* -------------------------------------------------------------------------- */ 1765 EXTERN_C_BEGIN 1766 #undef __FUNCT__ 1767 #define __FUNCT__ "SNESLineSearchSetMonitor_VI" 1768 PetscErrorCode PETSCSNES_DLLEXPORT SNESLineSearchSetMonitor_VI(SNES snes,PetscBool flg) 1769 { 1770 SNES_VI *vi = (SNES_VI*)snes->data; 1771 PetscErrorCode ierr; 1772 1773 PetscFunctionBegin; 1774 if (flg && !vi->lsmonitor) { 1775 ierr = PetscViewerASCIIMonitorCreate(((PetscObject)snes)->comm,"stdout",((PetscObject)snes)->tablevel,&vi->lsmonitor);CHKERRQ(ierr); 1776 } else if (!flg && vi->lsmonitor) { 1777 ierr = PetscViewerASCIIMonitorDestroy(vi->lsmonitor);CHKERRQ(ierr); 1778 } 1779 PetscFunctionReturn(0); 1780 } 1781 EXTERN_C_END 1782 1783 /* 1784 SNESView_VI - Prints info from the SNESVI data structure. 1785 1786 Input Parameters: 1787 . SNES - the SNES context 1788 . viewer - visualization context 1789 1790 Application Interface Routine: SNESView() 1791 */ 1792 #undef __FUNCT__ 1793 #define __FUNCT__ "SNESView_VI" 1794 static PetscErrorCode SNESView_VI(SNES snes,PetscViewer viewer) 1795 { 1796 SNES_VI *vi = (SNES_VI *)snes->data; 1797 const char *cstr,*tstr; 1798 PetscErrorCode ierr; 1799 PetscBool iascii; 1800 1801 PetscFunctionBegin; 1802 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1803 if (iascii) { 1804 if (vi->LineSearch == SNESLineSearchNo_VI) cstr = "SNESLineSearchNo"; 1805 else if (vi->LineSearch == SNESLineSearchQuadratic_VI) cstr = "SNESLineSearchQuadratic"; 1806 else if (vi->LineSearch == SNESLineSearchCubic_VI) cstr = "SNESLineSearchCubic"; 1807 else cstr = "unknown"; 1808 if (snes->ops->solve == SNESSolveVI_SS) tstr = "Semismooth"; 1809 else if (snes->ops->solve == SNESSolveVI_AS) tstr = "Active Set"; 1810 else if (snes->ops->solve == SNESSolveVI_RS) tstr = "Reduced Space"; 1811 else tstr = "unknown"; 1812 ierr = PetscViewerASCIIPrintf(viewer," VI algorithm: %s\n",tstr);CHKERRQ(ierr); 1813 ierr = PetscViewerASCIIPrintf(viewer," line search variant: %s\n",cstr);CHKERRQ(ierr); 1814 ierr = PetscViewerASCIIPrintf(viewer," alpha=%G, maxstep=%G, minlambda=%G\n",vi->alpha,vi->maxstep,vi->minlambda);CHKERRQ(ierr); 1815 } else { 1816 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for SNES EQ VI",((PetscObject)viewer)->type_name); 1817 } 1818 PetscFunctionReturn(0); 1819 } 1820 1821 /* 1822 SNESVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu. 1823 1824 Input Parameters: 1825 . snes - the SNES context. 1826 . xl - lower bound. 1827 . xu - upper bound. 1828 1829 Notes: 1830 If this routine is not called then the lower and upper bounds are set to 1831 PETSC_VI_INF and PETSC_VI_NINF respectively during SNESSetUp(). 1832 1833 */ 1834 1835 #undef __FUNCT__ 1836 #define __FUNCT__ "SNESVISetVariableBounds" 1837 PetscErrorCode SNESVISetVariableBounds(SNES snes, Vec xl, Vec xu) 1838 { 1839 SNES_VI *vi = (SNES_VI*)snes->data; 1840 1841 PetscFunctionBegin; 1842 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1843 PetscValidHeaderSpecific(xl,VEC_CLASSID,2); 1844 PetscValidHeaderSpecific(xu,VEC_CLASSID,3); 1845 1846 /* Check if SNESSetFunction is called */ 1847 if(!snes->vec_func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call SNESSetFunction() first"); 1848 1849 /* Check if the vector sizes are compatible for lower and upper bounds */ 1850 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); 1851 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); 1852 vi->usersetxbounds = PETSC_TRUE; 1853 vi->xl = xl; 1854 vi->xu = xu; 1855 1856 PetscFunctionReturn(0); 1857 } 1858 /* -------------------------------------------------------------------------- */ 1859 /* 1860 SNESSetFromOptions_VI - Sets various parameters for the SNESVI method. 1861 1862 Input Parameter: 1863 . snes - the SNES context 1864 1865 Application Interface Routine: SNESSetFromOptions() 1866 */ 1867 #undef __FUNCT__ 1868 #define __FUNCT__ "SNESSetFromOptions_VI" 1869 static PetscErrorCode SNESSetFromOptions_VI(SNES snes) 1870 { 1871 SNES_VI *vi = (SNES_VI *)snes->data; 1872 const char *lses[] = {"basic","basicnonorms","quadratic","cubic"}; 1873 const char *vies[] = {"ss","as","rs"}; 1874 PetscErrorCode ierr; 1875 PetscInt indx; 1876 PetscBool flg,set,flg2; 1877 1878 PetscFunctionBegin; 1879 ierr = PetscOptionsHead("SNES semismooth method options");CHKERRQ(ierr); 1880 ierr = PetscOptionsReal("-snes_vi_alpha","Function norm must decrease by","None",vi->alpha,&vi->alpha,0);CHKERRQ(ierr); 1881 ierr = PetscOptionsReal("-snes_vi_maxstep","Step must be less than","None",vi->maxstep,&vi->maxstep,0);CHKERRQ(ierr); 1882 ierr = PetscOptionsReal("-snes_vi_minlambda","Minimum lambda allowed","None",vi->minlambda,&vi->minlambda,0);CHKERRQ(ierr); 1883 ierr = PetscOptionsReal("-snes_vi_const_tol","constraint tolerance","None",vi->const_tol,&vi->const_tol,0);CHKERRQ(ierr); 1884 ierr = PetscOptionsBool("-snes_vi_lsmonitor","Print progress of line searches","SNESLineSearchSetMonitor",vi->lsmonitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);CHKERRQ(ierr); 1885 if (set) {ierr = SNESLineSearchSetMonitor(snes,flg);CHKERRQ(ierr);} 1886 ierr = PetscOptionsEList("-snes_vi_type","Semismooth algorithm used","",vies,3,"ss",&indx,&flg2);CHKERRQ(ierr); 1887 if (flg2) { 1888 switch (indx) { 1889 case 0: 1890 snes->ops->solve = SNESSolveVI_SS; 1891 break; 1892 case 1: 1893 snes->ops->solve = SNESSolveVI_AS; 1894 break; 1895 case 2: 1896 snes->ops->solve = SNESSolveVI_RS; 1897 break; 1898 } 1899 } 1900 ierr = PetscOptionsEList("-snes_vi_ls","Line search used","SNESLineSearchSet",lses,4,"cubic",&indx,&flg);CHKERRQ(ierr); 1901 if (flg) { 1902 switch (indx) { 1903 case 0: 1904 ierr = SNESLineSearchSet(snes,SNESLineSearchNo_VI,PETSC_NULL);CHKERRQ(ierr); 1905 break; 1906 case 1: 1907 ierr = SNESLineSearchSet(snes,SNESLineSearchNoNorms_VI,PETSC_NULL);CHKERRQ(ierr); 1908 break; 1909 case 2: 1910 ierr = SNESLineSearchSet(snes,SNESLineSearchQuadratic_VI,PETSC_NULL);CHKERRQ(ierr); 1911 break; 1912 case 3: 1913 ierr = SNESLineSearchSet(snes,SNESLineSearchCubic_VI,PETSC_NULL);CHKERRQ(ierr); 1914 break; 1915 } 1916 } 1917 ierr = PetscOptionsTail();CHKERRQ(ierr); 1918 PetscFunctionReturn(0); 1919 } 1920 /* -------------------------------------------------------------------------- */ 1921 /*MC 1922 SNESVI - Semismooth newton method based nonlinear solver that uses a line search 1923 1924 Options Database: 1925 + -snes_vi [cubic,quadratic,basic,basicnonorms] - Selects line search 1926 . -snes_vi_alpha <alpha> - Sets alpha 1927 . -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) 1928 . -snes_vi_minlambda <minlambda> - Sets the minimum lambda the line search will use minlambda / max_i ( y[i]/x[i] ) 1929 - -snes_vi_monitor - print information about progress of line searches 1930 1931 1932 Level: beginner 1933 1934 .seealso: SNESCreate(), SNES, SNESSetType(), SNESTR, SNESLineSearchSet(), 1935 SNESLineSearchSetPostCheck(), SNESLineSearchNo(), SNESLineSearchCubic(), SNESLineSearchQuadratic(), 1936 SNESLineSearchSet(), SNESLineSearchNoNorms(), SNESLineSearchSetPreCheck(), SNESLineSearchSetParams(), SNESLineSearchGetParams() 1937 1938 M*/ 1939 EXTERN_C_BEGIN 1940 #undef __FUNCT__ 1941 #define __FUNCT__ "SNESCreate_VI" 1942 PetscErrorCode PETSCSNES_DLLEXPORT SNESCreate_VI(SNES snes) 1943 { 1944 PetscErrorCode ierr; 1945 SNES_VI *vi; 1946 1947 PetscFunctionBegin; 1948 snes->ops->setup = SNESSetUp_VI; 1949 snes->ops->solve = SNESSolveVI_SS; 1950 snes->ops->destroy = SNESDestroy_VI; 1951 snes->ops->setfromoptions = SNESSetFromOptions_VI; 1952 snes->ops->view = SNESView_VI; 1953 snes->ops->converged = SNESDefaultConverged_VI; 1954 1955 ierr = PetscNewLog(snes,SNES_VI,&vi);CHKERRQ(ierr); 1956 snes->data = (void*)vi; 1957 vi->alpha = 1.e-4; 1958 vi->maxstep = 1.e8; 1959 vi->minlambda = 1.e-12; 1960 vi->LineSearch = SNESLineSearchCubic_VI; 1961 vi->lsP = PETSC_NULL; 1962 vi->postcheckstep = PETSC_NULL; 1963 vi->postcheck = PETSC_NULL; 1964 vi->precheckstep = PETSC_NULL; 1965 vi->precheck = PETSC_NULL; 1966 vi->const_tol = 2.2204460492503131e-16; 1967 vi->computessfunction = ComputeFischerFunction; 1968 1969 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetMonitor_C","SNESLineSearchSetMonitor_VI",SNESLineSearchSetMonitor_VI);CHKERRQ(ierr); 1970 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSet_C","SNESLineSearchSet_VI",SNESLineSearchSet_VI);CHKERRQ(ierr); 1971 1972 PetscFunctionReturn(0); 1973 } 1974 EXTERN_C_END 1975