Actual source code: slp.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: */
 10: /*
 11:    SLEPc nonlinear eigensolver: "slp"

 13:    Method: Succesive linear problems

 15:    Algorithm:

 17:        Newton-type iteration based on first order Taylor approximation.

 19:    References:

 21:        [1] A. Ruhe, "Algorithms for the nonlinear eigenvalue problem", SIAM J.
 22:            Numer. Anal. 10(4):674-689, 1973.
 23: */

 25: #include <slepc/private/nepimpl.h>         /*I "slepcnep.h" I*/
 26: #include <../src/nep/impls/nepdefl.h>

 28: typedef struct {
 29:   EPS        eps;      /* linear eigensolver for T*z = mu*Tp*z */
 30:   KSP        ksp;
 31: } NEP_SLP;

 33: typedef struct {
 34:   NEP_EXT_OP extop;
 35:   Vec        w;
 36: } NEP_SLP_EPS_MSHELL;

 38: PetscErrorCode NEPSetUp_SLP(NEP nep)
 39: {
 41:   NEP_SLP        *ctx = (NEP_SLP*)nep->data;
 42:   PetscBool      istrivial,flg;
 43:   ST             st;

 46:   if (nep->ncv) { PetscInfo(nep,"Setting ncv = nev, ignoring user-provided value\n"); }
 47:   nep->ncv = nep->nev;
 48:   if (nep->mpd) { PetscInfo(nep,"Setting mpd = nev, ignoring user-provided value\n"); }
 49:   nep->mpd = nep->nev;
 50:   if (nep->ncv>nep->nev+nep->mpd) SETERRQ(PetscObjectComm((PetscObject)nep),1,"The value of ncv must not be larger than nev+mpd");
 51:   if (!nep->max_it) nep->max_it = PetscMax(5000,2*nep->n/nep->ncv);
 52:   if (nep->which && nep->which!=NEP_TARGET_MAGNITUDE) SETERRQ(PetscObjectComm((PetscObject)nep),1,"Wrong value of which");

 54:   RGIsTrivial(nep->rg,&istrivial);
 55:   if (!istrivial) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"This solver does not support region filtering");

 57:   if (!ctx->eps) { NEPSLPGetEPS(nep,&ctx->eps); }
 58:   EPSGetST(ctx->eps,&st);
 59:   PetscObjectTypeCompareAny((PetscObject)st,&flg,STSINVERT,STCAYLEY,"");
 60:   if (flg) SETERRQ(PetscObjectComm((PetscObject)nep),1,"SLP does not support spectral transformation");
 61:   EPSSetDimensions(ctx->eps,1,PETSC_DECIDE,PETSC_DECIDE);
 62:   EPSSetWhichEigenpairs(ctx->eps,EPS_LARGEST_MAGNITUDE);
 63:   EPSSetTolerances(ctx->eps,nep->tol==PETSC_DEFAULT?SLEPC_DEFAULT_TOL/10.0:nep->tol/10.0,nep->max_it?nep->max_it:PETSC_DEFAULT);

 65:   NEPAllocateSolution(nep,0);
 66:   return(0);
 67: }

 69: PetscErrorCode NEPSLPEPSMatShell_MatMult(Mat M,Vec x,Vec y)
 70: {
 71:   PetscErrorCode     ierr;
 72:   NEP_SLP_EPS_MSHELL *ctx;

 75:   MatShellGetContext(M,(void**)&ctx);
 76:   MatMult(ctx->extop->MJ,x,ctx->w);
 77:   NEPDeflationFunctionSolve(ctx->extop,ctx->w,y);
 78:   return(0);
 79: }

 81: PetscErrorCode NEPSLPEPSMatShell_Destroy(Mat M)
 82: {
 83:   PetscErrorCode     ierr;
 84:   NEP_SLP_EPS_MSHELL *ctx;

 87:   MatShellGetContext(M,(void**)&ctx);
 88:   VecDestroy(&ctx->w);
 89:   PetscFree(ctx);
 90:   return(0);
 91: }

 93: PetscErrorCode NEPSLPEPSMatShell_CreateVecs(Mat M,Vec *left,Vec *right)
 94: {
 95:   PetscErrorCode     ierr;
 96:   NEP_SLP_EPS_MSHELL *ctx;

 99:   MatShellGetContext(M,(void**)&ctx);
100:   if (right) {
101:     VecDuplicate(ctx->w,right);
102:   }
103:   if (left) {
104:     VecDuplicate(ctx->w,left);
105:   }
106:   return(0);
107: }

