Actual source code: virs.c

petsc-3.13.6 2020-09-29
Report Typos and Errors

  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:   return(0);
 62: }

 64: /*
 65:      DMCreateInterpolation_SNESVI - Modifieds the interpolation obtained from the DM by removing all rows and columns associated with active constraints.

 67: */
 68: PetscErrorCode  DMCreateInterpolation_SNESVI(DM dm1,DM dm2,Mat *mat,Vec *vec)
 69: {
 71:   PetscContainer isnes;
 72:   DM_SNESVI      *dmsnesvi1,*dmsnesvi2;
 73:   Mat            interp;

 76:   PetscObjectQuery((PetscObject)dm1,"VI",(PetscObject*)&isnes);
 77:   if (!isnes) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_PLIB,"Composed VI data structure is missing");
 78:   PetscContainerGetPointer(isnes,(void**)&dmsnesvi1);
 79:   PetscObjectQuery((PetscObject)dm2,"VI",(PetscObject*)&isnes);
 80:   if (!isnes) SETERRQ(PetscObjectComm((PetscObject)dm2),PETSC_ERR_PLIB,"Composed VI data structure is missing");
 81:   PetscContainerGetPointer(isnes,(void**)&dmsnesvi2);

 83:   (*dmsnesvi1->createinterpolation)(dm1,dm2,&interp,NULL);
 84:   MatCreateSubMatrix(interp,dmsnesvi2->inactive,dmsnesvi1->inactive,MAT_INITIAL_MATRIX,mat);
 85:   MatDestroy(&interp);
 86:   *vec = 0;
 87:   return(0);
 88: }

 90: static PetscErrorCode DMSetVI(DM,IS);
 91: static PetscErrorCode DMDestroyVI(DM);

 93: /*
 94:      DMCoarsen_SNESVI - Computes the regular coarsened DM then computes additional information about its inactive set

 96: */
 97: PetscErrorCode  DMCoarsen_SNESVI(DM dm1,MPI_Comm comm,DM *dm2)
 98: {
100:   PetscContainer isnes;
101:   DM_SNESVI      *dmsnesvi1;
102:   Vec            finemarked,coarsemarked;
103:   IS             inactive;
104:   Mat            inject;
105:   const PetscInt *index;
106:   PetscInt       n,k,cnt = 0,rstart,*coarseindex;
107:   PetscScalar    *marked;

110:   PetscObjectQuery((PetscObject)dm1,"VI",(PetscObject*)&isnes);
111:   if (!isnes) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_PLIB,"Composed VI data structure is missing");
112:   PetscContainerGetPointer(isnes,(void**)&dmsnesvi1);

114:   /* get the original coarsen */
115:   (*dmsnesvi1->coarsen)(dm1,comm,dm2);

117:   /* not sure why this extra reference is needed, but without the dm2 disappears too early */
118:   /* Updating the KSPCreateVecs() to avoid using DMGetGlobalVector() when matrix is available removes the need for this reference? */
119:   /*  PetscObjectReference((PetscObject)*dm2);*/

121:   /* need to set back global vectors in order to use the original injection */
122:   DMClearGlobalVectors(dm1);

124:   dm1->ops->createglobalvector = dmsnesvi1->createglobalvector;

126:   DMCreateGlobalVector(dm1,&finemarked);
127:   DMCreateGlobalVector(*dm2,&coarsemarked);

129:   /*
130:      fill finemarked with locations of inactive points
131:   */
132:   ISGetIndices(dmsnesvi1->inactive,&index);
133:   ISGetLocalSize(dmsnesvi1->inactive,&n);
134:   VecSet(finemarked,0.0);
135:   for (k=0; k<n; k++) {
136:     VecSetValue(finemarked,index[k],1.0,INSERT_VALUES);
137:   }
138:   VecAssemblyBegin(finemarked);
139:   VecAssemblyEnd(finemarked);

141:   DMCreateInjection(*dm2,dm1,&inject);
142:   MatRestrict(inject,finemarked,coarsemarked);
143:   MatDestroy(&inject);

