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