Actual source code: viss.c


  2: #include <../src/snes/impls/vi/ss/vissimpl.h>

  4: /*@
  5:   SNESVIComputeMeritFunction - Evaluates the merit function for the mixed complementarity problem.

  7:    Input Parameter:
  8: .  phi - the `Vec` holding the evaluation of the semismooth function

 10:    Output Parameters:
 11: +  merit - the merit function 1/2 ||phi||^2
 12: -  phinorm -  the two-norm of the vector, ||phi||

 14:    Level: developer

 16: .seealso: `SNESVINEWTONSSLS`, `SNESVIComputeFunction()`
 17: @*/
 18: PetscErrorCode SNESVIComputeMeritFunction(Vec phi, PetscReal *merit, PetscReal *phinorm)
 19: {
 20:   VecNormBegin(phi, NORM_2, phinorm);
 21:   VecNormEnd(phi, NORM_2, phinorm);
 22:   *merit = 0.5 * (*phinorm) * (*phinorm);
 23:   return 0;
 24: }

 26: static inline PetscScalar Phi(PetscScalar a, PetscScalar b)
 27: {
 28:   return a + b - PetscSqrtScalar(a * a + b * b);
 29: }

 31: static inline PetscScalar DPhi(PetscScalar a, PetscScalar b)
 32: {
 33:   if ((PetscAbsScalar(a) >= 1.e-6) || (PetscAbsScalar(b) >= 1.e-6)) return 1.0 - a / PetscSqrtScalar(a * a + b * b);
 34:   else return .5;
 35: }

 37: /*@
 38:    SNESVIComputeFunction - Provides the function that reformulates a system of nonlinear equations in mixed complementarity form to a system of nonlinear
 39:    equations in semismooth form.

 41:    Input Parameters:
 42: +  snes - the SNES context
 43: .  X - current iterate
 44: -  functx - user defined function context

 46:    Output Parameter:
 47: .  phi - the evaluation of Semismooth function at X

 49:    Level: developer

 51: .seealso: `SNESVINEWTONSSLS`, `SNESVIComputeMeritFunction()`
 52: @*/
 53: PetscErrorCode SNESVIComputeFunction(SNES snes, Vec X, Vec phi, void *functx)
 54: {
 55:   SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS *)snes->data;
 56:   Vec                Xl = snes->xl, Xu = snes->xu, F = snes->vec_func;
 57:   PetscScalar       *phi_arr, *f_arr, *l, *u;
 58:   const PetscScalar *x_arr;
 59:   PetscInt           i, nlocal;

 61:   (*vi->computeuserfunction)(snes, X, F, functx);
 62:   VecGetLocalSize(X, &nlocal);
 63:   VecGetArrayRead(X, &x_arr);
 64:   VecGetArray(F, &f_arr);
 65:   VecGetArray(Xl, &l);
 66:   VecGetArray(Xu, &u);
 67:   VecGetArray(phi, &phi_arr);

 69:   for (i = 0; i < nlocal; i++) {
 70:     if ((PetscRealPart(l[i]) <= PETSC_NINFINITY) && (PetscRealPart(u[i]) >= PETSC_INFINITY)) { /* no constraints on variable */
 71:       phi_arr[i] = f_arr[i];
 72:     } else if (PetscRealPart(l[i]) <= PETSC_NINFINITY) { /* upper bound on variable only */
 73:       phi_arr[i] = -Phi(u[i] - x_arr[i], -f_arr[i]);
 74:     } else if (PetscRealPart(u[i]) >= PETSC_INFINITY) { /* lower bound on variable only */
 75:       phi_arr[i] = Phi(x_arr[i] - l[i], f_arr[i]);
 76:     } else if (l[i] == u[i]) {
 77:       phi_arr[i] = l[i] - x_arr[i];
 78:     } else { /* both bounds on variable */
 79:       phi_arr[i] = Phi(x_arr[i] - l[i], -Phi(u[i] - x_arr[i], -f_arr[i]));
 80:     }
 81:   }

 83:   VecRestoreArrayRead(X, &x_arr);
 84:   VecRestoreArray(F, &f_arr);
 85:   VecRestoreArray(Xl, &l);
 86:   VecRestoreArray(Xu, &u);
 87:   VecRestoreArray(phi, &phi_arr);
 88:   return 0;
 89: }

 91: /*
 92:    SNESVIComputeBsubdifferentialVectors - Computes the diagonal shift (Da) and row scaling (Db) vectors needed for the
 93:                                           the semismooth jacobian.
 94: */
 95: PetscErrorCode SNESVIComputeBsubdifferentialVectors(SNES snes, Vec X, Vec F, Mat jac, Vec Da, Vec Db)
 96: {
 97:   PetscScalar *l, *u, *x, *f, *da, *db, da1, da2, db1, db2;
 98:   PetscInt     i, nlocal;

100:   VecGetArray(X, &x);
101:   VecGetArray(F, &f);
102:   VecGetArray(snes->xl, &l);
103:   VecGetArray(snes->xu, &u);
104:   VecGetArray(Da, &da);
105:   VecGetArray(Db, &db);
106:   VecGetLocalSize(X, &nlocal);

108:   for (i = 0; i < nlocal; i++) {
109:     if ((PetscRealPart(l[i]) <= PETSC_NINFINITY) && (PetscRealPart(u[i]) >= PETSC_INFINITY)) { /* no constraints on variable */
110:       da[i] = 0;
111:       db[i] = 1;
112:     } else if (PetscRealPart(l[i]) <= PETSC_NINFINITY) { /* upper bound on variable only */
113:       da[i] = DPhi(u[i] - x[i], -f[i]);
114:       db[i] = DPhi(-f[i], u[i] - x[i]);
115:     } else if (PetscRealPart(u[i]) >= PETSC_INFINITY) { /* lower bound on variable only */
116:       da[i] = DPhi(x[i] - l[i], f[i]);
117:       db[i] = DPhi(f[i], x[i] - l[i]);
118:     } else if (l[i] == u[i]) { /* fixed variable */
119:       da[i] = 1;
120:       db[i] = 0;
121:     } else { /* upper and lower bounds on variable */
122:       da1   = DPhi(x[i] - l[i], -Phi(u[i] - x[i], -f[i]));
123:       db1   = DPhi(-Phi(u[i] - x[i], -f[i]), x[i] - l[i]);
124:       da2   = DPhi(u[i] - x[i], -f[i]);
125:       db2   = DPhi(-f[i], u[i] - x[i]);
126:       da[i] = da1 + db1 * da2;
127:       db[i] = db1 * db2;
128:     }
129:   }

131:   VecRestoreArray(X, &x);
132:   VecRestoreArray(F, &f);
133:   VecRestoreArray(snes->xl, &l);
134:   VecRestoreArray(snes->xu, &u);
135:   VecRestoreArray(Da, &da);
136:   VecRestoreArray(Db, &db);
137:   return 0;
138: }

