Actual source code: virs.c

petsc-3.9.0 2018-04-07
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 (*getinjection)(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: static PetscErrorCode DMHasCreateInjection_SNESVI(DM dm, PetscBool *flg)
 65: {
 69:   *flg = PETSC_FALSE;
 70:   return(0);
 71: }

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

 76: */
 77: PetscErrorCode  DMCreateInterpolation_SNESVI(DM dm1,DM dm2,Mat *mat,Vec *vec)
 78: {
 80:   PetscContainer isnes;
 81:   DM_SNESVI      *dmsnesvi1,*dmsnesvi2;
 82:   Mat            interp;

 85:   PetscObjectQuery((PetscObject)dm1,"VI",(PetscObject*)&isnes);
 86:   if (!isnes) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_PLIB,"Composed VI data structure is missing");
 87:   PetscContainerGetPointer(isnes,(void**)&dmsnesvi1);
 88:   PetscObjectQuery((PetscObject)dm2,"VI",(PetscObject*)&isnes);
 89:   if (!isnes) SETERRQ(PetscObjectComm((PetscObject)dm2),PETSC_ERR_PLIB,"Composed VI data structure is missing");
 90:   PetscContainerGetPointer(isnes,(void**)&dmsnesvi2);

 92:   (*dmsnesvi1->createinterpolation)(dm1,dm2,&interp,NULL);
 93:   MatCreateSubMatrix(interp,dmsnesvi2->inactive,dmsnesvi1->inactive,MAT_INITIAL_MATRIX,mat);
 94:   MatDestroy(&interp);
 95:   *vec = 0;
 96:   return(0);
 97: }

 99: extern PetscErrorCode  DMSetVI(DM,IS);

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

104: */
105: PetscErrorCode  DMCoarsen_SNESVI(DM dm1,MPI_Comm comm,DM *dm2)
106: {
108:   PetscContainer isnes;
109:   DM_SNESVI      *dmsnesvi1;
110:   Vec            finemarked,coarsemarked;
111:   IS             inactive;
112:   Mat            inject;
113:   const PetscInt *index;
114:   PetscInt       n,k,cnt = 0,rstart,*coarseindex;
115:   PetscScalar    *marked;

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

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

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

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

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

134:   DMCreateGlobalVector(dm1,&finemarked);
135:   DMCreateGlobalVector(*dm2,&coarsemarked);

137:   /*
138:      fill finemarked with locations of inactive points
139:   */
140:   ISGetIndices(dmsnesvi1->inactive,&index);
141:   ISGetLocalSize(dmsnesvi1->inactive,&n);
142:   VecSet(finemarked,0.0);
143:   for (k=0; k<n; k++) {
144:     VecSetValue(finemarked,index[k],1.0,INSERT_VALUES);
145:   }
146:   VecAssemblyBegin(finemarked);
147:   VecAssemblyEnd(finemarked);

149:   DMCreateInjection(*dm2,dm1,&inject);
150:   MatRestrict(inject,finemarked,coarsemarked);
151:   MatDestroy(&inject);

153:   /*
154:      create index set list of coarse inactive points from coarsemarked
155:   */
156:   VecGetLocalSize(coarsemarked,&n);
157:   VecGetOwnershipRange(coarsemarked,&rstart,NULL);
158:   VecGetArray(coarsemarked,&marked);
159:   for (k=0; k<n; k++) {
160:     if (marked[k] != 0.0) cnt++;
161:   }
162:   PetscMalloc1(cnt,&coarseindex);
163:   cnt  = 0;
164:   for (k=0; k<n; k++) {
165:     if (marked[k] != 0.0) coarseindex[cnt++] = k + rstart;
166:   }
167:   VecRestoreArray(coarsemarked,&marked);
168:   ISCreateGeneral(PetscObjectComm((PetscObject)coarsemarked),cnt,coarseindex,PETSC_OWN_POINTER,&inactive);

170:   DMClearGlobalVectors(dm1);

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

174:   DMSetVI(*dm2,inactive);

176:   VecDestroy(&finemarked);
177:   VecDestroy(&coarsemarked);
178:   ISDestroy(&inactive);
179:   return(0);
180: }

