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 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: {
 20:   VecNormBegin(phi,NORM_2,phinorm);
 21:   VecNormEnd(phi,NORM_2,phinorm);

 23:   *merit = 0.5*(*phinorm)*(*phinorm);
 24:   return 0;
 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 - Reformulates a system of nonlinear equations in mixed complementarity form to a system of nonlinear equations in semismooth form.

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

 46:    Output Parameters:
 47: .  phi - Semismooth function

 49: */
 50: static PetscErrorCode SNESVIComputeFunction(SNES snes,Vec X,Vec phi,void *functx)
 51: {
 52:   SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS*)snes->data;
 53:   Vec               Xl  = snes->xl,Xu = snes->xu,F = snes->vec_func;
 54:   PetscScalar       *phi_arr,*f_arr,*l,*u;
 55:   const PetscScalar *x_arr;
 56:   PetscInt          i,nlocal;

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

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

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

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

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

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

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

137: /*
138:    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.

140:    Input Parameters:
141: .  Da       - Diagonal shift vector for the semismooth jacobian.
142: .  Db       - Row scaling vector for the semismooth jacobian.

144:    Output Parameters:
145: .  jac      - semismooth jacobian
146: .  jac_pre  - optional preconditioning matrix

148:    Notes:
149:    The semismooth jacobian matrix is given by
150:    jac = Da + Db*jacfun
151:    where Db is the row scaling matrix stored as a vector,
152:          Da is the diagonal perturbation matrix stored as a vector
153:    and   jacfun is the jacobian of the original nonlinear function.
154: */
155: PetscErrorCode SNESVIComputeJacobian(Mat jac, Mat jac_pre,Vec Da, Vec Db)
156: {

158:   /* Do row scaling  and add diagonal perturbation */
159:   MatDiagonalScale(jac,Db,NULL);
160:   MatDiagonalSet(jac,Da,ADD_VALUES);
161:   if (jac != jac_pre) { /* If jac and jac_pre are different */
162:     MatDiagonalScale(jac_pre,Db,NULL);
163:     MatDiagonalSet(jac_pre,Da,ADD_VALUES);
164:   }
165:   return 0;
166: }

168: /*
169:    SNESVIComputeMeritFunctionGradient - Computes the gradient of the merit function psi.

171:    Input Parameters:
172:    phi - semismooth function.
173:    H   - semismooth jacobian

175:    Output Parameters:
176:    dpsi - merit function gradient

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

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

192:    Input Parameters:
193: .  snes - the SNES context

195:    Application Interface Routine: SNESSolve()

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

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

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

216:   SNESGetDM(snes,&dm);
217:   DMGetDMSNES(dm,&sdm);

219:   vi->computeuserfunction   = sdm->ops->computefunction;
220:   sdm->ops->computefunction = SNESVIComputeFunction;

222:   snes->numFailures            = 0;
223:   snes->numLinearSolveFailures = 0;
224:   snes->reason                 = SNES_CONVERGED_ITERATING;

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

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

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

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

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

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

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

265:     /* Call general purpose update function */
266:     if (snes->ops->update) {
267:       (*snes->ops->update)(snes, snes->iter);
268:     }

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

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

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

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

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

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

363:    Input Parameter:
364: .  snes - the SNES context

366:    Application Interface Routine: SNESSetUp()

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

377:   SNESSetUp_VI(snes);
378:   VecDuplicate(snes->vec_sol, &vi->dpsi);
379:   VecDuplicate(snes->vec_sol, &vi->phi);
380:   VecDuplicate(snes->vec_sol, &vi->Da);
381:   VecDuplicate(snes->vec_sol, &vi->Db);
382:   VecDuplicate(snes->vec_sol, &vi->z);
383:   VecDuplicate(snes->vec_sol, &vi->t);
384:   return 0;
385: }
386: /* -------------------------------------------------------------------------- */
387: PetscErrorCode SNESReset_VINEWTONSSLS(SNES snes)
388: {
389:   SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS*) snes->data;

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

401: /* -------------------------------------------------------------------------- */
402: /*
403:    SNESSetFromOptions_VINEWTONSSLS - Sets various parameters for the SNESVI method.

405:    Input Parameter:
406: .  snes - the SNES context

408:    Application Interface Routine: SNESSetFromOptions()
409: */
410: static PetscErrorCode SNESSetFromOptions_VINEWTONSSLS(PetscOptionItems *PetscOptionsObject,SNES snes)
411: {
412:   SNESSetFromOptions_VI(PetscOptionsObject,snes);
413:   PetscOptionsHead(PetscOptionsObject,"SNES semismooth method options");
414:   PetscOptionsTail();
415:   return 0;
416: }

418: /* -------------------------------------------------------------------------- */
419: /*MC
420:       SNESVINEWTONSSLS - Semi-smooth solver for variational inequalities based on Newton's method

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

426:    Level: beginner

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

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

436: M*/
437: PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONSSLS(SNES snes)
438: {
439:   SNES_VINEWTONSSLS *vi;
440:   SNESLineSearch    linesearch;

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

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

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

458:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

460:   PetscNewLog(snes,&vi);
461:   snes->data = (void*)vi;

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