Actual source code: virs.c
petsc-3.4.5 2014-06-29
2: #include <../src/snes/impls/vi/rs/virsimpl.h> /*I "petscsnes.h" I*/
3: #include <petsc-private/kspimpl.h>
4: #include <petsc-private/matimpl.h>
5: #include <petsc-private/dmimpl.h>
6: #include <petsc-private/vecimpl.h>
10: /*
11: SNESVIGetInactiveSet - Gets the global indices for the inactive set variables (these correspond to the degrees of freedom the linear
12: system is solved on)
14: Input parameter
15: . snes - the SNES context
17: Output parameter
18: . ISact - active set index set
20: */
21: PetscErrorCode SNESVIGetInactiveSet(SNES snes,IS *inact)
22: {
23: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data;
26: *inact = vi->IS_inact_prev;
27: return(0);
28: }
30: /*
31: 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
32: defined by the reduced space method.
34: Simple calls the regular DM interpolation and restricts it to operation on the variables not associated with active constraints.
36: <*/
37: typedef struct {
38: PetscInt n; /* size of vectors in the reduced DM space */
39: IS inactive;
41: PetscErrorCode (*createinterpolation)(DM,DM,Mat*,Vec*); /* DM's original routines */
42: PetscErrorCode (*coarsen)(DM, MPI_Comm, DM*);
43: PetscErrorCode (*createglobalvector)(DM,Vec*);
45: DM dm; /* when destroying this object we need to reset the above function into the base DM */
46: } DM_SNESVI;
50: /*
51: DMCreateGlobalVector_SNESVI - Creates global vector of the size of the reduced space
53: */
54: PetscErrorCode DMCreateGlobalVector_SNESVI(DM dm,Vec *vec)
55: {
57: PetscContainer isnes;
58: DM_SNESVI *dmsnesvi;
61: PetscObjectQuery((PetscObject)dm,"VI",(PetscObject*)&isnes);
62: if (!isnes) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Composed SNES is missing");
63: PetscContainerGetPointer(isnes,(void**)&dmsnesvi);
64: VecCreateMPI(PetscObjectComm((PetscObject)dm),dmsnesvi->n,PETSC_DETERMINE,vec);
65: return(0);
66: }
70: /*
71: DMCreateInterpolation_SNESVI - Modifieds the interpolation obtained from the DM by removing all rows and columns associated with active constraints.
73: */
74: PetscErrorCode DMCreateInterpolation_SNESVI(DM dm1,DM dm2,Mat *mat,Vec *vec)
75: {
77: PetscContainer isnes;
78: DM_SNESVI *dmsnesvi1,*dmsnesvi2;
79: Mat interp;
82: PetscObjectQuery((PetscObject)dm1,"VI",(PetscObject*)&isnes);
83: if (!isnes) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_PLIB,"Composed VI data structure is missing");
84: PetscContainerGetPointer(isnes,(void**)&dmsnesvi1);
85: PetscObjectQuery((PetscObject)dm2,"VI",(PetscObject*)&isnes);
86: if (!isnes) SETERRQ(PetscObjectComm((PetscObject)dm2),PETSC_ERR_PLIB,"Composed VI data structure is missing");
87: PetscContainerGetPointer(isnes,(void**)&dmsnesvi2);
89: (*dmsnesvi1->createinterpolation)(dm1,dm2,&interp,NULL);
90: MatGetSubMatrix(interp,dmsnesvi2->inactive,dmsnesvi1->inactive,MAT_INITIAL_MATRIX,mat);
91: MatDestroy(&interp);
92: *vec = 0;
93: return(0);
94: }
96: extern PetscErrorCode DMSetVI(DM,IS);
100: /*
101: DMCoarsen_SNESVI - Computes the regular coarsened DM then computes additional information about its inactive set
103: */
104: PetscErrorCode DMCoarsen_SNESVI(DM dm1,MPI_Comm comm,DM *dm2)
105: {
107: PetscContainer isnes;
108: DM_SNESVI *dmsnesvi1;
109: Vec finemarked,coarsemarked;
110: IS inactive;
111: VecScatter inject;
112: const PetscInt *index;
113: PetscInt n,k,cnt = 0,rstart,*coarseindex;
114: PetscScalar *marked;
117: PetscObjectQuery((PetscObject)dm1,"VI",(PetscObject*)&isnes);
118: if (!isnes) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_PLIB,"Composed VI data structure is missing");
119: PetscContainerGetPointer(isnes,(void**)&dmsnesvi1);
121: /* get the original coarsen */
122: (*dmsnesvi1->coarsen)(dm1,comm,dm2);
124: /* not sure why this extra reference is needed, but without the dm2 disappears too early */
125: PetscObjectReference((PetscObject)*dm2);
127: /* need to set back global vectors in order to use the original injection */
128: DMClearGlobalVectors(dm1);
130: dm1->ops->createglobalvector = dmsnesvi1->createglobalvector;
132: DMCreateGlobalVector(dm1,&finemarked);
133: DMCreateGlobalVector(*dm2,&coarsemarked);
135: /*
136: fill finemarked with locations of inactive points
137: */
138: ISGetIndices(dmsnesvi1->inactive,&index);
139: ISGetLocalSize(dmsnesvi1->inactive,&n);
140: VecSet(finemarked,0.0);
141: for (k=0; k<n; k++) {
142: VecSetValue(finemarked,index[k],1.0,INSERT_VALUES);
143: }
144: VecAssemblyBegin(finemarked);
145: VecAssemblyEnd(finemarked);
147: DMCreateInjection(*dm2,dm1,&inject);
148: VecScatterBegin(inject,finemarked,coarsemarked,INSERT_VALUES,SCATTER_FORWARD);
149: VecScatterEnd(inject,finemarked,coarsemarked,INSERT_VALUES,SCATTER_FORWARD);
150: VecScatterDestroy(&inject);
152: /*
153: create index set list of coarse inactive points from coarsemarked
154: */
155: VecGetLocalSize(coarsemarked,&n);
156: VecGetOwnershipRange(coarsemarked,&rstart,NULL);
157: VecGetArray(coarsemarked,&marked);
158: for (k=0; k<n; k++) {
159: if (marked[k] != 0.0) cnt++;
160: }
161: PetscMalloc(cnt*sizeof(PetscInt),&coarseindex);
162: cnt = 0;
163: for (k=0; k<n; k++) {
164: if (marked[k] != 0.0) coarseindex[cnt++] = k + rstart;
165: }
166: VecRestoreArray(coarsemarked,&marked);
167: ISCreateGeneral(PetscObjectComm((PetscObject)coarsemarked),cnt,coarseindex,PETSC_OWN_POINTER,&inactive);
169: DMClearGlobalVectors(dm1);
171: dm1->ops->createglobalvector = DMCreateGlobalVector_SNESVI;
173: DMSetVI(*dm2,inactive);
175: VecDestroy(&finemarked);
176: VecDestroy(&coarsemarked);
177: ISDestroy(&inactive);
178: return(0);
179: }
183: PetscErrorCode DMDestroy_SNESVI(DM_SNESVI *dmsnesvi)
184: {
188: /* reset the base methods in the DM object that were changed when the DM_SNESVI was reset */
189: dmsnesvi->dm->ops->createinterpolation = dmsnesvi->createinterpolation;
190: dmsnesvi->dm->ops->coarsen = dmsnesvi->coarsen;
191: dmsnesvi->dm->ops->createglobalvector = dmsnesvi->createglobalvector;
192: /* need to clear out this vectors because some of them may not have a reference to the DM
193: but they are counted as having references to the DM in DMDestroy() */
194: DMClearGlobalVectors(dmsnesvi->dm);
196: ISDestroy(&dmsnesvi->inactive);
197: PetscFree(dmsnesvi);
198: return(0);
199: }
203: /*
204: DMSetVI - Marks a DM as associated with a VI problem. This causes the interpolation/restriction operators to
205: be restricted to only those variables NOT associated with active constraints.
207: */
208: PetscErrorCode DMSetVI(DM dm,IS inactive)
209: {
211: PetscContainer isnes;
212: DM_SNESVI *dmsnesvi;
215: if (!dm) return(0);
217: PetscObjectReference((PetscObject)inactive);
219: PetscObjectQuery((PetscObject)dm,"VI",(PetscObject*)&isnes);
220: if (!isnes) {
221: PetscContainerCreate(PetscObjectComm((PetscObject)dm),&isnes);
222: PetscContainerSetUserDestroy(isnes,(PetscErrorCode (*)(void*))DMDestroy_SNESVI);
223: PetscNew(DM_SNESVI,&dmsnesvi);
224: PetscContainerSetPointer(isnes,(void*)dmsnesvi);
225: PetscObjectCompose((PetscObject)dm,"VI",(PetscObject)isnes);
226: PetscContainerDestroy(&isnes);
228: dmsnesvi->createinterpolation = dm->ops->createinterpolation;
229: dm->ops->createinterpolation = DMCreateInterpolation_SNESVI;
230: dmsnesvi->coarsen = dm->ops->coarsen;
231: dm->ops->coarsen = DMCoarsen_SNESVI;
232: dmsnesvi->createglobalvector = dm->ops->createglobalvector;
233: dm->ops->createglobalvector = DMCreateGlobalVector_SNESVI;
234: } else {
235: PetscContainerGetPointer(isnes,(void**)&dmsnesvi);
236: ISDestroy(&dmsnesvi->inactive);
237: }
238: DMClearGlobalVectors(dm);
239: ISGetLocalSize(inactive,&dmsnesvi->n);
241: dmsnesvi->inactive = inactive;
242: dmsnesvi->dm = dm;
243: return(0);
244: }
248: /*
249: DMDestroyVI - Frees the DM_SNESVI object contained in the DM
250: - also resets the function pointers in the DM for createinterpolation() etc to use the original DM
251: */
252: PetscErrorCode DMDestroyVI(DM dm)
253: {
257: if (!dm) return(0);
258: PetscObjectCompose((PetscObject)dm,"VI",(PetscObject)NULL);
259: return(0);
260: }
262: /* --------------------------------------------------------------------------------------------------------*/
269: PetscErrorCode SNESCreateIndexSets_VINEWTONRSLS(SNES snes,Vec X,Vec F,IS *ISact,IS *ISinact)
270: {
274: SNESVIGetActiveSetIS(snes,X,F,ISact);
275: ISComplement(*ISact,X->map->rstart,X->map->rend,ISinact);
276: return(0);
277: }
279: /* Create active and inactive set vectors. The local size of this vector is set and petsc computes the global size */
282: PetscErrorCode SNESCreateSubVectors_VINEWTONRSLS(SNES snes,PetscInt n,Vec *newv)
283: {
285: Vec v;
288: VecCreate(PetscObjectComm((PetscObject)snes),&v);
289: VecSetSizes(v,n,PETSC_DECIDE);
290: VecSetFromOptions(v);
291: *newv = v;
292: return(0);
293: }
295: /* Resets the snes PC and KSP when the active set sizes change */
298: PetscErrorCode SNESVIResetPCandKSP(SNES snes,Mat Amat,Mat Pmat)
299: {
301: KSP snesksp;
304: SNESGetKSP(snes,&snesksp);
305: KSPReset(snesksp);
307: /*
308: KSP kspnew;
309: PC pcnew;
310: const MatSolverPackage stype;
313: KSPCreate(PetscObjectComm((PetscObject)snes),&kspnew);
314: kspnew->pc_side = snesksp->pc_side;
315: kspnew->rtol = snesksp->rtol;
316: kspnew->abstol = snesksp->abstol;
317: kspnew->max_it = snesksp->max_it;
318: KSPSetType(kspnew,((PetscObject)snesksp)->type_name);
319: KSPGetPC(kspnew,&pcnew);
320: PCSetType(kspnew->pc,((PetscObject)snesksp->pc)->type_name);
321: PCSetOperators(kspnew->pc,Amat,Pmat,DIFFERENT_NONZERO_PATTERN);
322: PCFactorGetMatSolverPackage(snesksp->pc,&stype);
323: PCFactorSetMatSolverPackage(kspnew->pc,stype);
324: KSPDestroy(&snesksp);
325: snes->ksp = kspnew;
326: PetscLogObjectParent(snes,kspnew);
327: KSPSetFromOptions(kspnew);*/
328: return(0);
329: }
331: /* Variational Inequality solver using reduce space method. No semismooth algorithm is
332: implemented in this algorithm. It basically identifies the active constraints and does
333: a linear solve on the other variables (those not associated with the active constraints). */
336: PetscErrorCode SNESSolve_VINEWTONRSLS(SNES snes)
337: {
338: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data;
339: PetscErrorCode ierr;
340: PetscInt maxits,i,lits;
341: PetscBool lssucceed;
342: MatStructure flg = DIFFERENT_NONZERO_PATTERN;
343: PetscReal fnorm,gnorm,xnorm=0,ynorm;
344: Vec Y,X,F;
345: KSPConvergedReason kspreason;
348: snes->numFailures = 0;
349: snes->numLinearSolveFailures = 0;
350: snes->reason = SNES_CONVERGED_ITERATING;
352: maxits = snes->max_its; /* maximum number of iterations */
353: X = snes->vec_sol; /* solution vector */
354: F = snes->vec_func; /* residual vector */
355: Y = snes->work[0]; /* work vectors */
357: SNESLineSearchSetVIFunctions(snes->linesearch, SNESVIProjectOntoBounds, SNESVIComputeInactiveSetFnorm);
358: SNESLineSearchSetVecs(snes->linesearch, X, NULL, NULL, NULL, NULL);
359: SNESLineSearchSetUp(snes->linesearch);
361: PetscObjectAMSTakeAccess((PetscObject)snes);
362: snes->iter = 0;
363: snes->norm = 0.0;
364: PetscObjectAMSGrantAccess((PetscObject)snes);
366: SNESVIProjectOntoBounds(snes,X);
367: SNESComputeFunction(snes,X,F);
368: if (snes->domainerror) {
369: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
370: return(0);
371: }
372: SNESVIComputeInactiveSetFnorm(snes,F,X,&fnorm);
373: VecNormBegin(X,NORM_2,&xnorm); /* xnorm <- ||x|| */
374: VecNormEnd(X,NORM_2,&xnorm);
375: if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
377: PetscObjectAMSTakeAccess((PetscObject)snes);
378: snes->norm = fnorm;
379: PetscObjectAMSGrantAccess((PetscObject)snes);
380: SNESLogConvergenceHistory(snes,fnorm,0);
381: SNESMonitor(snes,0,fnorm);
383: /* set parameter for default relative tolerance convergence test */
384: snes->ttol = fnorm*snes->rtol;
385: /* test convergence */
386: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
387: if (snes->reason) return(0);
390: for (i=0; i<maxits; i++) {
392: IS IS_act,IS_inact; /* _act -> active set _inact -> inactive set */
393: IS IS_redact; /* redundant active set */
394: VecScatter scat_act,scat_inact;
395: PetscInt nis_act,nis_inact;
396: Vec Y_act,Y_inact,F_inact;
397: Mat jac_inact_inact,prejac_inact_inact;
398: PetscBool isequal;
400: /* Call general purpose update function */
401: if (snes->ops->update) {
402: (*snes->ops->update)(snes, snes->iter);
403: }
404: SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);
407: /* Create active and inactive index sets */
409: /*original
410: SNESVICreateIndexSets_RS(snes,X,F,&IS_act,&IS_inact);
411: */
412: SNESVIGetActiveSetIS(snes,X,F,&IS_act);
414: if (vi->checkredundancy) {
415: (*vi->checkredundancy)(snes,IS_act,&IS_redact,vi->ctxP);
416: if (IS_redact) {
417: ISSort(IS_redact);
418: ISComplement(IS_redact,X->map->rstart,X->map->rend,&IS_inact);
419: ISDestroy(&IS_redact);
420: } else {
421: ISComplement(IS_act,X->map->rstart,X->map->rend,&IS_inact);
422: }
423: } else {
424: ISComplement(IS_act,X->map->rstart,X->map->rend,&IS_inact);
425: }
428: /* Create inactive set submatrix */
429: MatGetSubMatrix(snes->jacobian,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);
431: if (0) { /* Dead code (temporary developer hack) */
432: IS keptrows;
433: MatFindNonzeroRows(jac_inact_inact,&keptrows);
434: if (keptrows) {
435: PetscInt cnt,*nrows,k;
436: const PetscInt *krows,*inact;
437: PetscInt rstart=jac_inact_inact->rmap->rstart;
439: MatDestroy(&jac_inact_inact);
440: ISDestroy(&IS_act);
442: ISGetLocalSize(keptrows,&cnt);
443: ISGetIndices(keptrows,&krows);
444: ISGetIndices(IS_inact,&inact);
445: PetscMalloc(cnt*sizeof(PetscInt),&nrows);
446: for (k=0; k<cnt; k++) nrows[k] = inact[krows[k]-rstart];
447: ISRestoreIndices(keptrows,&krows);
448: ISRestoreIndices(IS_inact,&inact);
449: ISDestroy(&keptrows);
450: ISDestroy(&IS_inact);
452: ISCreateGeneral(PetscObjectComm((PetscObject)snes),cnt,nrows,PETSC_OWN_POINTER,&IS_inact);
453: ISComplement(IS_inact,F->map->rstart,F->map->rend,&IS_act);
454: MatGetSubMatrix(snes->jacobian,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);
455: }
456: }
457: DMSetVI(snes->dm,IS_inact);
458: /* remove later */
460: /*
461: VecView(vi->xu,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xu))->comm));
462: VecView(vi->xl,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xl))->comm));
463: VecView(X,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)X)));
464: VecView(F,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)F)));
465: ISView(IS_inact,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)IS_inact)));
466: */
468: /* Get sizes of active and inactive sets */
469: ISGetLocalSize(IS_act,&nis_act);
470: ISGetLocalSize(IS_inact,&nis_inact);
472: /* Create active and inactive set vectors */
473: SNESCreateSubVectors_VINEWTONRSLS(snes,nis_inact,&F_inact);
474: SNESCreateSubVectors_VINEWTONRSLS(snes,nis_act,&Y_act);
475: SNESCreateSubVectors_VINEWTONRSLS(snes,nis_inact,&Y_inact);
477: /* Create scatter contexts */
478: VecScatterCreate(Y,IS_act,Y_act,NULL,&scat_act);
479: VecScatterCreate(Y,IS_inact,Y_inact,NULL,&scat_inact);
481: /* Do a vec scatter to active and inactive set vectors */
482: VecScatterBegin(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);
483: VecScatterEnd(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);
485: VecScatterBegin(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);
486: VecScatterEnd(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);
488: VecScatterBegin(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);
489: VecScatterEnd(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);
491: /* Active set direction = 0 */
492: VecSet(Y_act,0);
493: if (snes->jacobian != snes->jacobian_pre) {
494: MatGetSubMatrix(snes->jacobian_pre,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&prejac_inact_inact);
495: } else prejac_inact_inact = jac_inact_inact;
497: ISEqual(vi->IS_inact_prev,IS_inact,&isequal);
498: if (!isequal) {
499: SNESVIResetPCandKSP(snes,jac_inact_inact,prejac_inact_inact);
500: flg = DIFFERENT_NONZERO_PATTERN;
501: }
503: /* ISView(IS_inact,0); */
504: /* ISView(IS_act,0);*/
505: /* MatView(snes->jacobian_pre,0); */
509: KSPSetOperators(snes->ksp,jac_inact_inact,prejac_inact_inact,flg);
510: KSPSetUp(snes->ksp);
511: {
512: PC pc;
513: PetscBool flg;
514: KSPGetPC(snes->ksp,&pc);
515: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&flg);
516: if (flg) {
517: KSP *subksps;
518: PCFieldSplitGetSubKSP(pc,NULL,&subksps);
519: KSPGetPC(subksps[0],&pc);
520: PetscFree(subksps);
521: PetscObjectTypeCompare((PetscObject)pc,PCBJACOBI,&flg);
522: if (flg) {
523: PetscInt n,N = 101*101,j,cnts[3] = {0,0,0};
524: const PetscInt *ii;
526: ISGetSize(IS_inact,&n);
527: ISGetIndices(IS_inact,&ii);
528: for (j=0; j<n; j++) {
529: if (ii[j] < N) cnts[0]++;
530: else if (ii[j] < 2*N) cnts[1]++;
531: else if (ii[j] < 3*N) cnts[2]++;
532: }
533: ISRestoreIndices(IS_inact,&ii);
535: PCBJacobiSetTotalBlocks(pc,3,cnts);
536: }
537: }
538: }
540: KSPSolve(snes->ksp,F_inact,Y_inact);
541: KSPGetConvergedReason(snes->ksp,&kspreason);
542: if (kspreason < 0) {
543: if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
544: PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);
545: snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
546: break;
547: }
548: }
550: VecScatterBegin(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);
551: VecScatterEnd(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);
552: VecScatterBegin(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);
553: VecScatterEnd(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);
555: VecDestroy(&F_inact);
556: VecDestroy(&Y_act);
557: VecDestroy(&Y_inact);
558: VecScatterDestroy(&scat_act);
559: VecScatterDestroy(&scat_inact);
560: ISDestroy(&IS_act);
561: if (!isequal) {
562: ISDestroy(&vi->IS_inact_prev);
563: ISDuplicate(IS_inact,&vi->IS_inact_prev);
564: }
565: ISDestroy(&IS_inact);
566: MatDestroy(&jac_inact_inact);
567: if (snes->jacobian != snes->jacobian_pre) {
568: MatDestroy(&prejac_inact_inact);
569: }
570: KSPGetIterationNumber(snes->ksp,&lits);
571: snes->linear_its += lits;
572: PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);
573: /*
574: if (snes->ops->precheck) {
575: PetscBool changed_y = PETSC_FALSE;
576: (*snes->ops->precheck)(snes,X,Y,snes->precheck,&changed_y);
577: }
579: if (PetscLogPrintInfo) {
580: SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);
581: }
582: */
583: /* Compute a (scaled) negative update in the line search routine:
584: Y <- X - lambda*Y
585: and evaluate G = function(Y) (depends on the line search).
586: */
587: VecCopy(Y,snes->vec_sol_update);
588: ynorm = 1; gnorm = fnorm;
589: SNESLineSearchApply(snes->linesearch, X, F, &gnorm, Y);
590: SNESLineSearchGetSuccess(snes->linesearch, &lssucceed);
591: SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm);
592: PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",(double)fnorm,(double)gnorm,(double)ynorm,(int)lssucceed);
593: if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
594: if (snes->domainerror) {
595: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
596: DMDestroyVI(snes->dm);
597: return(0);
598: }
599: if (!lssucceed) {
600: if (++snes->numFailures >= snes->maxFailures) {
601: PetscBool ismin;
602: snes->reason = SNES_DIVERGED_LINE_SEARCH;
603: SNESVICheckLocalMin_Private(snes,snes->jacobian,F,X,gnorm,&ismin);
604: if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
605: break;
606: }
607: }
608: /* Update function and solution vectors */
609: fnorm = gnorm;
610: /* Monitor convergence */
611: PetscObjectAMSTakeAccess((PetscObject)snes);
612: snes->iter = i+1;
613: snes->norm = fnorm;
614: PetscObjectAMSGrantAccess((PetscObject)snes);
615: SNESLogConvergenceHistory(snes,snes->norm,lits);
616: SNESMonitor(snes,snes->iter,snes->norm);
617: /* Test for convergence, xnorm = || X || */
618: if (snes->ops->converged != SNESSkipConverged) { VecNorm(X,NORM_2,&xnorm); }
619: (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);
620: if (snes->reason) break;
621: }
622: DMDestroyVI(snes->dm);
623: if (i == maxits) {
624: PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);
625: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
626: }
627: return(0);
628: }
632: PetscErrorCode SNESVISetRedundancyCheck(SNES snes,PetscErrorCode (*func)(SNES,IS,IS*,void*),void *ctx)
633: {
634: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data;
638: vi->checkredundancy = func;
639: vi->ctxP = ctx;
640: return(0);
641: }
643: #if defined(PETSC_HAVE_MATLAB_ENGINE)
644: #include <engine.h>
645: #include <mex.h>
646: typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext;
650: PetscErrorCode SNESVIRedundancyCheck_Matlab(SNES snes,IS is_act,IS *is_redact,void *ctx)
651: {
652: PetscErrorCode ierr;
653: SNESMatlabContext *sctx = (SNESMatlabContext*)ctx;
654: int nlhs = 1, nrhs = 5;
655: mxArray *plhs[1], *prhs[5];
656: long long int l1 = 0, l2 = 0, ls = 0;
657: PetscInt *indices=NULL;
665: /* Create IS for reduced active set of size 0, its size and indices will
666: bet set by the Matlab function */
667: ISCreateGeneral(PetscObjectComm((PetscObject)snes),0,indices,PETSC_OWN_POINTER,is_redact);
668: /* call Matlab function in ctx */
669: PetscMemcpy(&ls,&snes,sizeof(snes));
670: PetscMemcpy(&l1,&is_act,sizeof(is_act));
671: PetscMemcpy(&l2,is_redact,sizeof(is_act));
672: prhs[0] = mxCreateDoubleScalar((double)ls);
673: prhs[1] = mxCreateDoubleScalar((double)l1);
674: prhs[2] = mxCreateDoubleScalar((double)l2);
675: prhs[3] = mxCreateString(sctx->funcname);
676: prhs[4] = sctx->ctx;
677: mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESVIRedundancyCheckInternal");
678: mxGetScalar(plhs[0]);
679: mxDestroyArray(prhs[0]);
680: mxDestroyArray(prhs[1]);
681: mxDestroyArray(prhs[2]);
682: mxDestroyArray(prhs[3]);
683: mxDestroyArray(plhs[0]);
684: return(0);
685: }
689: PetscErrorCode SNESVISetRedundancyCheckMatlab(SNES snes,const char *func,mxArray *ctx)
690: {
691: PetscErrorCode ierr;
692: SNESMatlabContext *sctx;
695: /* currently sctx is memory bleed */
696: PetscMalloc(sizeof(SNESMatlabContext),&sctx);
697: PetscStrallocpy(func,&sctx->funcname);
698: sctx->ctx = mxDuplicateArray(ctx);
699: SNESVISetRedundancyCheck(snes,SNESVIRedundancyCheck_Matlab,sctx);
700: return(0);
701: }
703: #endif
705: /* -------------------------------------------------------------------------- */
706: /*
707: SNESSetUp_VINEWTONRSLS - Sets up the internal data structures for the later use
708: of the SNESVI nonlinear solver.
710: Input Parameter:
711: . snes - the SNES context
712: . x - the solution vector
714: Application Interface Routine: SNESSetUp()
716: Notes:
717: For basic use of the SNES solvers, the user need not explicitly call
718: SNESSetUp(), since these actions will automatically occur during
719: the call to SNESSolve().
720: */
723: PetscErrorCode SNESSetUp_VINEWTONRSLS(SNES snes)
724: {
725: PetscErrorCode ierr;
726: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*) snes->data;
727: PetscInt *indices;
728: PetscInt i,n,rstart,rend;
729: SNESLineSearch linesearch;
732: SNESSetUp_VI(snes);
734: /* Set up previous active index set for the first snes solve
735: vi->IS_inact_prev = 0,1,2,....N */
737: VecGetOwnershipRange(snes->vec_sol,&rstart,&rend);
738: VecGetLocalSize(snes->vec_sol,&n);
739: PetscMalloc(n*sizeof(PetscInt),&indices);
740: for (i=0; i < n; i++) indices[i] = rstart + i;
741: ISCreateGeneral(PetscObjectComm((PetscObject)snes),n,indices,PETSC_OWN_POINTER,&vi->IS_inact_prev);
743: /* set the line search functions */
744: if (!snes->linesearch) {
745: SNESGetLineSearch(snes, &linesearch);
746: SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
747: }
748: return(0);
749: }
750: /* -------------------------------------------------------------------------- */
753: PetscErrorCode SNESReset_VINEWTONRSLS(SNES snes)
754: {
755: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*) snes->data;
756: PetscErrorCode ierr;
759: SNESReset_VI(snes);
760: ISDestroy(&vi->IS_inact_prev);
761: return(0);
762: }
764: /* -------------------------------------------------------------------------- */
765: /*MC
766: SNESVINEWTONRSLS - Reduced space active set solvers for variational inequalities based on Newton's method
768: Options Database:
769: + -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
770: - -snes_vi_monitor - prints the number of active constraints at each iteration.
772: Level: beginner
774: References:
775: - T. S. Munson, and S. Benson. Flexible Complementarity Solvers for Large-Scale
776: Applications, Optimization Methods and Software, 21 (2006).
778: .seealso: SNESVISetVariableBounds(), SNESVISetComputeVariableBounds(), SNESCreate(), SNES, SNESSetType(), SNESVINEWTONRSLS, SNESVINEWTONSSLS, SNESNEWTONTR, SNESLineSearchSet(),
779: SNESLineSearchSetPostCheck(), SNESLineSearchNo(), SNESLineSearchCubic(), SNESLineSearchQuadratic(),
780: SNESLineSearchSet(), SNESLineSearchNoNorms(), SNESLineSearchSetPreCheck(), SNESLineSearchSetParams(), SNESLineSearchGetParams()
782: M*/
785: PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONRSLS(SNES snes)
786: {
787: PetscErrorCode ierr;
788: SNES_VINEWTONRSLS *vi;
791: snes->ops->reset = SNESReset_VINEWTONRSLS;
792: snes->ops->setup = SNESSetUp_VINEWTONRSLS;
793: snes->ops->solve = SNESSolve_VINEWTONRSLS;
794: snes->ops->destroy = SNESDestroy_VI;
795: snes->ops->setfromoptions = SNESSetFromOptions_VI;
796: snes->ops->view = NULL;
797: snes->ops->converged = SNESConvergedDefault_VI;
799: snes->usesksp = PETSC_TRUE;
800: snes->usespc = PETSC_FALSE;
802: PetscNewLog(snes,SNES_VINEWTONRSLS,&vi);
803: snes->data = (void*)vi;
804: vi->checkredundancy = NULL;
806: PetscObjectComposeFunction((PetscObject)snes,"SNESVISetVariableBounds_C",SNESVISetVariableBounds_VI);
807: PetscObjectComposeFunction((PetscObject)snes,"SNESVISetComputeVariableBounds_C",SNESVISetComputeVariableBounds_VI);
808: return(0);
809: }