109: static PetscErrorCode NEPSLPSetUpLinearEP(NEP nep,NEP_EXT_OP extop,PetscScalar lambda,Vec u,PetscBool ini)
110: {
111:   PetscErrorCode     ierr;
112:   NEP_SLP            *slpctx = (NEP_SLP*)nep->data;
113:   Mat                Mshell;
114:   PetscInt           nloc,mloc;
115:   NEP_SLP_EPS_MSHELL *shellctx;

118:   if (ini) {
119:     /* Create mat shell */
120:     PetscNew(&shellctx);
121:     shellctx->extop = extop;
122:     NEPDeflationCreateVec(extop,&shellctx->w);
123:     MatGetLocalSize(nep->function,&mloc,&nloc);
124:     nloc += extop->szd; mloc += extop->szd;
125:     MatCreateShell(PetscObjectComm((PetscObject)nep),nloc,mloc,PETSC_DETERMINE,PETSC_DETERMINE,shellctx,&Mshell);
126:     MatShellSetOperation(Mshell,MATOP_MULT,(void(*)(void))NEPSLPEPSMatShell_MatMult);
127:     MatShellSetOperation(Mshell,MATOP_DESTROY,(void(*)(void))NEPSLPEPSMatShell_Destroy);
128:     MatShellSetOperation(Mshell,MATOP_CREATE_VECS,(void(*)(void))NEPSLPEPSMatShell_CreateVecs);
129:     EPSSetOperators(slpctx->eps,Mshell,NULL);
130:     MatDestroy(&Mshell);
131:   }
132:   NEPDeflationSolveSetUp(extop,lambda);
133:   NEPDeflationComputeJacobian(extop,lambda,NULL);
134:   EPSSetInitialSpace(slpctx->eps,1,&u);
135:   return(0);
136: }

