Actual source code: viss.c
petsc-3.14.6 2021-03-30
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: {
23: VecNormBegin(phi,NORM_2,phinorm);
24: VecNormEnd(phi,NORM_2,phinorm);
26: *merit = 0.5*(*phinorm)*(*phinorm);
27: return(0);
28: }
30: PETSC_STATIC_INLINE PetscScalar Phi(PetscScalar a,PetscScalar b)
31: {
32: return a + b - PetscSqrtScalar(a*a + b*b);
33: }
35: PETSC_STATIC_INLINE PetscScalar DPhi(PetscScalar a,PetscScalar b)
36: {
37: if ((PetscAbsScalar(a) >= 1.e-6) || (PetscAbsScalar(b) >= 1.e-6)) return 1.0 - a/ PetscSqrtScalar(a*a + b*b);
38: else return .5;
39: }
41: /*
42: SNESVIComputeFunction - Reformulates a system of nonlinear equations in mixed complementarity form to a system of nonlinear equations in semismooth form.
44: Input Parameters:
45: . snes - the SNES context
46: . X - current iterate
47: . functx - user defined function context
49: Output Parameters:
50: . phi - Semismooth function
52: */
53: static PetscErrorCode SNESVIComputeFunction(SNES snes,Vec X,Vec phi,void *functx)
54: {
55: PetscErrorCode ierr;
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;
63: (*vi->computeuserfunction)(snes,X,F,functx);
64: VecGetLocalSize(X,&nlocal);
65: VecGetArrayRead(X,&x_arr);
66: VecGetArray(F,&f_arr);
67: VecGetArray(Xl,&l);
68: VecGetArray(Xu,&u);
69: 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: VecRestoreArrayRead(X,&x_arr);
86: VecRestoreArray(F,&f_arr);
87: VecRestoreArray(Xl,&l);
88: VecRestoreArray(Xu,&u);
89: VecRestoreArray(phi,&phi_arr);
90: return(0);
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: {
100: PetscScalar *l,*u,*x,*f,*da,*db,da1,da2,db1,db2;
101: PetscInt i,nlocal;
104: VecGetArray(X,&x);
105: VecGetArray(F,&f);
106: VecGetArray(snes->xl,&l);
107: VecGetArray(snes->xu,&u);
108: VecGetArray(Da,&da);
109: VecGetArray(Db,&db);
110: VecGetLocalSize(X,&nlocal);
112: for (i=0; i< nlocal; i++) {
113: if ((PetscRealPart(l[i]) <= PETSC_NINFINITY) && (PetscRealPart(u[i]) >= PETSC_INFINITY)) { /* no constraints on variable */
114: da[i] = 0;
115: db[i] = 1;
116: } else if (PetscRealPart(l[i]) <= PETSC_NINFINITY) { /* upper bound on variable only */
117: da[i] = DPhi(u[i] - x[i], -f[i]);
118: db[i] = DPhi(-f[i],u[i] - x[i]);
119: } else if (PetscRealPart(u[i]) >= PETSC_INFINITY) { /* lower bound on variable only */
120: da[i] = DPhi(x[i] - l[i], f[i]);
121: db[i] = DPhi(f[i],x[i] - l[i]);
122: } else if (l[i] == u[i]) { /* fixed variable */
123: da[i] = 1;
124: db[i] = 0;
125: } else { /* upper and lower bounds on variable */
126: da1 = DPhi(x[i] - l[i], -Phi(u[i] - x[i], -f[i]));
127: db1 = DPhi(-Phi(u[i] - x[i], -f[i]),x[i] - l[i]);
128: da2 = DPhi(u[i] - x[i], -f[i]);
129: db2 = DPhi(-f[i],u[i] - x[i]);
130: da[i] = da1 + db1*da2;
131: db[i] = db1*db2;
132: }
133: }
135: VecRestoreArray(X,&x);
136: VecRestoreArray(F,&f);
137: VecRestoreArray(snes->xl,&l);
138: VecRestoreArray(snes->xu,&u);
139: VecRestoreArray(Da,&da);
140: VecRestoreArray(Db,&db);
141: return(0);
142: }
144: /*
145: 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.
147: Input Parameters:
148: . Da - Diagonal shift vector for the semismooth jacobian.
149: . Db - Row scaling vector for the semismooth jacobian.
151: Output Parameters:
152: . jac - semismooth jacobian
153: . jac_pre - optional preconditioning matrix
155: Notes:
156: The semismooth jacobian matrix is given by
157: jac = Da + Db*jacfun
158: where Db is the row scaling matrix stored as a vector,
159: Da is the diagonal perturbation matrix stored as a vector
160: and jacfun is the jacobian of the original nonlinear function.
161: */
162: PetscErrorCode SNESVIComputeJacobian(Mat jac, Mat jac_pre,Vec Da, Vec Db)
163: {
166: /* Do row scaling and add diagonal perturbation */
167: MatDiagonalScale(jac,Db,NULL);
168: MatDiagonalSet(jac,Da,ADD_VALUES);
169: if (jac != jac_pre) { /* If jac and jac_pre are different */
170: MatDiagonalScale(jac_pre,Db,NULL);
171: MatDiagonalSet(jac_pre,Da,ADD_VALUES);
172: }
173: return(0);
174: }
176: /*
177: SNESVIComputeMeritFunctionGradient - Computes the gradient of the merit function psi.
179: Input Parameters:
180: phi - semismooth function.
181: H - semismooth jacobian
183: Output Parameters:
184: dpsi - merit function gradient
186: Notes:
187: The merit function gradient is computed as follows
188: dpsi = H^T*phi
189: */
190: PetscErrorCode SNESVIComputeMeritFunctionGradient(Mat H, Vec phi, Vec dpsi)
191: {
195: MatMultTranspose(H,phi,dpsi);
196: return(0);
197: }
201: /*
202: SNESSolve_VINEWTONSSLS - Solves the complementarity problem with a semismooth Newton
203: method using a line search.
205: Input Parameters:
206: . snes - the SNES context
208: Application Interface Routine: SNESSolve()
210: Notes:
211: This implements essentially a semismooth Newton method with a
212: line search. The default line search does not do any line search
213: but rather takes a full Newton step.
215: 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
217: */
218: PetscErrorCode SNESSolve_VINEWTONSSLS(SNES snes)
219: {
220: SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS*)snes->data;
221: PetscErrorCode ierr;
222: PetscInt maxits,i,lits;
223: SNESLineSearchReason lssucceed;
224: PetscReal gnorm,xnorm=0,ynorm;
225: Vec Y,X,F;
226: KSPConvergedReason kspreason;
227: DM dm;
228: DMSNES sdm;
231: SNESGetDM(snes,&dm);
232: DMGetDMSNES(dm,&sdm);
234: vi->computeuserfunction = sdm->ops->computefunction;
235: sdm->ops->computefunction = SNESVIComputeFunction;
237: snes->numFailures = 0;
238: snes->numLinearSolveFailures = 0;
239: snes->reason = SNES_CONVERGED_ITERATING;
241: maxits = snes->max_its; /* maximum number of iterations */
242: X = snes->vec_sol; /* solution vector */
243: F = snes->vec_func; /* residual vector */
244: Y = snes->work[0]; /* work vectors */
246: PetscObjectSAWsTakeAccess((PetscObject)snes);
247: snes->iter = 0;
248: snes->norm = 0.0;
249: PetscObjectSAWsGrantAccess((PetscObject)snes);
251: SNESVIProjectOntoBounds(snes,X);
252: SNESComputeFunction(snes,X,vi->phi);
253: if (snes->domainerror) {
254: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
255: sdm->ops->computefunction = vi->computeuserfunction;
256: return(0);
257: }
258: /* Compute Merit function */
259: SNESVIComputeMeritFunction(vi->phi,&vi->merit,&vi->phinorm);
261: VecNormBegin(X,NORM_2,&xnorm); /* xnorm <- ||x|| */
262: VecNormEnd(X,NORM_2,&xnorm);
263: SNESCheckFunctionNorm(snes,vi->merit);
265: PetscObjectSAWsTakeAccess((PetscObject)snes);
266: snes->norm = vi->phinorm;
267: PetscObjectSAWsGrantAccess((PetscObject)snes);
268: SNESLogConvergenceHistory(snes,vi->phinorm,0);
269: SNESMonitor(snes,0,vi->phinorm);
271: /* test convergence */
272: (*snes->ops->converged)(snes,0,0.0,0.0,vi->phinorm,&snes->reason,snes->cnvP);
273: if (snes->reason) {
274: sdm->ops->computefunction = vi->computeuserfunction;
275: return(0);
276: }
278: for (i=0; i<maxits; i++) {
280: /* Call general purpose update function */
281: if (snes->ops->update) {
282: (*snes->ops->update)(snes, snes->iter);
283: }
285: /* Solve J Y = Phi, where J is the semismooth jacobian */
287: /* Get the jacobian -- note that the function must be the original function for snes_fd and snes_fd_color to work for this*/
288: sdm->ops->computefunction = vi->computeuserfunction;
289: SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
290: SNESCheckJacobianDomainerror(snes);
291: sdm->ops->computefunction = SNESVIComputeFunction;
293: /* Get the diagonal shift and row scaling vectors */
294: SNESVIComputeBsubdifferentialVectors(snes,X,F,snes->jacobian,vi->Da,vi->Db);
295: /* Compute the semismooth jacobian */
296: SNESVIComputeJacobian(snes->jacobian,snes->jacobian_pre,vi->Da,vi->Db);
297: /* Compute the merit function gradient */
298: SNESVIComputeMeritFunctionGradient(snes->jacobian,vi->phi,vi->dpsi);
299: KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);
300: KSPSolve(snes->ksp,vi->phi,Y);
301: KSPGetConvergedReason(snes->ksp,&kspreason);
303: if (kspreason < 0) {
304: if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
305: PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);
306: snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
307: break;
308: }
309: }
310: KSPGetIterationNumber(snes->ksp,&lits);
311: snes->linear_its += lits;
312: PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);
313: /*
314: if (snes->ops->precheck) {
315: PetscBool changed_y = PETSC_FALSE;
316: (*snes->ops->precheck)(snes,X,Y,snes->precheck,&changed_y);
317: }
319: if (PetscLogPrintInfo) {
320: SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);
321: }
322: */
323: /* Compute a (scaled) negative update in the line search routine:
324: Y <- X - lambda*Y
325: and evaluate G = function(Y) (depends on the line search).
326: */
327: VecCopy(Y,snes->vec_sol_update);
328: ynorm = 1; gnorm = vi->phinorm;
329: SNESLineSearchApply(snes->linesearch, X, vi->phi, &gnorm, Y);
330: SNESLineSearchGetReason(snes->linesearch, &lssucceed);
331: SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm);
332: PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",(double)vi->phinorm,(double)gnorm,(double)ynorm,(int)lssucceed);
333: if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
334: if (snes->domainerror) {
335: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
336: sdm->ops->computefunction = vi->computeuserfunction;
337: return(0);
338: }
339: if (lssucceed) {
340: if (++snes->numFailures >= snes->maxFailures) {
341: PetscBool ismin;
342: snes->reason = SNES_DIVERGED_LINE_SEARCH;
343: SNESVICheckLocalMin_Private(snes,snes->jacobian,vi->phi,X,gnorm,&ismin);
344: if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
345: break;
346: }
347: }
348: /* Update function and solution vectors */
349: vi->phinorm = gnorm;
350: vi->merit = 0.5*vi->phinorm*vi->phinorm;
351: /* Monitor convergence */
352: PetscObjectSAWsTakeAccess((PetscObject)snes);
353: snes->iter = i+1;
354: snes->norm = vi->phinorm;
355: snes->xnorm = xnorm;
356: snes->ynorm = ynorm;
357: PetscObjectSAWsGrantAccess((PetscObject)snes);
358: SNESLogConvergenceHistory(snes,snes->norm,lits);
359: SNESMonitor(snes,snes->iter,snes->norm);
360: /* Test for convergence, xnorm = || X || */
361: if (snes->ops->converged != SNESConvergedSkip) { VecNorm(X,NORM_2,&xnorm); }
362: (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,vi->phinorm,&snes->reason,snes->cnvP);
363: if (snes->reason) break;
364: }
365: if (i == maxits) {
366: PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);
367: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
368: }
369: sdm->ops->computefunction = vi->computeuserfunction;
370: return(0);
371: }
373: /* -------------------------------------------------------------------------- */
374: /*
375: SNESSetUp_VINEWTONSSLS - Sets up the internal data structures for the later use
376: of the SNES nonlinear solver.
378: Input Parameter:
379: . snes - the SNES context
381: Application Interface Routine: SNESSetUp()
383: Notes:
384: For basic use of the SNES solvers, the user need not explicitly call
385: SNESSetUp(), since these actions will automatically occur during
386: the call to SNESSolve().
387: */
388: PetscErrorCode SNESSetUp_VINEWTONSSLS(SNES snes)
389: {
390: PetscErrorCode ierr;
391: SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS*) snes->data;
394: SNESSetUp_VI(snes);
395: VecDuplicate(snes->vec_sol, &vi->dpsi);
396: VecDuplicate(snes->vec_sol, &vi->phi);
397: VecDuplicate(snes->vec_sol, &vi->Da);
398: VecDuplicate(snes->vec_sol, &vi->Db);
399: VecDuplicate(snes->vec_sol, &vi->z);
400: VecDuplicate(snes->vec_sol, &vi->t);
401: return(0);
402: }
403: /* -------------------------------------------------------------------------- */
404: PetscErrorCode SNESReset_VINEWTONSSLS(SNES snes)
405: {
406: SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS*) snes->data;
407: PetscErrorCode ierr;
410: SNESReset_VI(snes);
411: VecDestroy(&vi->dpsi);
412: VecDestroy(&vi->phi);
413: VecDestroy(&vi->Da);
414: VecDestroy(&vi->Db);
415: VecDestroy(&vi->z);
416: VecDestroy(&vi->t);
417: return(0);
418: }
420: /* -------------------------------------------------------------------------- */
421: /*
422: SNESSetFromOptions_VINEWTONSSLS - Sets various parameters for the SNESVI method.
424: Input Parameter:
425: . snes - the SNES context
427: Application Interface Routine: SNESSetFromOptions()
428: */
429: static PetscErrorCode SNESSetFromOptions_VINEWTONSSLS(PetscOptionItems *PetscOptionsObject,SNES snes)
430: {
434: SNESSetFromOptions_VI(PetscOptionsObject,snes);
435: PetscOptionsHead(PetscOptionsObject,"SNES semismooth method options");
436: PetscOptionsTail();
437: return(0);
438: }
441: /* -------------------------------------------------------------------------- */
442: /*MC
443: SNESVINEWTONSSLS - Semi-smooth solver for variational inequalities based on Newton's method
445: Options Database:
446: + -snes_type <vinewtonssls,vinewtonrsls> a semi-smooth solver, a reduced space active set method
447: - -snes_vi_monitor - prints the number of active constraints at each iteration.
449: Level: beginner
451: References:
452: + 1. - T. S. Munson, F. Facchinei, M. C. Ferris, A. Fischer, and C. Kanzow. The semismooth
453: algorithm for large scale complementarity problems. INFORMS Journal on Computing, 13 (2001).
454: - 2. - T. S. Munson, and S. Benson. Flexible Complementarity Solvers for Large Scale
455: Applications, Optimization Methods and Software, 21 (2006).
457: .seealso: SNESVISetVariableBounds(), SNESVISetComputeVariableBounds(), SNESCreate(), SNES, SNESSetType(), SNESVINEWTONRSLS, SNESNEWTONTR, SNESLineSearchSetType(),SNESLineSearchSetPostCheck(), SNESLineSearchSetPreCheck()
459: M*/
460: PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONSSLS(SNES snes)
461: {
462: PetscErrorCode ierr;
463: SNES_VINEWTONSSLS *vi;
464: SNESLineSearch linesearch;
467: snes->ops->reset = SNESReset_VINEWTONSSLS;
468: snes->ops->setup = SNESSetUp_VINEWTONSSLS;
469: snes->ops->solve = SNESSolve_VINEWTONSSLS;
470: snes->ops->destroy = SNESDestroy_VI;
471: snes->ops->setfromoptions = SNESSetFromOptions_VINEWTONSSLS;
472: snes->ops->view = NULL;
474: snes->usesksp = PETSC_TRUE;
475: snes->usesnpc = PETSC_FALSE;
477: SNESGetLineSearch(snes, &linesearch);
478: if (!((PetscObject)linesearch)->type_name) {
479: SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
480: SNESLineSearchBTSetAlpha(linesearch, 0.0);
481: }
483: snes->alwayscomputesfinalresidual = PETSC_FALSE;
485: PetscNewLog(snes,&vi);
486: snes->data = (void*)vi;
488: PetscObjectComposeFunction((PetscObject)snes,"SNESVISetVariableBounds_C",SNESVISetVariableBounds_VI);
489: PetscObjectComposeFunction((PetscObject)snes,"SNESVISetComputeVariableBounds_C",SNESVISetComputeVariableBounds_VI);
490: return(0);
491: }