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: {
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 */
168: MatDiagonalScale(jac,Db,NULL);
169: MatDiagonalSet(jac,Da,ADD_VALUES);
170: if (jac != jac_pre) { /* If jac and jac_pre are different */
171: MatDiagonalScale(jac_pre,Db,NULL);
172: MatDiagonalSet(jac_pre,Da,ADD_VALUES);
173: }
174: return(0);
175: }
177: /*
178: SNESVIComputeMeritFunctionGradient - Computes the gradient of the merit function psi.
180: Input Parameters:
181: phi - semismooth function.
182: H - semismooth jacobian
184: Output Parameters:
185: dpsi - merit function gradient
187: Notes:
188: The merit function gradient is computed as follows
189: dpsi = H^T*phi
190: */
191: PetscErrorCode SNESVIComputeMeritFunctionGradient(Mat H, Vec phi, Vec dpsi)
192: {
196: MatMultTranspose(H,phi,dpsi);
197: return(0);
198: }
200: /*
201: SNESSolve_VINEWTONSSLS - Solves the complementarity problem with a semismooth Newton
202: method using a line search.
204: Input Parameters:
205: . snes - the SNES context
207: Application Interface Routine: SNESSolve()
209: Notes:
210: This implements essentially a semismooth Newton method with a
211: line search. The default line search does not do any line search
212: but rather takes a full Newton step.
214: 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
216: */
217: PetscErrorCode SNESSolve_VINEWTONSSLS(SNES snes)
218: {
219: SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS*)snes->data;
220: PetscErrorCode ierr;
221: PetscInt maxits,i,lits;
222: SNESLineSearchReason lssucceed;
223: PetscReal gnorm,xnorm=0,ynorm;
224: Vec Y,X,F;
225: KSPConvergedReason kspreason;
226: DM dm;
227: DMSNES sdm;
230: SNESGetDM(snes,&dm);
231: DMGetDMSNES(dm,&sdm);
233: vi->computeuserfunction = sdm->ops->computefunction;
234: sdm->ops->computefunction = SNESVIComputeFunction;
236: snes->numFailures = 0;
237: snes->numLinearSolveFailures = 0;
238: snes->reason = SNES_CONVERGED_ITERATING;
240: maxits = snes->max_its; /* maximum number of iterations */
241: X = snes->vec_sol; /* solution vector */
242: F = snes->vec_func; /* residual vector */
243: Y = snes->work[0]; /* work vectors */
245: PetscObjectSAWsTakeAccess((PetscObject)snes);
246: snes->iter = 0;
247: snes->norm = 0.0;
248: PetscObjectSAWsGrantAccess((PetscObject)snes);
250: SNESVIProjectOntoBounds(snes,X);
251: SNESComputeFunction(snes,X,vi->phi);
252: if (snes->domainerror) {
253: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
254: sdm->ops->computefunction = vi->computeuserfunction;
255: return(0);
256: }
257: /* Compute Merit function */
258: SNESVIComputeMeritFunction(vi->phi,&vi->merit,&vi->phinorm);
260: VecNormBegin(X,NORM_2,&xnorm); /* xnorm <- ||x|| */
261: VecNormEnd(X,NORM_2,&xnorm);
262: SNESCheckFunctionNorm(snes,vi->merit);
264: PetscObjectSAWsTakeAccess((PetscObject)snes);
265: snes->norm = vi->phinorm;
266: PetscObjectSAWsGrantAccess((PetscObject)snes);
267: SNESLogConvergenceHistory(snes,vi->phinorm,0);
268: SNESMonitor(snes,0,vi->phinorm);
270: /* test convergence */
271: (*snes->ops->converged)(snes,0,0.0,0.0,vi->phinorm,&snes->reason,snes->cnvP);
272: if (snes->reason) {
273: sdm->ops->computefunction = vi->computeuserfunction;
274: return(0);
275: }
277: for (i=0; i<maxits; i++) {
279: /* Call general purpose update function */
280: if (snes->ops->update) {
281: (*snes->ops->update)(snes, snes->iter);
282: }
284: /* Solve J Y = Phi, where J is the semismooth jacobian */
286: /* Get the jacobian -- note that the function must be the original function for snes_fd and snes_fd_color to work for this*/
287: sdm->ops->computefunction = vi->computeuserfunction;
288: SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
289: SNESCheckJacobianDomainerror(snes);
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: snes->xnorm = xnorm;
355: snes->ynorm = ynorm;
356: PetscObjectSAWsGrantAccess((PetscObject)snes);
357: SNESLogConvergenceHistory(snes,snes->norm,lits);
358: SNESMonitor(snes,snes->iter,snes->norm);
359: /* Test for convergence, xnorm = || X || */
360: if (snes->ops->converged != SNESConvergedSkip) { VecNorm(X,NORM_2,&xnorm); }
361: (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,vi->phinorm,&snes->reason,snes->cnvP);
362: if (snes->reason) break;
363: }
364: if (i == maxits) {
365: PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);
366: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
367: }
368: sdm->ops->computefunction = vi->computeuserfunction;
369: return(0);
370: }
372: /* -------------------------------------------------------------------------- */
373: /*
374: SNESSetUp_VINEWTONSSLS - Sets up the internal data structures for the later use
375: of the SNES nonlinear solver.
377: Input Parameter:
378: . snes - the SNES context
380: Application Interface Routine: SNESSetUp()
382: Notes:
383: For basic use of the SNES solvers, the user need not explicitly call
384: SNESSetUp(), since these actions will automatically occur during
385: the call to SNESSolve().
386: */
387: PetscErrorCode SNESSetUp_VINEWTONSSLS(SNES snes)
388: {
389: PetscErrorCode ierr;
390: SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS*) snes->data;
393: SNESSetUp_VI(snes);
394: VecDuplicate(snes->vec_sol, &vi->dpsi);
395: VecDuplicate(snes->vec_sol, &vi->phi);
396: VecDuplicate(snes->vec_sol, &vi->Da);
397: VecDuplicate(snes->vec_sol, &vi->Db);
398: VecDuplicate(snes->vec_sol, &vi->z);
399: VecDuplicate(snes->vec_sol, &vi->t);
400: return(0);
401: }
402: /* -------------------------------------------------------------------------- */
403: PetscErrorCode SNESReset_VINEWTONSSLS(SNES snes)
404: {
405: SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS*) snes->data;
406: PetscErrorCode ierr;
409: SNESReset_VI(snes);
410: VecDestroy(&vi->dpsi);
411: VecDestroy(&vi->phi);
412: VecDestroy(&vi->Da);
413: VecDestroy(&vi->Db);
414: VecDestroy(&vi->z);
415: VecDestroy(&vi->t);
416: return(0);
417: }
419: /* -------------------------------------------------------------------------- */
420: /*
421: SNESSetFromOptions_VINEWTONSSLS - Sets various parameters for the SNESVI method.
423: Input Parameter:
424: . snes - the SNES context
426: Application Interface Routine: SNESSetFromOptions()
427: */
428: static PetscErrorCode SNESSetFromOptions_VINEWTONSSLS(PetscOptionItems *PetscOptionsObject,SNES snes)
429: {
433: SNESSetFromOptions_VI(PetscOptionsObject,snes);
434: PetscOptionsHead(PetscOptionsObject,"SNES semismooth method options");
435: PetscOptionsTail();
436: return(0);
437: }
439: /* -------------------------------------------------------------------------- */
440: /*MC
441: SNESVINEWTONSSLS - Semi-smooth solver for variational inequalities based on Newton's method
443: Options Database:
444: + -snes_type <vinewtonssls,vinewtonrsls> a semi-smooth solver, a reduced space active set method
445: - -snes_vi_monitor - prints the number of active constraints at each iteration.
447: Level: beginner
449: References:
450: + 1. - T. S. Munson, F. Facchinei, M. C. Ferris, A. Fischer, and C. Kanzow. The semismooth
451: algorithm for large scale complementarity problems. INFORMS Journal on Computing, 13 (2001).
452: - 2. - T. S. Munson, and S. Benson. Flexible Complementarity Solvers for Large Scale
453: Applications, Optimization Methods and Software, 21 (2006).
455: .seealso: SNESVISetVariableBounds(), SNESVISetComputeVariableBounds(), SNESCreate(), SNES, SNESSetType(), SNESVINEWTONRSLS, SNESNEWTONTR, SNESLineSearchSetType(),SNESLineSearchSetPostCheck(), SNESLineSearchSetPreCheck()
457: M*/
458: PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONSSLS(SNES snes)
459: {
460: PetscErrorCode ierr;
461: SNES_VINEWTONSSLS *vi;
462: SNESLineSearch linesearch;
465: snes->ops->reset = SNESReset_VINEWTONSSLS;
466: snes->ops->setup = SNESSetUp_VINEWTONSSLS;
467: snes->ops->solve = SNESSolve_VINEWTONSSLS;
468: snes->ops->destroy = SNESDestroy_VI;
469: snes->ops->setfromoptions = SNESSetFromOptions_VINEWTONSSLS;
470: snes->ops->view = NULL;
472: snes->usesksp = PETSC_TRUE;
473: snes->usesnpc = PETSC_FALSE;
475: SNESGetLineSearch(snes, &linesearch);
476: if (!((PetscObject)linesearch)->type_name) {
477: SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
478: SNESLineSearchBTSetAlpha(linesearch, 0.0);
479: }
481: snes->alwayscomputesfinalresidual = PETSC_FALSE;
483: PetscNewLog(snes,&vi);
484: snes->data = (void*)vi;
486: PetscObjectComposeFunction((PetscObject)snes,"SNESVISetVariableBounds_C",SNESVISetVariableBounds_VI);
487: PetscObjectComposeFunction((PetscObject)snes,"SNESVISetComputeVariableBounds_C",SNESVISetComputeVariableBounds_VI);
488: return(0);
489: }