Actual source code: viss.c
petsc-3.8.4 2018-03-24
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: sdm->ops->computefunction = SNESVIComputeFunction;
292: /* Get the diagonal shift and row scaling vectors */
293: SNESVIComputeBsubdifferentialVectors(snes,X,F,snes->jacobian,vi->Da,vi->Db);
294: /* Compute the semismooth jacobian */
295: SNESVIComputeJacobian(snes->jacobian,snes->jacobian_pre,vi->Da,vi->Db);
296: /* Compute the merit function gradient */
297: SNESVIComputeMeritFunctionGradient(snes->jacobian,vi->phi,vi->dpsi);
298: KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);
299: KSPSolve(snes->ksp,vi->phi,Y);
300: KSPGetConvergedReason(snes->ksp,&kspreason);
302: if (kspreason < 0) {
303: if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
304: PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);
305: snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
306: break;
307: }
308: }
309: KSPGetIterationNumber(snes->ksp,&lits);
310: snes->linear_its += lits;
311: PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);
312: /*
313: if (snes->ops->precheck) {
314: PetscBool changed_y = PETSC_FALSE;
315: (*snes->ops->precheck)(snes,X,Y,snes->precheck,&changed_y);
316: }
318: if (PetscLogPrintInfo) {
319: SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);
320: }
321: */
322: /* Compute a (scaled) negative update in the line search routine:
323: Y <- X - lambda*Y
324: and evaluate G = function(Y) (depends on the line search).
325: */
326: VecCopy(Y,snes->vec_sol_update);
327: ynorm = 1; gnorm = vi->phinorm;
328: SNESLineSearchApply(snes->linesearch, X, vi->phi, &gnorm, Y);
329: SNESLineSearchGetReason(snes->linesearch, &lssucceed);
330: SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm);
331: PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",(double)vi->phinorm,(double)gnorm,(double)ynorm,(int)lssucceed);
332: if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
333: if (snes->domainerror) {
334: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
335: sdm->ops->computefunction = vi->computeuserfunction;
336: return(0);
337: }
338: if (lssucceed) {
339: if (++snes->numFailures >= snes->maxFailures) {
340: PetscBool ismin;
341: snes->reason = SNES_DIVERGED_LINE_SEARCH;
342: SNESVICheckLocalMin_Private(snes,snes->jacobian,vi->phi,X,gnorm,&ismin);
343: if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
344: break;
345: }
346: }
347: /* Update function and solution vectors */
348: vi->phinorm = gnorm;
349: vi->merit = 0.5*vi->phinorm*vi->phinorm;
350: /* Monitor convergence */
351: PetscObjectSAWsTakeAccess((PetscObject)snes);
352: snes->iter = i+1;
353: snes->norm = vi->phinorm;
354: PetscObjectSAWsGrantAccess((PetscObject)snes);
355: SNESLogConvergenceHistory(snes,snes->norm,lits);
356: SNESMonitor(snes,snes->iter,snes->norm);
357: /* Test for convergence, xnorm = || X || */
358: if (snes->ops->converged != SNESConvergedSkip) { VecNorm(X,NORM_2,&xnorm); }
359: (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,vi->phinorm,&snes->reason,snes->cnvP);
360: if (snes->reason) break;
361: }
362: if (i == maxits) {
363: PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);
364: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
365: }
366: sdm->ops->computefunction = vi->computeuserfunction;
367: return(0);
368: }
370: /* -------------------------------------------------------------------------- */
371: /*
372: SNESSetUp_VINEWTONSSLS - Sets up the internal data structures for the later use
373: of the SNES nonlinear solver.
375: Input Parameter:
376: . snes - the SNES context
378: Application Interface Routine: SNESSetUp()
380: Notes:
381: For basic use of the SNES solvers, the user need not explicitly call
382: SNESSetUp(), since these actions will automatically occur during
383: the call to SNESSolve().
384: */
385: PetscErrorCode SNESSetUp_VINEWTONSSLS(SNES snes)
386: {
387: PetscErrorCode ierr;
388: SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS*) snes->data;
391: SNESSetUp_VI(snes);
392: VecDuplicate(snes->vec_sol, &vi->dpsi);
393: VecDuplicate(snes->vec_sol, &vi->phi);
394: VecDuplicate(snes->vec_sol, &vi->Da);
395: VecDuplicate(snes->vec_sol, &vi->Db);
396: VecDuplicate(snes->vec_sol, &vi->z);
397: VecDuplicate(snes->vec_sol, &vi->t);
398: return(0);
399: }
400: /* -------------------------------------------------------------------------- */
401: PetscErrorCode SNESReset_VINEWTONSSLS(SNES snes)
402: {
403: SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS*) snes->data;
404: PetscErrorCode ierr;
407: SNESReset_VI(snes);
408: VecDestroy(&vi->dpsi);
409: VecDestroy(&vi->phi);
410: VecDestroy(&vi->Da);
411: VecDestroy(&vi->Db);
412: VecDestroy(&vi->z);
413: VecDestroy(&vi->t);
414: return(0);
415: }
417: /* -------------------------------------------------------------------------- */
418: /*
419: SNESSetFromOptions_VINEWTONSSLS - Sets various parameters for the SNESVI method.
421: Input Parameter:
422: . snes - the SNES context
424: Application Interface Routine: SNESSetFromOptions()
425: */
426: static PetscErrorCode SNESSetFromOptions_VINEWTONSSLS(PetscOptionItems *PetscOptionsObject,SNES snes)
427: {
429: SNESLineSearch linesearch;
432: SNESSetFromOptions_VI(PetscOptionsObject,snes);
433: PetscOptionsHead(PetscOptionsObject,"SNES semismooth method options");
434: PetscOptionsTail();
435: /* set up the default line search */
436: if (!snes->linesearch) {
437: SNESGetLineSearch(snes, &linesearch);
438: SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
439: SNESLineSearchBTSetAlpha(linesearch, 0.0);
440: }
441: return(0);
442: }
445: /* -------------------------------------------------------------------------- */
446: /*MC
447: SNESVINEWTONSSLS - Semi-smooth solver for variational inequalities based on Newton's method
449: Options Database:
450: + -snes_type <vinewtonssls,vinewtonrsls> a semi-smooth solver, a reduced space active set method
451: - -snes_vi_monitor - prints the number of active constraints at each iteration.
453: Level: beginner
455: References:
456: + 1. - T. S. Munson, F. Facchinei, M. C. Ferris, A. Fischer, and C. Kanzow. The semismooth
457: algorithm for large scale complementarity problems. INFORMS Journal on Computing, 13 (2001).
458: - 2. - T. S. Munson, and S. Benson. Flexible Complementarity Solvers for Large Scale
459: Applications, Optimization Methods and Software, 21 (2006).
461: .seealso: SNESVISetVariableBounds(), SNESVISetComputeVariableBounds(), SNESCreate(), SNES, SNESSetType(), SNESVINEWTONRSLS, SNESNEWTONTR, SNESLineSearchSetType(),SNESLineSearchSetPostCheck(), SNESLineSearchSetPreCheck()
463: M*/
464: PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONSSLS(SNES snes)
465: {
466: PetscErrorCode ierr;
467: SNES_VINEWTONSSLS *vi;
470: snes->ops->reset = SNESReset_VINEWTONSSLS;
471: snes->ops->setup = SNESSetUp_VINEWTONSSLS;
472: snes->ops->solve = SNESSolve_VINEWTONSSLS;
473: snes->ops->destroy = SNESDestroy_VI;
474: snes->ops->setfromoptions = SNESSetFromOptions_VINEWTONSSLS;
475: snes->ops->view = NULL;
477: snes->usesksp = PETSC_TRUE;
478: snes->usesnpc = PETSC_FALSE;
480: snes->alwayscomputesfinalresidual = PETSC_FALSE;
482: PetscNewLog(snes,&vi);
483: snes->data = (void*)vi;
485: PetscObjectComposeFunction((PetscObject)snes,"SNESVISetVariableBounds_C",SNESVISetVariableBounds_VI);
486: PetscObjectComposeFunction((PetscObject)snes,"SNESVISetComputeVariableBounds_C",SNESVISetComputeVariableBounds_VI);
487: return(0);
488: }