182: PetscErrorCode DMDestroy_SNESVI(DM_SNESVI *dmsnesvi)
183: {

187:   /* reset the base methods in the DM object that were changed when the DM_SNESVI was reset */
188:   dmsnesvi->dm->ops->createinterpolation = dmsnesvi->createinterpolation;
189:   dmsnesvi->dm->ops->coarsen             = dmsnesvi->coarsen;
190:   dmsnesvi->dm->ops->createglobalvector  = dmsnesvi->createglobalvector;
191:   dmsnesvi->dm->ops->getinjection        = dmsnesvi->getinjection;
192:   dmsnesvi->dm->ops->hascreateinjection  = dmsnesvi->hascreateinjection;
193:   /* need to clear out this vectors because some of them may not have a reference to the DM
194:     but they are counted as having references to the DM in DMDestroy() */
195:   DMClearGlobalVectors(dmsnesvi->dm);

197:   ISDestroy(&dmsnesvi->inactive);
198:   PetscFree(dmsnesvi);
199:   return(0);
200: }

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

206: */
207: PetscErrorCode  DMSetVI(DM dm,IS inactive)
208: {
210:   PetscContainer isnes;
211:   DM_SNESVI      *dmsnesvi;

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

216:   PetscObjectReference((PetscObject)inactive);

218:   PetscObjectQuery((PetscObject)dm,"VI",(PetscObject*)&isnes);
219:   if (!isnes) {
220:     PetscContainerCreate(PetscObjectComm((PetscObject)dm),&isnes);
221:     PetscContainerSetUserDestroy(isnes,(PetscErrorCode (*)(void*))DMDestroy_SNESVI);
222:     PetscNew(&dmsnesvi);
223:     PetscContainerSetPointer(isnes,(void*)dmsnesvi);
224:     PetscObjectCompose((PetscObject)dm,"VI",(PetscObject)isnes);
225:     PetscContainerDestroy(&isnes);

227:     dmsnesvi->createinterpolation = dm->ops->createinterpolation;
228:     dm->ops->createinterpolation  = DMCreateInterpolation_SNESVI;
229:     dmsnesvi->coarsen             = dm->ops->coarsen;
230:     dm->ops->coarsen              = DMCoarsen_SNESVI;
231:     dmsnesvi->createglobalvector  = dm->ops->createglobalvector;
232:     dm->ops->createglobalvector   = DMCreateGlobalVector_SNESVI;
233:     dmsnesvi->getinjection        = dm->ops->getinjection;
234:     dm->ops->getinjection         = NULL;
235:     dmsnesvi->hascreateinjection  = dm->ops->hascreateinjection;
236:     dm->ops->hascreateinjection   = DMHasCreateInjection_SNESVI;
237:   } else {
238:     PetscContainerGetPointer(isnes,(void**)&dmsnesvi);
239:     ISDestroy(&dmsnesvi->inactive);
240:   }
241:   DMClearGlobalVectors(dm);
242:   ISGetLocalSize(inactive,&dmsnesvi->n);

244:   dmsnesvi->inactive = inactive;
245:   dmsnesvi->dm       = dm;
246:   return(0);
247: }

249: /*
250:      DMDestroyVI - Frees the DM_SNESVI object contained in the DM
251:          - also resets the function pointers in the DM for createinterpolation() etc to use the original DM
252: */
253: PetscErrorCode  DMDestroyVI(DM dm)
254: {

258:   if (!dm) return(0);
259:   PetscObjectCompose((PetscObject)dm,"VI",(PetscObject)NULL);
260:   return(0);
261: }

263: /* --------------------------------------------------------------------------------------------------------*/


266: PetscErrorCode SNESCreateIndexSets_VINEWTONRSLS(SNES snes,Vec X,Vec F,IS *ISact,IS *ISinact)
267: {

271:   SNESVIGetActiveSetIS(snes,X,F,ISact);
272:   ISComplement(*ISact,X->map->rstart,X->map->rend,ISinact);
273:   return(0);
274: }

276: /* Create active and inactive set vectors. The local size of this vector is set and petsc computes the global size */
277: PetscErrorCode SNESCreateSubVectors_VINEWTONRSLS(SNES snes,PetscInt n,Vec *newv)
278: {
280:   Vec            v;

283:   VecCreate(PetscObjectComm((PetscObject)snes),&v);
284:   VecSetSizes(v,n,PETSC_DECIDE);
285:   VecSetType(v,VECSTANDARD);
286:   *newv = v;
287:   return(0);
288: }