145:   /*
146:      create index set list of coarse inactive points from coarsemarked
147:   */
148:   VecGetLocalSize(coarsemarked,&n);
149:   VecGetOwnershipRange(coarsemarked,&rstart,NULL);
150:   VecGetArray(coarsemarked,&marked);
151:   for (k=0; k<n; k++) {
152:     if (marked[k] != 0.0) cnt++;
153:   }
154:   PetscMalloc1(cnt,&coarseindex);
155:   cnt  = 0;
156:   for (k=0; k<n; k++) {
157:     if (marked[k] != 0.0) coarseindex[cnt++] = k + rstart;
158:   }
159:   VecRestoreArray(coarsemarked,&marked);
160:   ISCreateGeneral(PetscObjectComm((PetscObject)coarsemarked),cnt,coarseindex,PETSC_OWN_POINTER,&inactive);

162:   DMClearGlobalVectors(dm1);

164:   dm1->ops->createglobalvector = DMCreateGlobalVector_SNESVI;

166:   DMSetVI(*dm2,inactive);

168:   VecDestroy(&finemarked);
169:   VecDestroy(&coarsemarked);
170:   ISDestroy(&inactive);
171:   return(0);
172: }

174: PetscErrorCode DMDestroy_SNESVI(DM_SNESVI *dmsnesvi)
175: {

179:   /* reset the base methods in the DM object that were changed when the DM_SNESVI was reset */
180:   dmsnesvi->dm->ops->createinterpolation = dmsnesvi->createinterpolation;
181:   dmsnesvi->dm->ops->coarsen             = dmsnesvi->coarsen;
182:   dmsnesvi->dm->ops->createglobalvector  = dmsnesvi->createglobalvector;
183:   dmsnesvi->dm->ops->createinjection     = dmsnesvi->createinjection;
184:   dmsnesvi->dm->ops->hascreateinjection  = dmsnesvi->hascreateinjection;
185:   /* need to clear out this vectors because some of them may not have a reference to the DM
186:     but they are counted as having references to the DM in DMDestroy() */
187:   DMClearGlobalVectors(dmsnesvi->dm);

189:   ISDestroy(&dmsnesvi->inactive);
190:   PetscFree(dmsnesvi);
191:   return(0);
192: }

194: /*
195:      DMSetVI - Marks a DM as associated with a VI problem. This causes the interpolation/restriction operators to
196:                be restricted to only those variables NOT associated with active constraints.

198: */
199: static PetscErrorCode DMSetVI(DM dm,IS inactive)
200: {
202:   PetscContainer isnes;
203:   DM_SNESVI      *dmsnesvi;

206:   if (!dm) return(0);

208:   PetscObjectReference((PetscObject)inactive);

210:   PetscObjectQuery((PetscObject)dm,"VI",(PetscObject*)&isnes);
211:   if (!isnes) {
212:     PetscContainerCreate(PetscObjectComm((PetscObject)dm),&isnes);
213:     PetscContainerSetUserDestroy(isnes,(PetscErrorCode (*)(void*))DMDestroy_SNESVI);
214:     PetscNew(&dmsnesvi);
215:     PetscContainerSetPointer(isnes,(void*)dmsnesvi);
216:     PetscObjectCompose((PetscObject)dm,"VI",(PetscObject)isnes);
217:     PetscContainerDestroy(&isnes);

219:     dmsnesvi->createinterpolation = dm->ops->createinterpolation;
220:     dm->ops->createinterpolation  = DMCreateInterpolation_SNESVI;
221:     dmsnesvi->coarsen             = dm->ops->coarsen;
222:     dm->ops->coarsen              = DMCoarsen_SNESVI;
223:     dmsnesvi->createglobalvector  = dm->ops->createglobalvector;
224:     dm->ops->createglobalvector   = DMCreateGlobalVector_SNESVI;
225:     dmsnesvi->createinjection     = dm->ops->createinjection;
226:     dm->ops->createinjection      = NULL;
227:     dmsnesvi->hascreateinjection  = dm->ops->hascreateinjection;
228:     dm->ops->hascreateinjection   = NULL;
229:   } else {
230:     PetscContainerGetPointer(isnes,(void**)&dmsnesvi);
231:     ISDestroy(&dmsnesvi->inactive);
232:   }
233:   DMClearGlobalVectors(dm);
234:   ISGetLocalSize(inactive,&dmsnesvi->n);

236:   dmsnesvi->inactive = inactive;
237:   dmsnesvi->dm       = dm;
238:   return(0);
239: }

