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:   PetscFunctionBegin;
 21:   PetscCall(VecNormBegin(phi, NORM_2, phinorm));
 22:   PetscCall(VecNormEnd(phi, NORM_2, phinorm));
 23:   *merit = 0.5 * (*phinorm) * (*phinorm);
 24:   PetscFunctionReturn(PETSC_SUCCESS);
 25: }

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

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

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

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

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

 50:    Level: developer

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

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

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

 85:   PetscCall(VecRestoreArrayRead(X, &x_arr));
 86:   PetscCall(VecRestoreArray(F, &f_arr));
 87:   PetscCall(VecRestoreArray(Xl, &l));
 88:   PetscCall(VecRestoreArray(Xu, &u));
 89:   PetscCall(VecRestoreArray(phi, &phi_arr));
 90:   PetscFunctionReturn(PETSC_SUCCESS);
 91: }

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

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

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

134:   PetscCall(VecRestoreArray(X, &x));
135:   PetscCall(VecRestoreArray(F, &f));
136:   PetscCall(VecRestoreArray(snes->xl, &l));
137:   PetscCall(VecRestoreArray(snes->xu, &u));
138:   PetscCall(VecRestoreArray(Da, &da));
139:   PetscCall(VecRestoreArray(Db, &db));
140:   PetscFunctionReturn(PETSC_SUCCESS);
141: }

143: /*
144:    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.

146:    Input Parameters:
147: .  Da       - Diagonal shift vector for the semismooth jacobian.
148: .  Db       - Row scaling vector for the semismooth jacobian.

150:    Output Parameters:
151: .  jac      - semismooth jacobian
152: .  jac_pre  - optional preconditioning matrix

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

174: /*
175:    SNESVIComputeMeritFunctionGradient - Computes the gradient of the merit function psi.

177:    Input Parameters:
178:    phi - semismooth function.
179:    H   - semismooth jacobian

181:    Output Parameter:
182:    dpsi - merit function gradient

184:    Note:
185:   The merit function gradient is computed as follows
186:         dpsi = H^T*phi
187: */
188: PetscErrorCode SNESVIComputeMeritFunctionGradient(Mat H, Vec phi, Vec dpsi)
189: {
190:   PetscFunctionBegin;
191:   PetscCall(MatMultTranspose(H, phi, dpsi));
192:   PetscFunctionReturn(PETSC_SUCCESS);
193: }

195: /*
196:    SNESSolve_VINEWTONSSLS - Solves the complementarity problem with a semismooth Newton
197:    method using a line search.

199:    Input Parameter:
200: .  snes - the SNES context

202:    Application Interface Routine: SNESSolve()

204:    Note:
205:    This implements essentially a semismooth Newton method with a
206:    line search. The default line search does not do any line search
207:    but rather takes a full Newton step.

209:    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

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

223:   PetscFunctionBegin;
224:   PetscCall(SNESGetDM(snes, &dm));
225:   PetscCall(DMGetDMSNES(dm, &sdm));

227:   vi->computeuserfunction   = sdm->ops->computefunction;
228:   sdm->ops->computefunction = SNESVIComputeFunction;

230:   snes->numFailures            = 0;
231:   snes->numLinearSolveFailures = 0;
232:   snes->reason                 = SNES_CONVERGED_ITERATING;

234:   maxits = snes->max_its;  /* maximum number of iterations */
235:   X      = snes->vec_sol;  /* solution vector */
236:   F      = snes->vec_func; /* residual vector */
237:   Y      = snes->work[0];  /* work vectors */

239:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
240:   snes->iter = 0;
241:   snes->norm = 0.0;
242:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));

244:   PetscCall(SNESVIProjectOntoBounds(snes, X));
245:   PetscCall(SNESComputeFunction(snes, X, vi->phi));
246:   if (snes->domainerror) {
247:     snes->reason              = SNES_DIVERGED_FUNCTION_DOMAIN;
248:     sdm->ops->computefunction = vi->computeuserfunction;
249:     PetscFunctionReturn(PETSC_SUCCESS);
250:   }
251:   /* Compute Merit function */
252:   PetscCall(SNESVIComputeMeritFunction(vi->phi, &vi->merit, &vi->phinorm));

254:   PetscCall(VecNormBegin(X, NORM_2, &xnorm)); /* xnorm <- ||x||  */
255:   PetscCall(VecNormEnd(X, NORM_2, &xnorm));
256:   SNESCheckFunctionNorm(snes, vi->merit);

258:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
259:   snes->norm = vi->phinorm;
260:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
261:   PetscCall(SNESLogConvergenceHistory(snes, vi->phinorm, 0));
262:   PetscCall(SNESMonitor(snes, 0, vi->phinorm));

264:   /* test convergence */
265:   PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, vi->phinorm, &snes->reason, snes->cnvP);
266:   if (snes->reason) {
267:     sdm->ops->computefunction = vi->computeuserfunction;
268:     PetscFunctionReturn(PETSC_SUCCESS);
269:   }

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

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

277:     /* Get the jacobian -- note that the function must be the original function for snes_fd and snes_fd_color to work for this*/
278:     sdm->ops->computefunction = vi->computeuserfunction;
279:     PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));
280:     SNESCheckJacobianDomainerror(snes);
281:     sdm->ops->computefunction = SNESVIComputeFunction;

