Actual source code: virs.c
2: #include <../src/snes/impls/vi/rs/virsimpl.h>
3: #include <petsc/private/dmimpl.h>
4: #include <petsc/private/vecimpl.h>
6: /*
7: SNESVIGetInactiveSet - Gets the global indices for the inactive set variables (these correspond to the degrees of freedom the linear
8: system is solved on)
10: Input parameter:
11: . snes - the SNES context
13: Output parameter:
14: . inact - inactive set index set
16: */
17: PetscErrorCode SNESVIGetInactiveSet(SNES snes,IS *inact)
18: {
19: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data;
21: *inact = vi->IS_inact;
22: return 0;
23: }
25: /*
26: Provides a wrapper to a DM to allow it to be used to generated the interpolation/restriction from the DM for the smaller matrices and vectors
27: defined by the reduced space method.
29: Simple calls the regular DM interpolation and restricts it to operation on the variables not associated with active constraints.
31: */
32: typedef struct {
33: PetscInt n; /* size of vectors in the reduced DM space */
34: IS inactive;
36: PetscErrorCode (*createinterpolation)(DM,DM,Mat*,Vec*); /* DM's original routines */
37: PetscErrorCode (*coarsen)(DM, MPI_Comm, DM*);
38: PetscErrorCode (*createglobalvector)(DM,Vec*);
39: PetscErrorCode (*createinjection)(DM,DM,Mat*);
40: PetscErrorCode (*hascreateinjection)(DM,PetscBool*);
42: DM dm; /* when destroying this object we need to reset the above function into the base DM */
43: } DM_SNESVI;
45: /*
46: DMCreateGlobalVector_SNESVI - Creates global vector of the size of the reduced space
48: */
49: PetscErrorCode DMCreateGlobalVector_SNESVI(DM dm,Vec *vec)
50: {
51: PetscContainer isnes;
52: DM_SNESVI *dmsnesvi;
54: PetscObjectQuery((PetscObject)dm,"VI",(PetscObject*)&isnes);
56: PetscContainerGetPointer(isnes,(void**)&dmsnesvi);
57: VecCreateMPI(PetscObjectComm((PetscObject)dm),dmsnesvi->n,PETSC_DETERMINE,vec);
58: VecSetDM(*vec, dm);
59: return 0;
60: }
62: /*
63: DMCreateInterpolation_SNESVI - Modifieds the interpolation obtained from the DM by removing all rows and columns associated with active constraints.
65: */
66: PetscErrorCode DMCreateInterpolation_SNESVI(DM dm1,DM dm2,Mat *mat,Vec *vec)
67: {
68: PetscContainer isnes;
69: DM_SNESVI *dmsnesvi1,*dmsnesvi2;
70: Mat interp;
72: PetscObjectQuery((PetscObject)dm1,"VI",(PetscObject*)&isnes);
74: PetscContainerGetPointer(isnes,(void**)&dmsnesvi1);
75: PetscObjectQuery((PetscObject)dm2,"VI",(PetscObject*)&isnes);
77: PetscContainerGetPointer(isnes,(void**)&dmsnesvi2);
79: (*dmsnesvi1->createinterpolation)(dm1,dm2,&interp,NULL);
80: MatCreateSubMatrix(interp,dmsnesvi2->inactive,dmsnesvi1->inactive,MAT_INITIAL_MATRIX,mat);
81: MatDestroy(&interp);
82: *vec = NULL;
83: return 0;
84: }
86: static PetscErrorCode DMSetVI(DM,IS);
87: static PetscErrorCode DMDestroyVI(DM);
89: /*
90: DMCoarsen_SNESVI - Computes the regular coarsened DM then computes additional information about its inactive set
92: */
93: PetscErrorCode DMCoarsen_SNESVI(DM dm1,MPI_Comm comm,DM *dm2)
94: {
95: PetscContainer isnes;
96: DM_SNESVI *dmsnesvi1;
97: Vec finemarked,coarsemarked;
98: IS inactive;
99: Mat inject;
100: const PetscInt *index;
101: PetscInt n,k,cnt = 0,rstart,*coarseindex;
102: PetscScalar *marked;
104: PetscObjectQuery((PetscObject)dm1,"VI",(PetscObject*)&isnes);
106: PetscContainerGetPointer(isnes,(void**)&dmsnesvi1);
108: /* get the original coarsen */
109: (*dmsnesvi1->coarsen)(dm1,comm,dm2);
111: /* not sure why this extra reference is needed, but without the dm2 disappears too early */
112: /* Updating the KSPCreateVecs() to avoid using DMGetGlobalVector() when matrix is available removes the need for this reference? */
113: /* PetscObjectReference((PetscObject)*dm2);*/
115: /* need to set back global vectors in order to use the original injection */
116: DMClearGlobalVectors(dm1);
118: dm1->ops->createglobalvector = dmsnesvi1->createglobalvector;
120: DMCreateGlobalVector(dm1,&finemarked);
121: DMCreateGlobalVector(*dm2,&coarsemarked);
123: /*
124: fill finemarked with locations of inactive points
125: */
126: ISGetIndices(dmsnesvi1->inactive,&index);
127: ISGetLocalSize(dmsnesvi1->inactive,&n);
128: VecSet(finemarked,0.0);
129: for (k=0; k<n; k++) {
130: VecSetValue(finemarked,index[k],1.0,INSERT_VALUES);
131: }
132: VecAssemblyBegin(finemarked);
133: VecAssemblyEnd(finemarked);
135: DMCreateInjection(*dm2,dm1,&inject);
136: MatRestrict(inject,finemarked,coarsemarked);
137: MatDestroy(&inject);
139: /*
140: create index set list of coarse inactive points from coarsemarked
141: */
142: VecGetLocalSize(coarsemarked,&n);
143: VecGetOwnershipRange(coarsemarked,&rstart,NULL);
144: VecGetArray(coarsemarked,&marked);
145: for (k=0; k<n; k++) {
146: if (marked[k] != 0.0) cnt++;
147: }
148: PetscMalloc1(cnt,&coarseindex);
149: cnt = 0;
150: for (k=0; k<n; k++) {
151: if (marked[k] != 0.0) coarseindex[cnt++] = k + rstart;
152: }
153: VecRestoreArray(coarsemarked,&marked);
154: ISCreateGeneral(PetscObjectComm((PetscObject)coarsemarked),cnt,coarseindex,PETSC_OWN_POINTER,&inactive);
156: DMClearGlobalVectors(dm1);
158: dm1->ops->createglobalvector = DMCreateGlobalVector_SNESVI;
160: DMSetVI(*dm2,inactive);
162: VecDestroy(&finemarked);
163: VecDestroy(&coarsemarked);
164: ISDestroy(&inactive);
165: return 0;
166: }
168: PetscErrorCode DMDestroy_SNESVI(DM_SNESVI *dmsnesvi)
169: {
170: /* reset the base methods in the DM object that were changed when the DM_SNESVI was reset */
171: dmsnesvi->dm->ops->createinterpolation = dmsnesvi->createinterpolation;
172: dmsnesvi->dm->ops->coarsen = dmsnesvi->coarsen;
173: dmsnesvi->dm->ops->createglobalvector = dmsnesvi->createglobalvector;
174: dmsnesvi->dm->ops->createinjection = dmsnesvi->createinjection;
175: dmsnesvi->dm->ops->hascreateinjection = dmsnesvi->hascreateinjection;
176: /* need to clear out this vectors because some of them may not have a reference to the DM
177: but they are counted as having references to the DM in DMDestroy() */
178: DMClearGlobalVectors(dmsnesvi->dm);
180: ISDestroy(&dmsnesvi->inactive);
181: PetscFree(dmsnesvi);
182: return 0;
183: }
185: /*
186: DMSetVI - Marks a DM as associated with a VI problem. This causes the interpolation/restriction operators to
187: be restricted to only those variables NOT associated with active constraints.
189: */
190: static PetscErrorCode DMSetVI(DM dm,IS inactive)
191: {
192: PetscContainer isnes;
193: DM_SNESVI *dmsnesvi;
195: if (!dm) return 0;
197: PetscObjectReference((PetscObject)inactive);
199: PetscObjectQuery((PetscObject)dm,"VI",(PetscObject*)&isnes);
200: if (!isnes) {
201: PetscContainerCreate(PetscObjectComm((PetscObject)dm),&isnes);
202: PetscContainerSetUserDestroy(isnes,(PetscErrorCode (*)(void*))DMDestroy_SNESVI);
203: PetscNew(&dmsnesvi);
204: PetscContainerSetPointer(isnes,(void*)dmsnesvi);
205: PetscObjectCompose((PetscObject)dm,"VI",(PetscObject)isnes);
206: PetscContainerDestroy(&isnes);
208: dmsnesvi->createinterpolation = dm->ops->createinterpolation;
209: dm->ops->createinterpolation = DMCreateInterpolation_SNESVI;
210: dmsnesvi->coarsen = dm->ops->coarsen;
211: dm->ops->coarsen = DMCoarsen_SNESVI;
212: dmsnesvi->createglobalvector = dm->ops->createglobalvector;
213: dm->ops->createglobalvector = DMCreateGlobalVector_SNESVI;
214: dmsnesvi->createinjection = dm->ops->createinjection;
215: dm->ops->createinjection = NULL;
216: dmsnesvi->hascreateinjection = dm->ops->hascreateinjection;
217: dm->ops->hascreateinjection = NULL;
218: } else {
219: PetscContainerGetPointer(isnes,(void**)&dmsnesvi);
220: ISDestroy(&dmsnesvi->inactive);
221: }
222: DMClearGlobalVectors(dm);
223: ISGetLocalSize(inactive,&dmsnesvi->n);
225: dmsnesvi->inactive = inactive;
226: dmsnesvi->dm = dm;
227: return 0;
228: }
230: /*
231: DMDestroyVI - Frees the DM_SNESVI object contained in the DM
232: - also resets the function pointers in the DM for createinterpolation() etc to use the original DM
233: */
234: static PetscErrorCode DMDestroyVI(DM dm)
235: {
236: if (!dm) return 0;
237: PetscObjectCompose((PetscObject)dm,"VI",(PetscObject)NULL);
238: return 0;
239: }
241: /* --------------------------------------------------------------------------------------------------------*/
243: PetscErrorCode SNESCreateIndexSets_VINEWTONRSLS(SNES snes,Vec X,Vec F,IS *ISact,IS *ISinact)
244: {
245: SNESVIGetActiveSetIS(snes,X,F,ISact);
246: ISComplement(*ISact,X->map->rstart,X->map->rend,ISinact);
247: return 0;
248: }
250: /* Create active and inactive set vectors. The local size of this vector is set and petsc computes the global size */
251: PetscErrorCode SNESCreateSubVectors_VINEWTONRSLS(SNES snes,PetscInt n,Vec *newv)
252: {
253: Vec v;
255: VecCreate(PetscObjectComm((PetscObject)snes),&v);
256: VecSetSizes(v,n,PETSC_DECIDE);
257: VecSetType(v,VECSTANDARD);
258: *newv = v;
259: return 0;
260: }
262: /* Resets the snes PC and KSP when the active set sizes change */
263: PetscErrorCode SNESVIResetPCandKSP(SNES snes,Mat Amat,Mat Pmat)
264: {
265: KSP snesksp;
267: SNESGetKSP(snes,&snesksp);
268: KSPReset(snesksp);
269: KSPResetFromOptions(snesksp);
271: /*
272: KSP kspnew;
273: PC pcnew;
274: MatSolverType stype;
276: KSPCreate(PetscObjectComm((PetscObject)snes),&kspnew);
277: kspnew->pc_side = snesksp->pc_side;
278: kspnew->rtol = snesksp->rtol;
279: kspnew->abstol = snesksp->abstol;
280: kspnew->max_it = snesksp->max_it;
281: KSPSetType(kspnew,((PetscObject)snesksp)->type_name);
282: KSPGetPC(kspnew,&pcnew);
283: PCSetType(kspnew->pc,((PetscObject)snesksp->pc)->type_name);
284: PCSetOperators(kspnew->pc,Amat,Pmat);
285: PCFactorGetMatSolverType(snesksp->pc,&stype);
286: PCFactorSetMatSolverType(kspnew->pc,stype);
287: KSPDestroy(&snesksp);
288: snes->ksp = kspnew;
289: PetscLogObjectParent((PetscObject)snes,(PetscObject)kspnew);
290: KSPSetFromOptions(kspnew);*/
291: return 0;
292: }
294: /* Variational Inequality solver using reduce space method. No semismooth algorithm is
295: implemented in this algorithm. It basically identifies the active constraints and does
296: a linear solve on the other variables (those not associated with the active constraints). */
297: PetscErrorCode SNESSolve_VINEWTONRSLS(SNES snes)
298: {
299: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data;
300: PetscInt maxits,i,lits;
301: SNESLineSearchReason lssucceed;
302: PetscReal fnorm,gnorm,xnorm=0,ynorm;
303: Vec Y,X,F;
304: KSPConvergedReason kspreason;
305: KSP ksp;
306: PC pc;
308: /* Multigrid must use Galerkin for coarse grids with active set/reduced space methods; cannot rediscretize on coarser grids*/
309: SNESGetKSP(snes,&ksp);
310: KSPGetPC(ksp,&pc);
311: PCMGSetGalerkin(pc,PC_MG_GALERKIN_BOTH);
313: snes->numFailures = 0;
314: snes->numLinearSolveFailures = 0;
315: snes->reason = SNES_CONVERGED_ITERATING;
317: maxits = snes->max_its; /* maximum number of iterations */
318: X = snes->vec_sol; /* solution vector */
319: F = snes->vec_func; /* residual vector */
320: Y = snes->work[0]; /* work vectors */
322: SNESLineSearchSetVIFunctions(snes->linesearch, SNESVIProjectOntoBounds, SNESVIComputeInactiveSetFnorm);
323: SNESLineSearchSetVecs(snes->linesearch, X, NULL, NULL, NULL, NULL);
324: SNESLineSearchSetUp(snes->linesearch);
326: PetscObjectSAWsTakeAccess((PetscObject)snes);
327: snes->iter = 0;
328: snes->norm = 0.0;
329: PetscObjectSAWsGrantAccess((PetscObject)snes);
331: SNESVIProjectOntoBounds(snes,X);
332: SNESComputeFunction(snes,X,F);
333: SNESVIComputeInactiveSetFnorm(snes,F,X,&fnorm);
334: VecNorm(X,NORM_2,&xnorm); /* xnorm <- ||x|| */
335: SNESCheckFunctionNorm(snes,fnorm);
336: PetscObjectSAWsTakeAccess((PetscObject)snes);
337: snes->norm = fnorm;
338: PetscObjectSAWsGrantAccess((PetscObject)snes);
339: SNESLogConvergenceHistory(snes,fnorm,0);
340: SNESMonitor(snes,0,fnorm);
342: /* test convergence */
343: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
344: if (snes->reason) return 0;
346: for (i=0; i<maxits; i++) {
348: IS IS_act; /* _act -> active set _inact -> inactive set */
349: IS IS_redact; /* redundant active set */
350: VecScatter scat_act,scat_inact;
351: PetscInt nis_act,nis_inact;
352: Vec Y_act,Y_inact,F_inact;
353: Mat jac_inact_inact,prejac_inact_inact;
354: PetscBool isequal;
356: /* Call general purpose update function */
357: if (snes->ops->update) {
358: (*snes->ops->update)(snes, snes->iter);
359: }
360: SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
361: SNESCheckJacobianDomainerror(snes);
363: /* Create active and inactive index sets */
365: /*original
366: SNESVICreateIndexSets_RS(snes,X,F,&IS_act,&vi->IS_inact);
367: */
368: SNESVIGetActiveSetIS(snes,X,F,&IS_act);
370: if (vi->checkredundancy) {
371: (*vi->checkredundancy)(snes,IS_act,&IS_redact,vi->ctxP);
372: if (IS_redact) {
373: ISSort(IS_redact);
374: ISComplement(IS_redact,X->map->rstart,X->map->rend,&vi->IS_inact);
375: ISDestroy(&IS_redact);
376: } else {
377: ISComplement(IS_act,X->map->rstart,X->map->rend,&vi->IS_inact);
378: }
379: } else {
380: ISComplement(IS_act,X->map->rstart,X->map->rend,&vi->IS_inact);
381: }
383: /* Create inactive set submatrix */
384: MatCreateSubMatrix(snes->jacobian,vi->IS_inact,vi->IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);
386: if (0) { /* Dead code (temporary developer hack) */
387: IS keptrows;
388: MatFindNonzeroRows(jac_inact_inact,&keptrows);
389: if (keptrows) {
390: PetscInt cnt,*nrows,k;
391: const PetscInt *krows,*inact;
392: PetscInt rstart;
394: MatGetOwnershipRange(jac_inact_inact,&rstart,NULL);
395: MatDestroy(&jac_inact_inact);
396: ISDestroy(&IS_act);
398: ISGetLocalSize(keptrows,&cnt);
399: ISGetIndices(keptrows,&krows);
400: ISGetIndices(vi->IS_inact,&inact);
401: PetscMalloc1(cnt,&nrows);
402: for (k=0; k<cnt; k++) nrows[k] = inact[krows[k]-rstart];
403: ISRestoreIndices(keptrows,&krows);
404: ISRestoreIndices(vi->IS_inact,&inact);
405: ISDestroy(&keptrows);
406: ISDestroy(&vi->IS_inact);
408: ISCreateGeneral(PetscObjectComm((PetscObject)snes),cnt,nrows,PETSC_OWN_POINTER,&vi->IS_inact);
409: ISComplement(vi->IS_inact,F->map->rstart,F->map->rend,&IS_act);
410: MatCreateSubMatrix(snes->jacobian,vi->IS_inact,vi->IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);
411: }
412: }
413: DMSetVI(snes->dm,vi->IS_inact);
414: /* remove later */
416: /*
417: VecView(vi->xu,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xu))->comm));
418: VecView(vi->xl,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xl))->comm));
419: VecView(X,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)X)));
420: VecView(F,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)F)));
421: ISView(vi->IS_inact,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)vi->IS_inact)));
422: */
424: /* Get sizes of active and inactive sets */
425: ISGetLocalSize(IS_act,&nis_act);
426: ISGetLocalSize(vi->IS_inact,&nis_inact);
428: /* Create active and inactive set vectors */
429: SNESCreateSubVectors_VINEWTONRSLS(snes,nis_inact,&F_inact);
430: SNESCreateSubVectors_VINEWTONRSLS(snes,nis_act,&Y_act);
431: SNESCreateSubVectors_VINEWTONRSLS(snes,nis_inact,&Y_inact);
433: /* Create scatter contexts */
434: VecScatterCreate(Y,IS_act,Y_act,NULL,&scat_act);
435: VecScatterCreate(Y,vi->IS_inact,Y_inact,NULL,&scat_inact);
437: /* Do a vec scatter to active and inactive set vectors */
438: VecScatterBegin(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);
439: VecScatterEnd(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);
441: VecScatterBegin(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);
442: VecScatterEnd(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);
444: VecScatterBegin(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);
445: VecScatterEnd(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);
447: /* Active set direction = 0 */
448: VecSet(Y_act,0);
449: if (snes->jacobian != snes->jacobian_pre) {
450: MatCreateSubMatrix(snes->jacobian_pre,vi->IS_inact,vi->IS_inact,MAT_INITIAL_MATRIX,&prejac_inact_inact);
451: } else prejac_inact_inact = jac_inact_inact;
453: ISEqual(vi->IS_inact_prev,vi->IS_inact,&isequal);
454: if (!isequal) {
455: SNESVIResetPCandKSP(snes,jac_inact_inact,prejac_inact_inact);
456: PCFieldSplitRestrictIS(pc,vi->IS_inact);
457: }
459: /* ISView(vi->IS_inact,0); */
460: /* ISView(IS_act,0);*/
461: /* MatView(snes->jacobian_pre,0); */
463: KSPSetOperators(snes->ksp,jac_inact_inact,prejac_inact_inact);
464: KSPSetUp(snes->ksp);
465: {
466: PC pc;
467: PetscBool flg;
468: KSPGetPC(snes->ksp,&pc);
469: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&flg);
470: if (flg) {
471: KSP *subksps;
472: PCFieldSplitGetSubKSP(pc,NULL,&subksps);
473: KSPGetPC(subksps[0],&pc);
474: PetscFree(subksps);
475: PetscObjectTypeCompare((PetscObject)pc,PCBJACOBI,&flg);
476: if (flg) {
477: PetscInt n,N = 101*101,j,cnts[3] = {0,0,0};
478: const PetscInt *ii;
480: ISGetSize(vi->IS_inact,&n);
481: ISGetIndices(vi->IS_inact,&ii);
482: for (j=0; j<n; j++) {
483: if (ii[j] < N) cnts[0]++;
484: else if (ii[j] < 2*N) cnts[1]++;
485: else if (ii[j] < 3*N) cnts[2]++;
486: }
487: ISRestoreIndices(vi->IS_inact,&ii);
489: PCBJacobiSetTotalBlocks(pc,3,cnts);
490: }
491: }
492: }
494: KSPSolve(snes->ksp,F_inact,Y_inact);
495: VecScatterBegin(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);
496: VecScatterEnd(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);
497: VecScatterBegin(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);
498: VecScatterEnd(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);
500: VecDestroy(&F_inact);
501: VecDestroy(&Y_act);
502: VecDestroy(&Y_inact);
503: VecScatterDestroy(&scat_act);
504: VecScatterDestroy(&scat_inact);
505: ISDestroy(&IS_act);
506: if (!isequal) {
507: ISDestroy(&vi->IS_inact_prev);
508: ISDuplicate(vi->IS_inact,&vi->IS_inact_prev);
509: }
510: ISDestroy(&vi->IS_inact);
511: MatDestroy(&jac_inact_inact);
512: if (snes->jacobian != snes->jacobian_pre) {
513: MatDestroy(&prejac_inact_inact);
514: }
516: KSPGetConvergedReason(snes->ksp,&kspreason);
517: if (kspreason < 0) {
518: if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
519: PetscInfo(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);
520: snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
521: break;
522: }
523: }
525: KSPGetIterationNumber(snes->ksp,&lits);
526: snes->linear_its += lits;
527: PetscInfo(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);
528: /*
529: if (snes->ops->precheck) {
530: PetscBool changed_y = PETSC_FALSE;
531: (*snes->ops->precheck)(snes,X,Y,snes->precheck,&changed_y);
532: }
534: if (PetscLogPrintInfo) {
535: SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);
536: }
537: */
538: /* Compute a (scaled) negative update in the line search routine:
539: Y <- X - lambda*Y
540: and evaluate G = function(Y) (depends on the line search).
541: */
542: VecCopy(Y,snes->vec_sol_update);
543: ynorm = 1; gnorm = fnorm;
544: SNESLineSearchApply(snes->linesearch, X, F, &gnorm, Y);
545: SNESLineSearchGetReason(snes->linesearch, &lssucceed);
546: SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm);
547: PetscInfo(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",(double)fnorm,(double)gnorm,(double)ynorm,(int)lssucceed);
548: if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
549: if (snes->domainerror) {
550: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
551: DMDestroyVI(snes->dm);
552: return 0;
553: }
554: if (lssucceed) {
555: if (++snes->numFailures >= snes->maxFailures) {
556: PetscBool ismin;
557: snes->reason = SNES_DIVERGED_LINE_SEARCH;
558: SNESVICheckLocalMin_Private(snes,snes->jacobian,F,X,gnorm,&ismin);
559: if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
560: break;
561: }
562: }
563: DMDestroyVI(snes->dm);
564: /* Update function and solution vectors */
565: fnorm = gnorm;
566: /* Monitor convergence */
567: PetscObjectSAWsTakeAccess((PetscObject)snes);
568: snes->iter = i+1;
569: snes->norm = fnorm;
570: snes->xnorm = xnorm;
571: snes->ynorm = ynorm;
572: PetscObjectSAWsGrantAccess((PetscObject)snes);
573: SNESLogConvergenceHistory(snes,snes->norm,lits);
574: SNESMonitor(snes,snes->iter,snes->norm);
575: /* Test for convergence, xnorm = || X || */
576: if (snes->ops->converged != SNESConvergedSkip) VecNorm(X,NORM_2,&xnorm);
577: (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);
578: if (snes->reason) break;
579: }
580: /* make sure that the VI information attached to the DM is removed if the for loop above was broken early due to some exceptional conditional */
581: DMDestroyVI(snes->dm);
582: if (i == maxits) {
583: PetscInfo(snes,"Maximum number of iterations has been reached: %D\n",maxits);
584: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
585: }
586: return 0;
587: }
589: PetscErrorCode SNESVISetRedundancyCheck(SNES snes,PetscErrorCode (*func)(SNES,IS,IS*,void*),void *ctx)
590: {
591: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data;
594: vi->checkredundancy = func;
595: vi->ctxP = ctx;
596: return 0;
597: }
599: #if defined(PETSC_HAVE_MATLAB_ENGINE)
600: #include <engine.h>
601: #include <mex.h>
602: typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext;
604: PetscErrorCode SNESVIRedundancyCheck_Matlab(SNES snes,IS is_act,IS *is_redact,void *ctx)
605: {
606: SNESMatlabContext *sctx = (SNESMatlabContext*)ctx;
607: int nlhs = 1, nrhs = 5;
608: mxArray *plhs[1], *prhs[5];
609: long long int l1 = 0, l2 = 0, ls = 0;
610: PetscInt *indices=NULL;
617: /* Create IS for reduced active set of size 0, its size and indices will
618: bet set by the Matlab function */
619: ISCreateGeneral(PetscObjectComm((PetscObject)snes),0,indices,PETSC_OWN_POINTER,is_redact);
620: /* call Matlab function in ctx */
621: PetscArraycpy(&ls,&snes,1);
622: PetscArraycpy(&l1,&is_act,1);
623: PetscArraycpy(&l2,is_redact,1);
624: prhs[0] = mxCreateDoubleScalar((double)ls);
625: prhs[1] = mxCreateDoubleScalar((double)l1);
626: prhs[2] = mxCreateDoubleScalar((double)l2);
627: prhs[3] = mxCreateString(sctx->funcname);
628: prhs[4] = sctx->ctx;
629: mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESVIRedundancyCheckInternal");
630: mxGetScalar(plhs[0]);
631: mxDestroyArray(prhs[0]);
632: mxDestroyArray(prhs[1]);
633: mxDestroyArray(prhs[2]);
634: mxDestroyArray(prhs[3]);
635: mxDestroyArray(plhs[0]);
636: return 0;
637: }
639: PetscErrorCode SNESVISetRedundancyCheckMatlab(SNES snes,const char *func,mxArray *ctx)
640: {
641: SNESMatlabContext *sctx;
643: /* currently sctx is memory bleed */
644: PetscNew(&sctx);
645: PetscStrallocpy(func,&sctx->funcname);
646: sctx->ctx = mxDuplicateArray(ctx);
647: SNESVISetRedundancyCheck(snes,SNESVIRedundancyCheck_Matlab,sctx);
648: return 0;
649: }
651: #endif
653: /* -------------------------------------------------------------------------- */
654: /*
655: SNESSetUp_VINEWTONRSLS - Sets up the internal data structures for the later use
656: of the SNESVI nonlinear solver.
658: Input Parameter:
659: . snes - the SNES context
661: Application Interface Routine: SNESSetUp()
663: Notes:
664: For basic use of the SNES solvers, the user need not explicitly call
665: SNESSetUp(), since these actions will automatically occur during
666: the call to SNESSolve().
667: */
668: PetscErrorCode SNESSetUp_VINEWTONRSLS(SNES snes)
669: {
670: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*) snes->data;
671: PetscInt *indices;
672: PetscInt i,n,rstart,rend;
673: SNESLineSearch linesearch;
675: SNESSetUp_VI(snes);
677: /* Set up previous active index set for the first snes solve
678: vi->IS_inact_prev = 0,1,2,....N */
680: VecGetOwnershipRange(snes->vec_sol,&rstart,&rend);
681: VecGetLocalSize(snes->vec_sol,&n);
682: PetscMalloc1(n,&indices);
683: for (i=0; i < n; i++) indices[i] = rstart + i;
684: ISCreateGeneral(PetscObjectComm((PetscObject)snes),n,indices,PETSC_OWN_POINTER,&vi->IS_inact_prev);
686: /* set the line search functions */
687: if (!snes->linesearch) {
688: SNESGetLineSearch(snes, &linesearch);
689: SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
690: }
691: return 0;
692: }
693: /* -------------------------------------------------------------------------- */
694: PetscErrorCode SNESReset_VINEWTONRSLS(SNES snes)
695: {
696: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*) snes->data;
698: SNESReset_VI(snes);
699: ISDestroy(&vi->IS_inact_prev);
700: return 0;
701: }
703: /* -------------------------------------------------------------------------- */
704: /*MC
705: SNESVINEWTONRSLS - Reduced space active set solvers for variational inequalities based on Newton's method
707: Options Database:
708: + -snes_type <vinewtonssls,vinewtonrsls> - a semi-smooth solver, a reduced space active set method
709: - -snes_vi_monitor - prints the number of active constraints at each iteration.
711: Level: beginner
713: References:
714: . * - T. S. Munson, and S. Benson. Flexible Complementarity Solvers for Large Scale
715: Applications, Optimization Methods and Software, 21 (2006).
717: .seealso: SNESVISetVariableBounds(), SNESVISetComputeVariableBounds(), SNESCreate(), SNES, SNESSetType(), SNESVINEWTONSSLS, SNESNEWTONTR, SNESLineSearchSetType(),SNESLineSearchSetPostCheck(), SNESLineSearchSetPreCheck()
719: M*/
720: PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONRSLS(SNES snes)
721: {
722: SNES_VINEWTONRSLS *vi;
723: SNESLineSearch linesearch;
725: snes->ops->reset = SNESReset_VINEWTONRSLS;
726: snes->ops->setup = SNESSetUp_VINEWTONRSLS;
727: snes->ops->solve = SNESSolve_VINEWTONRSLS;
728: snes->ops->destroy = SNESDestroy_VI;
729: snes->ops->setfromoptions = SNESSetFromOptions_VI;
730: snes->ops->view = NULL;
731: snes->ops->converged = SNESConvergedDefault_VI;
733: snes->usesksp = PETSC_TRUE;
734: snes->usesnpc = PETSC_FALSE;
736: SNESGetLineSearch(snes, &linesearch);
737: if (!((PetscObject)linesearch)->type_name) {
738: SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
739: }
740: SNESLineSearchBTSetAlpha(linesearch, 0.0);
742: snes->alwayscomputesfinalresidual = PETSC_TRUE;
744: PetscNewLog(snes,&vi);
745: snes->data = (void*)vi;
746: vi->checkredundancy = NULL;
748: PetscObjectComposeFunction((PetscObject)snes,"SNESVISetVariableBounds_C",SNESVISetVariableBounds_VI);
749: PetscObjectComposeFunction((PetscObject)snes,"SNESVISetComputeVariableBounds_C",SNESVISetComputeVariableBounds_VI);
750: return 0;
751: }