290: /* Resets the snes PC and KSP when the active set sizes change */
291: PetscErrorCode SNESVIResetPCandKSP(SNES snes,Mat Amat,Mat Pmat)
292: {
294:   KSP            snesksp;

297:   SNESGetKSP(snes,&snesksp);
298:   KSPReset(snesksp);

300:   /*
301:   KSP                    kspnew;
302:   PC                     pcnew;
303:   MatSolverType          stype;


306:   KSPCreate(PetscObjectComm((PetscObject)snes),&kspnew);
307:   kspnew->pc_side = snesksp->pc_side;
308:   kspnew->rtol    = snesksp->rtol;
309:   kspnew->abstol    = snesksp->abstol;
310:   kspnew->max_it  = snesksp->max_it;
311:   KSPSetType(kspnew,((PetscObject)snesksp)->type_name);
312:   KSPGetPC(kspnew,&pcnew);
313:   PCSetType(kspnew->pc,((PetscObject)snesksp->pc)->type_name);
314:   PCSetOperators(kspnew->pc,Amat,Pmat);
315:   PCFactorGetMatSolverType(snesksp->pc,&stype);
316:   PCFactorSetMatSolverType(kspnew->pc,stype);
317:   KSPDestroy(&snesksp);
318:   snes->ksp = kspnew;
319:   PetscLogObjectParent((PetscObject)snes,(PetscObject)kspnew);
320:    KSPSetFromOptions(kspnew);*/
321:   return(0);
322: }

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

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

345:   snes->numFailures            = 0;
346:   snes->numLinearSolveFailures = 0;
347:   snes->reason                 = SNES_CONVERGED_ITERATING;

349:   maxits = snes->max_its;               /* maximum number of iterations */
350:   X      = snes->vec_sol;               /* solution vector */
351:   F      = snes->vec_func;              /* residual vector */
352:   Y      = snes->work[0];               /* work vectors */

354:   SNESLineSearchSetVIFunctions(snes->linesearch, SNESVIProjectOntoBounds, SNESVIComputeInactiveSetFnorm);
355:   SNESLineSearchSetVecs(snes->linesearch, X, NULL, NULL, NULL, NULL);
356:   SNESLineSearchSetUp(snes->linesearch);

358:   PetscObjectSAWsTakeAccess((PetscObject)snes);
359:   snes->iter = 0;
360:   snes->norm = 0.0;
361:   PetscObjectSAWsGrantAccess((PetscObject)snes);

363:   SNESVIProjectOntoBounds(snes,X);
364:   SNESComputeFunction(snes,X,F);
365:   SNESVIComputeInactiveSetFnorm(snes,F,X,&fnorm);
366:   VecNorm(X,NORM_2,&xnorm);        /* xnorm <- ||x||  */
367:   SNESCheckFunctionNorm(snes,fnorm);
368:   PetscObjectSAWsTakeAccess((PetscObject)snes);
369:   snes->norm = fnorm;
370:   PetscObjectSAWsGrantAccess((PetscObject)snes);
371:   SNESLogConvergenceHistory(snes,fnorm,0);
372:   SNESMonitor(snes,0,fnorm);

374:   /* test convergence */
375:   (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
376:   if (snes->reason) return(0);


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

381:     IS         IS_act; /* _act -> active set _inact -> inactive set */
382:     IS         IS_redact; /* redundant active set */
383:     VecScatter scat_act,scat_inact;
384:     PetscInt   nis_act,nis_inact;
385:     Vec        Y_act,Y_inact,F_inact;
386:     Mat        jac_inact_inact,prejac_inact_inact;
387:     PetscBool  isequal;

389:     /* Call general purpose update function */
390:     if (snes->ops->update) {
391:       (*snes->ops->update)(snes, snes->iter);
392:     }
393:     SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);


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

398:     /*original
399:     SNESVICreateIndexSets_RS(snes,X,F,&IS_act,&vi->IS_inact);
400:      */
401:     SNESVIGetActiveSetIS(snes,X,F,&IS_act);

403:     if (vi->checkredundancy) {
404:       (*vi->checkredundancy)(snes,IS_act,&IS_redact,vi->ctxP);
405:       if (IS_redact) {
406:         ISSort(IS_redact);
407:         ISComplement(IS_redact,X->map->rstart,X->map->rend,&vi->IS_inact);
408:         ISDestroy(&IS_redact);
409:       } else {
410:         ISComplement(IS_act,X->map->rstart,X->map->rend,&vi->IS_inact);
411:       }
412:     } else {
413:       ISComplement(IS_act,X->map->rstart,X->map->rend,&vi->IS_inact);
414:     }


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

