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;
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*);
40: PetscErrorCode (*createinjection)(DM,DM,Mat*);
41: PetscErrorCode (*hascreateinjection)(DM,PetscBool*);
43: DM dm; /* when destroying this object we need to reset the above function into the base DM */
44: } DM_SNESVI;
46: /*
47: DMCreateGlobalVector_SNESVI - Creates global vector of the size of the reduced space
49: */
50: PetscErrorCode DMCreateGlobalVector_SNESVI(DM dm,Vec *vec)
51: {
53: PetscContainer isnes;
54: DM_SNESVI *dmsnesvi;
57: PetscObjectQuery((PetscObject)dm,"VI",(PetscObject*)&isnes);
58: if (!isnes) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Composed SNES is missing");
59: PetscContainerGetPointer(isnes,(void**)&dmsnesvi);
60: VecCreateMPI(PetscObjectComm((PetscObject)dm),dmsnesvi->n,PETSC_DETERMINE,vec);
61: VecSetDM(*vec, dm);
62: return(0);
63: }
65: /*
66: DMCreateInterpolation_SNESVI - Modifieds the interpolation obtained from the DM by removing all rows and columns associated with active constraints.
68: */
69: PetscErrorCode DMCreateInterpolation_SNESVI(DM dm1,DM dm2,Mat *mat,Vec *vec)
70: {
72: PetscContainer isnes;
73: DM_SNESVI *dmsnesvi1,*dmsnesvi2;
74: Mat interp;
77: PetscObjectQuery((PetscObject)dm1,"VI",(PetscObject*)&isnes);
78: if (!isnes) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_PLIB,"Composed VI data structure is missing");
79: PetscContainerGetPointer(isnes,(void**)&dmsnesvi1);
80: PetscObjectQuery((PetscObject)dm2,"VI",(PetscObject*)&isnes);
81: if (!isnes) SETERRQ(PetscObjectComm((PetscObject)dm2),PETSC_ERR_PLIB,"Composed VI data structure is missing");
82: PetscContainerGetPointer(isnes,(void**)&dmsnesvi2);
84: (*dmsnesvi1->createinterpolation)(dm1,dm2,&interp,NULL);
85: MatCreateSubMatrix(interp,dmsnesvi2->inactive,dmsnesvi1->inactive,MAT_INITIAL_MATRIX,mat);
86: MatDestroy(&interp);
87: *vec = NULL;
88: return(0);
89: }
91: static PetscErrorCode DMSetVI(DM,IS);
92: static PetscErrorCode DMDestroyVI(DM);
94: /*
95: DMCoarsen_SNESVI - Computes the regular coarsened DM then computes additional information about its inactive set
97: */
98: PetscErrorCode DMCoarsen_SNESVI(DM dm1,MPI_Comm comm,DM *dm2)
99: {
101: PetscContainer isnes;
102: DM_SNESVI *dmsnesvi1;
103: Vec finemarked,coarsemarked;
104: IS inactive;
105: Mat inject;
106: const PetscInt *index;
107: PetscInt n,k,cnt = 0,rstart,*coarseindex;
108: PetscScalar *marked;
111: PetscObjectQuery((PetscObject)dm1,"VI",(PetscObject*)&isnes);
112: if (!isnes) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_PLIB,"Composed VI data structure is missing");
113: PetscContainerGetPointer(isnes,(void**)&dmsnesvi1);
115: /* get the original coarsen */
116: (*dmsnesvi1->coarsen)(dm1,comm,dm2);
118: /* not sure why this extra reference is needed, but without the dm2 disappears too early */
119: /* Updating the KSPCreateVecs() to avoid using DMGetGlobalVector() when matrix is available removes the need for this reference? */
120: /* PetscObjectReference((PetscObject)*dm2);*/
122: /* need to set back global vectors in order to use the original injection */
123: DMClearGlobalVectors(dm1);
125: dm1->ops->createglobalvector = dmsnesvi1->createglobalvector;
127: DMCreateGlobalVector(dm1,&finemarked);
128: DMCreateGlobalVector(*dm2,&coarsemarked);
130: /*
131: fill finemarked with locations of inactive points
132: */
133: ISGetIndices(dmsnesvi1->inactive,&index);
134: ISGetLocalSize(dmsnesvi1->inactive,&n);
135: VecSet(finemarked,0.0);
136: for (k=0; k<n; k++) {
137: VecSetValue(finemarked,index[k],1.0,INSERT_VALUES);
138: }
139: VecAssemblyBegin(finemarked);
140: VecAssemblyEnd(finemarked);
142: DMCreateInjection(*dm2,dm1,&inject);
143: MatRestrict(inject,finemarked,coarsemarked);
144: MatDestroy(&inject);
146: /*
147: create index set list of coarse inactive points from coarsemarked
148: */
149: VecGetLocalSize(coarsemarked,&n);
150: VecGetOwnershipRange(coarsemarked,&rstart,NULL);
151: VecGetArray(coarsemarked,&marked);
152: for (k=0; k<n; k++) {
153: if (marked[k] != 0.0) cnt++;
154: }
155: PetscMalloc1(cnt,&coarseindex);
156: cnt = 0;
157: for (k=0; k<n; k++) {
158: if (marked[k] != 0.0) coarseindex[cnt++] = k + rstart;
159: }
160: VecRestoreArray(coarsemarked,&marked);
161: ISCreateGeneral(PetscObjectComm((PetscObject)coarsemarked),cnt,coarseindex,PETSC_OWN_POINTER,&inactive);
163: DMClearGlobalVectors(dm1);
165: dm1->ops->createglobalvector = DMCreateGlobalVector_SNESVI;
167: DMSetVI(*dm2,inactive);
169: VecDestroy(&finemarked);
170: VecDestroy(&coarsemarked);
171: ISDestroy(&inactive);
172: return(0);
173: }
175: PetscErrorCode DMDestroy_SNESVI(DM_SNESVI *dmsnesvi)
176: {
180: /* reset the base methods in the DM object that were changed when the DM_SNESVI was reset */
181: dmsnesvi->dm->ops->createinterpolation = dmsnesvi->createinterpolation;
182: dmsnesvi->dm->ops->coarsen = dmsnesvi->coarsen;
183: dmsnesvi->dm->ops->createglobalvector = dmsnesvi->createglobalvector;
184: dmsnesvi->dm->ops->createinjection = dmsnesvi->createinjection;
185: dmsnesvi->dm->ops->hascreateinjection = dmsnesvi->hascreateinjection;
186: /* need to clear out this vectors because some of them may not have a reference to the DM
187: but they are counted as having references to the DM in DMDestroy() */
188: DMClearGlobalVectors(dmsnesvi->dm);
190: ISDestroy(&dmsnesvi->inactive);
191: PetscFree(dmsnesvi);
192: return(0);
193: }
195: /*
196: DMSetVI - Marks a DM as associated with a VI problem. This causes the interpolation/restriction operators to
197: be restricted to only those variables NOT associated with active constraints.
199: */
200: static PetscErrorCode DMSetVI(DM dm,IS inactive)
201: {
203: PetscContainer isnes;
204: DM_SNESVI *dmsnesvi;
207: if (!dm) return(0);
209: PetscObjectReference((PetscObject)inactive);
211: PetscObjectQuery((PetscObject)dm,"VI",(PetscObject*)&isnes);
212: if (!isnes) {
213: PetscContainerCreate(PetscObjectComm((PetscObject)dm),&isnes);
214: PetscContainerSetUserDestroy(isnes,(PetscErrorCode (*)(void*))DMDestroy_SNESVI);
215: PetscNew(&dmsnesvi);
216: PetscContainerSetPointer(isnes,(void*)dmsnesvi);
217: PetscObjectCompose((PetscObject)dm,"VI",(PetscObject)isnes);
218: PetscContainerDestroy(&isnes);
220: dmsnesvi->createinterpolation = dm->ops->createinterpolation;
221: dm->ops->createinterpolation = DMCreateInterpolation_SNESVI;
222: dmsnesvi->coarsen = dm->ops->coarsen;
223: dm->ops->coarsen = DMCoarsen_SNESVI;
224: dmsnesvi->createglobalvector = dm->ops->createglobalvector;
225: dm->ops->createglobalvector = DMCreateGlobalVector_SNESVI;
226: dmsnesvi->createinjection = dm->ops->createinjection;
227: dm->ops->createinjection = NULL;
228: dmsnesvi->hascreateinjection = dm->ops->hascreateinjection;
229: dm->ops->hascreateinjection = NULL;
230: } else {
231: PetscContainerGetPointer(isnes,(void**)&dmsnesvi);
232: ISDestroy(&dmsnesvi->inactive);
233: }
234: DMClearGlobalVectors(dm);
235: ISGetLocalSize(inactive,&dmsnesvi->n);
237: dmsnesvi->inactive = inactive;
238: dmsnesvi->dm = dm;
239: return(0);
240: }
242: /*
243: DMDestroyVI - Frees the DM_SNESVI object contained in the DM
244: - also resets the function pointers in the DM for createinterpolation() etc to use the original DM
245: */
246: static PetscErrorCode DMDestroyVI(DM dm)
247: {
251: if (!dm) return(0);
252: PetscObjectCompose((PetscObject)dm,"VI",(PetscObject)NULL);
253: return(0);
254: }
256: /* --------------------------------------------------------------------------------------------------------*/
258: PetscErrorCode SNESCreateIndexSets_VINEWTONRSLS(SNES snes,Vec X,Vec F,IS *ISact,IS *ISinact)
259: {
263: SNESVIGetActiveSetIS(snes,X,F,ISact);
264: ISComplement(*ISact,X->map->rstart,X->map->rend,ISinact);
265: return(0);
266: }
268: /* Create active and inactive set vectors. The local size of this vector is set and petsc computes the global size */
269: PetscErrorCode SNESCreateSubVectors_VINEWTONRSLS(SNES snes,PetscInt n,Vec *newv)
270: {
272: Vec v;
275: VecCreate(PetscObjectComm((PetscObject)snes),&v);
276: VecSetSizes(v,n,PETSC_DECIDE);
277: VecSetType(v,VECSTANDARD);
278: *newv = v;
279: return(0);
280: }
282: /* Resets the snes PC and KSP when the active set sizes change */
283: PetscErrorCode SNESVIResetPCandKSP(SNES snes,Mat Amat,Mat Pmat)
284: {
286: KSP snesksp;
289: SNESGetKSP(snes,&snesksp);
290: KSPReset(snesksp);
291: KSPResetFromOptions(snesksp);
293: /*
294: KSP kspnew;
295: PC pcnew;
296: MatSolverType stype;
298: KSPCreate(PetscObjectComm((PetscObject)snes),&kspnew);
299: kspnew->pc_side = snesksp->pc_side;
300: kspnew->rtol = snesksp->rtol;
301: kspnew->abstol = snesksp->abstol;
302: kspnew->max_it = snesksp->max_it;
303: KSPSetType(kspnew,((PetscObject)snesksp)->type_name);
304: KSPGetPC(kspnew,&pcnew);
305: PCSetType(kspnew->pc,((PetscObject)snesksp->pc)->type_name);
306: PCSetOperators(kspnew->pc,Amat,Pmat);
307: PCFactorGetMatSolverType(snesksp->pc,&stype);
308: PCFactorSetMatSolverType(kspnew->pc,stype);
309: KSPDestroy(&snesksp);
310: snes->ksp = kspnew;
311: PetscLogObjectParent((PetscObject)snes,(PetscObject)kspnew);
312: KSPSetFromOptions(kspnew);*/
313: return(0);
314: }
316: /* Variational Inequality solver using reduce space method. No semismooth algorithm is
317: implemented in this algorithm. It basically identifies the active constraints and does
318: a linear solve on the other variables (those not associated with the active constraints). */
319: PetscErrorCode SNESSolve_VINEWTONRSLS(SNES snes)
320: {
321: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data;
322: PetscErrorCode ierr;
323: PetscInt maxits,i,lits;
324: SNESLineSearchReason lssucceed;
325: PetscReal fnorm,gnorm,xnorm=0,ynorm;
326: Vec Y,X,F;
327: KSPConvergedReason kspreason;
328: KSP ksp;
329: PC pc;
332: /* Multigrid must use Galerkin for coarse grids with active set/reduced space methods; cannot rediscretize on coarser grids*/
333: SNESGetKSP(snes,&ksp);
334: KSPGetPC(ksp,&pc);
335: PCMGSetGalerkin(pc,PC_MG_GALERKIN_BOTH);
337: snes->numFailures = 0;
338: snes->numLinearSolveFailures = 0;
339: snes->reason = SNES_CONVERGED_ITERATING;
341: maxits = snes->max_its; /* maximum number of iterations */
342: X = snes->vec_sol; /* solution vector */
343: F = snes->vec_func; /* residual vector */
344: Y = snes->work[0]; /* work vectors */
346: SNESLineSearchSetVIFunctions(snes->linesearch, SNESVIProjectOntoBounds, SNESVIComputeInactiveSetFnorm);
347: SNESLineSearchSetVecs(snes->linesearch, X, NULL, NULL, NULL, NULL);
348: SNESLineSearchSetUp(snes->linesearch);
350: PetscObjectSAWsTakeAccess((PetscObject)snes);
351: snes->iter = 0;
352: snes->norm = 0.0;
353: PetscObjectSAWsGrantAccess((PetscObject)snes);
355: SNESVIProjectOntoBounds(snes,X);
356: SNESComputeFunction(snes,X,F);
357: SNESVIComputeInactiveSetFnorm(snes,F,X,&fnorm);
358: VecNorm(X,NORM_2,&xnorm); /* xnorm <- ||x|| */
359: SNESCheckFunctionNorm(snes,fnorm);
360: PetscObjectSAWsTakeAccess((PetscObject)snes);
361: snes->norm = fnorm;
362: PetscObjectSAWsGrantAccess((PetscObject)snes);
363: SNESLogConvergenceHistory(snes,fnorm,0);
364: SNESMonitor(snes,0,fnorm);
366: /* test convergence */
367: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
368: if (snes->reason) return(0);
370: for (i=0; i<maxits; i++) {
372: IS IS_act; /* _act -> active set _inact -> inactive set */
373: IS IS_redact; /* redundant active set */
374: VecScatter scat_act,scat_inact;
375: PetscInt nis_act,nis_inact;
376: Vec Y_act,Y_inact,F_inact;
377: Mat jac_inact_inact,prejac_inact_inact;
378: PetscBool isequal;
380: /* Call general purpose update function */
381: if (snes->ops->update) {
382: (*snes->ops->update)(snes, snes->iter);
383: }
384: SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
385: SNESCheckJacobianDomainerror(snes);
387: /* Create active and inactive index sets */
389: /*original
390: SNESVICreateIndexSets_RS(snes,X,F,&IS_act,&vi->IS_inact);
391: */
392: SNESVIGetActiveSetIS(snes,X,F,&IS_act);
394: if (vi->checkredundancy) {
395: (*vi->checkredundancy)(snes,IS_act,&IS_redact,vi->ctxP);
396: if (IS_redact) {
397: ISSort(IS_redact);
398: ISComplement(IS_redact,X->map->rstart,X->map->rend,&vi->IS_inact);
399: ISDestroy(&IS_redact);
400: } else {
401: ISComplement(IS_act,X->map->rstart,X->map->rend,&vi->IS_inact);
402: }
403: } else {
404: ISComplement(IS_act,X->map->rstart,X->map->rend,&vi->IS_inact);
405: }
407: /* Create inactive set submatrix */
408: MatCreateSubMatrix(snes->jacobian,vi->IS_inact,vi->IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);
410: if (0) { /* Dead code (temporary developer hack) */
411: IS keptrows;
412: MatFindNonzeroRows(jac_inact_inact,&keptrows);
413: if (keptrows) {
414: PetscInt cnt,*nrows,k;
415: const PetscInt *krows,*inact;
416: PetscInt rstart;
418: MatGetOwnershipRange(jac_inact_inact,&rstart,NULL);
419: MatDestroy(&jac_inact_inact);
420: ISDestroy(&IS_act);
422: ISGetLocalSize(keptrows,&cnt);
423: ISGetIndices(keptrows,&krows);
424: ISGetIndices(vi->IS_inact,&inact);
425: PetscMalloc1(cnt,&nrows);
426: for (k=0; k<cnt; k++) nrows[k] = inact[krows[k]-rstart];
427: ISRestoreIndices(keptrows,&krows);
428: ISRestoreIndices(vi->IS_inact,&inact);
429: ISDestroy(&keptrows);
430: ISDestroy(&vi->IS_inact);
432: ISCreateGeneral(PetscObjectComm((PetscObject)snes),cnt,nrows,PETSC_OWN_POINTER,&vi->IS_inact);
433: ISComplement(vi->IS_inact,F->map->rstart,F->map->rend,&IS_act);
434: MatCreateSubMatrix(snes->jacobian,vi->IS_inact,vi->IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);
435: }
436: }
437: DMSetVI(snes->dm,vi->IS_inact);
438: /* remove later */
440: /*
441: VecView(vi->xu,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xu))->comm));
442: VecView(vi->xl,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xl))->comm));
443: VecView(X,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)X)));
444: VecView(F,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)F)));
445: ISView(vi->IS_inact,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)vi->IS_inact)));
446: */
448: /* Get sizes of active and inactive sets */
449: ISGetLocalSize(IS_act,&nis_act);
450: ISGetLocalSize(vi->IS_inact,&nis_inact);
452: /* Create active and inactive set vectors */
453: SNESCreateSubVectors_VINEWTONRSLS(snes,nis_inact,&F_inact);
454: SNESCreateSubVectors_VINEWTONRSLS(snes,nis_act,&Y_act);
455: SNESCreateSubVectors_VINEWTONRSLS(snes,nis_inact,&Y_inact);
457: /* Create scatter contexts */
458: VecScatterCreate(Y,IS_act,Y_act,NULL,&scat_act);
459: VecScatterCreate(Y,vi->IS_inact,Y_inact,NULL,&scat_inact);
461: /* Do a vec scatter to active and inactive set vectors */
462: VecScatterBegin(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);
463: VecScatterEnd(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);
465: VecScatterBegin(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);
466: VecScatterEnd(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);
468: VecScatterBegin(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);
469: VecScatterEnd(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);
471: /* Active set direction = 0 */
472: VecSet(Y_act,0);
473: if (snes->jacobian != snes->jacobian_pre) {
474: MatCreateSubMatrix(snes->jacobian_pre,vi->IS_inact,vi->IS_inact,MAT_INITIAL_MATRIX,&prejac_inact_inact);
475: } else prejac_inact_inact = jac_inact_inact;
477: ISEqual(vi->IS_inact_prev,vi->IS_inact,&isequal);
478: if (!isequal) {
479: SNESVIResetPCandKSP(snes,jac_inact_inact,prejac_inact_inact);
480: PCFieldSplitRestrictIS(pc,vi->IS_inact);
481: }
483: /* ISView(vi->IS_inact,0); */
484: /* ISView(IS_act,0);*/
485: /* MatView(snes->jacobian_pre,0); */
487: KSPSetOperators(snes->ksp,jac_inact_inact,prejac_inact_inact);
488: KSPSetUp(snes->ksp);
489: {
490: PC pc;
491: PetscBool flg;
492: KSPGetPC(snes->ksp,&pc);
493: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&flg);
494: if (flg) {
495: KSP *subksps;
496: PCFieldSplitGetSubKSP(pc,NULL,&subksps);
497: KSPGetPC(subksps[0],&pc);
498: PetscFree(subksps);
499: PetscObjectTypeCompare((PetscObject)pc,PCBJACOBI,&flg);
500: if (flg) {
501: PetscInt n,N = 101*101,j,cnts[3] = {0,0,0};
502: const PetscInt *ii;
504: ISGetSize(vi->IS_inact,&n);
505: ISGetIndices(vi->IS_inact,&ii);
506: for (j=0; j<n; j++) {
507: if (ii[j] < N) cnts[0]++;
508: else if (ii[j] < 2*N) cnts[1]++;
509: else if (ii[j] < 3*N) cnts[2]++;
510: }
511: ISRestoreIndices(vi->IS_inact,&ii);
513: PCBJacobiSetTotalBlocks(pc,3,cnts);
514: }
515: }
516: }
518: KSPSolve(snes->ksp,F_inact,Y_inact);
519: VecScatterBegin(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);
520: VecScatterEnd(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);
521: VecScatterBegin(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);
522: VecScatterEnd(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);
524: VecDestroy(&F_inact);
525: VecDestroy(&Y_act);
526: VecDestroy(&Y_inact);
527: VecScatterDestroy(&scat_act);
528: VecScatterDestroy(&scat_inact);
529: ISDestroy(&IS_act);
530: if (!isequal) {
531: ISDestroy(&vi->IS_inact_prev);
532: ISDuplicate(vi->IS_inact,&vi->IS_inact_prev);
533: }
534: ISDestroy(&vi->IS_inact);
535: MatDestroy(&jac_inact_inact);
536: if (snes->jacobian != snes->jacobian_pre) {
537: MatDestroy(&prejac_inact_inact);
538: }
540: KSPGetConvergedReason(snes->ksp,&kspreason);
541: if (kspreason < 0) {
542: if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
543: PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);
544: snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
545: break;
546: }
547: }
549: KSPGetIterationNumber(snes->ksp,&lits);
550: snes->linear_its += lits;
551: PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);
552: /*
553: if (snes->ops->precheck) {
554: PetscBool changed_y = PETSC_FALSE;
555: (*snes->ops->precheck)(snes,X,Y,snes->precheck,&changed_y);
556: }
558: if (PetscLogPrintInfo) {
559: SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);
560: }
561: */
562: /* Compute a (scaled) negative update in the line search routine:
563: Y <- X - lambda*Y
564: and evaluate G = function(Y) (depends on the line search).
565: */
566: VecCopy(Y,snes->vec_sol_update);
567: ynorm = 1; gnorm = fnorm;
568: SNESLineSearchApply(snes->linesearch, X, F, &gnorm, Y);
569: SNESLineSearchGetReason(snes->linesearch, &lssucceed);
570: SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm);
571: PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",(double)fnorm,(double)gnorm,(double)ynorm,(int)lssucceed);
572: if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
573: if (snes->domainerror) {
574: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
575: DMDestroyVI(snes->dm);
576: return(0);
577: }
578: if (lssucceed) {
579: if (++snes->numFailures >= snes->maxFailures) {
580: PetscBool ismin;
581: snes->reason = SNES_DIVERGED_LINE_SEARCH;
582: SNESVICheckLocalMin_Private(snes,snes->jacobian,F,X,gnorm,&ismin);
583: if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
584: break;
585: }
586: }
587: DMDestroyVI(snes->dm);
588: /* Update function and solution vectors */
589: fnorm = gnorm;
590: /* Monitor convergence */
591: PetscObjectSAWsTakeAccess((PetscObject)snes);
592: snes->iter = i+1;
593: snes->norm = fnorm;
594: snes->xnorm = xnorm;
595: snes->ynorm = ynorm;
596: PetscObjectSAWsGrantAccess((PetscObject)snes);
597: SNESLogConvergenceHistory(snes,snes->norm,lits);
598: SNESMonitor(snes,snes->iter,snes->norm);
599: /* Test for convergence, xnorm = || X || */
600: if (snes->ops->converged != SNESConvergedSkip) { VecNorm(X,NORM_2,&xnorm); }
601: (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);
602: if (snes->reason) break;
603: }
604: /* 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 */
605: DMDestroyVI(snes->dm);
606: if (i == maxits) {
607: PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);
608: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
609: }
610: return(0);
611: }
613: PetscErrorCode SNESVISetRedundancyCheck(SNES snes,PetscErrorCode (*func)(SNES,IS,IS*,void*),void *ctx)
614: {
615: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data;
619: vi->checkredundancy = func;
620: vi->ctxP = ctx;
621: return(0);
622: }
624: #if defined(PETSC_HAVE_MATLAB_ENGINE)
625: #include <engine.h>
626: #include <mex.h>
627: typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext;
629: PetscErrorCode SNESVIRedundancyCheck_Matlab(SNES snes,IS is_act,IS *is_redact,void *ctx)
630: {
631: PetscErrorCode ierr;
632: SNESMatlabContext *sctx = (SNESMatlabContext*)ctx;
633: int nlhs = 1, nrhs = 5;
634: mxArray *plhs[1], *prhs[5];
635: long long int l1 = 0, l2 = 0, ls = 0;
636: PetscInt *indices=NULL;
644: /* Create IS for reduced active set of size 0, its size and indices will
645: bet set by the Matlab function */
646: ISCreateGeneral(PetscObjectComm((PetscObject)snes),0,indices,PETSC_OWN_POINTER,is_redact);
647: /* call Matlab function in ctx */
648: PetscArraycpy(&ls,&snes,1);
649: PetscArraycpy(&l1,&is_act,1);
650: PetscArraycpy(&l2,is_redact,1);
651: prhs[0] = mxCreateDoubleScalar((double)ls);
652: prhs[1] = mxCreateDoubleScalar((double)l1);
653: prhs[2] = mxCreateDoubleScalar((double)l2);
654: prhs[3] = mxCreateString(sctx->funcname);
655: prhs[4] = sctx->ctx;
656: mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESVIRedundancyCheckInternal");
657: mxGetScalar(plhs[0]);
658: mxDestroyArray(prhs[0]);
659: mxDestroyArray(prhs[1]);
660: mxDestroyArray(prhs[2]);
661: mxDestroyArray(prhs[3]);
662: mxDestroyArray(plhs[0]);
663: return(0);
664: }
666: PetscErrorCode SNESVISetRedundancyCheckMatlab(SNES snes,const char *func,mxArray *ctx)
667: {
668: PetscErrorCode ierr;
669: SNESMatlabContext *sctx;
672: /* currently sctx is memory bleed */
673: PetscNew(&sctx);
674: PetscStrallocpy(func,&sctx->funcname);
675: sctx->ctx = mxDuplicateArray(ctx);
676: SNESVISetRedundancyCheck(snes,SNESVIRedundancyCheck_Matlab,sctx);
677: return(0);
678: }
680: #endif
682: /* -------------------------------------------------------------------------- */
683: /*
684: SNESSetUp_VINEWTONRSLS - Sets up the internal data structures for the later use
685: of the SNESVI nonlinear solver.
687: Input Parameter:
688: . snes - the SNES context
690: Application Interface Routine: SNESSetUp()
692: Notes:
693: For basic use of the SNES solvers, the user need not explicitly call
694: SNESSetUp(), since these actions will automatically occur during
695: the call to SNESSolve().
696: */
697: PetscErrorCode SNESSetUp_VINEWTONRSLS(SNES snes)
698: {
699: PetscErrorCode ierr;
700: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*) snes->data;
701: PetscInt *indices;
702: PetscInt i,n,rstart,rend;
703: SNESLineSearch linesearch;
706: SNESSetUp_VI(snes);
708: /* Set up previous active index set for the first snes solve
709: vi->IS_inact_prev = 0,1,2,....N */
711: VecGetOwnershipRange(snes->vec_sol,&rstart,&rend);
712: VecGetLocalSize(snes->vec_sol,&n);
713: PetscMalloc1(n,&indices);
714: for (i=0; i < n; i++) indices[i] = rstart + i;
715: ISCreateGeneral(PetscObjectComm((PetscObject)snes),n,indices,PETSC_OWN_POINTER,&vi->IS_inact_prev);
717: /* set the line search functions */
718: if (!snes->linesearch) {
719: SNESGetLineSearch(snes, &linesearch);
720: SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
721: }
722: return(0);
723: }
724: /* -------------------------------------------------------------------------- */
725: PetscErrorCode SNESReset_VINEWTONRSLS(SNES snes)
726: {
727: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*) snes->data;
728: PetscErrorCode ierr;
731: SNESReset_VI(snes);
732: ISDestroy(&vi->IS_inact_prev);
733: return(0);
734: }
736: /* -------------------------------------------------------------------------- */
737: /*MC
738: SNESVINEWTONRSLS - Reduced space active set solvers for variational inequalities based on Newton's method
740: Options Database:
741: + -snes_type <vinewtonssls,vinewtonrsls> - a semi-smooth solver, a reduced space active set method
742: - -snes_vi_monitor - prints the number of active constraints at each iteration.
744: Level: beginner
746: References:
747: . 1. - T. S. Munson, and S. Benson. Flexible Complementarity Solvers for Large Scale
748: Applications, Optimization Methods and Software, 21 (2006).
750: .seealso: SNESVISetVariableBounds(), SNESVISetComputeVariableBounds(), SNESCreate(), SNES, SNESSetType(), SNESVINEWTONSSLS, SNESNEWTONTR, SNESLineSearchSetType(),SNESLineSearchSetPostCheck(), SNESLineSearchSetPreCheck()
752: M*/
753: PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONRSLS(SNES snes)
754: {
755: PetscErrorCode ierr;
756: SNES_VINEWTONRSLS *vi;
757: SNESLineSearch linesearch;
760: snes->ops->reset = SNESReset_VINEWTONRSLS;
761: snes->ops->setup = SNESSetUp_VINEWTONRSLS;
762: snes->ops->solve = SNESSolve_VINEWTONRSLS;
763: snes->ops->destroy = SNESDestroy_VI;
764: snes->ops->setfromoptions = SNESSetFromOptions_VI;
765: snes->ops->view = NULL;
766: snes->ops->converged = SNESConvergedDefault_VI;
768: snes->usesksp = PETSC_TRUE;
769: snes->usesnpc = PETSC_FALSE;
771: SNESGetLineSearch(snes, &linesearch);
772: if (!((PetscObject)linesearch)->type_name) {
773: SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
774: }
775: SNESLineSearchBTSetAlpha(linesearch, 0.0);
777: snes->alwayscomputesfinalresidual = PETSC_TRUE;
779: PetscNewLog(snes,&vi);
780: snes->data = (void*)vi;
781: vi->checkredundancy = NULL;
783: PetscObjectComposeFunction((PetscObject)snes,"SNESVISetVariableBounds_C",SNESVISetVariableBounds_VI);
784: PetscObjectComposeFunction((PetscObject)snes,"SNESVISetComputeVariableBounds_C",SNESVISetComputeVariableBounds_VI);
785: return(0);
786: }