241: /*
242:      DMDestroyVI - Frees the DM_SNESVI object contained in the DM
243:          - also resets the function pointers in the DM for createinterpolation() etc to use the original DM
244: */
245: static PetscErrorCode DMDestroyVI(DM dm)
246: {

250:   if (!dm) return(0);
251:   PetscObjectCompose((PetscObject)dm,"VI",(PetscObject)NULL);
252:   return(0);
253: }

255: /* --------------------------------------------------------------------------------------------------------*/


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;


299:   KSPCreate(PetscObjectComm((PetscObject)snes),&kspnew);
300:   kspnew->pc_side = snesksp->pc_side;
301:   kspnew->rtol    = snesksp->rtol;
302:   kspnew->abstol    = snesksp->abstol;
303:   kspnew->max_it  = snesksp->max_it;
304:   KSPSetType(kspnew,((PetscObject)snesksp)->type_name);
305:   KSPGetPC(kspnew,&pcnew);
306:   PCSetType(kspnew->pc,((PetscObject)snesksp->pc)->type_name);
307:   PCSetOperators(kspnew->pc,Amat,Pmat);
308:   PCFactorGetMatSolverType(snesksp->pc,&stype);
309:   PCFactorSetMatSolverType(kspnew->pc,stype);
310:   KSPDestroy(&snesksp);
311:   snes->ksp = kspnew;
312:   PetscLogObjectParent((PetscObject)snes,(PetscObject)kspnew);
313:    KSPSetFromOptions(kspnew);*/
314:   return(0);
315: }