420:     if (0) {                    /* Dead code (temporary developer hack) */
421:       IS keptrows;
422:       MatFindNonzeroRows(jac_inact_inact,&keptrows);
423:       if (keptrows) {
424:         PetscInt       cnt,*nrows,k;
425:         const PetscInt *krows,*inact;
426:         PetscInt       rstart;

428:         MatGetOwnershipRange(jac_inact_inact,&rstart,NULL);
429:         MatDestroy(&jac_inact_inact);
430:         ISDestroy(&IS_act);

432:         ISGetLocalSize(keptrows,&cnt);
433:         ISGetIndices(keptrows,&krows);
434:         ISGetIndices(vi->IS_inact,&inact);
435:         PetscMalloc1(cnt,&nrows);
436:         for (k=0; k<cnt; k++) nrows[k] = inact[krows[k]-rstart];
437:         ISRestoreIndices(keptrows,&krows);
438:         ISRestoreIndices(vi->IS_inact,&inact);
439:         ISDestroy(&keptrows);
440:         ISDestroy(&vi->IS_inact);

442:         ISCreateGeneral(PetscObjectComm((PetscObject)snes),cnt,nrows,PETSC_OWN_POINTER,&vi->IS_inact);
443:         ISComplement(vi->IS_inact,F->map->rstart,F->map->rend,&IS_act);
444:         MatCreateSubMatrix(snes->jacobian,vi->IS_inact,vi->IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);
445:       }
446:     }
447:     DMSetVI(snes->dm,vi->IS_inact);
448:     /* remove later */

450:     /*
451:     VecView(vi->xu,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xu))->comm));
452:     VecView(vi->xl,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xl))->comm));
453:     VecView(X,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)X)));
454:     VecView(F,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)F)));
455:     ISView(vi->IS_inact,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)vi->IS_inact)));
456:      */

458:     /* Get sizes of active and inactive sets */
459:     ISGetLocalSize(IS_act,&nis_act);
460:     ISGetLocalSize(vi->IS_inact,&nis_inact);

462:     /* Create active and inactive set vectors */
463:     SNESCreateSubVectors_VINEWTONRSLS(snes,nis_inact,&F_inact);
464:     SNESCreateSubVectors_VINEWTONRSLS(snes,nis_act,&Y_act);
465:     SNESCreateSubVectors_VINEWTONRSLS(snes,nis_inact,&Y_inact);

467:     /* Create scatter contexts */
468:     VecScatterCreate(Y,IS_act,Y_act,NULL,&scat_act);
469:     VecScatterCreate(Y,vi->IS_inact,Y_inact,NULL,&scat_inact);

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

475:     VecScatterBegin(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);
476:     VecScatterEnd(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);

478:     VecScatterBegin(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);
479:     VecScatterEnd(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);

481:     /* Active set direction = 0 */
482:     VecSet(Y_act,0);
483:     if (snes->jacobian != snes->jacobian_pre) {
484:       MatCreateSubMatrix(snes->jacobian_pre,vi->IS_inact,vi->IS_inact,MAT_INITIAL_MATRIX,&prejac_inact_inact);
485:     } else prejac_inact_inact = jac_inact_inact;

487:     ISEqual(vi->IS_inact_prev,vi->IS_inact,&isequal);
488:     if (!isequal) {
489:       SNESVIResetPCandKSP(snes,jac_inact_inact,prejac_inact_inact);
490:       PCFieldSplitRestrictIS(pc,vi->IS_inact);
491:     }

493:     /*      ISView(vi->IS_inact,0); */
494:     /*      ISView(IS_act,0);*/
495:     /*      MatView(snes->jacobian_pre,0); */



499:     KSPSetOperators(snes->ksp,jac_inact_inact,prejac_inact_inact);
500:     KSPSetUp(snes->ksp);
501:     {
502:       PC        pc;
503:       PetscBool flg;
504:       KSPGetPC(snes->ksp,&pc);
505:       PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&flg);
506:       if (flg) {
507:         KSP *subksps;
508:         PCFieldSplitGetSubKSP(pc,NULL,&subksps);
509:         KSPGetPC(subksps[0],&pc);
510:         PetscFree(subksps);
511:         PetscObjectTypeCompare((PetscObject)pc,PCBJACOBI,&flg);
512:         if (flg) {
513:           PetscInt       n,N = 101*101,j,cnts[3] = {0,0,0};
514:           const PetscInt *ii;

516:           ISGetSize(vi->IS_inact,&n);
517:           ISGetIndices(vi->IS_inact,&ii);
518:           for (j=0; j<n; j++) {
519:             if (ii[j] < N) cnts[0]++;
520:             else if (ii[j] < 2*N) cnts[1]++;
521:             else if (ii[j] < 3*N) cnts[2]++;
522:           }
523:           ISRestoreIndices(vi->IS_inact,&ii);

525:           PCBJacobiSetTotalBlocks(pc,3,cnts);
526:         }
527:       }
528:     }