140: /*
141:    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.

143:    Input Parameters:
144: .  Da       - Diagonal shift vector for the semismooth jacobian.
145: .  Db       - Row scaling vector for the semismooth jacobian.

147:    Output Parameters:
148: .  jac      - semismooth jacobian
149: .  jac_pre  - optional preconditioning matrix

151:    Note:
152:    The semismooth jacobian matrix is given by
153:    jac = Da + Db*jacfun
154:    where Db is the row scaling matrix stored as a vector,
155:          Da is the diagonal perturbation matrix stored as a vector
156:    and   jacfun is the jacobian of the original nonlinear function.
157: */
158: PetscErrorCode SNESVIComputeJacobian(Mat jac, Mat jac_pre, Vec Da, Vec Db)
159: {
160:   /* Do row scaling  and add diagonal perturbation */
161:   MatDiagonalScale(jac, Db, NULL);
162:   MatDiagonalSet(jac, Da, ADD_VALUES);
163:   if (jac != jac_pre) { /* If jac and jac_pre are different */
164:     MatDiagonalScale(jac_pre, Db, NULL);
165:     MatDiagonalSet(jac_pre, Da, ADD_VALUES);
166:   }
167:   return 0;
168: }

170: /*
171:    SNESVIComputeMeritFunctionGradient - Computes the gradient of the merit function psi.

173:    Input Parameters:
174:    phi - semismooth function.
175:    H   - semismooth jacobian

177:    Output Parameter:
178:    dpsi - merit function gradient

180:    Note:
181:   The merit function gradient is computed as follows
182:         dpsi = H^T*phi
183: */
184: PetscErrorCode SNESVIComputeMeritFunctionGradient(Mat H, Vec phi, Vec dpsi)
185: {
186:   MatMultTranspose(H, phi, dpsi);
187:   return 0;
188: }

