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: VecNormBegin(phi, NORM_2, phinorm);
21: VecNormEnd(phi, NORM_2, phinorm);
22: *merit = 0.5 * (*phinorm) * (*phinorm);
23: return 0;
24: }
26: static inline PetscScalar Phi(PetscScalar a, PetscScalar b)
27: {
28: return a + b - PetscSqrtScalar(a * a + b * b);
29: }
31: static inline PetscScalar DPhi(PetscScalar a, PetscScalar b)
32: {
33: if ((PetscAbsScalar(a) >= 1.e-6) || (PetscAbsScalar(b) >= 1.e-6)) return 1.0 - a / PetscSqrtScalar(a * a + b * b);
34: else return .5;
35: }
37: /*@
38: SNESVIComputeFunction - Provides the function that reformulates a system of nonlinear equations in mixed complementarity form to a system of nonlinear
39: equations in semismooth form.
41: Input Parameters:
42: + snes - the SNES context
43: . X - current iterate
44: - functx - user defined function context
46: Output Parameter:
47: . phi - the evaluation of Semismooth function at X
49: Level: developer
51: .seealso: `SNESVINEWTONSSLS`, `SNESVIComputeMeritFunction()`
52: @*/
53: PetscErrorCode SNESVIComputeFunction(SNES snes, Vec X, Vec phi, void *functx)
54: {
55: SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS *)snes->data;
56: Vec Xl = snes->xl, Xu = snes->xu, F = snes->vec_func;
57: PetscScalar *phi_arr, *f_arr, *l, *u;
58: const PetscScalar *x_arr;
59: PetscInt i, nlocal;
61: (*vi->computeuserfunction)(snes, X, F, functx);
62: VecGetLocalSize(X, &nlocal);
63: VecGetArrayRead(X, &x_arr);
64: VecGetArray(F, &f_arr);
65: VecGetArray(Xl, &l);
66: VecGetArray(Xu, &u);
67: VecGetArray(phi, &phi_arr);
69: for (i = 0; i < nlocal; i++) {
70: if ((PetscRealPart(l[i]) <= PETSC_NINFINITY) && (PetscRealPart(u[i]) >= PETSC_INFINITY)) { /* no constraints on variable */
71: phi_arr[i] = f_arr[i];
72: } else if (PetscRealPart(l[i]) <= PETSC_NINFINITY) { /* upper bound on variable only */
73: phi_arr[i] = -Phi(u[i] - x_arr[i], -f_arr[i]);
74: } else if (PetscRealPart(u[i]) >= PETSC_INFINITY) { /* lower bound on variable only */
75: phi_arr[i] = Phi(x_arr[i] - l[i], f_arr[i]);
76: } else if (l[i] == u[i]) {
77: phi_arr[i] = l[i] - x_arr[i];
78: } else { /* both bounds on variable */
79: phi_arr[i] = Phi(x_arr[i] - l[i], -Phi(u[i] - x_arr[i], -f_arr[i]));
80: }
81: }
83: VecRestoreArrayRead(X, &x_arr);
84: VecRestoreArray(F, &f_arr);
85: VecRestoreArray(Xl, &l);
86: VecRestoreArray(Xu, &u);
87: VecRestoreArray(phi, &phi_arr);
88: return 0;
89: }
91: /*
92: SNESVIComputeBsubdifferentialVectors - Computes the diagonal shift (Da) and row scaling (Db) vectors needed for the
93: the semismooth jacobian.
94: */
95: PetscErrorCode SNESVIComputeBsubdifferentialVectors(SNES snes, Vec X, Vec F, Mat jac, Vec Da, Vec Db)
96: {
97: PetscScalar *l, *u, *x, *f, *da, *db, da1, da2, db1, db2;
98: PetscInt i, nlocal;
100: VecGetArray(X, &x);
101: VecGetArray(F, &f);
102: VecGetArray(snes->xl, &l);
103: VecGetArray(snes->xu, &u);
104: VecGetArray(Da, &da);
105: VecGetArray(Db, &db);
106: VecGetLocalSize(X, &nlocal);
108: for (i = 0; i < nlocal; i++) {
109: if ((PetscRealPart(l[i]) <= PETSC_NINFINITY) && (PetscRealPart(u[i]) >= PETSC_INFINITY)) { /* no constraints on variable */
110: da[i] = 0;
111: db[i] = 1;
112: } else if (PetscRealPart(l[i]) <= PETSC_NINFINITY) { /* upper bound on variable only */
113: da[i] = DPhi(u[i] - x[i], -f[i]);
114: db[i] = DPhi(-f[i], u[i] - x[i]);
115: } else if (PetscRealPart(u[i]) >= PETSC_INFINITY) { /* lower bound on variable only */
116: da[i] = DPhi(x[i] - l[i], f[i]);
117: db[i] = DPhi(f[i], x[i] - l[i]);
118: } else if (l[i] == u[i]) { /* fixed variable */
119: da[i] = 1;
120: db[i] = 0;
121: } else { /* upper and lower bounds on variable */
122: da1 = DPhi(x[i] - l[i], -Phi(u[i] - x[i], -f[i]));
123: db1 = DPhi(-Phi(u[i] - x[i], -f[i]), x[i] - l[i]);
124: da2 = DPhi(u[i] - x[i], -f[i]);
125: db2 = DPhi(-f[i], u[i] - x[i]);
126: da[i] = da1 + db1 * da2;
127: db[i] = db1 * db2;
128: }
129: }
131: VecRestoreArray(X, &x);
132: VecRestoreArray(F, &f);
133: VecRestoreArray(snes->xl, &l);
134: VecRestoreArray(snes->xu, &u);
135: VecRestoreArray(Da, &da);
136: VecRestoreArray(Db, &db);
137: return 0;
138: }
140: /*
141: 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.
143: Input Parameters:
144: . Da - Diagonal shift vector for the semismooth jacobian.
145: . Db - Row scaling vector for the semismooth jacobian.
147: Output Parameters:
148: . jac - semismooth jacobian
149: . jac_pre - optional preconditioning matrix
151: Note:
152: The semismooth jacobian matrix is given by
153: jac = Da + Db*jacfun
154: where Db is the row scaling matrix stored as a vector,
155: Da is the diagonal perturbation matrix stored as a vector
156: and jacfun is the jacobian of the original nonlinear function.
157: */
158: PetscErrorCode SNESVIComputeJacobian(Mat jac, Mat jac_pre, Vec Da, Vec Db)
159: {
160: /* Do row scaling and add diagonal perturbation */
161: MatDiagonalScale(jac, Db, NULL);
162: MatDiagonalSet(jac, Da, ADD_VALUES);
163: if (jac != jac_pre) { /* If jac and jac_pre are different */
164: MatDiagonalScale(jac_pre, Db, NULL);
165: MatDiagonalSet(jac_pre, Da, ADD_VALUES);
166: }
167: return 0;
168: }
170: /*
171: SNESVIComputeMeritFunctionGradient - Computes the gradient of the merit function psi.
173: Input Parameters:
174: phi - semismooth function.
175: H - semismooth jacobian
177: Output Parameter:
178: dpsi - merit function gradient
180: Note:
181: The merit function gradient is computed as follows
182: dpsi = H^T*phi
183: */
184: PetscErrorCode SNESVIComputeMeritFunctionGradient(Mat H, Vec phi, Vec dpsi)
185: {
186: MatMultTranspose(H, phi, dpsi);
187: return 0;
188: }
190: /*
191: SNESSolve_VINEWTONSSLS - Solves the complementarity problem with a semismooth Newton
192: method using a line search.
194: Input Parameter:
195: . snes - the SNES context
197: Application Interface Routine: SNESSolve()
199: Note:
200: This implements essentially a semismooth Newton method with a
201: line search. The default line search does not do any line search
202: but rather takes a full Newton step.
204: 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
206: */
207: PetscErrorCode SNESSolve_VINEWTONSSLS(SNES snes)
208: {
209: SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS *)snes->data;
210: PetscInt maxits, i, lits;
211: SNESLineSearchReason lssucceed;
212: PetscReal gnorm, xnorm = 0, ynorm;
213: Vec Y, X, F;
214: KSPConvergedReason kspreason;
215: DM dm;
216: DMSNES sdm;
218: SNESGetDM(snes, &dm);
219: DMGetDMSNES(dm, &sdm);
221: vi->computeuserfunction = sdm->ops->computefunction;
222: sdm->ops->computefunction = SNESVIComputeFunction;
224: snes->numFailures = 0;
225: snes->numLinearSolveFailures = 0;
226: snes->reason = SNES_CONVERGED_ITERATING;
228: maxits = snes->max_its; /* maximum number of iterations */
229: X = snes->vec_sol; /* solution vector */
230: F = snes->vec_func; /* residual vector */
231: Y = snes->work[0]; /* work vectors */
233: PetscObjectSAWsTakeAccess((PetscObject)snes);
234: snes->iter = 0;
235: snes->norm = 0.0;
236: PetscObjectSAWsGrantAccess((PetscObject)snes);
238: SNESVIProjectOntoBounds(snes, X);
239: SNESComputeFunction(snes, X, vi->phi);
240: if (snes->domainerror) {
241: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
242: sdm->ops->computefunction = vi->computeuserfunction;
243: return 0;
244: }
245: /* Compute Merit function */
246: SNESVIComputeMeritFunction(vi->phi, &vi->merit, &vi->phinorm);
248: VecNormBegin(X, NORM_2, &xnorm); /* xnorm <- ||x|| */
249: VecNormEnd(X, NORM_2, &xnorm);
250: SNESCheckFunctionNorm(snes, vi->merit);
252: PetscObjectSAWsTakeAccess((PetscObject)snes);
253: snes->norm = vi->phinorm;
254: PetscObjectSAWsGrantAccess((PetscObject)snes);
255: SNESLogConvergenceHistory(snes, vi->phinorm, 0);
256: SNESMonitor(snes, 0, vi->phinorm);
258: /* test convergence */
259: PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, vi->phinorm, &snes->reason, snes->cnvP);
260: if (snes->reason) {
261: sdm->ops->computefunction = vi->computeuserfunction;
262: return 0;
263: }
265: for (i = 0; i < maxits; i++) {
266: /* Call general purpose update function */
267: PetscTryTypeMethod(snes, update, snes->iter);
269: /* Solve J Y = Phi, where J is the semismooth jacobian */
271: /* Get the jacobian -- note that the function must be the original function for snes_fd and snes_fd_color to work for this*/
272: sdm->ops->computefunction = vi->computeuserfunction;
273: SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre);
274: SNESCheckJacobianDomainerror(snes);
275: sdm->ops->computefunction = SNESVIComputeFunction;
277: /* Get the diagonal shift and row scaling vectors */
278: SNESVIComputeBsubdifferentialVectors(snes, X, F, snes->jacobian, vi->Da, vi->Db);
279: /* Compute the semismooth jacobian */
280: SNESVIComputeJacobian(snes->jacobian, snes->jacobian_pre, vi->Da, vi->Db);
281: /* Compute the merit function gradient */
282: SNESVIComputeMeritFunctionGradient(snes->jacobian, vi->phi, vi->dpsi);
283: KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre);
284: KSPSolve(snes->ksp, vi->phi, Y);
285: KSPGetConvergedReason(snes->ksp, &kspreason);
287: if (kspreason < 0) {
288: if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
289: PetscInfo(snes, "iter=%" PetscInt_FMT ", number linear solve failures %" PetscInt_FMT " greater than current SNES allowed, stopping solve\n", snes->iter, snes->numLinearSolveFailures);
290: snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
291: break;
292: }
293: }
294: KSPGetIterationNumber(snes->ksp, &lits);
295: snes->linear_its += lits;
296: PetscInfo(snes, "iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits);
297: /*
298: if (snes->ops->precheck) {
299: PetscBool changed_y = PETSC_FALSE;
300: PetscUseTypeMethod(snes,precheck ,X,Y,snes->precheck,&changed_y);
301: }
303: if (PetscLogPrintInfo) SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);
304: */
305: /* Compute a (scaled) negative update in the line search routine:
306: Y <- X - lambda*Y
307: and evaluate G = function(Y) (depends on the line search).
308: */
309: VecCopy(Y, snes->vec_sol_update);
310: ynorm = 1;
311: gnorm = vi->phinorm;
312: SNESLineSearchApply(snes->linesearch, X, vi->phi, &gnorm, Y);
313: SNESLineSearchGetReason(snes->linesearch, &lssucceed);
314: SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm);
315: PetscInfo(snes, "fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n", (double)vi->phinorm, (double)gnorm, (double)ynorm, (int)lssucceed);
316: if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
317: if (snes->domainerror) {
318: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
319: sdm->ops->computefunction = vi->computeuserfunction;
320: return 0;
321: }
322: if (lssucceed) {
323: if (++snes->numFailures >= snes->maxFailures) {
324: PetscBool ismin;
325: snes->reason = SNES_DIVERGED_LINE_SEARCH;
326: SNESVICheckLocalMin_Private(snes, snes->jacobian, vi->phi, X, gnorm, &ismin);
327: if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
328: break;
329: }
330: }
331: /* Update function and solution vectors */
332: vi->phinorm = gnorm;
333: vi->merit = 0.5 * vi->phinorm * vi->phinorm;
334: /* Monitor convergence */
335: PetscObjectSAWsTakeAccess((PetscObject)snes);
336: snes->iter = i + 1;
337: snes->norm = vi->phinorm;
338: snes->xnorm = xnorm;
339: snes->ynorm = ynorm;
340: PetscObjectSAWsGrantAccess((PetscObject)snes);
341: SNESLogConvergenceHistory(snes, snes->norm, lits);
342: SNESMonitor(snes, snes->iter, snes->norm);
343: /* Test for convergence, xnorm = || X || */
344: if (snes->ops->converged != SNESConvergedSkip) VecNorm(X, NORM_2, &xnorm);
345: PetscUseTypeMethod(snes, converged, snes->iter, xnorm, ynorm, vi->phinorm, &snes->reason, snes->cnvP);
346: if (snes->reason) break;
347: }
348: if (i == maxits) {
349: PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", maxits);
350: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
351: }
352: sdm->ops->computefunction = vi->computeuserfunction;
353: return 0;
354: }
356: /*
357: SNESSetUp_VINEWTONSSLS - Sets up the internal data structures for the later use
358: of the SNES nonlinear solver.
360: Input Parameter:
361: . snes - the SNES context
363: Application Interface Routine: SNESSetUp()
365: Note:
366: For basic use of the SNES solvers, the user need not explicitly call
367: SNESSetUp(), since these actions will automatically occur during
368: the call to SNESSolve().
369: */
370: PetscErrorCode SNESSetUp_VINEWTONSSLS(SNES snes)
371: {
372: SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS *)snes->data;
374: SNESSetUp_VI(snes);
375: VecDuplicate(snes->vec_sol, &vi->dpsi);
376: VecDuplicate(snes->vec_sol, &vi->phi);
377: VecDuplicate(snes->vec_sol, &vi->Da);
378: VecDuplicate(snes->vec_sol, &vi->Db);
379: VecDuplicate(snes->vec_sol, &vi->z);
380: VecDuplicate(snes->vec_sol, &vi->t);
381: return 0;
382: }
384: PetscErrorCode SNESReset_VINEWTONSSLS(SNES snes)
385: {
386: SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS *)snes->data;
388: SNESReset_VI(snes);
389: VecDestroy(&vi->dpsi);
390: VecDestroy(&vi->phi);
391: VecDestroy(&vi->Da);
392: VecDestroy(&vi->Db);
393: VecDestroy(&vi->z);
394: VecDestroy(&vi->t);
395: return 0;
396: }
398: /*
399: SNESSetFromOptions_VINEWTONSSLS - Sets various parameters for the SNESVI method.
401: Input Parameter:
402: . snes - the SNES context
404: Application Interface Routine: SNESSetFromOptions()
405: */
406: static PetscErrorCode SNESSetFromOptions_VINEWTONSSLS(SNES snes, PetscOptionItems *PetscOptionsObject)
407: {
408: SNESSetFromOptions_VI(snes, PetscOptionsObject);
409: PetscOptionsHeadBegin(PetscOptionsObject, "SNES semismooth method options");
410: PetscOptionsHeadEnd();
411: return 0;
412: }
414: /*MC
415: SNESVINEWTONSSLS - Semi-smooth solver for variational inequalities based on Newton's method
417: Options Database Keys:
418: + -snes_type <vinewtonssls,vinewtonrsls> a semi-smooth solver, a reduced space active set method
419: - -snes_vi_monitor - prints the number of active constraints at each iteration.
421: Level: beginner
423: References:
424: + * - T. S. Munson, F. Facchinei, M. C. Ferris, A. Fischer, and C. Kanzow. The semismooth
425: algorithm for large scale complementarity problems. INFORMS Journal on Computing, 13 (2001).
426: - * - T. S. Munson, and S. Benson. Flexible Complementarity Solvers for Large Scale
427: Applications, Optimization Methods and Software, 21 (2006).
429: Notes:
430: This family of algorithm is much like an interior point method.
432: The reduced space active set solvers `SNESVINEWTONRSLS` provide an alternative approach that does not result in extremely ill-conditioned linear systems
434: .seealso: `SNESVINEWTONRSLS`, `SNESVISetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESVINEWTONRSLS`, `SNESNEWTONTR`, `SNESLineSearchSetType()`, `SNESLineSearchSetPostCheck()`, `SNESLineSearchSetPreCheck()`
435: M*/
436: PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONSSLS(SNES snes)
437: {
438: SNES_VINEWTONSSLS *vi;
439: SNESLineSearch linesearch;
441: snes->ops->reset = SNESReset_VINEWTONSSLS;
442: snes->ops->setup = SNESSetUp_VINEWTONSSLS;
443: snes->ops->solve = SNESSolve_VINEWTONSSLS;
444: snes->ops->destroy = SNESDestroy_VI;
445: snes->ops->setfromoptions = SNESSetFromOptions_VINEWTONSSLS;
446: snes->ops->view = NULL;
448: snes->usesksp = PETSC_TRUE;
449: snes->usesnpc = PETSC_FALSE;
451: SNESGetLineSearch(snes, &linesearch);
452: if (!((PetscObject)linesearch)->type_name) {
453: SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
454: SNESLineSearchBTSetAlpha(linesearch, 0.0);
455: }
457: snes->alwayscomputesfinalresidual = PETSC_FALSE;
459: PetscNew(&vi);
460: snes->data = (void *)vi;
462: PetscObjectComposeFunction((PetscObject)snes, "SNESVISetVariableBounds_C", SNESVISetVariableBounds_VI);
463: PetscObjectComposeFunction((PetscObject)snes, "SNESVISetComputeVariableBounds_C", SNESVISetComputeVariableBounds_VI);
464: return 0;
465: }