530:     KSPSolve(snes->ksp,F_inact,Y_inact);
531:     VecScatterBegin(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);
532:     VecScatterEnd(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);
533:     VecScatterBegin(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);
534:     VecScatterEnd(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);

536:     VecDestroy(&F_inact);
537:     VecDestroy(&Y_act);
538:     VecDestroy(&Y_inact);
539:     VecScatterDestroy(&scat_act);
540:     VecScatterDestroy(&scat_inact);
541:     ISDestroy(&IS_act);
542:     if (!isequal) {
543:       ISDestroy(&vi->IS_inact_prev);
544:       ISDuplicate(vi->IS_inact,&vi->IS_inact_prev);
545:     }
546:     ISDestroy(&vi->IS_inact);
547:     MatDestroy(&jac_inact_inact);
548:     if (snes->jacobian != snes->jacobian_pre) {
549:       MatDestroy(&prejac_inact_inact);
550:     }

552:     KSPGetConvergedReason(snes->ksp,&kspreason);
553:     if (kspreason < 0) {
554:       if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
555:         PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);
556:         snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
557:         break;
558:       }
559:     }

561:     KSPGetIterationNumber(snes->ksp,&lits);
562:     snes->linear_its += lits;
563:     PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);
564:     /*
565:     if (snes->ops->precheck) {
566:       PetscBool changed_y = PETSC_FALSE;
567:       (*snes->ops->precheck)(snes,X,Y,snes->precheck,&changed_y);
568:     }

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

623: PetscErrorCode SNESVISetRedundancyCheck(SNES snes,PetscErrorCode (*func)(SNES,IS,IS*,void*),void *ctx)
624: {
625:   SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data;

629:   vi->checkredundancy = func;
630:   vi->ctxP            = ctx;
631:   return(0);
632: }

634: #if defined(PETSC_HAVE_MATLAB_ENGINE)
635: #include <engine.h>
636: #include <mex.h>
637: typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext;

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


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

676: PetscErrorCode SNESVISetRedundancyCheckMatlab(SNES snes,const char *func,mxArray *ctx)
677: {
678:   PetscErrorCode    ierr;
679:   SNESMatlabContext *sctx;

682:   /* currently sctx is memory bleed */
683:   PetscNew(&sctx);
684:   PetscStrallocpy(func,&sctx->funcname);
685:   sctx->ctx = mxDuplicateArray(ctx);
686:   SNESVISetRedundancyCheck(snes,SNESVIRedundancyCheck_Matlab,sctx);
687:   return(0);
688: }

690: #endif

692: /* -------------------------------------------------------------------------- */
693: /*
694:    SNESSetUp_VINEWTONRSLS - Sets up the internal data structures for the later use
695:    of the SNESVI nonlinear solver.

697:    Input Parameter:
698: .  snes - the SNES context

700:    Application Interface Routine: SNESSetUp()

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

716:   SNESSetUp_VI(snes);

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

721:   VecGetOwnershipRange(snes->vec_sol,&rstart,&rend);
722:   VecGetLocalSize(snes->vec_sol,&n);
723:   PetscMalloc1(n,&indices);
724:   for (i=0; i < n; i++) indices[i] = rstart + i;
725:   ISCreateGeneral(PetscObjectComm((PetscObject)snes),n,indices,PETSC_OWN_POINTER,&vi->IS_inact_prev);

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

741:   SNESReset_VI(snes);
742:   ISDestroy(&vi->IS_inact_prev);
743:   return(0);
744: }

746: /* -------------------------------------------------------------------------- */
747: /*MC
748:       SNESVINEWTONRSLS - Reduced space active set solvers for variational inequalities based on Newton's method

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

754:    Level: beginner

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

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

762: M*/
763: PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONRSLS(SNES snes)
764: {
765:   PetscErrorCode    ierr;
766:   SNES_VINEWTONRSLS *vi;

769:   snes->ops->reset          = SNESReset_VINEWTONRSLS;
770:   snes->ops->setup          = SNESSetUp_VINEWTONRSLS;
771:   snes->ops->solve          = SNESSolve_VINEWTONRSLS;
772:   snes->ops->destroy        = SNESDestroy_VI;
773:   snes->ops->setfromoptions = SNESSetFromOptions_VI;
774:   snes->ops->view           = NULL;
775:   snes->ops->converged      = SNESConvergedDefault_VI;

777:   snes->usesksp = PETSC_TRUE;
778:   snes->usesnpc = PETSC_FALSE;

780:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

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

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