283:     /* Get the diagonal shift and row scaling vectors */
284:     PetscCall(SNESVIComputeBsubdifferentialVectors(snes, X, F, snes->jacobian, vi->Da, vi->Db));
285:     /* Compute the semismooth jacobian */
286:     PetscCall(SNESVIComputeJacobian(snes->jacobian, snes->jacobian_pre, vi->Da, vi->Db));
287:     /* Compute the merit function gradient */
288:     PetscCall(SNESVIComputeMeritFunctionGradient(snes->jacobian, vi->phi, vi->dpsi));
289:     PetscCall(KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre));
290:     PetscCall(KSPSolve(snes->ksp, vi->phi, Y));
291:     PetscCall(KSPGetConvergedReason(snes->ksp, &kspreason));

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

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

362: /*
363:    SNESSetUp_VINEWTONSSLS - Sets up the internal data structures for the later use
364:    of the SNES nonlinear solver.

366:    Input Parameter:
367: .  snes - the SNES context

369:    Application Interface Routine: SNESSetUp()

371:    Note:
372:    For basic use of the SNES solvers, the user need not explicitly call
373:    SNESSetUp(), since these actions will automatically occur during
374:    the call to SNESSolve().
375:  */
376: PetscErrorCode SNESSetUp_VINEWTONSSLS(SNES snes)
377: {
378:   SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS *)snes->data;

380:   PetscFunctionBegin;
381:   PetscCall(SNESSetUp_VI(snes));
382:   PetscCall(VecDuplicate(snes->vec_sol, &vi->dpsi));
383:   PetscCall(VecDuplicate(snes->vec_sol, &vi->phi));
384:   PetscCall(VecDuplicate(snes->vec_sol, &vi->Da));
385:   PetscCall(VecDuplicate(snes->vec_sol, &vi->Db));
386:   PetscCall(VecDuplicate(snes->vec_sol, &vi->z));
387:   PetscCall(VecDuplicate(snes->vec_sol, &vi->t));
388:   PetscFunctionReturn(PETSC_SUCCESS);
389: }

391: PetscErrorCode SNESReset_VINEWTONSSLS(SNES snes)
392: {
393:   SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS *)snes->data;

395:   PetscFunctionBegin;
396:   PetscCall(SNESReset_VI(snes));
397:   PetscCall(VecDestroy(&vi->dpsi));
398:   PetscCall(VecDestroy(&vi->phi));
399:   PetscCall(VecDestroy(&vi->Da));
400:   PetscCall(VecDestroy(&vi->Db));
401:   PetscCall(VecDestroy(&vi->z));
402:   PetscCall(VecDestroy(&vi->t));
403:   PetscFunctionReturn(PETSC_SUCCESS);
404: }

406: /*
407:    SNESSetFromOptions_VINEWTONSSLS - Sets various parameters for the SNESVI method.

409:    Input Parameter:
410: .  snes - the SNES context

412:    Application Interface Routine: SNESSetFromOptions()
413: */
414: static PetscErrorCode SNESSetFromOptions_VINEWTONSSLS(SNES snes, PetscOptionItems *PetscOptionsObject)
415: {
416:   PetscFunctionBegin;
417:   PetscCall(SNESSetFromOptions_VI(snes, PetscOptionsObject));
418:   PetscOptionsHeadBegin(PetscOptionsObject, "SNES semismooth method options");
419:   PetscOptionsHeadEnd();
420:   PetscFunctionReturn(PETSC_SUCCESS);
421: }

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

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

430:    Level: beginner

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

438:    Notes:
439:    This family of algorithm is much like an interior point method.

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

443: .seealso: `SNESVINEWTONRSLS`, `SNESVISetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESVINEWTONRSLS`, `SNESNEWTONTR`, `SNESLineSearchSetType()`, `SNESLineSearchSetPostCheck()`, `SNESLineSearchSetPreCheck()`
444: M*/
445: PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONSSLS(SNES snes)
446: {
447:   SNES_VINEWTONSSLS *vi;
448:   SNESLineSearch     linesearch;

450:   PetscFunctionBegin;
451:   snes->ops->reset          = SNESReset_VINEWTONSSLS;
452:   snes->ops->setup          = SNESSetUp_VINEWTONSSLS;
453:   snes->ops->solve          = SNESSolve_VINEWTONSSLS;
454:   snes->ops->destroy        = SNESDestroy_VI;
455:   snes->ops->setfromoptions = SNESSetFromOptions_VINEWTONSSLS;
456:   snes->ops->view           = NULL;

458:   snes->usesksp = PETSC_TRUE;
459:   snes->usesnpc = PETSC_FALSE;

461:   PetscCall(SNESGetLineSearch(snes, &linesearch));
462:   if (!((PetscObject)linesearch)->type_name) {
463:     PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT));
464:     PetscCall(SNESLineSearchBTSetAlpha(linesearch, 0.0));
465:   }

467:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

469:   PetscCall(PetscNew(&vi));
470:   snes->data = (void *)vi;

472:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESVISetVariableBounds_C", SNESVISetVariableBounds_VI));
473:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESVISetComputeVariableBounds_C", SNESVISetComputeVariableBounds_VI));
474:   PetscFunctionReturn(PETSC_SUCCESS);
475: }