138: PetscErrorCode NEPSolve_SLP(NEP nep)
139: {
141:   NEP_SLP        *ctx = (NEP_SLP*)nep->data;
142:   Mat            F,H;
143:   Vec            uu,u,r;
144:   PetscScalar    sigma,lambda,mu,im,*Hp,*Ap;
145:   PetscReal      resnorm;
146:   PetscInt       nconv,ldh,ldds,i,j;
147:   PetscBool      skip=PETSC_FALSE;
148:   NEP_EXT_OP     extop=NULL;    /* Extended operator for deflation */

151:   /* get initial approximation of eigenvalue and eigenvector */
152:   NEPGetDefaultShift(nep,&sigma);
153:   if (!nep->nini) {
154:     BVSetRandomColumn(nep->V,0);
155:   }
156:   lambda = sigma;
157:   if (!ctx->ksp) { NEPSLPGetKSP(nep,&ctx->ksp); }
158:   NEPDeflationInitialize(nep,nep->V,ctx->ksp,PETSC_TRUE,nep->nev,&extop);
159:   NEPDeflationCreateVec(extop,&u);
160:   VecDuplicate(u,&r);
161:   BVGetColumn(nep->V,0,&uu);
162:   NEPDeflationCopyToExtendedVec(extop,uu,NULL,u,PETSC_FALSE);
163:   BVRestoreColumn(nep->V,0,&uu);

165:   /* Restart loop */
166:   while (nep->reason == NEP_CONVERGED_ITERATING) {
167:     nep->its++;

169:     /* form residual,  r = T(lambda)*u (used in convergence test only) */
170:     NEPDeflationComputeFunction(extop,lambda,&F);
171:     MatMult(F,u,r);

173:     /* convergence test */
174:     VecNorm(r,NORM_2,&resnorm);
175:     (*nep->converged)(nep,lambda,0,resnorm,&nep->errest[nep->nconv],nep->convergedctx);
176:     nep->eigr[nep->nconv] = lambda;
177:     if (nep->errest[nep->nconv]<=nep->tol) {
178:       nep->nconv = nep->nconv + 1;
179:       skip = PETSC_TRUE;     
180:       NEPDeflationLocking(extop,u,lambda);
181:     }
182:     (*nep->stopping)(nep,nep->its,nep->max_it,nep->nconv,nep->nev,&nep->reason,nep->stoppingctx);
183:     if (!skip || nep->reason>0) {
184:       NEPMonitor(nep,nep->its,nep->nconv,nep->eigr,nep->eigi,nep->errest,(nep->reason>0)?nep->nconv:nep->nconv+1);
185:     }

187:     if (nep->reason == NEP_CONVERGED_ITERATING) {
188:       if (!skip) {
189:         /* evaluate T(lambda) and T'(lambda) */
190:         NEPSLPSetUpLinearEP(nep,extop,lambda,u,nep->its==1?PETSC_TRUE:PETSC_FALSE);
191:         /* compute new eigenvalue correction mu and eigenvector approximation u */
192:         EPSSolve(ctx->eps);
193:         EPSGetConverged(ctx->eps,&nconv);
194:         if (!nconv) {
195:           PetscInfo1(nep,"iter=%D, inner iteration failed, stopping solve\n",nep->its);
196:           nep->reason = NEP_DIVERGED_LINEAR_SOLVE;
197:           break;
198:         }
199:         EPSGetEigenpair(ctx->eps,0,&mu,&im,u,NULL);
200:         mu = 1.0/mu;
201:         if (PetscAbsScalar(im)>PETSC_MACHINE_EPSILON) SETERRQ(PetscObjectComm((PetscObject)nep),1,"Complex eigenvalue approximation - not implemented in real scalars");
202:       } else {
203:         nep->its--;  /* do not count this as a full iteration */
204:         /* use second eigenpair computed in previous iteration */
205:         EPSGetConverged(ctx->eps,&nconv);
206:         if (nconv>=2) {
207:           EPSGetEigenpair(ctx->eps,1,&mu,&im,u,NULL);
208:           mu = 1.0/mu;
209:         } else {
210:           NEPDeflationSetRandomVec(extop,u);
211:           mu = lambda-sigma;
212:         }
213:         skip = PETSC_FALSE;
214:       }
215:       /* correct eigenvalue */
216:       lambda = lambda - mu;
217:     }
218:   }
219:   NEPDeflationGetInvariantPair(extop,NULL,&H);
220:   MatGetSize(H,NULL,&ldh);
221:   DSSetType(nep->ds,DSNHEP);
222:   DSAllocate(nep->ds,PetscMax(nep->nconv,1));
223:   DSGetLeadingDimension(nep->ds,&ldds);
224:   MatDenseGetArray(H,&Hp);
225:   DSGetArray(nep->ds,DS_MAT_A,&Ap);
226:   for (j=0;j<nep->nconv;j++)
227:     for (i=0;i<nep->nconv;i++) Ap[j*ldds+i] = Hp[j*ldh+i];
228:   DSRestoreArray(nep->ds,DS_MAT_A,&Ap);
229:   MatDenseRestoreArray(H,&Hp);
230:   MatDestroy(&H);
231:   DSSetDimensions(nep->ds,nep->nconv,0,0,nep->nconv);
232:   DSSolve(nep->ds,nep->eigr,nep->eigi);
233:   NEPDeflationReset(extop);
234:   VecDestroy(&u);
235:   VecDestroy(&r);
236:   return(0);
237: }