317: /* Variational Inequality solver using reduce space method. No semismooth algorithm is
318:    implemented in this algorithm. It basically identifies the active constraints and does
319:    a linear solve on the other variables (those not associated with the active constraints). */
320: PetscErrorCode SNESSolve_VINEWTONRSLS(SNES snes)
321: {
322:   SNES_VINEWTONRSLS    *vi = (SNES_VINEWTONRSLS*)snes->data;
323:   PetscErrorCode       ierr;
324:   PetscInt             maxits,i,lits;
325:   SNESLineSearchReason lssucceed;
326:   PetscReal            fnorm,gnorm,xnorm=0,ynorm;
327:   Vec                  Y,X,F;
328:   KSPConvergedReason   kspreason;
329:   KSP                  ksp;
330:   PC                   pc;

333:   /* Multigrid must use Galerkin for coarse grids with active set/reduced space methods; cannot rediscretize on coarser grids*/
334:   SNESGetKSP(snes,&ksp);
335:   KSPGetPC(ksp,&pc);
336:   PCMGSetGalerkin(pc,PC_MG_GALERKIN_BOTH);

338:   snes->numFailures            = 0;
339:   snes->numLinearSolveFailures = 0;
340:   snes->reason                 = SNES_CONVERGED_ITERATING;

342:   maxits = snes->max_its;               /* maximum number of iterations */
343:   X      = snes->vec_sol;               /* solution vector */
344:   F      = snes->vec_func;              /* residual vector */
345:   Y      = snes->work[0];               /* work vectors */

347:   SNESLineSearchSetVIFunctions(snes->linesearch, SNESVIProjectOntoBounds, SNESVIComputeInactiveSetFnorm);
348:   SNESLineSearchSetVecs(snes->linesearch, X, NULL, NULL, NULL, NULL);
349:   SNESLineSearchSetUp(snes->linesearch);

351:   PetscObjectSAWsTakeAccess((PetscObject)snes);
352:   snes->iter = 0;
353:   snes->norm = 0.0;
354:   PetscObjectSAWsGrantAccess((PetscObject)snes);

356:   SNESVIProjectOntoBounds(snes,X);
357:   SNESComputeFunction(snes,X,F);
358:   SNESVIComputeInactiveSetFnorm(snes,F,X,&fnorm);
359:   VecNorm(X,NORM_2,&xnorm);        /* xnorm <- ||x||  */
360:   SNESCheckFunctionNorm(snes,fnorm);
361:   PetscObjectSAWsTakeAccess((PetscObject)snes);
362:   snes->norm = fnorm;
363:   PetscObjectSAWsGrantAccess((PetscObject)snes);
364:   SNESLogConvergenceHistory(snes,fnorm,0);
365:   SNESMonitor(snes,0,fnorm);

367:   /* test convergence */
368:   (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
369:   if (snes->reason) return(0);


372:   for (i=0; i<maxits; i++) {

374:     IS         IS_act; /* _act -> active set _inact -> inactive set */
375:     IS         IS_redact; /* redundant active set */
376:     VecScatter scat_act,scat_inact;
377:     PetscInt   nis_act,nis_inact;
378:     Vec        Y_act,Y_inact,F_inact;
379:     Mat        jac_inact_inact,prejac_inact_inact;
380:     PetscBool  isequal;

382:     /* Call general purpose update function */
383:     if (snes->ops->update) {
384:       (*snes->ops->update)(snes, snes->iter);
385:     }
386:     SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
387:     SNESCheckJacobianDomainerror(snes);

389:     /* Create active and inactive index sets */

391:     /*original
392:     SNESVICreateIndexSets_RS(snes,X,F,&IS_act,&vi->IS_inact);
393:      */
394:     SNESVIGetActiveSetIS(snes,X,F,&IS_act);

396:     if (vi->checkredundancy) {
397:       (*vi->checkredundancy)(snes,IS_act,&IS_redact,vi->ctxP);
398:       if (IS_redact) {
399:         ISSort(IS_redact);
400:         ISComplement(IS_redact,X->map->rstart,X->map->rend,&vi->IS_inact);
401:         ISDestroy(&IS_redact);
402:       } else {
403:         ISComplement(IS_act,X->map->rstart,X->map->rend,&vi->IS_inact);
404:       }
405:     } else {
406:       ISComplement(IS_act,X->map->rstart,X->map->rend,&vi->IS_inact);
407:     }


410:     /* Create inactive set submatrix */
411:     MatCreateSubMatrix(snes->jacobian,vi->IS_inact,vi->IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);

413:     if (0) {                    /* Dead code (temporary developer hack) */
414:       IS keptrows;
415:       MatFindNonzeroRows(jac_inact_inact,&keptrows);
416:       if (keptrows) {
417:         PetscInt       cnt,*nrows,k;
418:         const PetscInt *krows,*inact;
419:         PetscInt       rstart;

421:         MatGetOwnershipRange(jac_inact_inact,&rstart,NULL);
422:         MatDestroy(&jac_inact_inact);
423:         ISDestroy(&IS_act);

425:         ISGetLocalSize(keptrows,&cnt);
426:         ISGetIndices(keptrows,&krows);
427:         ISGetIndices(vi->IS_inact,&inact);
428:         PetscMalloc1(cnt,&nrows);
429:         for (k=0; k<cnt; k++) nrows[k] = inact[krows[k]-rstart];
430:         ISRestoreIndices(keptrows,&krows);
431:         ISRestoreIndices(vi->IS_inact,&inact);
432:         ISDestroy(&keptrows);
433:         ISDestroy(&vi->IS_inact);

435:         ISCreateGeneral(PetscObjectComm((PetscObject)snes),cnt,nrows,PETSC_OWN_POINTER,&vi->IS_inact);
436:         ISComplement(vi->IS_inact,F->map->rstart,F->map->rend,&IS_act);
437:         MatCreateSubMatrix(snes->jacobian,vi->IS_inact,vi->IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);
438:       }
439:     }
440:     DMSetVI(snes->dm,vi->IS_inact);
441:     /* remove later */

443:     /*
444:     VecView(vi->xu,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xu))->comm));
445:     VecView(vi->xl,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xl))->comm));
446:     VecView(X,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)X)));
447:     VecView(F,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)F)));
448:     ISView(vi->IS_inact,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)vi->IS_inact)));
449:      */

451:     /* Get sizes of active and inactive sets */
452:     ISGetLocalSize(IS_act,&nis_act);
453:     ISGetLocalSize(vi->IS_inact,&nis_inact);

455:     /* Create active and inactive set vectors */
456:     SNESCreateSubVectors_VINEWTONRSLS(snes,nis_inact,&F_inact);
457:     SNESCreateSubVectors_VINEWTONRSLS(snes,nis_act,&Y_act);
458:     SNESCreateSubVectors_VINEWTONRSLS(snes,nis_inact,&Y_inact);

460:     /* Create scatter contexts */
461:     VecScatterCreate(Y,IS_act,Y_act,NULL,&scat_act);
462:     VecScatterCreate(Y,vi->IS_inact,Y_inact,NULL,&scat_inact);

464:     /* Do a vec scatter to active and inactive set vectors */
465:     VecScatterBegin(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);
466:     VecScatterEnd(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);

468:     VecScatterBegin(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);
469:     VecScatterEnd(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);

471:     VecScatterBegin(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);
472:     VecScatterEnd(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);

474:     /* Active set direction = 0 */
475:     VecSet(Y_act,0);
476:     if (snes->jacobian != snes->jacobian_pre) {
477:       MatCreateSubMatrix(snes->jacobian_pre,vi->IS_inact,vi->IS_inact,MAT_INITIAL_MATRIX,&prejac_inact_inact);
478:     } else prejac_inact_inact = jac_inact_inact;

480:     ISEqual(vi->IS_inact_prev,vi->IS_inact,&isequal);
481:     if (!isequal) {
482:       SNESVIResetPCandKSP(snes,jac_inact_inact,prejac_inact_inact);
483:       PCFieldSplitRestrictIS(pc,vi->IS_inact);
484:     }

486:     /*      ISView(vi->IS_inact,0); */
487:     /*      ISView(IS_act,0);*/
488:     /*      MatView(snes->jacobian_pre,0); */



492:     KSPSetOperators(snes->ksp,jac_inact_inact,prejac_inact_inact);
493:     KSPSetUp(snes->ksp);
494:     {
495:       PC        pc;
496:       PetscBool flg;
497:       KSPGetPC(snes->ksp,&pc);
498:       PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&flg);
499:       if (flg) {
500:         KSP *subksps;
501:         PCFieldSplitGetSubKSP(pc,NULL,&subksps);
502:         KSPGetPC(subksps[0],&pc);
503:         PetscFree(subksps);
504:         PetscObjectTypeCompare((PetscObject)pc,PCBJACOBI,&flg);
505:         if (flg) {
506:           PetscInt       n,N = 101*101,j,cnts[3] = {0,0,0};
507:           const PetscInt *ii;

509:           ISGetSize(vi->IS_inact,&n);
510:           ISGetIndices(vi->IS_inact,&ii);
511:           for (j=0; j<n; j++) {
512:             if (ii[j] < N) cnts[0]++;
513:             else if (ii[j] < 2*N) cnts[1]++;
514:             else if (ii[j] < 3*N) cnts[2]++;
515:           }
516:           ISRestoreIndices(vi->IS_inact,&ii);

518:           PCBJacobiSetTotalBlocks(pc,3,cnts);
519:         }
520:       }
521:     }

523:     KSPSolve(snes->ksp,F_inact,Y_inact);
524:     VecScatterBegin(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);
525:     VecScatterEnd(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);
526:     VecScatterBegin(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);
527:     VecScatterEnd(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);

529:     VecDestroy(&F_inact);
530:     VecDestroy(&Y_act);
531:     VecDestroy(&Y_inact);
532:     VecScatterDestroy(&scat_act);
533:     VecScatterDestroy(&scat_inact);
534:     ISDestroy(&IS_act);
535:     if (!isequal) {
536:       ISDestroy(&vi->IS_inact_prev);
537:       ISDuplicate(vi->IS_inact,&vi->IS_inact_prev);
538:     }
539:     ISDestroy(&vi->IS_inact);
540:     MatDestroy(&jac_inact_inact);
541:     if (snes->jacobian != snes->jacobian_pre) {
542:       MatDestroy(&prejac_inact_inact);
543:     }

545:     KSPGetConvergedReason(snes->ksp,&kspreason);
546:     if (kspreason < 0) {
547:       if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
548:         PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);
549:         snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
550:         break;
551:       }
552:     }

