Actual source code: virs.c
petsc-3.8.4 2018-03-24
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;
22: *inact = vi->IS_inact;
23: return(0);
24: }
26: /*
27: 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
28: defined by the reduced space method.
30: Simple calls the regular DM interpolation and restricts it to operation on the variables not associated with active constraints.
32: <*/
33: typedef struct {
34: PetscInt n; /* size of vectors in the reduced DM space */
35: IS inactive;
37: PetscErrorCode (*createinterpolation)(DM,DM,Mat*,Vec*); /* DM's original routines */
38: PetscErrorCode (*coarsen)(DM, MPI_Comm, DM*);
39: PetscErrorCode (*createglobalvector)(DM,Vec*);
41: DM dm; /* when destroying this object we need to reset the above function into the base DM */
42: } DM_SNESVI;
44: /*
45: DMCreateGlobalVector_SNESVI - Creates global vector of the size of the reduced space
47: */
48: PetscErrorCode DMCreateGlobalVector_SNESVI(DM dm,Vec *vec)
49: {
51: PetscContainer isnes;
52: DM_SNESVI *dmsnesvi;
55: PetscObjectQuery((PetscObject)dm,"VI",(PetscObject*)&isnes);
56: if (!isnes) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Composed SNES is missing");
57: PetscContainerGetPointer(isnes,(void**)&dmsnesvi);
58: VecCreateMPI(PetscObjectComm((PetscObject)dm),dmsnesvi->n,PETSC_DETERMINE,vec);
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: {
69: PetscContainer isnes;
70: DM_SNESVI *dmsnesvi1,*dmsnesvi2;
71: Mat interp;
74: PetscObjectQuery((PetscObject)dm1,"VI",(PetscObject*)&isnes);
75: if (!isnes) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_PLIB,"Composed VI data structure is missing");
76: PetscContainerGetPointer(isnes,(void**)&dmsnesvi1);
77: PetscObjectQuery((PetscObject)dm2,"VI",(PetscObject*)&isnes);
78: if (!isnes) SETERRQ(PetscObjectComm((PetscObject)dm2),PETSC_ERR_PLIB,"Composed VI data structure is missing");
79: PetscContainerGetPointer(isnes,(void**)&dmsnesvi2);
81: (*dmsnesvi1->createinterpolation)(dm1,dm2,&interp,NULL);
82: MatCreateSubMatrix(interp,dmsnesvi2->inactive,dmsnesvi1->inactive,MAT_INITIAL_MATRIX,mat);
83: MatDestroy(&interp);
84: *vec = 0;
85: return(0);
86: }
88: extern PetscErrorCode DMSetVI(DM,IS);
90: /*
91: DMCoarsen_SNESVI - Computes the regular coarsened DM then computes additional information about its inactive set
93: */
94: PetscErrorCode DMCoarsen_SNESVI(DM dm1,MPI_Comm comm,DM *dm2)
95: {
97: PetscContainer isnes;
98: DM_SNESVI *dmsnesvi1;
99: Vec finemarked,coarsemarked;
100: IS inactive;
101: Mat inject;
102: const PetscInt *index;
103: PetscInt n,k,cnt = 0,rstart,*coarseindex;
104: PetscScalar *marked;
107: PetscObjectQuery((PetscObject)dm1,"VI",(PetscObject*)&isnes);
108: if (!isnes) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_PLIB,"Composed VI data structure is missing");
109: PetscContainerGetPointer(isnes,(void**)&dmsnesvi1);
111: /* get the original coarsen */
112: (*dmsnesvi1->coarsen)(dm1,comm,dm2);
114: /* not sure why this extra reference is needed, but without the dm2 disappears too early */
115: /* Updating the KSPCreateVecs() to avoid using DMGetGlobalVector() when matrix is available removes the need for this reference? */
116: /* PetscObjectReference((PetscObject)*dm2);*/
118: /* need to set back global vectors in order to use the original injection */
119: DMClearGlobalVectors(dm1);
121: dm1->ops->createglobalvector = dmsnesvi1->createglobalvector;
123: DMCreateGlobalVector(dm1,&finemarked);
124: DMCreateGlobalVector(*dm2,&coarsemarked);
126: /*
127: fill finemarked with locations of inactive points
128: */
129: ISGetIndices(dmsnesvi1->inactive,&index);
130: ISGetLocalSize(dmsnesvi1->inactive,&n);
131: VecSet(finemarked,0.0);
132: for (k=0; k<n; k++) {
133: VecSetValue(finemarked,index[k],1.0,INSERT_VALUES);
134: }
135: VecAssemblyBegin(finemarked);
136: VecAssemblyEnd(finemarked);
138: DMCreateInjection(*dm2,dm1,&inject);
139: MatRestrict(inject,finemarked,coarsemarked);
140: MatDestroy(&inject);
142: /*
143: create index set list of coarse inactive points from coarsemarked
144: */
145: VecGetLocalSize(coarsemarked,&n);
146: VecGetOwnershipRange(coarsemarked,&rstart,NULL);
147: VecGetArray(coarsemarked,&marked);
148: for (k=0; k<n; k++) {
149: if (marked[k] != 0.0) cnt++;
150: }
151: PetscMalloc1(cnt,&coarseindex);
152: cnt = 0;
153: for (k=0; k<n; k++) {
154: if (marked[k] != 0.0) coarseindex[cnt++] = k + rstart;
155: }
156: VecRestoreArray(coarsemarked,&marked);
157: ISCreateGeneral(PetscObjectComm((PetscObject)coarsemarked),cnt,coarseindex,PETSC_OWN_POINTER,&inactive);
159: DMClearGlobalVectors(dm1);
161: dm1->ops->createglobalvector = DMCreateGlobalVector_SNESVI;
163: DMSetVI(*dm2,inactive);
165: VecDestroy(&finemarked);
166: VecDestroy(&coarsemarked);
167: ISDestroy(&inactive);
168: return(0);
169: }
171: PetscErrorCode DMDestroy_SNESVI(DM_SNESVI *dmsnesvi)
172: {
176: /* reset the base methods in the DM object that were changed when the DM_SNESVI was reset */
177: dmsnesvi->dm->ops->createinterpolation = dmsnesvi->createinterpolation;
178: dmsnesvi->dm->ops->coarsen = dmsnesvi->coarsen;
179: dmsnesvi->dm->ops->createglobalvector = dmsnesvi->createglobalvector;
180: /* need to clear out this vectors because some of them may not have a reference to the DM
181: but they are counted as having references to the DM in DMDestroy() */
182: DMClearGlobalVectors(dmsnesvi->dm);
184: ISDestroy(&dmsnesvi->inactive);
185: PetscFree(dmsnesvi);
186: return(0);
187: }
189: /*
190: DMSetVI - Marks a DM as associated with a VI problem. This causes the interpolation/restriction operators to
191: be restricted to only those variables NOT associated with active constraints.
193: */
194: PetscErrorCode DMSetVI(DM dm,IS inactive)
195: {
197: PetscContainer isnes;
198: DM_SNESVI *dmsnesvi;
201: if (!dm) return(0);
203: PetscObjectReference((PetscObject)inactive);
205: PetscObjectQuery((PetscObject)dm,"VI",(PetscObject*)&isnes);
206: if (!isnes) {
207: PetscContainerCreate(PetscObjectComm((PetscObject)dm),&isnes);
208: PetscContainerSetUserDestroy(isnes,(PetscErrorCode (*)(void*))DMDestroy_SNESVI);
209: PetscNew(&dmsnesvi);
210: PetscContainerSetPointer(isnes,(void*)dmsnesvi);
211: PetscObjectCompose((PetscObject)dm,"VI",(PetscObject)isnes);
212: PetscContainerDestroy(&isnes);
214: dmsnesvi->createinterpolation = dm->ops->createinterpolation;
215: dm->ops->createinterpolation = DMCreateInterpolation_SNESVI;
216: dmsnesvi->coarsen = dm->ops->coarsen;
217: dm->ops->coarsen = DMCoarsen_SNESVI;
218: dmsnesvi->createglobalvector = dm->ops->createglobalvector;
219: dm->ops->createglobalvector = DMCreateGlobalVector_SNESVI;
220: } else {
221: PetscContainerGetPointer(isnes,(void**)&dmsnesvi);
222: ISDestroy(&dmsnesvi->inactive);
223: }
224: DMClearGlobalVectors(dm);
225: ISGetLocalSize(inactive,&dmsnesvi->n);
227: dmsnesvi->inactive = inactive;
228: dmsnesvi->dm = dm;
229: return(0);
230: }
232: /*
233: DMDestroyVI - Frees the DM_SNESVI object contained in the DM
234: - also resets the function pointers in the DM for createinterpolation() etc to use the original DM
235: */
236: PetscErrorCode DMDestroyVI(DM dm)
237: {
241: if (!dm) return(0);
242: PetscObjectCompose((PetscObject)dm,"VI",(PetscObject)NULL);
243: return(0);
244: }
246: /* --------------------------------------------------------------------------------------------------------*/
249: PetscErrorCode SNESCreateIndexSets_VINEWTONRSLS(SNES snes,Vec X,Vec F,IS *ISact,IS *ISinact)
250: {
254: SNESVIGetActiveSetIS(snes,X,F,ISact);
255: ISComplement(*ISact,X->map->rstart,X->map->rend,ISinact);
256: return(0);
257: }
259: /* Create active and inactive set vectors. The local size of this vector is set and petsc computes the global size */
260: PetscErrorCode SNESCreateSubVectors_VINEWTONRSLS(SNES snes,PetscInt n,Vec *newv)
261: {
263: Vec v;
266: VecCreate(PetscObjectComm((PetscObject)snes),&v);
267: VecSetSizes(v,n,PETSC_DECIDE);
268: VecSetType(v,VECSTANDARD);
269: *newv = v;
270: return(0);
271: }
273: /* Resets the snes PC and KSP when the active set sizes change */
274: PetscErrorCode SNESVIResetPCandKSP(SNES snes,Mat Amat,Mat Pmat)
275: {
277: KSP snesksp;
280: SNESGetKSP(snes,&snesksp);
281: KSPReset(snesksp);
283: /*
284: KSP kspnew;
285: PC pcnew;
286: const MatSolverPackage stype;
289: KSPCreate(PetscObjectComm((PetscObject)snes),&kspnew);
290: kspnew->pc_side = snesksp->pc_side;
291: kspnew->rtol = snesksp->rtol;
292: kspnew->abstol = snesksp->abstol;
293: kspnew->max_it = snesksp->max_it;
294: KSPSetType(kspnew,((PetscObject)snesksp)->type_name);
295: KSPGetPC(kspnew,&pcnew);
296: PCSetType(kspnew->pc,((PetscObject)snesksp->pc)->type_name);
297: PCSetOperators(kspnew->pc,Amat,Pmat);
298: PCFactorGetMatSolverPackage(snesksp->pc,&stype);
299: PCFactorSetMatSolverPackage(kspnew->pc,stype);
300: KSPDestroy(&snesksp);
301: snes->ksp = kspnew;
302: PetscLogObjectParent((PetscObject)snes,(PetscObject)kspnew);
303: KSPSetFromOptions(kspnew);*/
304: return(0);
305: }
307: /* Variational Inequality solver using reduce space method. No semismooth algorithm is
308: implemented in this algorithm. It basically identifies the active constraints and does
309: a linear solve on the other variables (those not associated with the active constraints). */
310: PetscErrorCode SNESSolve_VINEWTONRSLS(SNES snes)
311: {
312: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data;
313: PetscErrorCode ierr;
314: PetscInt maxits,i,lits;
315: SNESLineSearchReason lssucceed;
316: PetscReal fnorm,gnorm,xnorm=0,ynorm;
317: Vec Y,X,F;
318: KSPConvergedReason kspreason;
319: KSP ksp;
320: PC pc;
323: /* Multigrid must use Galerkin for coarse grids with active set/reduced space methods; cannot rediscretize on coarser grids*/
324: SNESGetKSP(snes,&ksp);
325: KSPGetPC(ksp,&pc);
326: PCMGSetGalerkin(pc,PC_MG_GALERKIN_BOTH);
328: snes->numFailures = 0;
329: snes->numLinearSolveFailures = 0;
330: snes->reason = SNES_CONVERGED_ITERATING;
332: maxits = snes->max_its; /* maximum number of iterations */
333: X = snes->vec_sol; /* solution vector */
334: F = snes->vec_func; /* residual vector */
335: Y = snes->work[0]; /* work vectors */
337: SNESLineSearchSetVIFunctions(snes->linesearch, SNESVIProjectOntoBounds, SNESVIComputeInactiveSetFnorm);
338: SNESLineSearchSetVecs(snes->linesearch, X, NULL, NULL, NULL, NULL);
339: SNESLineSearchSetUp(snes->linesearch);
341: PetscObjectSAWsTakeAccess((PetscObject)snes);
342: snes->iter = 0;
343: snes->norm = 0.0;
344: PetscObjectSAWsGrantAccess((PetscObject)snes);
346: SNESVIProjectOntoBounds(snes,X);
347: SNESComputeFunction(snes,X,F);
348: SNESVIComputeInactiveSetFnorm(snes,F,X,&fnorm);
349: VecNormBegin(X,NORM_2,&xnorm); /* xnorm <- ||x|| */
350: VecNormEnd(X,NORM_2,&xnorm);
351: SNESCheckFunctionNorm(snes,fnorm);
352: PetscObjectSAWsTakeAccess((PetscObject)snes);
353: snes->norm = fnorm;
354: PetscObjectSAWsGrantAccess((PetscObject)snes);
355: SNESLogConvergenceHistory(snes,fnorm,0);
356: SNESMonitor(snes,0,fnorm);
358: /* test convergence */
359: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
360: if (snes->reason) return(0);
363: for (i=0; i<maxits; i++) {
365: IS IS_act; /* _act -> active set _inact -> inactive set */
366: IS IS_redact; /* redundant active set */
367: VecScatter scat_act,scat_inact;
368: PetscInt nis_act,nis_inact;
369: Vec Y_act,Y_inact,F_inact;
370: Mat jac_inact_inact,prejac_inact_inact;
371: PetscBool isequal;
373: /* Call general purpose update function */
374: if (snes->ops->update) {
375: (*snes->ops->update)(snes, snes->iter);
376: }
377: SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
380: /* Create active and inactive index sets */
382: /*original
383: SNESVICreateIndexSets_RS(snes,X,F,&IS_act,&vi->IS_inact);
384: */
385: SNESVIGetActiveSetIS(snes,X,F,&IS_act);
387: if (vi->checkredundancy) {
388: (*vi->checkredundancy)(snes,IS_act,&IS_redact,vi->ctxP);
389: if (IS_redact) {
390: ISSort(IS_redact);
391: ISComplement(IS_redact,X->map->rstart,X->map->rend,&vi->IS_inact);
392: ISDestroy(&IS_redact);
393: } else {
394: ISComplement(IS_act,X->map->rstart,X->map->rend,&vi->IS_inact);
395: }
396: } else {
397: ISComplement(IS_act,X->map->rstart,X->map->rend,&vi->IS_inact);
398: }
401: /* Create inactive set submatrix */
402: MatCreateSubMatrix(snes->jacobian,vi->IS_inact,vi->IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);
404: if (0) { /* Dead code (temporary developer hack) */
405: IS keptrows;
406: MatFindNonzeroRows(jac_inact_inact,&keptrows);
407: if (keptrows) {
408: PetscInt cnt,*nrows,k;
409: const PetscInt *krows,*inact;
410: PetscInt rstart;
412: MatGetOwnershipRange(jac_inact_inact,&rstart,NULL);
413: MatDestroy(&jac_inact_inact);
414: ISDestroy(&IS_act);
416: ISGetLocalSize(keptrows,&cnt);
417: ISGetIndices(keptrows,&krows);
418: ISGetIndices(vi->IS_inact,&inact);
419: PetscMalloc1(cnt,&nrows);
420: for (k=0; k<cnt; k++) nrows[k] = inact[krows[k]-rstart];
421: ISRestoreIndices(keptrows,&krows);
422: ISRestoreIndices(vi->IS_inact,&inact);
423: ISDestroy(&keptrows);
424: ISDestroy(&vi->IS_inact);
426: ISCreateGeneral(PetscObjectComm((PetscObject)snes),cnt,nrows,PETSC_OWN_POINTER,&vi->IS_inact);
427: ISComplement(vi->IS_inact,F->map->rstart,F->map->rend,&IS_act);
428: MatCreateSubMatrix(snes->jacobian,vi->IS_inact,vi->IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);
429: }
430: }
431: DMSetVI(snes->dm,vi->IS_inact);
432: /* remove later */
434: /*
435: VecView(vi->xu,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xu))->comm));
436: VecView(vi->xl,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xl))->comm));
437: VecView(X,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)X)));
438: VecView(F,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)F)));
439: ISView(vi->IS_inact,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)vi->IS_inact)));
440: */
442: /* Get sizes of active and inactive sets */
443: ISGetLocalSize(IS_act,&nis_act);
444: ISGetLocalSize(vi->IS_inact,&nis_inact);
446: /* Create active and inactive set vectors */
447: SNESCreateSubVectors_VINEWTONRSLS(snes,nis_inact,&F_inact);
448: SNESCreateSubVectors_VINEWTONRSLS(snes,nis_act,&Y_act);
449: SNESCreateSubVectors_VINEWTONRSLS(snes,nis_inact,&Y_inact);
451: /* Create scatter contexts */
452: VecScatterCreate(Y,IS_act,Y_act,NULL,&scat_act);
453: VecScatterCreate(Y,vi->IS_inact,Y_inact,NULL,&scat_inact);
455: /* Do a vec scatter to active and inactive set vectors */
456: VecScatterBegin(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);
457: VecScatterEnd(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);
459: VecScatterBegin(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);
460: VecScatterEnd(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);
462: VecScatterBegin(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);
463: VecScatterEnd(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);
465: /* Active set direction = 0 */
466: VecSet(Y_act,0);
467: if (snes->jacobian != snes->jacobian_pre) {
468: MatCreateSubMatrix(snes->jacobian_pre,vi->IS_inact,vi->IS_inact,MAT_INITIAL_MATRIX,&prejac_inact_inact);
469: } else prejac_inact_inact = jac_inact_inact;
471: ISEqual(vi->IS_inact_prev,vi->IS_inact,&isequal);
472: if (!isequal) {
473: SNESVIResetPCandKSP(snes,jac_inact_inact,prejac_inact_inact);
474: PCFieldSplitRestrictIS(pc,vi->IS_inact);
475: }
477: /* ISView(vi->IS_inact,0); */
478: /* ISView(IS_act,0);*/
479: /* MatView(snes->jacobian_pre,0); */
483: KSPSetOperators(snes->ksp,jac_inact_inact,prejac_inact_inact);
484: KSPSetUp(snes->ksp);
485: {
486: PC pc;
487: PetscBool flg;
488: KSPGetPC(snes->ksp,&pc);
489: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&flg);
490: if (flg) {
491: KSP *subksps;
492: PCFieldSplitGetSubKSP(pc,NULL,&subksps);
493: KSPGetPC(subksps[0],&pc);
494: PetscFree(subksps);
495: PetscObjectTypeCompare((PetscObject)pc,PCBJACOBI,&flg);
496: if (flg) {
497: PetscInt n,N = 101*101,j,cnts[3] = {0,0,0};
498: const PetscInt *ii;
500: ISGetSize(vi->IS_inact,&n);
501: ISGetIndices(vi->IS_inact,&ii);
502: for (j=0; j<n; j++) {
503: if (ii[j] < N) cnts[0]++;
504: else if (ii[j] < 2*N) cnts[1]++;
505: else if (ii[j] < 3*N) cnts[2]++;
506: }
507: ISRestoreIndices(vi->IS_inact,&ii);
509: PCBJacobiSetTotalBlocks(pc,3,cnts);
510: }
511: }
512: }
514: KSPSolve(snes->ksp,F_inact,Y_inact);
515: VecScatterBegin(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);
516: VecScatterEnd(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);
517: VecScatterBegin(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);
518: VecScatterEnd(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);
520: VecDestroy(&F_inact);
521: VecDestroy(&Y_act);
522: VecDestroy(&Y_inact);
523: VecScatterDestroy(&scat_act);
524: VecScatterDestroy(&scat_inact);
525: ISDestroy(&IS_act);
526: if (!isequal) {
527: ISDestroy(&vi->IS_inact_prev);
528: ISDuplicate(vi->IS_inact,&vi->IS_inact_prev);
529: }
530: ISDestroy(&vi->IS_inact);
531: MatDestroy(&jac_inact_inact);
532: if (snes->jacobian != snes->jacobian_pre) {
533: MatDestroy(&prejac_inact_inact);
534: }
536: KSPGetConvergedReason(snes->ksp,&kspreason);
537: if (kspreason < 0) {
538: if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
539: PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);
540: snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
541: break;
542: }
543: }
545: KSPGetIterationNumber(snes->ksp,&lits);
546: snes->linear_its += lits;
547: PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);
548: /*
549: if (snes->ops->precheck) {
550: PetscBool changed_y = PETSC_FALSE;
551: (*snes->ops->precheck)(snes,X,Y,snes->precheck,&changed_y);
552: }
554: if (PetscLogPrintInfo) {
555: SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);
556: }
557: */
558: /* Compute a (scaled) negative update in the line search routine:
559: Y <- X - lambda*Y
560: and evaluate G = function(Y) (depends on the line search).
561: */
562: VecCopy(Y,snes->vec_sol_update);
563: ynorm = 1; gnorm = fnorm;
564: SNESLineSearchApply(snes->linesearch, X, F, &gnorm, Y);
565: SNESLineSearchGetReason(snes->linesearch, &lssucceed);
566: SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm);
567: PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",(double)fnorm,(double)gnorm,(double)ynorm,(int)lssucceed);
568: if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
569: if (snes->domainerror) {
570: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
571: DMDestroyVI(snes->dm);
572: return(0);
573: }
574: if (lssucceed) {
575: if (++snes->numFailures >= snes->maxFailures) {
576: PetscBool ismin;
577: snes->reason = SNES_DIVERGED_LINE_SEARCH;
578: SNESVICheckLocalMin_Private(snes,snes->jacobian,F,X,gnorm,&ismin);
579: if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
580: break;
581: }
582: }
583: DMDestroyVI(snes->dm);
584: /* Update function and solution vectors */
585: fnorm = gnorm;
586: /* Monitor convergence */
587: PetscObjectSAWsTakeAccess((PetscObject)snes);
588: snes->iter = i+1;
589: snes->norm = fnorm;
590: PetscObjectSAWsGrantAccess((PetscObject)snes);
591: SNESLogConvergenceHistory(snes,snes->norm,lits);
592: SNESMonitor(snes,snes->iter,snes->norm);
593: /* Test for convergence, xnorm = || X || */
594: if (snes->ops->converged != SNESConvergedSkip) { VecNorm(X,NORM_2,&xnorm); }
595: (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);
596: if (snes->reason) break;
597: }
598: /* 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 */
599: DMDestroyVI(snes->dm);
600: if (i == maxits) {
601: PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);
602: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
603: }
604: return(0);
605: }
607: PetscErrorCode SNESVISetRedundancyCheck(SNES snes,PetscErrorCode (*func)(SNES,IS,IS*,void*),void *ctx)
608: {
609: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data;
613: vi->checkredundancy = func;
614: vi->ctxP = ctx;
615: return(0);
616: }
618: #if defined(PETSC_HAVE_MATLAB_ENGINE)
619: #include <engine.h>
620: #include <mex.h>
621: typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext;
623: PetscErrorCode SNESVIRedundancyCheck_Matlab(SNES snes,IS is_act,IS *is_redact,void *ctx)
624: {
625: PetscErrorCode ierr;
626: SNESMatlabContext *sctx = (SNESMatlabContext*)ctx;
627: int nlhs = 1, nrhs = 5;
628: mxArray *plhs[1], *prhs[5];
629: long long int l1 = 0, l2 = 0, ls = 0;
630: PetscInt *indices=NULL;
638: /* Create IS for reduced active set of size 0, its size and indices will
639: bet set by the Matlab function */
640: ISCreateGeneral(PetscObjectComm((PetscObject)snes),0,indices,PETSC_OWN_POINTER,is_redact);
641: /* call Matlab function in ctx */
642: PetscMemcpy(&ls,&snes,sizeof(snes));
643: PetscMemcpy(&l1,&is_act,sizeof(is_act));
644: PetscMemcpy(&l2,is_redact,sizeof(is_act));
645: prhs[0] = mxCreateDoubleScalar((double)ls);
646: prhs[1] = mxCreateDoubleScalar((double)l1);
647: prhs[2] = mxCreateDoubleScalar((double)l2);
648: prhs[3] = mxCreateString(sctx->funcname);
649: prhs[4] = sctx->ctx;
650: mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESVIRedundancyCheckInternal");
651: mxGetScalar(plhs[0]);
652: mxDestroyArray(prhs[0]);
653: mxDestroyArray(prhs[1]);
654: mxDestroyArray(prhs[2]);
655: mxDestroyArray(prhs[3]);
656: mxDestroyArray(plhs[0]);
657: return(0);
658: }
660: PetscErrorCode SNESVISetRedundancyCheckMatlab(SNES snes,const char *func,mxArray *ctx)
661: {
662: PetscErrorCode ierr;
663: SNESMatlabContext *sctx;
666: /* currently sctx is memory bleed */
667: PetscNew(&sctx);
668: PetscStrallocpy(func,&sctx->funcname);
669: sctx->ctx = mxDuplicateArray(ctx);
670: SNESVISetRedundancyCheck(snes,SNESVIRedundancyCheck_Matlab,sctx);
671: return(0);
672: }
674: #endif
676: /* -------------------------------------------------------------------------- */
677: /*
678: SNESSetUp_VINEWTONRSLS - Sets up the internal data structures for the later use
679: of the SNESVI nonlinear solver.
681: Input Parameter:
682: . snes - the SNES context
684: Application Interface Routine: SNESSetUp()
686: Notes:
687: For basic use of the SNES solvers, the user need not explicitly call
688: SNESSetUp(), since these actions will automatically occur during
689: the call to SNESSolve().
690: */
691: PetscErrorCode SNESSetUp_VINEWTONRSLS(SNES snes)
692: {
693: PetscErrorCode ierr;
694: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*) snes->data;
695: PetscInt *indices;
696: PetscInt i,n,rstart,rend;
697: SNESLineSearch linesearch;
700: SNESSetUp_VI(snes);
702: /* Set up previous active index set for the first snes solve
703: vi->IS_inact_prev = 0,1,2,....N */
705: VecGetOwnershipRange(snes->vec_sol,&rstart,&rend);
706: VecGetLocalSize(snes->vec_sol,&n);
707: PetscMalloc1(n,&indices);
708: for (i=0; i < n; i++) indices[i] = rstart + i;
709: ISCreateGeneral(PetscObjectComm((PetscObject)snes),n,indices,PETSC_OWN_POINTER,&vi->IS_inact_prev);
711: /* set the line search functions */
712: if (!snes->linesearch) {
713: SNESGetLineSearch(snes, &linesearch);
714: SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
715: }
716: return(0);
717: }
718: /* -------------------------------------------------------------------------- */
719: PetscErrorCode SNESReset_VINEWTONRSLS(SNES snes)
720: {
721: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*) snes->data;
722: PetscErrorCode ierr;
725: SNESReset_VI(snes);
726: ISDestroy(&vi->IS_inact_prev);
727: return(0);
728: }
730: /* -------------------------------------------------------------------------- */
731: /*MC
732: SNESVINEWTONRSLS - Reduced space active set solvers for variational inequalities based on Newton's method
734: Options Database:
735: + -snes_type <vinewtonssls,vinewtonrsls> - a semi-smooth solver, a reduced space active set method
736: - -snes_vi_monitor - prints the number of active constraints at each iteration.
738: Level: beginner
740: References:
741: . 1. - T. S. Munson, and S. Benson. Flexible Complementarity Solvers for Large Scale
742: Applications, Optimization Methods and Software, 21 (2006).
744: .seealso: SNESVISetVariableBounds(), SNESVISetComputeVariableBounds(), SNESCreate(), SNES, SNESSetType(), SNESVINEWTONSSLS, SNESNEWTONTR, SNESLineSearchSetType(),SNESLineSearchSetPostCheck(), SNESLineSearchSetPreCheck()
746: M*/
747: PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONRSLS(SNES snes)
748: {
749: PetscErrorCode ierr;
750: SNES_VINEWTONRSLS *vi;
753: snes->ops->reset = SNESReset_VINEWTONRSLS;
754: snes->ops->setup = SNESSetUp_VINEWTONRSLS;
755: snes->ops->solve = SNESSolve_VINEWTONRSLS;
756: snes->ops->destroy = SNESDestroy_VI;
757: snes->ops->setfromoptions = SNESSetFromOptions_VI;
758: snes->ops->view = NULL;
759: snes->ops->converged = SNESConvergedDefault_VI;
761: snes->usesksp = PETSC_TRUE;
762: snes->usesnpc = PETSC_FALSE;
764: snes->alwayscomputesfinalresidual = PETSC_TRUE;
766: PetscNewLog(snes,&vi);
767: snes->data = (void*)vi;
768: vi->checkredundancy = NULL;
770: PetscObjectComposeFunction((PetscObject)snes,"SNESVISetVariableBounds_C",SNESVISetVariableBounds_VI);
771: PetscObjectComposeFunction((PetscObject)snes,"SNESVISetComputeVariableBounds_C",SNESVISetComputeVariableBounds_VI);
772: return(0);
773: }