239: PetscErrorCode NEPSetFromOptions_SLP(PetscOptionItems *PetscOptionsObject,NEP nep)
240: {
242:   NEP_SLP        *ctx = (NEP_SLP*)nep->data;

245:   PetscOptionsHead(PetscOptionsObject,"NEP SLP Options");
246:   PetscOptionsTail();

248:   if (!ctx->eps) { NEPSLPGetEPS(nep,&ctx->eps); }
249:   EPSSetFromOptions(ctx->eps);
250:   if (!ctx->ksp) { NEPSLPGetKSP(nep,&ctx->ksp); }
251:   KSPSetFromOptions(ctx->ksp);
252:   return(0);
253: }

255: static PetscErrorCode NEPSLPSetEPS_SLP(NEP nep,EPS eps)
256: {
258:   NEP_SLP        *ctx = (NEP_SLP*)nep->data;

261:   PetscObjectReference((PetscObject)eps);
262:   EPSDestroy(&ctx->eps);
263:   ctx->eps = eps;
264:   PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->eps);
265:   nep->state = NEP_STATE_INITIAL;
266:   return(0);
267: }

269: /*@
270:    NEPSLPSetEPS - Associate a linear eigensolver object (EPS) to the
271:    nonlinear eigenvalue solver.

273:    Collective on NEP

275:    Input Parameters:
276: +  nep - nonlinear eigenvalue solver
277: -  eps - the eigensolver object

279:    Level: advanced

281: .seealso: NEPSLPGetEPS()
282: @*/
283: PetscErrorCode NEPSLPSetEPS(NEP nep,EPS eps)
284: {

291:   PetscTryMethod(nep,"NEPSLPSetEPS_C",(NEP,EPS),(nep,eps));
292:   return(0);
293: }

295: static PetscErrorCode NEPSLPGetEPS_SLP(NEP nep,EPS *eps)
296: {
298:   NEP_SLP        *ctx = (NEP_SLP*)nep->data;

301:   if (!ctx->eps) {
302:     EPSCreate(PetscObjectComm((PetscObject)nep),&ctx->eps);
303:     PetscObjectIncrementTabLevel((PetscObject)ctx->eps,(PetscObject)nep,1);
304:     EPSSetOptionsPrefix(ctx->eps,((PetscObject)nep)->prefix);
305:     EPSAppendOptionsPrefix(ctx->eps,"nep_slp_");
306:     PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->eps);
307:     PetscObjectSetOptions((PetscObject)ctx->eps,((PetscObject)nep)->options);
308:   }
309:   *eps = ctx->eps;
310:   return(0);
311: }

313: /*@
314:    NEPSLPGetEPS - Retrieve the linear eigensolver object (EPS) associated
315:    to the nonlinear eigenvalue solver.

317:    Not Collective

319:    Input Parameter:
320: .  nep - nonlinear eigenvalue solver

322:    Output Parameter:
323: .  eps - the eigensolver object

325:    Level: advanced

327: .seealso: NEPSLPSetEPS()
328: @*/
329: PetscErrorCode NEPSLPGetEPS(NEP nep,EPS *eps)
330: {

336:   PetscUseMethod(nep,"NEPSLPGetEPS_C",(NEP,EPS*),(nep,eps));
337:   return(0);
338: }

340: static PetscErrorCode NEPSLPSetKSP_SLP(NEP nep,KSP ksp)
341: {
343:   NEP_SLP        *ctx = (NEP_SLP*)nep->data;

346:   PetscObjectReference((PetscObject)ksp);
347:   KSPDestroy(&ctx->ksp);
348:   ctx->ksp = ksp;
349:   PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->ksp);
350:   nep->state = NEP_STATE_INITIAL;
351:   return(0);
352: }

354: /*@
355:    NEPSLPSetKSP - Associate a linear solver object (KSP) to the nonlinear
356:    eigenvalue solver.

358:    Collective on NEP

360:    Input Parameters:
361: +  nep - eigenvalue solver
362: -  ksp - the linear solver object

364:    Level: advanced

366: .seealso: NEPSLPGetKSP()
367: @*/
368: PetscErrorCode NEPSLPSetKSP(NEP nep,KSP ksp)
369: {

376:   PetscTryMethod(nep,"NEPSLPSetKSP_C",(NEP,KSP),(nep,ksp));
377:   return(0);
378: }