554:     KSPGetIterationNumber(snes->ksp,&lits);
555:     snes->linear_its += lits;
556:     PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);
557:     /*
558:     if (snes->ops->precheck) {
559:       PetscBool changed_y = PETSC_FALSE;
560:       (*snes->ops->precheck)(snes,X,Y,snes->precheck,&changed_y);
561:     }

563:     if (PetscLogPrintInfo) {
564:       SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);
565:     }
566:     */
567:     /* Compute a (scaled) negative update in the line search routine:
568:          Y <- X - lambda*Y
569:        and evaluate G = function(Y) (depends on the line search).
570:     */
571:     VecCopy(Y,snes->vec_sol_update);
572:     ynorm = 1; gnorm = fnorm;
573:     SNESLineSearchApply(snes->linesearch, X, F, &gnorm, Y);
574:     SNESLineSearchGetReason(snes->linesearch, &lssucceed);
575:     SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm);
576:     PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",(double)fnorm,(double)gnorm,(double)ynorm,(int)lssucceed);
577:     if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
578:     if (snes->domainerror) {
579:       snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
580:       DMDestroyVI(snes->dm);
581:       return(0);
582:     }
583:     if (lssucceed) {
584:       if (++snes->numFailures >= snes->maxFailures) {
585:         PetscBool ismin;
586:         snes->reason = SNES_DIVERGED_LINE_SEARCH;
587:         SNESVICheckLocalMin_Private(snes,snes->jacobian,F,X,gnorm,&ismin);
588:         if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
589:         break;
590:       }
591:    }
592:    DMDestroyVI(snes->dm);
593:     /* Update function and solution vectors */
594:     fnorm = gnorm;
595:     /* Monitor convergence */
596:     PetscObjectSAWsTakeAccess((PetscObject)snes);
597:     snes->iter = i+1;
598:     snes->norm = fnorm;
599:     snes->xnorm = xnorm;
600:     snes->ynorm = ynorm;
601:     PetscObjectSAWsGrantAccess((PetscObject)snes);
602:     SNESLogConvergenceHistory(snes,snes->norm,lits);
603:     SNESMonitor(snes,snes->iter,snes->norm);
604:     /* Test for convergence, xnorm = || X || */
605:     if (snes->ops->converged != SNESConvergedSkip) { VecNorm(X,NORM_2,&xnorm); }
606:     (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);
607:     if (snes->reason) break;
608:   }
609:   /* 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 */
610:   DMDestroyVI(snes->dm);
611:   if (i == maxits) {
612:     PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);
613:     if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
614:   }
615:   return(0);
616: }

