Actual source code: viss.c

petsc-3.15.0 2021-04-05
Report Typos and Errors
```
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 semismooth function

10:   Output Parameter:
11: . merit - the merit function
12: . phinorm - ||phi||

14:   Notes:
15:   The merit function for the mixed complementarity problem is defined as
16:      merit = 0.5*phi^T*phi
17: */
18: static PetscErrorCode SNESVIComputeMeritFunction(Vec phi, PetscReal *merit,PetscReal *phinorm)
19: {

23:   VecNormBegin(phi,NORM_2,phinorm);
24:   VecNormEnd(phi,NORM_2,phinorm);

26:   *merit = 0.5*(*phinorm)*(*phinorm);
27:   return(0);
28: }

30: PETSC_STATIC_INLINE PetscScalar Phi(PetscScalar a,PetscScalar b)
31: {
32:   return a + b - PetscSqrtScalar(a*a + b*b);
33: }

35: PETSC_STATIC_INLINE PetscScalar DPhi(PetscScalar a,PetscScalar b)
36: {
37:   if ((PetscAbsScalar(a) >= 1.e-6) || (PetscAbsScalar(b) >= 1.e-6)) return 1.0 - a/ PetscSqrtScalar(a*a + b*b);
38:   else return .5;
39: }

41: /*
42:    SNESVIComputeFunction - Reformulates a system of nonlinear equations in mixed complementarity form to a system of nonlinear equations in semismooth form.

44:    Input Parameters:
45: .  snes - the SNES context
46: .  X - current iterate
47: .  functx - user defined function context

49:    Output Parameters:
50: .  phi - Semismooth function

52: */
53: static PetscErrorCode SNESVIComputeFunction(SNES snes,Vec X,Vec phi,void *functx)
54: {
55:   PetscErrorCode    ierr;
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;

63:   (*vi->computeuserfunction)(snes,X,F,functx);
64:   VecGetLocalSize(X,&nlocal);
66:   VecGetArray(F,&f_arr);
67:   VecGetArray(Xl,&l);
68:   VecGetArray(Xu,&u);
69:   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:   }

86:   VecRestoreArray(F,&f_arr);
87:   VecRestoreArray(Xl,&l);
88:   VecRestoreArray(Xu,&u);
89:   VecRestoreArray(phi,&phi_arr);
90:   return(0);
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: {
100:   PetscScalar    *l,*u,*x,*f,*da,*db,da1,da2,db1,db2;
101:   PetscInt       i,nlocal;

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

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

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

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

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

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

155:    Notes:
156:    The semismooth jacobian matrix is given by
157:    jac = Da + Db*jacfun
158:    where Db is the row scaling matrix stored as a vector,
159:          Da is the diagonal perturbation matrix stored as a vector
160:    and   jacfun is the jacobian of the original nonlinear function.
161: */
162: PetscErrorCode SNESVIComputeJacobian(Mat jac, Mat jac_pre,Vec Da, Vec Db)
163: {

166:   /* Do row scaling  and add diagonal perturbation */
167:   MatDiagonalScale(jac,Db,NULL);
169:   if (jac != jac_pre) { /* If jac and jac_pre are different */
170:     MatDiagonalScale(jac_pre,Db,NULL);
172:   }
173:   return(0);
174: }

176: /*

179:    Input Parameters:
180:    phi - semismooth function.
181:    H   - semismooth jacobian

183:    Output Parameters:
184:    dpsi - merit function gradient

186:    Notes:
187:   The merit function gradient is computed as follows
188:         dpsi = H^T*phi
189: */
190: PetscErrorCode SNESVIComputeMeritFunctionGradient(Mat H, Vec phi, Vec dpsi)
191: {

195:   MatMultTranspose(H,phi,dpsi);
196:   return(0);
197: }

201: /*
202:    SNESSolve_VINEWTONSSLS - Solves the complementarity problem with a semismooth Newton
203:    method using a line search.

205:    Input Parameters:
206: .  snes - the SNES context

208:    Application Interface Routine: SNESSolve()

210:    Notes:
211:    This implements essentially a semismooth Newton method with a
212:    line search. The default line search does not do any line search
213:    but rather takes a full Newton step.

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

217: */
218: PetscErrorCode SNESSolve_VINEWTONSSLS(SNES snes)
219: {
220:   SNES_VINEWTONSSLS    *vi = (SNES_VINEWTONSSLS*)snes->data;
221:   PetscErrorCode       ierr;
222:   PetscInt             maxits,i,lits;
223:   SNESLineSearchReason lssucceed;
224:   PetscReal            gnorm,xnorm=0,ynorm;
225:   Vec                  Y,X,F;
226:   KSPConvergedReason   kspreason;
227:   DM                   dm;
228:   DMSNES               sdm;

231:   SNESGetDM(snes,&dm);
232:   DMGetDMSNES(dm,&sdm);

234:   vi->computeuserfunction   = sdm->ops->computefunction;
235:   sdm->ops->computefunction = SNESVIComputeFunction;

237:   snes->numFailures            = 0;
238:   snes->numLinearSolveFailures = 0;
239:   snes->reason                 = SNES_CONVERGED_ITERATING;

241:   maxits = snes->max_its;               /* maximum number of iterations */
242:   X      = snes->vec_sol;               /* solution vector */
243:   F      = snes->vec_func;              /* residual vector */
244:   Y      = snes->work[0];               /* work vectors */

246:   PetscObjectSAWsTakeAccess((PetscObject)snes);
247:   snes->iter = 0;
248:   snes->norm = 0.0;
249:   PetscObjectSAWsGrantAccess((PetscObject)snes);

251:   SNESVIProjectOntoBounds(snes,X);
252:   SNESComputeFunction(snes,X,vi->phi);
253:   if (snes->domainerror) {
254:     snes->reason              = SNES_DIVERGED_FUNCTION_DOMAIN;
255:     sdm->ops->computefunction = vi->computeuserfunction;
256:     return(0);
257:   }
258:   /* Compute Merit function */
259:   SNESVIComputeMeritFunction(vi->phi,&vi->merit,&vi->phinorm);

261:   VecNormBegin(X,NORM_2,&xnorm);        /* xnorm <- ||x||  */
262:   VecNormEnd(X,NORM_2,&xnorm);
263:   SNESCheckFunctionNorm(snes,vi->merit);

265:   PetscObjectSAWsTakeAccess((PetscObject)snes);
266:   snes->norm = vi->phinorm;
267:   PetscObjectSAWsGrantAccess((PetscObject)snes);
268:   SNESLogConvergenceHistory(snes,vi->phinorm,0);
269:   SNESMonitor(snes,0,vi->phinorm);

271:   /* test convergence */
272:   (*snes->ops->converged)(snes,0,0.0,0.0,vi->phinorm,&snes->reason,snes->cnvP);
273:   if (snes->reason) {
274:     sdm->ops->computefunction = vi->computeuserfunction;
275:     return(0);
276:   }

278:   for (i=0; i<maxits; i++) {

280:     /* Call general purpose update function */
281:     if (snes->ops->update) {
282:       (*snes->ops->update)(snes, snes->iter);
283:     }

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

287:     /* Get the jacobian -- note that the function must be the original function for snes_fd and snes_fd_color to work for this*/
288:     sdm->ops->computefunction = vi->computeuserfunction;
289:     SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
290:     SNESCheckJacobianDomainerror(snes);
291:     sdm->ops->computefunction = SNESVIComputeFunction;

293:     /* Get the diagonal shift and row scaling vectors */
294:     SNESVIComputeBsubdifferentialVectors(snes,X,F,snes->jacobian,vi->Da,vi->Db);
295:     /* Compute the semismooth jacobian */
296:     SNESVIComputeJacobian(snes->jacobian,snes->jacobian_pre,vi->Da,vi->Db);
297:     /* Compute the merit function gradient */
299:     KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);
300:     KSPSolve(snes->ksp,vi->phi,Y);
301:     KSPGetConvergedReason(snes->ksp,&kspreason);

303:     if (kspreason < 0) {
304:       if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
305:         PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);
306:         snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
307:         break;
308:       }
309:     }
310:     KSPGetIterationNumber(snes->ksp,&lits);
311:     snes->linear_its += lits;
312:     PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);
313:     /*
314:     if (snes->ops->precheck) {
315:       PetscBool changed_y = PETSC_FALSE;
316:       (*snes->ops->precheck)(snes,X,Y,snes->precheck,&changed_y);
317:     }

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

373: /* -------------------------------------------------------------------------- */
374: /*
375:    SNESSetUp_VINEWTONSSLS - Sets up the internal data structures for the later use
376:    of the SNES nonlinear solver.

378:    Input Parameter:
379: .  snes - the SNES context

381:    Application Interface Routine: SNESSetUp()

383:    Notes:
384:    For basic use of the SNES solvers, the user need not explicitly call
385:    SNESSetUp(), since these actions will automatically occur during
386:    the call to SNESSolve().
387:  */
388: PetscErrorCode SNESSetUp_VINEWTONSSLS(SNES snes)
389: {
390:   PetscErrorCode    ierr;
391:   SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS*) snes->data;

394:   SNESSetUp_VI(snes);
395:   VecDuplicate(snes->vec_sol, &vi->dpsi);
396:   VecDuplicate(snes->vec_sol, &vi->phi);
397:   VecDuplicate(snes->vec_sol, &vi->Da);
398:   VecDuplicate(snes->vec_sol, &vi->Db);
399:   VecDuplicate(snes->vec_sol, &vi->z);
400:   VecDuplicate(snes->vec_sol, &vi->t);
401:   return(0);
402: }
403: /* -------------------------------------------------------------------------- */
404: PetscErrorCode SNESReset_VINEWTONSSLS(SNES snes)
405: {
406:   SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS*) snes->data;
407:   PetscErrorCode    ierr;

410:   SNESReset_VI(snes);
411:   VecDestroy(&vi->dpsi);
412:   VecDestroy(&vi->phi);
413:   VecDestroy(&vi->Da);
414:   VecDestroy(&vi->Db);
415:   VecDestroy(&vi->z);
416:   VecDestroy(&vi->t);
417:   return(0);
418: }

420: /* -------------------------------------------------------------------------- */
421: /*
422:    SNESSetFromOptions_VINEWTONSSLS - Sets various parameters for the SNESVI method.

424:    Input Parameter:
425: .  snes - the SNES context

427:    Application Interface Routine: SNESSetFromOptions()
428: */
429: static PetscErrorCode SNESSetFromOptions_VINEWTONSSLS(PetscOptionItems *PetscOptionsObject,SNES snes)
430: {

434:   SNESSetFromOptions_VI(PetscOptionsObject,snes);
436:   PetscOptionsTail();
437:   return(0);
438: }

441: /* -------------------------------------------------------------------------- */
442: /*MC
443:       SNESVINEWTONSSLS - Semi-smooth solver for variational inequalities based on Newton's method

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

449:    Level: beginner

451:    References:
452: +  1. -  T. S. Munson, F. Facchinei, M. C. Ferris, A. Fischer, and C. Kanzow. The semismooth
453:      algorithm for large scale complementarity problems. INFORMS Journal on Computing, 13 (2001).
454: -  2. -  T. S. Munson, and S. Benson. Flexible Complementarity Solvers for Large Scale
455:      Applications, Optimization Methods and Software, 21 (2006).

457: .seealso:  SNESVISetVariableBounds(), SNESVISetComputeVariableBounds(), SNESCreate(), SNES, SNESSetType(), SNESVINEWTONRSLS, SNESNEWTONTR, SNESLineSearchSetType(),SNESLineSearchSetPostCheck(), SNESLineSearchSetPreCheck()

459: M*/
460: PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONSSLS(SNES snes)
461: {
462:   PetscErrorCode    ierr;
463:   SNES_VINEWTONSSLS *vi;
464:   SNESLineSearch    linesearch;

467:   snes->ops->reset          = SNESReset_VINEWTONSSLS;
468:   snes->ops->setup          = SNESSetUp_VINEWTONSSLS;
469:   snes->ops->solve          = SNESSolve_VINEWTONSSLS;
470:   snes->ops->destroy        = SNESDestroy_VI;
471:   snes->ops->setfromoptions = SNESSetFromOptions_VINEWTONSSLS;
472:   snes->ops->view           = NULL;

474:   snes->usesksp = PETSC_TRUE;
475:   snes->usesnpc = PETSC_FALSE;

477:   SNESGetLineSearch(snes, &linesearch);
478:   if (!((PetscObject)linesearch)->type_name) {
479:     SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
480:     SNESLineSearchBTSetAlpha(linesearch, 0.0);
481:   }

483:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

485:   PetscNewLog(snes,&vi);
486:   snes->data = (void*)vi;

488:   PetscObjectComposeFunction((PetscObject)snes,"SNESVISetVariableBounds_C",SNESVISetVariableBounds_VI);
489:   PetscObjectComposeFunction((PetscObject)snes,"SNESVISetComputeVariableBounds_C",SNESVISetComputeVariableBounds_VI);
490:   return(0);
491: }

```