380: static PetscErrorCode NEPSLPGetKSP_SLP(NEP nep,KSP *ksp)
381: {
383:   NEP_SLP        *ctx = (NEP_SLP*)nep->data;

386:   if (!ctx->ksp) {
387:     KSPCreate(PetscObjectComm((PetscObject)nep),&ctx->ksp);
388:     PetscObjectIncrementTabLevel((PetscObject)ctx->ksp,(PetscObject)nep,1);
389:     KSPSetOptionsPrefix(ctx->ksp,((PetscObject)nep)->prefix);
390:     KSPAppendOptionsPrefix(ctx->ksp,"nep_slp_");
391:     PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->ksp);
392:     KSPSetErrorIfNotConverged(ctx->ksp,PETSC_TRUE);
393:     KSPSetTolerances(ctx->ksp,SLEPC_DEFAULT_TOL,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
394:   }
395:   *ksp = ctx->ksp;
396:   return(0);
397: }

399: /*@
400:    NEPSLPGetKSP - Retrieve the linear solver object (KSP) associated with
401:    the nonlinear eigenvalue solver.

403:    Not Collective

405:    Input Parameter:
406: .  nep - nonlinear eigenvalue solver

408:    Output Parameter:
409: .  ksp - the linear solver object

411:    Level: advanced

413: .seealso: NEPSLPSetKSP()
414: @*/
415: PetscErrorCode NEPSLPGetKSP(NEP nep,KSP *ksp)
416: {

422:   PetscUseMethod(nep,"NEPSLPGetKSP_C",(NEP,KSP*),(nep,ksp));
423:   return(0);
424: }

426: PetscErrorCode NEPView_SLP(NEP nep,PetscViewer viewer)
427: {
429:   NEP_SLP        *ctx = (NEP_SLP*)nep->data;
430:   PetscBool      isascii;

433:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
434:   if (isascii) {
435:     if (!ctx->eps) { NEPSLPGetEPS(nep,&ctx->eps); }
436:     PetscViewerASCIIPushTab(viewer);
437:     EPSView(ctx->eps,viewer);
438:     PetscViewerASCIIPopTab(viewer);
439:   }
440:   return(0);
441: }

443: PetscErrorCode NEPReset_SLP(NEP nep)
444: {
446:   NEP_SLP        *ctx = (NEP_SLP*)nep->data;

449:   EPSReset(ctx->eps);
450:   KSPReset(ctx->ksp);
451:   return(0);
452: }

454: PetscErrorCode NEPDestroy_SLP(NEP nep)
455: {
457:   NEP_SLP        *ctx = (NEP_SLP*)nep->data;

460:   KSPDestroy(&ctx->ksp);
461:   EPSDestroy(&ctx->eps);
462:   PetscFree(nep->data);
463:   PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetEPS_C",NULL);
464:   PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetEPS_C",NULL);
465:   PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetKSP_C",NULL);
466:   PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetKSP_C",NULL);
467:   return(0);
468: }

470: SLEPC_EXTERN PetscErrorCode NEPCreate_SLP(NEP nep)
471: {
473:   NEP_SLP        *ctx;

476:   PetscNewLog(nep,&ctx);
477:   nep->data = (void*)ctx;

479:   nep->useds = PETSC_TRUE;

481:   nep->ops->solve          = NEPSolve_SLP;
482:   nep->ops->setup          = NEPSetUp_SLP;
483:   nep->ops->setfromoptions = NEPSetFromOptions_SLP;
484:   nep->ops->reset          = NEPReset_SLP;
485:   nep->ops->destroy        = NEPDestroy_SLP;
486:   nep->ops->view           = NEPView_SLP;
487:   nep->ops->computevectors = NEPComputeVectors_Schur;

489:   PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetEPS_C",NEPSLPSetEPS_SLP);
490:   PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetEPS_C",NEPSLPGetEPS_SLP);
491:   PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetKSP_C",NEPSLPSetKSP_SLP);
492:   PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetKSP_C",NEPSLPGetKSP_SLP);
493:   return(0);
494: }