618: PetscErrorCode SNESVISetRedundancyCheck(SNES snes,PetscErrorCode (*func)(SNES,IS,IS*,void*),void *ctx)
619: {
620:   SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data;

624:   vi->checkredundancy = func;
625:   vi->ctxP            = ctx;
626:   return(0);
627: }

629: #if defined(PETSC_HAVE_MATLAB_ENGINE)
630: #include <engine.h>
631: #include <mex.h>
632: typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext;

634: PetscErrorCode SNESVIRedundancyCheck_Matlab(SNES snes,IS is_act,IS *is_redact,void *ctx)
635: {
636:   PetscErrorCode    ierr;
637:   SNESMatlabContext *sctx = (SNESMatlabContext*)ctx;
638:   int               nlhs  = 1, nrhs = 5;
639:   mxArray           *plhs[1], *prhs[5];
640:   long long int     l1      = 0, l2 = 0, ls = 0;
641:   PetscInt          *indices=NULL;


649:   /* Create IS for reduced active set of size 0, its size and indices will
650:    bet set by the Matlab function */
651:   ISCreateGeneral(PetscObjectComm((PetscObject)snes),0,indices,PETSC_OWN_POINTER,is_redact);
652:   /* call Matlab function in ctx */
653:   PetscArraycpy(&ls,&snes,1);
654:   PetscArraycpy(&l1,&is_act,1);
655:   PetscArraycpy(&l2,is_redact,1);
656:   prhs[0] = mxCreateDoubleScalar((double)ls);
657:   prhs[1] = mxCreateDoubleScalar((double)l1);
658:   prhs[2] = mxCreateDoubleScalar((double)l2);
659:   prhs[3] = mxCreateString(sctx->funcname);
660:   prhs[4] = sctx->ctx;
661:   mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESVIRedundancyCheckInternal");
662:   mxGetScalar(plhs[0]);
663:   mxDestroyArray(prhs[0]);
664:   mxDestroyArray(prhs[1]);
665:   mxDestroyArray(prhs[2]);
666:   mxDestroyArray(prhs[3]);
667:   mxDestroyArray(plhs[0]);
668:   return(0);
669: }

671: PetscErrorCode SNESVISetRedundancyCheckMatlab(SNES snes,const char *func,mxArray *ctx)
672: {
673:   PetscErrorCode    ierr;
674:   SNESMatlabContext *sctx;

677:   /* currently sctx is memory bleed */
678:   PetscNew(&sctx);
679:   PetscStrallocpy(func,&sctx->funcname);
680:   sctx->ctx = mxDuplicateArray(ctx);
681:   SNESVISetRedundancyCheck(snes,SNESVIRedundancyCheck_Matlab,sctx);
682:   return(0);
683: }

685: #endif

