Actual source code: viss.c
petsc-3.5.4 2015-05-23
2: #include <../src/snes/impls/vi/ss/vissimpl.h> /*I "petscsnes.h" I*/
3: #include <../include/petsc-private/kspimpl.h>
4: #include <../include/petsc-private/matimpl.h>
5: #include <../include/petsc-private/dmimpl.h>
8: /*
9: SNESVIComputeMeritFunction - Evaluates the merit function for the mixed complementarity problem.
11: Input Parameter:
12: . phi - the semismooth function
14: Output Parameter:
15: . merit - the merit function
16: . phinorm - ||phi||
18: Notes:
19: The merit function for the mixed complementarity problem is defined as
20: merit = 0.5*phi^T*phi
21: */
24: static PetscErrorCode SNESVIComputeMeritFunction(Vec phi, PetscReal *merit,PetscReal *phinorm)
25: {
29: VecNormBegin(phi,NORM_2,phinorm);
30: VecNormEnd(phi,NORM_2,phinorm);
32: *merit = 0.5*(*phinorm)*(*phinorm);
33: return(0);
34: }
36: PETSC_STATIC_INLINE PetscScalar Phi(PetscScalar a,PetscScalar b)
37: {
38: return a + b - PetscSqrtScalar(a*a + b*b);
39: }
41: PETSC_STATIC_INLINE PetscScalar DPhi(PetscScalar a,PetscScalar b)
42: {
43: if ((PetscAbsScalar(a) >= 1.e-6) || (PetscAbsScalar(b) >= 1.e-6)) return 1.0 - a/ PetscSqrtScalar(a*a + b*b);
44: else return .5;
45: }
47: /*
48: SNESVIComputeFunction - Reformulates a system of nonlinear equations in mixed complementarity form to a system of nonlinear equations in semismooth form.
50: Input Parameters:
51: . snes - the SNES context
52: . x - current iterate
53: . functx - user defined function context
55: Output Parameters:
56: . phi - Semismooth function
58: */
61: static PetscErrorCode SNESVIComputeFunction(SNES snes,Vec X,Vec phi,void *functx)
62: {
63: PetscErrorCode ierr;
64: SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS*)snes->data;
65: Vec Xl = snes->xl,Xu = snes->xu,F = snes->vec_func;
66: PetscScalar *phi_arr,*x_arr,*f_arr,*l,*u;
67: PetscInt i,nlocal;
70: (*vi->computeuserfunction)(snes,X,F,functx);
71: VecGetLocalSize(X,&nlocal);
72: VecGetArray(X,&x_arr);
73: VecGetArray(F,&f_arr);
74: VecGetArray(Xl,&l);
75: VecGetArray(Xu,&u);
76: VecGetArray(phi,&phi_arr);
78: for (i=0; i < nlocal; i++) {
79: if ((PetscRealPart(l[i]) <= PETSC_NINFINITY) && (PetscRealPart(u[i]) >= PETSC_INFINITY)) { /* no constraints on variable */
80: phi_arr[i] = f_arr[i];
81: } else if (PetscRealPart(l[i]) <= PETSC_NINFINITY) { /* upper bound on variable only */
82: phi_arr[i] = -Phi(u[i] - x_arr[i],-f_arr[i]);
83: } else if (PetscRealPart(u[i]) >= PETSC_INFINITY) { /* lower bound on variable only */
84: phi_arr[i] = Phi(x_arr[i] - l[i],f_arr[i]);
85: } else if (l[i] == u[i]) {
86: phi_arr[i] = l[i] - x_arr[i];
87: } else { /* both bounds on variable */
88: phi_arr[i] = Phi(x_arr[i] - l[i],-Phi(u[i] - x_arr[i],-f_arr[i]));
89: }
90: }
92: VecRestoreArray(X,&x_arr);
93: VecRestoreArray(F,&f_arr);
94: VecRestoreArray(Xl,&l);
95: VecRestoreArray(Xu,&u);
96: VecRestoreArray(phi,&phi_arr);
97: return(0);
98: }
100: /*
101: SNESVIComputeBsubdifferentialVectors - Computes the diagonal shift (Da) and row scaling (Db) vectors needed for the
102: the semismooth jacobian.
103: */
106: PetscErrorCode SNESVIComputeBsubdifferentialVectors(SNES snes,Vec X,Vec F,Mat jac,Vec Da,Vec Db)
107: {
109: PetscScalar *l,*u,*x,*f,*da,*db,da1,da2,db1,db2;
110: PetscInt i,nlocal;
113: VecGetArray(X,&x);
114: VecGetArray(F,&f);
115: VecGetArray(snes->xl,&l);
116: VecGetArray(snes->xu,&u);
117: VecGetArray(Da,&da);
118: VecGetArray(Db,&db);
119: VecGetLocalSize(X,&nlocal);
121: for (i=0; i< nlocal; i++) {
122: if ((PetscRealPart(l[i]) <= PETSC_NINFINITY) && (PetscRealPart(u[i]) >= PETSC_INFINITY)) { /* no constraints on variable */
123: da[i] = 0;
124: db[i] = 1;
125: } else if (PetscRealPart(l[i]) <= PETSC_NINFINITY) { /* upper bound on variable only */
126: da[i] = DPhi(u[i] - x[i], -f[i]);
127: db[i] = DPhi(-f[i],u[i] - x[i]);
128: } else if (PetscRealPart(u[i]) >= PETSC_INFINITY) { /* lower bound on variable only */
129: da[i] = DPhi(x[i] - l[i], f[i]);
130: db[i] = DPhi(f[i],x[i] - l[i]);
131: } else if (l[i] == u[i]) { /* fixed variable */
132: da[i] = 1;
133: db[i] = 0;
134: } else { /* upper and lower bounds on variable */
135: da1 = DPhi(x[i] - l[i], -Phi(u[i] - x[i], -f[i]));
136: db1 = DPhi(-Phi(u[i] - x[i], -f[i]),x[i] - l[i]);
137: da2 = DPhi(u[i] - x[i], -f[i]);
138: db2 = DPhi(-f[i],u[i] - x[i]);
139: da[i] = da1 + db1*da2;
140: db[i] = db1*db2;
141: }
142: }
144: VecRestoreArray(X,&x);
145: VecRestoreArray(F,&f);
146: VecRestoreArray(snes->xl,&l);
147: VecRestoreArray(snes->xu,&u);
148: VecRestoreArray(Da,&da);
149: VecRestoreArray(Db,&db);
150: return(0);
151: }
153: /*
154: 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.
156: Input Parameters:
157: . Da - Diagonal shift vector for the semismooth jacobian.
158: . Db - Row scaling vector for the semismooth jacobian.
160: Output Parameters:
161: . jac - semismooth jacobian
162: . jac_pre - optional preconditioning matrix
164: Notes:
165: The semismooth jacobian matrix is given by
166: jac = Da + Db*jacfun
167: where Db is the row scaling matrix stored as a vector,
168: Da is the diagonal perturbation matrix stored as a vector
169: and jacfun is the jacobian of the original nonlinear function.
170: */
173: PetscErrorCode SNESVIComputeJacobian(Mat jac, Mat jac_pre,Vec Da, Vec Db)
174: {
177: /* Do row scaling and add diagonal perturbation */
178: MatDiagonalScale(jac,Db,NULL);
179: MatDiagonalSet(jac,Da,ADD_VALUES);
180: if (jac != jac_pre) { /* If jac and jac_pre are different */
181: MatDiagonalScale(jac_pre,Db,NULL);
182: MatDiagonalSet(jac_pre,Da,ADD_VALUES);
183: }
184: return(0);
185: }
187: /*
188: SNESVIComputeMeritFunctionGradient - Computes the gradient of the merit function psi.
190: Input Parameters:
191: phi - semismooth function.
192: H - semismooth jacobian
194: Output Parameters:
195: dpsi - merit function gradient
197: Notes:
198: The merit function gradient is computed as follows
199: dpsi = H^T*phi
200: */
203: PetscErrorCode SNESVIComputeMeritFunctionGradient(Mat H, Vec phi, Vec dpsi)
204: {
208: MatMultTranspose(H,phi,dpsi);
209: return(0);
210: }
214: /*
215: SNESSolve_VINEWTONSSLS - Solves the complementarity problem with a semismooth Newton
216: method using a line search.
218: Input Parameters:
219: . snes - the SNES context
221: Output Parameter:
222: . outits - number of iterations until termination
224: Application Interface Routine: SNESSolve()
226: Notes:
227: This implements essentially a semismooth Newton method with a
228: line search. The default line search does not do any line seach
229: but rather takes a full newton step.
231: 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
233: */
236: PetscErrorCode SNESSolve_VINEWTONSSLS(SNES snes)
237: {
238: SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS*)snes->data;
239: PetscErrorCode ierr;
240: PetscInt maxits,i,lits;
241: PetscBool lssucceed;
242: PetscReal gnorm,xnorm=0,ynorm;
243: Vec Y,X,F;
244: KSPConvergedReason kspreason;
245: DM dm;
246: DMSNES sdm;
249: SNESGetDM(snes,&dm);
250: DMGetDMSNES(dm,&sdm);
252: vi->computeuserfunction = sdm->ops->computefunction;
253: sdm->ops->computefunction = SNESVIComputeFunction;
255: snes->numFailures = 0;
256: snes->numLinearSolveFailures = 0;
257: snes->reason = SNES_CONVERGED_ITERATING;
259: maxits = snes->max_its; /* maximum number of iterations */
260: X = snes->vec_sol; /* solution vector */
261: F = snes->vec_func; /* residual vector */
262: Y = snes->work[0]; /* work vectors */
264: PetscObjectSAWsTakeAccess((PetscObject)snes);
265: snes->iter = 0;
266: snes->norm = 0.0;
267: PetscObjectSAWsGrantAccess((PetscObject)snes);
269: SNESVIProjectOntoBounds(snes,X);
270: SNESComputeFunction(snes,X,vi->phi);
271: if (snes->domainerror) {
272: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
273: sdm->ops->computefunction = vi->computeuserfunction;
274: return(0);
275: }
276: /* Compute Merit function */
277: SNESVIComputeMeritFunction(vi->phi,&vi->merit,&vi->phinorm);
279: VecNormBegin(X,NORM_2,&xnorm); /* xnorm <- ||x|| */
280: VecNormEnd(X,NORM_2,&xnorm);
281: if (PetscIsInfOrNanReal(vi->merit)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
283: PetscObjectSAWsTakeAccess((PetscObject)snes);
284: snes->norm = vi->phinorm;
285: PetscObjectSAWsGrantAccess((PetscObject)snes);
286: SNESLogConvergenceHistory(snes,vi->phinorm,0);
287: SNESMonitor(snes,0,vi->phinorm);
289: /* test convergence */
290: (*snes->ops->converged)(snes,0,0.0,0.0,vi->phinorm,&snes->reason,snes->cnvP);
291: if (snes->reason) {
292: sdm->ops->computefunction = vi->computeuserfunction;
293: return(0);
294: }
296: for (i=0; i<maxits; i++) {
298: /* Call general purpose update function */
299: if (snes->ops->update) {
300: (*snes->ops->update)(snes, snes->iter);
301: }
303: /* Solve J Y = Phi, where J is the semismooth jacobian */
305: /* Get the jacobian -- note that the function must be the original function for snes_fd and snes_fd_color to work for this*/
306: sdm->ops->computefunction = vi->computeuserfunction;
307: SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
308: sdm->ops->computefunction = SNESVIComputeFunction;
310: /* Get the diagonal shift and row scaling vectors */
311: SNESVIComputeBsubdifferentialVectors(snes,X,F,snes->jacobian,vi->Da,vi->Db);
312: /* Compute the semismooth jacobian */
313: SNESVIComputeJacobian(snes->jacobian,snes->jacobian_pre,vi->Da,vi->Db);
314: /* Compute the merit function gradient */
315: SNESVIComputeMeritFunctionGradient(snes->jacobian,vi->phi,vi->dpsi);
316: KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);
317: KSPSolve(snes->ksp,vi->phi,Y);
318: KSPGetConvergedReason(snes->ksp,&kspreason);
320: if (kspreason < 0) {
321: if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
322: PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);
323: snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
324: break;
325: }
326: }
327: KSPGetIterationNumber(snes->ksp,&lits);
328: snes->linear_its += lits;
329: PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);
330: /*
331: if (snes->ops->precheck) {
332: PetscBool changed_y = PETSC_FALSE;
333: (*snes->ops->precheck)(snes,X,Y,snes->precheck,&changed_y);
334: }
336: if (PetscLogPrintInfo) {
337: SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);
338: }
339: */
340: /* Compute a (scaled) negative update in the line search routine:
341: Y <- X - lambda*Y
342: and evaluate G = function(Y) (depends on the line search).
343: */
344: VecCopy(Y,snes->vec_sol_update);
345: ynorm = 1; gnorm = vi->phinorm;
346: SNESLineSearchApply(snes->linesearch, X, vi->phi, &gnorm, Y);
347: SNESLineSearchGetSuccess(snes->linesearch, &lssucceed);
348: SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm);
349: PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",(double)vi->phinorm,(double)gnorm,(double)ynorm,(int)lssucceed);
350: if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
351: if (snes->domainerror) {
352: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
353: sdm->ops->computefunction = vi->computeuserfunction;
354: return(0);
355: }
356: if (!lssucceed) {
357: if (++snes->numFailures >= snes->maxFailures) {
358: PetscBool ismin;
359: snes->reason = SNES_DIVERGED_LINE_SEARCH;
360: SNESVICheckLocalMin_Private(snes,snes->jacobian,vi->phi,X,gnorm,&ismin);
361: if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
362: break;
363: }
364: }
365: /* Update function and solution vectors */
366: vi->phinorm = gnorm;
367: vi->merit = 0.5*vi->phinorm*vi->phinorm;
368: /* Monitor convergence */
369: PetscObjectSAWsTakeAccess((PetscObject)snes);
370: snes->iter = i+1;
371: snes->norm = vi->phinorm;
372: PetscObjectSAWsGrantAccess((PetscObject)snes);
373: SNESLogConvergenceHistory(snes,snes->norm,lits);
374: SNESMonitor(snes,snes->iter,snes->norm);
375: /* Test for convergence, xnorm = || X || */
376: if (snes->ops->converged != SNESConvergedSkip) { VecNorm(X,NORM_2,&xnorm); }
377: (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,vi->phinorm,&snes->reason,snes->cnvP);
378: if (snes->reason) break;
379: }
380: if (i == maxits) {
381: PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);
382: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
383: }
384: sdm->ops->computefunction = vi->computeuserfunction;
385: return(0);
386: }
388: /* -------------------------------------------------------------------------- */
389: /*
390: SNESSetUp_VINEWTONSSLS - Sets up the internal data structures for the later use
391: of the SNES nonlinear solver.
393: Input Parameter:
394: . snes - the SNES context
395: . x - the solution vector
397: Application Interface Routine: SNESSetUp()
399: Notes:
400: For basic use of the SNES solvers, the user need not explicitly call
401: SNESSetUp(), since these actions will automatically occur during
402: the call to SNESSolve().
403: */
406: PetscErrorCode SNESSetUp_VINEWTONSSLS(SNES snes)
407: {
408: PetscErrorCode ierr;
409: SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS*) snes->data;
412: SNESSetUp_VI(snes);
413: VecDuplicate(snes->vec_sol, &vi->dpsi);
414: VecDuplicate(snes->vec_sol, &vi->phi);
415: VecDuplicate(snes->vec_sol, &vi->Da);
416: VecDuplicate(snes->vec_sol, &vi->Db);
417: VecDuplicate(snes->vec_sol, &vi->z);
418: VecDuplicate(snes->vec_sol, &vi->t);
419: return(0);
420: }
421: /* -------------------------------------------------------------------------- */
424: PetscErrorCode SNESReset_VINEWTONSSLS(SNES snes)
425: {
426: SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS*) snes->data;
427: PetscErrorCode ierr;
430: SNESReset_VI(snes);
431: VecDestroy(&vi->dpsi);
432: VecDestroy(&vi->phi);
433: VecDestroy(&vi->Da);
434: VecDestroy(&vi->Db);
435: VecDestroy(&vi->z);
436: VecDestroy(&vi->t);
437: return(0);
438: }
440: /* -------------------------------------------------------------------------- */
441: /*
442: SNESSetFromOptions_VINEWTONSSLS - Sets various parameters for the SNESVI method.
444: Input Parameter:
445: . snes - the SNES context
447: Application Interface Routine: SNESSetFromOptions()
448: */
451: static PetscErrorCode SNESSetFromOptions_VINEWTONSSLS(SNES snes)
452: {
454: SNESLineSearch linesearch;
457: SNESSetFromOptions_VI(snes);
458: PetscOptionsHead("SNES semismooth method options");
459: PetscOptionsTail();
460: /* set up the default line search */
461: if (!snes->linesearch) {
462: SNESGetLineSearch(snes, &linesearch);
463: SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
464: SNESLineSearchBTSetAlpha(linesearch, 0.0);
465: }
466: return(0);
467: }
470: /* -------------------------------------------------------------------------- */
471: /*MC
472: SNESVINEWTONSSLS - Semi-smooth solver for variational inequalities based on Newton's method
474: Options Database:
475: + -snes_vi_type <ss,rs,rsaug> a semi-smooth solver, a reduced space active set method, and a reduced space active set method that does not eliminate the active constraints from the Jacobian instead augments the Jacobian with additional variables that enforce the constraints
476: - -snes_vi_monitor - prints the number of active constraints at each iteration.
478: Level: beginner
480: References:
481: - T. S. Munson, F. Facchinei, M. C. Ferris, A. Fischer, and C. Kanzow. The semismooth
482: algorithm for large scale complementarity problems. INFORMS Journal on Computing, 13 (2001).
484: .seealso: SNESVISetVariableBounds(), SNESVISetComputeVariableBounds(), SNESCreate(), SNES, SNESSetType(), SNESVINEWTONRSLS, SNESVINEWTONSSLS, SNESNEWTONTR, SNESLineSearchSet(),
485: SNESLineSearchSetPostCheck(), SNESLineSearchNo(), SNESLineSearchCubic(), SNESLineSearchQuadratic(),
486: SNESLineSearchSet(), SNESLineSearchNoNorms(), SNESLineSearchSetPreCheck(), SNESLineSearchSetParams(), SNESLineSearchGetParams()
488: M*/
491: PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONSSLS(SNES snes)
492: {
493: PetscErrorCode ierr;
494: SNES_VINEWTONSSLS *vi;
497: snes->ops->reset = SNESReset_VINEWTONSSLS;
498: snes->ops->setup = SNESSetUp_VINEWTONSSLS;
499: snes->ops->solve = SNESSolve_VINEWTONSSLS;
500: snes->ops->destroy = SNESDestroy_VI;
501: snes->ops->setfromoptions = SNESSetFromOptions_VINEWTONSSLS;
502: snes->ops->view = NULL;
504: snes->usesksp = PETSC_TRUE;
505: snes->usespc = PETSC_FALSE;
507: PetscNewLog(snes,&vi);
508: snes->data = (void*)vi;
510: PetscObjectComposeFunction((PetscObject)snes,"SNESVISetVariableBounds_C",SNESVISetVariableBounds_VI);
511: PetscObjectComposeFunction((PetscObject)snes,"SNESVISetComputeVariableBounds_C",SNESVISetComputeVariableBounds_VI);
512: return(0);
513: }