Actual source code: virs.c

petsc-3.4.4 2014-03-13
  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: }