687: /* -------------------------------------------------------------------------- */
688: /*
689:    SNESSetUp_VINEWTONRSLS - Sets up the internal data structures for the later use
690:    of the SNESVI nonlinear solver.

692:    Input Parameter:
693: .  snes - the SNES context

695:    Application Interface Routine: SNESSetUp()

697:    Notes:
698:    For basic use of the SNES solvers, the user need not explicitly call
699:    SNESSetUp(), since these actions will automatically occur during
700:    the call to SNESSolve().
701:  */
702: PetscErrorCode SNESSetUp_VINEWTONRSLS(SNES snes)
703: {
704:   PetscErrorCode    ierr;
705:   SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*) snes->data;
706:   PetscInt          *indices;
707:   PetscInt          i,n,rstart,rend;
708:   SNESLineSearch    linesearch;

711:   SNESSetUp_VI(snes);

713:   /* Set up previous active index set for the first snes solve
714:    vi->IS_inact_prev = 0,1,2,....N */

716:   VecGetOwnershipRange(snes->vec_sol,&rstart,&rend);
717:   VecGetLocalSize(snes->vec_sol,&n);
718:   PetscMalloc1(n,&indices);
719:   for (i=0; i < n; i++) indices[i] = rstart + i;
720:   ISCreateGeneral(PetscObjectComm((PetscObject)snes),n,indices,PETSC_OWN_POINTER,&vi->IS_inact_prev);

722:   /* set the line search functions */
723:   if (!snes->linesearch) {
724:     SNESGetLineSearch(snes, &linesearch);
725:     SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
726:   }
727:   return(0);
728: }
729: /* -------------------------------------------------------------------------- */
730: PetscErrorCode SNESReset_VINEWTONRSLS(SNES snes)
731: {
732:   SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*) snes->data;
733:   PetscErrorCode    ierr;

736:   SNESReset_VI(snes);
737:   ISDestroy(&vi->IS_inact_prev);
738:   return(0);
739: }

741: /* -------------------------------------------------------------------------- */
742: /*MC
743:       SNESVINEWTONRSLS - Reduced space active set solvers for variational inequalities based on Newton's method

745:    Options Database:
746: +   -snes_type <vinewtonssls,vinewtonrsls> - a semi-smooth solver, a reduced space active set method
747: -   -snes_vi_monitor - prints the number of active constraints at each iteration.

749:    Level: beginner

751:    References:
752: .  1. - T. S. Munson, and S. Benson. Flexible Complementarity Solvers for Large Scale
753:      Applications, Optimization Methods and Software, 21 (2006).

755: .seealso:  SNESVISetVariableBounds(), SNESVISetComputeVariableBounds(), SNESCreate(), SNES, SNESSetType(), SNESVINEWTONSSLS, SNESNEWTONTR, SNESLineSearchSetType(),SNESLineSearchSetPostCheck(), SNESLineSearchSetPreCheck()

757: M*/
758: PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONRSLS(SNES snes)
759: {
760:   PetscErrorCode    ierr;
761:   SNES_VINEWTONRSLS *vi;
762:   SNESLineSearch    linesearch;

765:   snes->ops->reset          = SNESReset_VINEWTONRSLS;
766:   snes->ops->setup          = SNESSetUp_VINEWTONRSLS;
767:   snes->ops->solve          = SNESSolve_VINEWTONRSLS;
768:   snes->ops->destroy        = SNESDestroy_VI;
769:   snes->ops->setfromoptions = SNESSetFromOptions_VI;
770:   snes->ops->view           = NULL;
771:   snes->ops->converged      = SNESConvergedDefault_VI;

773:   snes->usesksp = PETSC_TRUE;
774:   snes->usesnpc = PETSC_FALSE;

776:   SNESGetLineSearch(snes, &linesearch);
777:   if (!((PetscObject)linesearch)->type_name) {
778:     SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
779:   }
780:   SNESLineSearchBTSetAlpha(linesearch, 0.0);

782:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

784:   PetscNewLog(snes,&vi);
785:   snes->data          = (void*)vi;
786:   vi->checkredundancy = NULL;

788:   PetscObjectComposeFunction((PetscObject)snes,"SNESVISetVariableBounds_C",SNESVISetVariableBounds_VI);
789:   PetscObjectComposeFunction((PetscObject)snes,"SNESVISetComputeVariableBounds_C",SNESVISetComputeVariableBounds_VI);
790:   return(0);
791: }