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