Actual source code: viss.c
petsc-3.4.5 2014-06-29
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,*x_arr,*f_arr,*l,*u;
67: PetscInt i,nlocal;
70: (*vi->computeuserfunction)(snes,X,F,functx);
71: VecGetLocalSize(X,&nlocal);
72: VecGetArray(X,&x_arr);
73: VecGetArray(F,&f_arr);
74: VecGetArray(Xl,&l);
75: VecGetArray(Xu,&u);
76: VecGetArray(phi,&phi_arr);
78: for (i=0; i < nlocal; i++) {
79: if ((PetscRealPart(l[i]) <= SNES_VI_NINF) && (PetscRealPart(u[i]) >= SNES_VI_INF)) { /* no constraints on variable */
80: phi_arr[i] = f_arr[i];
81: } else if (PetscRealPart(l[i]) <= SNES_VI_NINF) { /* upper bound on variable only */
82: phi_arr[i] = -Phi(u[i] - x_arr[i],-f_arr[i]);
83: } else if (PetscRealPart(u[i]) >= SNES_VI_INF) { /* lower bound on variable only */
84: phi_arr[i] = Phi(x_arr[i] - l[i],f_arr[i]);
85: } else if (l[i] == u[i]) {
86: phi_arr[i] = l[i] - x_arr[i];
87: } else { /* both bounds on variable */
88: phi_arr[i] = Phi(x_arr[i] - l[i],-Phi(u[i] - x_arr[i],-f_arr[i]));
89: }
90: }
92: VecRestoreArray(X,&x_arr);
93: VecRestoreArray(F,&f_arr);
94: VecRestoreArray(Xl,&l);
95: VecRestoreArray(Xu,&u);
96: VecRestoreArray(phi,&phi_arr);
97: return(0);
98: }
100: /*
101: SNESVIComputeBsubdifferentialVectors - Computes the diagonal shift (Da) and row scaling (Db) vectors needed for the
102: the semismooth jacobian.
103: */
106: PetscErrorCode SNESVIComputeBsubdifferentialVectors(SNES snes,Vec X,Vec F,Mat jac,Vec Da,Vec Db)
107: {
109: PetscScalar *l,*u,*x,*f,*da,*db,da1,da2,db1,db2;
110: PetscInt i,nlocal;
113: VecGetArray(X,&x);
114: VecGetArray(F,&f);
115: VecGetArray(snes->xl,&l);
116: VecGetArray(snes->xu,&u);
117: VecGetArray(Da,&da);
118: VecGetArray(Db,&db);
119: VecGetLocalSize(X,&nlocal);
121: for (i=0; i< nlocal; i++) {
122: if ((PetscRealPart(l[i]) <= SNES_VI_NINF) && (PetscRealPart(u[i]) >= SNES_VI_INF)) { /* no constraints on variable */
123: da[i] = 0;
124: db[i] = 1;
125: } else if (PetscRealPart(l[i]) <= SNES_VI_NINF) { /* upper bound on variable only */
126: da[i] = DPhi(u[i] - x[i], -f[i]);
127: db[i] = DPhi(-f[i],u[i] - x[i]);
128: } else if (PetscRealPart(u[i]) >= SNES_VI_INF) { /* lower bound on variable only */
129: da[i] = DPhi(x[i] - l[i], f[i]);
130: db[i] = DPhi(f[i],x[i] - l[i]);
131: } else if (l[i] == u[i]) { /* fixed variable */
132: da[i] = 1;
133: db[i] = 0;
134: } else { /* upper and lower bounds on variable */
135: da1 = DPhi(x[i] - l[i], -Phi(u[i] - x[i], -f[i]));
136: db1 = DPhi(-Phi(u[i] - x[i], -f[i]),x[i] - l[i]);
137: da2 = DPhi(u[i] - x[i], -f[i]);
138: db2 = DPhi(-f[i],u[i] - x[i]);
139: da[i] = da1 + db1*da2;
140: db[i] = db1*db2;
141: }
142: }
144: VecRestoreArray(X,&x);
145: VecRestoreArray(F,&f);
146: VecRestoreArray(snes->xl,&l);
147: VecRestoreArray(snes->xu,&u);
148: VecRestoreArray(Da,&da);
149: VecRestoreArray(Db,&db);
150: return(0);
151: }
153: /*
154: 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.
156: Input Parameters:
157: . Da - Diagonal shift vector for the semismooth jacobian.
158: . Db - Row scaling vector for the semismooth jacobian.
160: Output Parameters:
161: . jac - semismooth jacobian
162: . jac_pre - optional preconditioning matrix
164: Notes:
165: The semismooth jacobian matrix is given by
166: jac = Da + Db*jacfun
167: where Db is the row scaling matrix stored as a vector,
168: Da is the diagonal perturbation matrix stored as a vector
169: and jacfun is the jacobian of the original nonlinear function.
170: */
173: PetscErrorCode SNESVIComputeJacobian(Mat jac, Mat jac_pre,Vec Da, Vec Db)
174: {
177: /* Do row scaling and add diagonal perturbation */
178: MatDiagonalScale(jac,Db,NULL);
179: MatDiagonalSet(jac,Da,ADD_VALUES);
180: if (jac != jac_pre) { /* If jac and jac_pre are different */
181: MatDiagonalScale(jac_pre,Db,NULL);
182: MatDiagonalSet(jac_pre,Da,ADD_VALUES);
183: }
184: return(0);
185: }
187: /*
188: SNESVIComputeMeritFunctionGradient - Computes the gradient of the merit function psi.
190: Input Parameters:
191: phi - semismooth function.
192: H - semismooth jacobian
194: Output Parameters:
195: dpsi - merit function gradient
197: Notes:
198: The merit function gradient is computed as follows
199: dpsi = H^T*phi
200: */
203: PetscErrorCode SNESVIComputeMeritFunctionGradient(Mat H, Vec phi, Vec dpsi)
204: {
208: MatMultTranspose(H,phi,dpsi);
209: return(0);
210: }
214: /*
215: SNESSolve_VINEWTONSSLS - Solves the complementarity problem with a semismooth Newton
216: method using a line search.
218: Input Parameters:
219: . snes - the SNES context
221: Output Parameter:
222: . outits - number of iterations until termination
224: Application Interface Routine: SNESSolve()
226: Notes:
227: This implements essentially a semismooth Newton method with a
228: line search. The default line search does not do any line seach
229: but rather takes a full newton step.
231: 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
233: */
236: PetscErrorCode SNESSolve_VINEWTONSSLS(SNES snes)
237: {
238: SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS*)snes->data;
239: PetscErrorCode ierr;
240: PetscInt maxits,i,lits;
241: PetscBool lssucceed;
242: MatStructure flg = DIFFERENT_NONZERO_PATTERN;
243: PetscReal gnorm,xnorm=0,ynorm;
244: Vec Y,X,F;
245: KSPConvergedReason kspreason;
246: DM dm;
247: DMSNES sdm;
250: SNESGetDM(snes,&dm);
251: DMGetDMSNES(dm,&sdm);
253: vi->computeuserfunction = sdm->ops->computefunction;
254: sdm->ops->computefunction = SNESVIComputeFunction;
256: snes->numFailures = 0;
257: snes->numLinearSolveFailures = 0;
258: snes->reason = SNES_CONVERGED_ITERATING;
260: maxits = snes->max_its; /* maximum number of iterations */
261: X = snes->vec_sol; /* solution vector */
262: F = snes->vec_func; /* residual vector */
263: Y = snes->work[0]; /* work vectors */
265: PetscObjectAMSTakeAccess((PetscObject)snes);
266: snes->iter = 0;
267: snes->norm = 0.0;
268: PetscObjectAMSGrantAccess((PetscObject)snes);
270: SNESVIProjectOntoBounds(snes,X);
271: SNESComputeFunction(snes,X,vi->phi);
272: if (snes->domainerror) {
273: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
274: sdm->ops->computefunction = vi->computeuserfunction;
275: return(0);
276: }
277: /* Compute Merit function */
278: SNESVIComputeMeritFunction(vi->phi,&vi->merit,&vi->phinorm);
280: VecNormBegin(X,NORM_2,&xnorm); /* xnorm <- ||x|| */
281: VecNormEnd(X,NORM_2,&xnorm);
282: if (PetscIsInfOrNanReal(vi->merit)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
284: PetscObjectAMSTakeAccess((PetscObject)snes);
285: snes->norm = vi->phinorm;
286: PetscObjectAMSGrantAccess((PetscObject)snes);
287: SNESLogConvergenceHistory(snes,vi->phinorm,0);
288: SNESMonitor(snes,0,vi->phinorm);
290: /* set parameter for default relative tolerance convergence test */
291: snes->ttol = vi->phinorm*snes->rtol;
292: /* test convergence */
293: (*snes->ops->converged)(snes,0,0.0,0.0,vi->phinorm,&snes->reason,snes->cnvP);
294: if (snes->reason) {
295: sdm->ops->computefunction = vi->computeuserfunction;
296: return(0);
297: }
299: for (i=0; i<maxits; i++) {
301: /* Call general purpose update function */
302: if (snes->ops->update) {
303: (*snes->ops->update)(snes, snes->iter);
304: }
306: /* Solve J Y = Phi, where J is the semismooth jacobian */
308: /* Get the jacobian -- note that the function must be the original function for snes_fd and snes_fd_color to work for this*/
309: sdm->ops->computefunction = vi->computeuserfunction;
310: SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);
311: sdm->ops->computefunction = SNESVIComputeFunction;
313: /* Get the diagonal shift and row scaling vectors */
314: SNESVIComputeBsubdifferentialVectors(snes,X,F,snes->jacobian,vi->Da,vi->Db);
315: /* Compute the semismooth jacobian */
316: SNESVIComputeJacobian(snes->jacobian,snes->jacobian_pre,vi->Da,vi->Db);
317: /* Compute the merit function gradient */
318: SNESVIComputeMeritFunctionGradient(snes->jacobian,vi->phi,vi->dpsi);
319: KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,flg);
320: KSPSolve(snes->ksp,vi->phi,Y);
321: KSPGetConvergedReason(snes->ksp,&kspreason);
323: if (kspreason < 0) {
324: if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
325: PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);
326: snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
327: break;
328: }
329: }
330: KSPGetIterationNumber(snes->ksp,&lits);
331: snes->linear_its += lits;
332: PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);
333: /*
334: if (snes->ops->precheck) {
335: PetscBool changed_y = PETSC_FALSE;
336: (*snes->ops->precheck)(snes,X,Y,snes->precheck,&changed_y);
337: }
339: if (PetscLogPrintInfo) {
340: SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);
341: }
342: */
343: /* Compute a (scaled) negative update in the line search routine:
344: Y <- X - lambda*Y
345: and evaluate G = function(Y) (depends on the line search).
346: */
347: VecCopy(Y,snes->vec_sol_update);
348: ynorm = 1; gnorm = vi->phinorm;
349: SNESLineSearchApply(snes->linesearch, X, vi->phi, &gnorm, Y);
350: SNESLineSearchGetSuccess(snes->linesearch, &lssucceed);
351: SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm);
352: PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",(double)vi->phinorm,(double)gnorm,(double)ynorm,(int)lssucceed);
353: if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
354: if (snes->domainerror) {
355: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
356: sdm->ops->computefunction = vi->computeuserfunction;
357: return(0);
358: }
359: if (!lssucceed) {
360: if (++snes->numFailures >= snes->maxFailures) {
361: PetscBool ismin;
362: snes->reason = SNES_DIVERGED_LINE_SEARCH;
363: SNESVICheckLocalMin_Private(snes,snes->jacobian,vi->phi,X,gnorm,&ismin);
364: if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
365: break;
366: }
367: }
368: /* Update function and solution vectors */
369: vi->phinorm = gnorm;
370: vi->merit = 0.5*vi->phinorm*vi->phinorm;
371: /* Monitor convergence */
372: PetscObjectAMSTakeAccess((PetscObject)snes);
373: snes->iter = i+1;
374: snes->norm = vi->phinorm;
375: PetscObjectAMSGrantAccess((PetscObject)snes);
376: SNESLogConvergenceHistory(snes,snes->norm,lits);
377: SNESMonitor(snes,snes->iter,snes->norm);
378: /* Test for convergence, xnorm = || X || */
379: if (snes->ops->converged != SNESSkipConverged) { VecNorm(X,NORM_2,&xnorm); }
380: (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,vi->phinorm,&snes->reason,snes->cnvP);
381: if (snes->reason) break;
382: }
383: if (i == maxits) {
384: PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);
385: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
386: }
387: sdm->ops->computefunction = vi->computeuserfunction;
388: return(0);
389: }
391: /* -------------------------------------------------------------------------- */
392: /*
393: SNESSetUp_VINEWTONSSLS - Sets up the internal data structures for the later use
394: of the SNES nonlinear solver.
396: Input Parameter:
397: . snes - the SNES context
398: . x - the solution vector
400: Application Interface Routine: SNESSetUp()
402: Notes:
403: For basic use of the SNES solvers, the user need not explicitly call
404: SNESSetUp(), since these actions will automatically occur during
405: the call to SNESSolve().
406: */
409: PetscErrorCode SNESSetUp_VINEWTONSSLS(SNES snes)
410: {
411: PetscErrorCode ierr;
412: SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS*) snes->data;
415: SNESSetUp_VI(snes);
416: VecDuplicate(snes->vec_sol, &vi->dpsi);
417: VecDuplicate(snes->vec_sol, &vi->phi);
418: VecDuplicate(snes->vec_sol, &vi->Da);
419: VecDuplicate(snes->vec_sol, &vi->Db);
420: VecDuplicate(snes->vec_sol, &vi->z);
421: VecDuplicate(snes->vec_sol, &vi->t);
422: return(0);
423: }
424: /* -------------------------------------------------------------------------- */
427: PetscErrorCode SNESReset_VINEWTONSSLS(SNES snes)
428: {
429: SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS*) snes->data;
430: PetscErrorCode ierr;
433: SNESReset_VI(snes);
434: VecDestroy(&vi->dpsi);
435: VecDestroy(&vi->phi);
436: VecDestroy(&vi->Da);
437: VecDestroy(&vi->Db);
438: VecDestroy(&vi->z);
439: VecDestroy(&vi->t);
440: return(0);
441: }
443: /* -------------------------------------------------------------------------- */
444: /*
445: SNESSetFromOptions_VINEWTONSSLS - Sets various parameters for the SNESVI method.
447: Input Parameter:
448: . snes - the SNES context
450: Application Interface Routine: SNESSetFromOptions()
451: */
454: static PetscErrorCode SNESSetFromOptions_VINEWTONSSLS(SNES snes)
455: {
457: SNESLineSearch linesearch;
460: SNESSetFromOptions_VI(snes);
461: PetscOptionsHead("SNES semismooth method options");
462: PetscOptionsTail();
463: /* set up the default line search */
464: if (!snes->linesearch) {
465: SNESGetLineSearch(snes, &linesearch);
466: SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
467: SNESLineSearchBTSetAlpha(linesearch, 0.0);
468: }
469: return(0);
470: }
473: /* -------------------------------------------------------------------------- */
474: /*MC
475: SNESVINEWTONSSLS - Semi-smooth solver for variational inequalities based on Newton's method
477: Options Database:
478: + -snes_vi_type <ss,rs,rsaug> a semi-smooth solver, a reduced space active set method, and a reduced space active set method that does not eliminate the active constraints from the Jacobian instead augments the Jacobian with additional variables that enforce the constraints
479: - -snes_vi_monitor - prints the number of active constraints at each iteration.
481: Level: beginner
483: References:
484: - T. S. Munson, F. Facchinei, M. C. Ferris, A. Fischer, and C. Kanzow. The semismooth
485: algorithm for large scale complementarity problems. INFORMS Journal on Computing, 13 (2001).
487: .seealso: SNESVISetVariableBounds(), SNESVISetComputeVariableBounds(), SNESCreate(), SNES, SNESSetType(), SNESVINEWTONRSLS, SNESVINEWTONSSLS, SNESNEWTONTR, SNESLineSearchSet(),
488: SNESLineSearchSetPostCheck(), SNESLineSearchNo(), SNESLineSearchCubic(), SNESLineSearchQuadratic(),
489: SNESLineSearchSet(), SNESLineSearchNoNorms(), SNESLineSearchSetPreCheck(), SNESLineSearchSetParams(), SNESLineSearchGetParams()
491: M*/
494: PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONSSLS(SNES snes)
495: {
496: PetscErrorCode ierr;
497: SNES_VINEWTONSSLS *vi;
500: snes->ops->reset = SNESReset_VINEWTONSSLS;
501: snes->ops->setup = SNESSetUp_VINEWTONSSLS;
502: snes->ops->solve = SNESSolve_VINEWTONSSLS;
503: snes->ops->destroy = SNESDestroy_VI;
504: snes->ops->setfromoptions = SNESSetFromOptions_VINEWTONSSLS;
505: snes->ops->view = NULL;
507: snes->usesksp = PETSC_TRUE;
508: snes->usespc = PETSC_FALSE;
510: PetscNewLog(snes,SNES_VINEWTONSSLS,&vi);
511: snes->data = (void*)vi;
513: PetscObjectComposeFunction((PetscObject)snes,"SNESVISetVariableBounds_C",SNESVISetVariableBounds_VI);
514: PetscObjectComposeFunction((PetscObject)snes,"SNESVISetComputeVariableBounds_C",SNESVISetComputeVariableBounds_VI);
515: return(0);
516: }