190: /*
191:    SNESSolve_VINEWTONSSLS - Solves the complementarity problem with a semismooth Newton
192:    method using a line search.

194:    Input Parameter:
195: .  snes - the SNES context

197:    Application Interface Routine: SNESSolve()

199:    Note:
200:    This implements essentially a semismooth Newton method with a
201:    line search. The default line search does not do any line search
202:    but rather takes a full Newton step.

204:    Developer Note: the code in this file should be slightly modified so that this routine need not exist and the SNESSolve_NEWTONLS() routine is called directly with the appropriate wrapped function and Jacobian evaluations

206: */
207: PetscErrorCode SNESSolve_VINEWTONSSLS(SNES snes)
208: {
209:   SNES_VINEWTONSSLS   *vi = (SNES_VINEWTONSSLS *)snes->data;
210:   PetscInt             maxits, i, lits;
211:   SNESLineSearchReason lssucceed;
212:   PetscReal            gnorm, xnorm = 0, ynorm;
213:   Vec                  Y, X, F;
214:   KSPConvergedReason   kspreason;
215:   DM                   dm;
216:   DMSNES               sdm;

218:   SNESGetDM(snes, &dm);
219:   DMGetDMSNES(dm, &sdm);

221:   vi->computeuserfunction   = sdm->ops->computefunction;
222:   sdm->ops->computefunction = SNESVIComputeFunction;

224:   snes->numFailures            = 0;
225:   snes->numLinearSolveFailures = 0;
226:   snes->reason                 = SNES_CONVERGED_ITERATING;

228:   maxits = snes->max_its;  /* maximum number of iterations */
229:   X      = snes->vec_sol;  /* solution vector */
230:   F      = snes->vec_func; /* residual vector */
231:   Y      = snes->work[0];  /* work vectors */

233:   PetscObjectSAWsTakeAccess((PetscObject)snes);
234:   snes->iter = 0;
235:   snes->norm = 0.0;
236:   PetscObjectSAWsGrantAccess((PetscObject)snes);

238:   SNESVIProjectOntoBounds(snes, X);
239:   SNESComputeFunction(snes, X, vi->phi);
240:   if (snes->domainerror) {
241:     snes->reason              = SNES_DIVERGED_FUNCTION_DOMAIN;
242:     sdm->ops->computefunction = vi->computeuserfunction;
243:     return 0;
244:   }
245:   /* Compute Merit function */
246:   SNESVIComputeMeritFunction(vi->phi, &vi->merit, &vi->phinorm);

248:   VecNormBegin(X, NORM_2, &xnorm); /* xnorm <- ||x||  */
249:   VecNormEnd(X, NORM_2, &xnorm);
250:   SNESCheckFunctionNorm(snes, vi->merit);

252:   PetscObjectSAWsTakeAccess((PetscObject)snes);
253:   snes->norm = vi->phinorm;
254:   PetscObjectSAWsGrantAccess((PetscObject)snes);
255:   SNESLogConvergenceHistory(snes, vi->phinorm, 0);
256:   SNESMonitor(snes, 0, vi->phinorm);

258:   /* test convergence */
259:   PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, vi->phinorm, &snes->reason, snes->cnvP);
260:   if (snes->reason) {
261:     sdm->ops->computefunction = vi->computeuserfunction;
262:     return 0;
263:   }

265:   for (i = 0; i < maxits; i++) {
266:     /* Call general purpose update function */
267:     PetscTryTypeMethod(snes, update, snes->iter);

269:     /* Solve J Y = Phi, where J is the semismooth jacobian */

271:     /* Get the jacobian -- note that the function must be the original function for snes_fd and snes_fd_color to work for this*/
272:     sdm->ops->computefunction = vi->computeuserfunction;
273:     SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre);
274:     SNESCheckJacobianDomainerror(snes);
275:     sdm->ops->computefunction = SNESVIComputeFunction;

