Actual source code: viss.c
petsc-3.6.1 2015-08-06
2: #include <../src/snes/impls/vi/ss/vissimpl.h> /*I "petscsnes.h" I*/
3: #include <../include/petsc/private/kspimpl.h>
4: #include <../include/petsc/private/matimpl.h>
5: #include <../include/petsc/private/dmimpl.h>
8: /*
9: SNESVIComputeMeritFunction - Evaluates the merit function for the mixed complementarity problem.
11: Input Parameter:
12: . phi - the semismooth function
14: Output Parameter:
15: . merit - the merit function
16: . phinorm - ||phi||
18: Notes:
19: The merit function for the mixed complementarity problem is defined as
20: merit = 0.5*phi^T*phi
21: */
24: static PetscErrorCode SNESVIComputeMeritFunction(Vec phi, PetscReal *merit,PetscReal *phinorm)
25: {
29: VecNormBegin(phi,NORM_2,phinorm);
30: VecNormEnd(phi,NORM_2,phinorm);
32: *merit = 0.5*(*phinorm)*(*phinorm);
33: return(0);
34: }
36: PETSC_STATIC_INLINE PetscScalar Phi(PetscScalar a,PetscScalar b)
37: {
38: return a + b - PetscSqrtScalar(a*a + b*b);
39: }
41: PETSC_STATIC_INLINE PetscScalar DPhi(PetscScalar a,PetscScalar b)
42: {
43: if ((PetscAbsScalar(a) >= 1.e-6) || (PetscAbsScalar(b) >= 1.e-6)) return 1.0 - a/ PetscSqrtScalar(a*a + b*b);
44: else return .5;
45: }
47: /*
48: SNESVIComputeFunction - Reformulates a system of nonlinear equations in mixed complementarity form to a system of nonlinear equations in semismooth form.
50: Input Parameters:
51: . snes - the SNES context
52: . X - current iterate
53: . functx - user defined function context
55: Output Parameters:
56: . phi - Semismooth function
58: */
61: static PetscErrorCode SNESVIComputeFunction(SNES snes,Vec X,Vec phi,void *functx)
62: {
63: PetscErrorCode ierr;
64: SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS*)snes->data;
65: Vec Xl = snes->xl,Xu = snes->xu,F = snes->vec_func;
66: PetscScalar *phi_arr,*f_arr,*l,*u;
67: const PetscScalar *x_arr;
68: PetscInt i,nlocal;
71: (*vi->computeuserfunction)(snes,X,F,functx);
72: VecGetLocalSize(X,&nlocal);
73: VecGetArrayRead(X,&x_arr);
74: VecGetArray(F,&f_arr);
75: VecGetArray(Xl,&l);
76: VecGetArray(Xu,&u);
77: VecGetArray(phi,&phi_arr);
79: for (i=0; i < nlocal; i++) {
80: if ((PetscRealPart(l[i]) <= PETSC_NINFINITY) && (PetscRealPart(u[i]) >= PETSC_INFINITY)) { /* no constraints on variable */
81: phi_arr[i] = f_arr[i];
82: } else if (PetscRealPart(l[i]) <= PETSC_NINFINITY) { /* upper bound on variable only */
83: phi_arr[i] = -Phi(u[i] - x_arr[i],-f_arr[i]);
84: } else if (PetscRealPart(u[i]) >= PETSC_INFINITY) { /* lower bound on variable only */
85: phi_arr[i] = Phi(x_arr[i] - l[i],f_arr[i]);
86: } else if (l[i] == u[i]) {
87: phi_arr[i] = l[i] - x_arr[i];
88: } else { /* both bounds on variable */
89: phi_arr[i] = Phi(x_arr[i] - l[i],-Phi(u[i] - x_arr[i],-f_arr[i]));
90: }
91: }
93: VecRestoreArrayRead(X,&x_arr);
94: VecRestoreArray(F,&f_arr);
95: VecRestoreArray(Xl,&l);
96: VecRestoreArray(Xu,&u);
97: VecRestoreArray(phi,&phi_arr);
98: return(0);
99: }
101: /*
102: SNESVIComputeBsubdifferentialVectors - Computes the diagonal shift (Da) and row scaling (Db) vectors needed for the
103: the semismooth jacobian.
104: */
107: PetscErrorCode SNESVIComputeBsubdifferentialVectors(SNES snes,Vec X,Vec F,Mat jac,Vec Da,Vec Db)
108: {
110: PetscScalar *l,*u,*x,*f,*da,*db,da1,da2,db1,db2;
111: PetscInt i,nlocal;
114: VecGetArray(X,&x);
115: VecGetArray(F,&f);
116: VecGetArray(snes->xl,&l);
117: VecGetArray(snes->xu,&u);
118: VecGetArray(Da,&da);
119: VecGetArray(Db,&db);
120: VecGetLocalSize(X,&nlocal);
122: for (i=0; i< nlocal; i++) {
123: if ((PetscRealPart(l[i]) <= PETSC_NINFINITY) && (PetscRealPart(u[i]) >= PETSC_INFINITY)) { /* no constraints on variable */
124: da[i] = 0;
125: db[i] = 1;
126: } else if (PetscRealPart(l[i]) <= PETSC_NINFINITY) { /* upper bound on variable only */
127: da[i] = DPhi(u[i] - x[i], -f[i]);
128: db[i] = DPhi(-f[i],u[i] - x[i]);
129: } else if (PetscRealPart(u[i]) >= PETSC_INFINITY) { /* lower bound on variable only */
130: da[i] = DPhi(x[i] - l[i], f[i]);
131: db[i] = DPhi(f[i],x[i] - l[i]);
132: } else if (l[i] == u[i]) { /* fixed variable */
133: da[i] = 1;
134: db[i] = 0;
135: } else { /* upper and lower bounds on variable */
136: da1 = DPhi(x[i] - l[i], -Phi(u[i] - x[i], -f[i]));
137: db1 = DPhi(-Phi(u[i] - x[i], -f[i]),x[i] - l[i]);
138: da2 = DPhi(u[i] - x[i], -f[i]);
139: db2 = DPhi(-f[i],u[i] - x[i]);
140: da[i] = da1 + db1*da2;
141: db[i] = db1*db2;
142: }
143: }
145: VecRestoreArray(X,&x);
146: VecRestoreArray(F,&f);
147: VecRestoreArray(snes->xl,&l);
148: VecRestoreArray(snes->xu,&u);
149: VecRestoreArray(Da,&da);
150: VecRestoreArray(Db,&db);
151: return(0);
152: }
154: /*
155: 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.
157: Input Parameters:
158: . Da - Diagonal shift vector for the semismooth jacobian.
159: . Db - Row scaling vector for the semismooth jacobian.
161: Output Parameters:
162: . jac - semismooth jacobian
163: . jac_pre - optional preconditioning matrix
165: Notes:
166: The semismooth jacobian matrix is given by
167: jac = Da + Db*jacfun
168: where Db is the row scaling matrix stored as a vector,
169: Da is the diagonal perturbation matrix stored as a vector
170: and jacfun is the jacobian of the original nonlinear function.
171: */
174: PetscErrorCode SNESVIComputeJacobian(Mat jac, Mat jac_pre,Vec Da, Vec Db)
175: {
178: /* Do row scaling and add diagonal perturbation */
179: MatDiagonalScale(jac,Db,NULL);
180: MatDiagonalSet(jac,Da,ADD_VALUES);
181: if (jac != jac_pre) { /* If jac and jac_pre are different */
182: MatDiagonalScale(jac_pre,Db,NULL);
183: MatDiagonalSet(jac_pre,Da,ADD_VALUES);
184: }
185: return(0);
186: }
188: /*
189: SNESVIComputeMeritFunctionGradient - Computes the gradient of the merit function psi.
191: Input Parameters:
192: phi - semismooth function.
193: H - semismooth jacobian
195: Output Parameters:
196: dpsi - merit function gradient
198: Notes:
199: The merit function gradient is computed as follows
200: dpsi = H^T*phi
201: */
204: PetscErrorCode SNESVIComputeMeritFunctionGradient(Mat H, Vec phi, Vec dpsi)
205: {
209: MatMultTranspose(H,phi,dpsi);
210: return(0);
211: }
215: /*
216: SNESSolve_VINEWTONSSLS - Solves the complementarity problem with a semismooth Newton
217: method using a line search.
219: Input Parameters:
220: . snes - the SNES context
222: Application Interface Routine: SNESSolve()
224: Notes:
225: This implements essentially a semismooth Newton method with a
226: line search. The default line search does not do any line search
227: but rather takes a full Newton step.
229: 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
231: */
234: PetscErrorCode SNESSolve_VINEWTONSSLS(SNES snes)
235: {
236: SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS*)snes->data;
237: PetscErrorCode ierr;
238: PetscInt maxits,i,lits;
239: SNESLineSearchReason lssucceed;
240: PetscReal gnorm,xnorm=0,ynorm;
241: Vec Y,X,F;
242: KSPConvergedReason kspreason;
243: DM dm;
244: DMSNES sdm;
247: SNESGetDM(snes,&dm);
248: DMGetDMSNES(dm,&sdm);
250: vi->computeuserfunction = sdm->ops->computefunction;
251: sdm->ops->computefunction = SNESVIComputeFunction;
253: snes->numFailures = 0;
254: snes->numLinearSolveFailures = 0;
255: snes->reason = SNES_CONVERGED_ITERATING;
257: maxits = snes->max_its; /* maximum number of iterations */
258: X = snes->vec_sol; /* solution vector */
259: F = snes->vec_func; /* residual vector */
260: Y = snes->work[0]; /* work vectors */
262: PetscObjectSAWsTakeAccess((PetscObject)snes);
263: snes->iter = 0;
264: snes->norm = 0.0;
265: PetscObjectSAWsGrantAccess((PetscObject)snes);
267: SNESVIProjectOntoBounds(snes,X);
268: SNESComputeFunction(snes,X,vi->phi);
269: if (snes->domainerror) {
270: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
271: sdm->ops->computefunction = vi->computeuserfunction;
272: return(0);
273: }
274: /* Compute Merit function */
275: SNESVIComputeMeritFunction(vi->phi,&vi->merit,&vi->phinorm);
277: VecNormBegin(X,NORM_2,&xnorm); /* xnorm <- ||x|| */
278: VecNormEnd(X,NORM_2,&xnorm);
279: SNESCheckFunctionNorm(snes,vi->merit);
281: PetscObjectSAWsTakeAccess((PetscObject)snes);
282: snes->norm = vi->phinorm;
283: PetscObjectSAWsGrantAccess((PetscObject)snes);
284: SNESLogConvergenceHistory(snes,vi->phinorm,0);
285: SNESMonitor(snes,0,vi->phinorm);
287: /* test convergence */
288: (*snes->ops->converged)(snes,0,0.0,0.0,vi->phinorm,&snes->reason,snes->cnvP);
289: if (snes->reason) {
290: sdm->ops->computefunction = vi->computeuserfunction;
291: return(0);
292: }
294: for (i=0; i<maxits; i++) {
296: /* Call general purpose update function */
297: if (snes->ops->update) {
298: (*snes->ops->update)(snes, snes->iter);
299: }
301: /* Solve J Y = Phi, where J is the semismooth jacobian */
303: /* Get the jacobian -- note that the function must be the original function for snes_fd and snes_fd_color to work for this*/
304: sdm->ops->computefunction = vi->computeuserfunction;
305: SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
306: sdm->ops->computefunction = SNESVIComputeFunction;
308: /* Get the diagonal shift and row scaling vectors */
309: SNESVIComputeBsubdifferentialVectors(snes,X,F,snes->jacobian,vi->Da,vi->Db);
310: /* Compute the semismooth jacobian */
311: SNESVIComputeJacobian(snes->jacobian,snes->jacobian_pre,vi->Da,vi->Db);
312: /* Compute the merit function gradient */
313: SNESVIComputeMeritFunctionGradient(snes->jacobian,vi->phi,vi->dpsi);
314: KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);
315: KSPSolve(snes->ksp,vi->phi,Y);
316: KSPGetConvergedReason(snes->ksp,&kspreason);
318: if (kspreason < 0) {
319: if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
320: PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);
321: snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
322: break;
323: }
324: }
325: KSPGetIterationNumber(snes->ksp,&lits);
326: snes->linear_its += lits;
327: PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);
328: /*
329: if (snes->ops->precheck) {
330: PetscBool changed_y = PETSC_FALSE;
331: (*snes->ops->precheck)(snes,X,Y,snes->precheck,&changed_y);
332: }
334: if (PetscLogPrintInfo) {
335: SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);
336: }
337: */
338: /* Compute a (scaled) negative update in the line search routine:
339: Y <- X - lambda*Y
340: and evaluate G = function(Y) (depends on the line search).
341: */
342: VecCopy(Y,snes->vec_sol_update);
343: ynorm = 1; gnorm = vi->phinorm;
344: SNESLineSearchApply(snes->linesearch, X, vi->phi, &gnorm, Y);
345: SNESLineSearchGetReason(snes->linesearch, &lssucceed);
346: SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm);
347: PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",(double)vi->phinorm,(double)gnorm,(double)ynorm,(int)lssucceed);
348: if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
349: if (snes->domainerror) {
350: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
351: sdm->ops->computefunction = vi->computeuserfunction;
352: return(0);
353: }
354: if (lssucceed) {
355: if (++snes->numFailures >= snes->maxFailures) {
356: PetscBool ismin;
357: snes->reason = SNES_DIVERGED_LINE_SEARCH;
358: SNESVICheckLocalMin_Private(snes,snes->jacobian,vi->phi,X,gnorm,&ismin);
359: if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
360: break;
361: }
362: }
363: /* Update function and solution vectors */
364: vi->phinorm = gnorm;
365: vi->merit = 0.5*vi->phinorm*vi->phinorm;
366: /* Monitor convergence */
367: PetscObjectSAWsTakeAccess((PetscObject)snes);
368: snes->iter = i+1;
369: snes->norm = vi->phinorm;
370: PetscObjectSAWsGrantAccess((PetscObject)snes);
371: SNESLogConvergenceHistory(snes,snes->norm,lits);
372: SNESMonitor(snes,snes->iter,snes->norm);
373: /* Test for convergence, xnorm = || X || */
374: if (snes->ops->converged != SNESConvergedSkip) { VecNorm(X,NORM_2,&xnorm); }
375: (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,vi->phinorm,&snes->reason,snes->cnvP);
376: if (snes->reason) break;
377: }
378: if (i == maxits) {
379: PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);
380: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
381: }
382: sdm->ops->computefunction = vi->computeuserfunction;
383: return(0);
384: }
386: /* -------------------------------------------------------------------------- */
387: /*
388: SNESSetUp_VINEWTONSSLS - Sets up the internal data structures for the later use
389: of the SNES nonlinear solver.
391: Input Parameter:
392: . snes - the SNES context
394: Application Interface Routine: SNESSetUp()
396: Notes:
397: For basic use of the SNES solvers, the user need not explicitly call
398: SNESSetUp(), since these actions will automatically occur during
399: the call to SNESSolve().
400: */
403: PetscErrorCode SNESSetUp_VINEWTONSSLS(SNES snes)
404: {
405: PetscErrorCode ierr;
406: SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS*) snes->data;
409: SNESSetUp_VI(snes);
410: VecDuplicate(snes->vec_sol, &vi->dpsi);
411: VecDuplicate(snes->vec_sol, &vi->phi);
412: VecDuplicate(snes->vec_sol, &vi->Da);
413: VecDuplicate(snes->vec_sol, &vi->Db);
414: VecDuplicate(snes->vec_sol, &vi->z);
415: VecDuplicate(snes->vec_sol, &vi->t);
416: return(0);
417: }
418: /* -------------------------------------------------------------------------- */
421: PetscErrorCode SNESReset_VINEWTONSSLS(SNES snes)
422: {
423: SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS*) snes->data;
424: PetscErrorCode ierr;
427: SNESReset_VI(snes);
428: VecDestroy(&vi->dpsi);
429: VecDestroy(&vi->phi);
430: VecDestroy(&vi->Da);
431: VecDestroy(&vi->Db);
432: VecDestroy(&vi->z);
433: VecDestroy(&vi->t);
434: return(0);
435: }
437: /* -------------------------------------------------------------------------- */
438: /*
439: SNESSetFromOptions_VINEWTONSSLS - Sets various parameters for the SNESVI method.
441: Input Parameter:
442: . snes - the SNES context
444: Application Interface Routine: SNESSetFromOptions()
445: */
448: static PetscErrorCode SNESSetFromOptions_VINEWTONSSLS(PetscOptions *PetscOptionsObject,SNES snes)
449: {
451: SNESLineSearch linesearch;
454: SNESSetFromOptions_VI(PetscOptionsObject,snes);
455: PetscOptionsHead(PetscOptionsObject,"SNES semismooth method options");
456: PetscOptionsTail();
457: /* set up the default line search */
458: if (!snes->linesearch) {
459: SNESGetLineSearch(snes, &linesearch);
460: SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
461: SNESLineSearchBTSetAlpha(linesearch, 0.0);
462: }
463: return(0);
464: }
467: /* -------------------------------------------------------------------------- */
468: /*MC
469: SNESVINEWTONSSLS - Semi-smooth solver for variational inequalities based on Newton's method
471: Options Database:
472: + -snes_type <vinewtonssls,vinewtonrsls> a semi-smooth solver, a reduced space active set method
473: - -snes_vi_monitor - prints the number of active constraints at each iteration.
475: Level: beginner
477: References:
478: - T. S. Munson, F. Facchinei, M. C. Ferris, A. Fischer, and C. Kanzow. The semismooth
479: algorithm for large scale complementarity problems. INFORMS Journal on Computing, 13 (2001).
480: - T. S. Munson, and S. Benson. Flexible Complementarity Solvers for Large-Scale
481: Applications, Optimization Methods and Software, 21 (2006).
483: .seealso: SNESVISetVariableBounds(), SNESVISetComputeVariableBounds(), SNESCreate(), SNES, SNESSetType(), SNESVINEWTONRSLS, SNESNEWTONTR, SNESLineSearchSet(),SNESLineSearchSetPostCheck(), SNESLineSearchSetPreCheck()
485: M*/
488: PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONSSLS(SNES snes)
489: {
490: PetscErrorCode ierr;
491: SNES_VINEWTONSSLS *vi;
494: snes->ops->reset = SNESReset_VINEWTONSSLS;
495: snes->ops->setup = SNESSetUp_VINEWTONSSLS;
496: snes->ops->solve = SNESSolve_VINEWTONSSLS;
497: snes->ops->destroy = SNESDestroy_VI;
498: snes->ops->setfromoptions = SNESSetFromOptions_VINEWTONSSLS;
499: snes->ops->view = NULL;
501: snes->usesksp = PETSC_TRUE;
502: snes->usespc = PETSC_FALSE;
504: PetscNewLog(snes,&vi);
505: snes->data = (void*)vi;
507: PetscObjectComposeFunction((PetscObject)snes,"SNESVISetVariableBounds_C",SNESVISetVariableBounds_VI);
508: PetscObjectComposeFunction((PetscObject)snes,"SNESVISetComputeVariableBounds_C",SNESVISetComputeVariableBounds_VI);
509: return(0);
510: }