277:     /* Get the diagonal shift and row scaling vectors */
278:     SNESVIComputeBsubdifferentialVectors(snes, X, F, snes->jacobian, vi->Da, vi->Db);
279:     /* Compute the semismooth jacobian */
280:     SNESVIComputeJacobian(snes->jacobian, snes->jacobian_pre, vi->Da, vi->Db);
281:     /* Compute the merit function gradient */
282:     SNESVIComputeMeritFunctionGradient(snes->jacobian, vi->phi, vi->dpsi);
283:     KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre);
284:     KSPSolve(snes->ksp, vi->phi, Y);
285:     KSPGetConvergedReason(snes->ksp, &kspreason);

287:     if (kspreason < 0) {
288:       if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
289:         PetscInfo(snes, "iter=%" PetscInt_FMT ", number linear solve failures %" PetscInt_FMT " greater than current SNES allowed, stopping solve\n", snes->iter, snes->numLinearSolveFailures);
290:         snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
291:         break;
292:       }
293:     }
294:     KSPGetIterationNumber(snes->ksp, &lits);
295:     snes->linear_its += lits;
296:     PetscInfo(snes, "iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits);
297:     /*
298:     if (snes->ops->precheck) {
299:       PetscBool changed_y = PETSC_FALSE;
300:       PetscUseTypeMethod(snes,precheck ,X,Y,snes->precheck,&changed_y);
301:     }

303:     if (PetscLogPrintInfo) SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);
304:     */
305:     /* Compute a (scaled) negative update in the line search routine:
306:          Y <- X - lambda*Y
307:        and evaluate G = function(Y) (depends on the line search).
308:     */
309:     VecCopy(Y, snes->vec_sol_update);
310:     ynorm = 1;
311:     gnorm = vi->phinorm;
312:     SNESLineSearchApply(snes->linesearch, X, vi->phi, &gnorm, Y);
313:     SNESLineSearchGetReason(snes->linesearch, &lssucceed);
314:     SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm);
315:     PetscInfo(snes, "fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n", (double)vi->phinorm, (double)gnorm, (double)ynorm, (int)lssucceed);
316:     if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
317:     if (snes->domainerror) {
318:       snes->reason              = SNES_DIVERGED_FUNCTION_DOMAIN;
319:       sdm->ops->computefunction = vi->computeuserfunction;
320:       return 0;
321:     }
322:     if (lssucceed) {
323:       if (++snes->numFailures >= snes->maxFailures) {
324:         PetscBool ismin;
325:         snes->reason = SNES_DIVERGED_LINE_SEARCH;
326:         SNESVICheckLocalMin_Private(snes, snes->jacobian, vi->phi, X, gnorm, &ismin);
327:         if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
328:         break;
329:       }
330:     }
331:     /* Update function and solution vectors */
332:     vi->phinorm = gnorm;
333:     vi->merit   = 0.5 * vi->phinorm * vi->phinorm;
334:     /* Monitor convergence */
335:     PetscObjectSAWsTakeAccess((PetscObject)snes);
336:     snes->iter  = i + 1;
337:     snes->norm  = vi->phinorm;
338:     snes->xnorm = xnorm;
339:     snes->ynorm = ynorm;
340:     PetscObjectSAWsGrantAccess((PetscObject)snes);
341:     SNESLogConvergenceHistory(snes, snes->norm, lits);
342:     SNESMonitor(snes, snes->iter, snes->norm);
343:     /* Test for convergence, xnorm = || X || */
344:     if (snes->ops->converged != SNESConvergedSkip) VecNorm(X, NORM_2, &xnorm);
345:     PetscUseTypeMethod(snes, converged, snes->iter, xnorm, ynorm, vi->phinorm, &snes->reason, snes->cnvP);
346:     if (snes->reason) break;
347:   }
348:   if (i == maxits) {
349:     PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", maxits);
350:     if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
351:   }
352:   sdm->ops->computefunction = vi->computeuserfunction;
353:   return 0;
354: }

356: /*
357:    SNESSetUp_VINEWTONSSLS - Sets up the internal data structures for the later use
358:    of the SNES nonlinear solver.

360:    Input Parameter:
361: .  snes - the SNES context

363:    Application Interface Routine: SNESSetUp()

365:    Note:
366:    For basic use of the SNES solvers, the user need not explicitly call
367:    SNESSetUp(), since these actions will automatically occur during
368:    the call to SNESSolve().
369:  */
370: PetscErrorCode SNESSetUp_VINEWTONSSLS(SNES snes)
371: {
372:   SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS *)snes->data;

374:   SNESSetUp_VI(snes);
375:   VecDuplicate(snes->vec_sol, &vi->dpsi);
376:   VecDuplicate(snes->vec_sol, &vi->phi);
377:   VecDuplicate(snes->vec_sol, &vi->Da);
378:   VecDuplicate(snes->vec_sol, &vi->Db);
379:   VecDuplicate(snes->vec_sol, &vi->z);
380:   VecDuplicate(snes->vec_sol, &vi->t);
381:   return 0;
382: }

384: PetscErrorCode SNESReset_VINEWTONSSLS(SNES snes)
385: {
386:   SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS *)snes->data;

388:   SNESReset_VI(snes);
389:   VecDestroy(&vi->dpsi);
390:   VecDestroy(&vi->phi);
391:   VecDestroy(&vi->Da);
392:   VecDestroy(&vi->Db);
393:   VecDestroy(&vi->z);
394:   VecDestroy(&vi->t);
395:   return 0;
396: }

398: /*
399:    SNESSetFromOptions_VINEWTONSSLS - Sets various parameters for the SNESVI method.

401:    Input Parameter:
402: .  snes - the SNES context

404:    Application Interface Routine: SNESSetFromOptions()
405: */
406: static PetscErrorCode SNESSetFromOptions_VINEWTONSSLS(SNES snes, PetscOptionItems *PetscOptionsObject)
407: {
408:   SNESSetFromOptions_VI(snes, PetscOptionsObject);
409:   PetscOptionsHeadBegin(PetscOptionsObject, "SNES semismooth method options");
410:   PetscOptionsHeadEnd();
411:   return 0;
412: }

414: /*MC
415:       SNESVINEWTONSSLS - Semi-smooth solver for variational inequalities based on Newton's method

417:    Options Database Keys:
418: +   -snes_type <vinewtonssls,vinewtonrsls> a semi-smooth solver, a reduced space active set method
419: -   -snes_vi_monitor - prints the number of active constraints at each iteration.

421:    Level: beginner

423:    References:
424: +  * -  T. S. Munson, F. Facchinei, M. C. Ferris, A. Fischer, and C. Kanzow. The semismooth
425:      algorithm for large scale complementarity problems. INFORMS Journal on Computing, 13 (2001).
426: -  * -  T. S. Munson, and S. Benson. Flexible Complementarity Solvers for Large Scale
427:      Applications, Optimization Methods and Software, 21 (2006).

429:    Notes:
430:    This family of algorithm is much like an interior point method.

432:    The reduced space active set solvers `SNESVINEWTONRSLS` provide an alternative approach that does not result in extremely ill-conditioned linear systems

434: .seealso: `SNESVINEWTONRSLS`, `SNESVISetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESVINEWTONRSLS`, `SNESNEWTONTR`, `SNESLineSearchSetType()`, `SNESLineSearchSetPostCheck()`, `SNESLineSearchSetPreCheck()`
435: M*/
436: PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONSSLS(SNES snes)
437: {
438:   SNES_VINEWTONSSLS *vi;
439:   SNESLineSearch     linesearch;

441:   snes->ops->reset          = SNESReset_VINEWTONSSLS;
442:   snes->ops->setup          = SNESSetUp_VINEWTONSSLS;
443:   snes->ops->solve          = SNESSolve_VINEWTONSSLS;
444:   snes->ops->destroy        = SNESDestroy_VI;
445:   snes->ops->setfromoptions = SNESSetFromOptions_VINEWTONSSLS;
446:   snes->ops->view           = NULL;

448:   snes->usesksp = PETSC_TRUE;
449:   snes->usesnpc = PETSC_FALSE;

451:   SNESGetLineSearch(snes, &linesearch);
452:   if (!((PetscObject)linesearch)->type_name) {
453:     SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
454:     SNESLineSearchBTSetAlpha(linesearch, 0.0);
455:   }

457:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

459:   PetscNew(&vi);
460:   snes->data = (void *)vi;

462:   PetscObjectComposeFunction((PetscObject)snes, "SNESVISetVariableBounds_C", SNESVISetVariableBounds_VI);
463:   PetscObjectComposeFunction((PetscObject)snes, "SNESVISetComputeVariableBounds_C", SNESVISetComputeVariableBounds_VI);
464:   return 0;
465: }