Actual source code: narnoldi.c

slepc-3.22.1 2024-10-28
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, 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: "narnoldi"

 13:    Method: Nonlinear Arnoldi

 15:    Algorithm:

 17:        Arnoldi for nonlinear eigenproblems.

 19:    References:

 21:        [1] H. Voss, "An Arnoldi method for nonlinear eigenvalue problems",
 22:            BIT 44:387-401, 2004.
 23: */

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

 28: typedef struct {
 29:   PetscInt lag;             /* interval to rebuild preconditioner */
 30:   KSP      ksp;             /* linear solver object */
 31: } NEP_NARNOLDI;

 33: static PetscErrorCode NEPSetUp_NArnoldi(NEP nep)
 34: {
 35:   PetscFunctionBegin;
 36:   PetscCall(NEPSetDimensions_Default(nep,nep->nev,&nep->ncv,&nep->mpd));
 37:   PetscCheck(nep->ncv<=nep->nev+nep->mpd,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER_INPUT,"The value of ncv must not be larger than nev+mpd");
 38:   if (nep->max_it==PETSC_DETERMINE) nep->max_it = nep->nev*nep->ncv;
 39:   if (!nep->which) nep->which = NEP_TARGET_MAGNITUDE;
 40:   PetscCheck(nep->which==NEP_TARGET_MAGNITUDE,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"This solver supports only target magnitude eigenvalues");
 41:   NEPCheckUnsupported(nep,NEP_FEATURE_CALLBACK | NEP_FEATURE_REGION | NEP_FEATURE_TWOSIDED);
 42:   PetscCall(NEPAllocateSolution(nep,0));
 43:   PetscCall(NEPSetWorkVecs(nep,3));
 44:   PetscFunctionReturn(PETSC_SUCCESS);
 45: }

 47: static PetscErrorCode NEPSolve_NArnoldi(NEP nep)
 48: {
 49:   NEP_NARNOLDI       *ctx = (NEP_NARNOLDI*)nep->data;
 50:   Mat                T,H,A;
 51:   Vec                f,r,u,uu;
 52:   PetscScalar        *X,lambda=0.0,lambda2=0.0,*eigr,sigma;
 53:   PetscReal          beta,resnorm=0.0,nrm,perr=0.0;
 54:   PetscInt           n;
 55:   PetscBool          breakdown,skip=PETSC_FALSE;
 56:   BV                 Vext;
 57:   DS                 ds;
 58:   NEP_EXT_OP         extop=NULL;
 59:   SlepcSC            sc;
 60:   KSPConvergedReason kspreason;

 62:   PetscFunctionBegin;
 63:   /* get initial space and shift */
 64:   PetscCall(NEPGetDefaultShift(nep,&sigma));
 65:   if (!nep->nini) {
 66:     PetscCall(BVSetRandomColumn(nep->V,0));
 67:     PetscCall(BVNormColumn(nep->V,0,NORM_2,&nrm));
 68:     PetscCall(BVScaleColumn(nep->V,0,1.0/nrm));
 69:     n = 1;
 70:   } else n = nep->nini;

 72:   if (!ctx->ksp) PetscCall(NEPNArnoldiGetKSP(nep,&ctx->ksp));
 73:   PetscCall(NEPDeflationInitialize(nep,nep->V,ctx->ksp,PETSC_FALSE,nep->nev,&extop));
 74:   PetscCall(NEPDeflationCreateBV(extop,nep->ncv,&Vext));

 76:   /* prepare linear solver */
 77:   PetscCall(NEPDeflationSolveSetUp(extop,sigma));

 79:   PetscCall(BVGetColumn(Vext,0,&f));
 80:   PetscCall(VecDuplicate(f,&r));
 81:   PetscCall(VecDuplicate(f,&u));
 82:   PetscCall(BVGetColumn(nep->V,0,&uu));
 83:   PetscCall(NEPDeflationCopyToExtendedVec(extop,uu,NULL,f,PETSC_FALSE));
 84:   PetscCall(BVRestoreColumn(nep->V,0,&uu));
 85:   PetscCall(VecCopy(f,r));
 86:   PetscCall(NEPDeflationFunctionSolve(extop,r,f));
 87:   PetscCall(VecNorm(f,NORM_2,&nrm));
 88:   PetscCall(VecScale(f,1.0/nrm));
 89:   PetscCall(BVRestoreColumn(Vext,0,&f));

 91:   PetscCall(DSCreate(PetscObjectComm((PetscObject)nep),&ds));
 92:   PetscCall(DSSetType(ds,DSNEP));
 93:   PetscCall(DSNEPSetFN(ds,nep->nt,nep->f));
 94:   PetscCall(DSAllocate(ds,nep->ncv));
 95:   PetscCall(DSGetSlepcSC(ds,&sc));
 96:   sc->comparison    = nep->sc->comparison;
 97:   sc->comparisonctx = nep->sc->comparisonctx;
 98:   PetscCall(DSSetFromOptions(ds));

100:   /* build projected matrices for initial space */
101:   PetscCall(DSSetDimensions(ds,n,0,0));
102:   PetscCall(NEPDeflationProjectOperator(extop,Vext,ds,0,n));

104:   PetscCall(PetscMalloc1(nep->ncv,&eigr));

106:   /* Restart loop */
107:   while (nep->reason == NEP_CONVERGED_ITERATING) {
108:     nep->its++;

110:     /* solve projected problem */
111:     PetscCall(DSSetDimensions(ds,n,0,0));
112:     PetscCall(DSSetState(ds,DS_STATE_RAW));
113:     PetscCall(DSSolve(ds,eigr,NULL));
114:     PetscCall(DSSynchronize(ds,eigr,NULL));
115:     if (nep->its>1) lambda2 = lambda;
116:     lambda = eigr[0];
117:     nep->eigr[nep->nconv] = lambda;

119:     /* compute Ritz vector, x = V*s */
120:     PetscCall(DSGetArray(ds,DS_MAT_X,&X));
121:     PetscCall(BVSetActiveColumns(Vext,0,n));
122:     PetscCall(BVMultVec(Vext,1.0,0.0,u,X));
123:     PetscCall(DSRestoreArray(ds,DS_MAT_X,&X));

125:     /* compute the residual, r = T(lambda)*x */
126:     PetscCall(NEPDeflationComputeFunction(extop,lambda,&T));
127:     PetscCall(MatMult(T,u,r));

129:     /* convergence test */
130:     PetscCall(VecNorm(r,NORM_2,&resnorm));
131:     if (nep->its>1) perr=nep->errest[nep->nconv];
132:     PetscCall((*nep->converged)(nep,lambda,0,resnorm,&nep->errest[nep->nconv],nep->convergedctx));
133:     if (nep->errest[nep->nconv]<=nep->tol) {
134:       nep->nconv = nep->nconv + 1;
135:       PetscCall(NEPDeflationLocking(extop,u,lambda));
136:       skip = PETSC_TRUE;
137:     }
138:     PetscCall((*nep->stopping)(nep,nep->its,nep->max_it,nep->nconv,nep->nev,&nep->reason,nep->stoppingctx));
139:     if (!skip || nep->reason>0) PetscCall(NEPMonitor(nep,nep->its,nep->nconv,nep->eigr,nep->eigi,nep->errest,(nep->reason>0)?nep->nconv:nep->nconv+1));

141:     if (nep->reason == NEP_CONVERGED_ITERATING) {
142:       if (!skip) {
143:         if (n>=nep->ncv) {
144:           nep->reason = NEP_DIVERGED_SUBSPACE_EXHAUSTED;
145:           break;
146:         }
147:         if (ctx->lag && !(nep->its%ctx->lag) && nep->its>=2*ctx->lag && perr && nep->errest[nep->nconv]>.5*perr) PetscCall(NEPDeflationSolveSetUp(extop,lambda2));

149:         /* continuation vector: f = T(sigma)\r */
150:         PetscCall(BVGetColumn(Vext,n,&f));
151:         PetscCall(NEPDeflationFunctionSolve(extop,r,f));
152:         PetscCall(BVRestoreColumn(Vext,n,&f));
153:         PetscCall(KSPGetConvergedReason(ctx->ksp,&kspreason));
154:         if (kspreason<0) {
155:           PetscCall(PetscInfo(nep,"iter=%" PetscInt_FMT ", linear solve failed, stopping solve\n",nep->its));
156:           nep->reason = NEP_DIVERGED_LINEAR_SOLVE;
157:           break;
158:         }

160:         /* orthonormalize */
161:         PetscCall(BVOrthonormalizeColumn(Vext,n,PETSC_FALSE,&beta,&breakdown));
162:         if (breakdown || beta==0.0) {
163:           PetscCall(PetscInfo(nep,"iter=%" PetscInt_FMT ", orthogonalization failed, stopping solve\n",nep->its));
164:           nep->reason = NEP_DIVERGED_BREAKDOWN;
165:           break;
166:         }

168:         /* update projected matrices */
169:         PetscCall(DSSetDimensions(ds,n+1,0,0));
170:         PetscCall(NEPDeflationProjectOperator(extop,Vext,ds,n,n+1));
171:         n++;
172:       } else {
173:         nep->its--;  /* do not count this as a full iteration */
174:         PetscCall(BVGetColumn(Vext,0,&f));
175:         PetscCall(NEPDeflationSetRandomVec(extop,f));
176:         PetscCall(NEPDeflationSolveSetUp(extop,sigma));
177:         PetscCall(VecCopy(f,r));
178:         PetscCall(NEPDeflationFunctionSolve(extop,r,f));
179:         PetscCall(VecNorm(f,NORM_2,&nrm));
180:         PetscCall(VecScale(f,1.0/nrm));
181:         PetscCall(BVRestoreColumn(Vext,0,&f));
182:         n = 1;
183:         PetscCall(DSSetDimensions(ds,n,0,0));
184:         PetscCall(NEPDeflationProjectOperator(extop,Vext,ds,n-1,n));
185:         skip = PETSC_FALSE;
186:       }
187:     }
188:   }

190:   PetscCall(NEPDeflationGetInvariantPair(extop,NULL,&H));
191:   PetscCall(DSSetType(nep->ds,DSNHEP));
192:   PetscCall(DSAllocate(nep->ds,PetscMax(nep->nconv,1)));
193:   PetscCall(DSSetDimensions(nep->ds,nep->nconv,0,nep->nconv));
194:   PetscCall(DSGetMat(nep->ds,DS_MAT_A,&A));
195:   PetscCall(MatCopy(H,A,SAME_NONZERO_PATTERN));
196:   PetscCall(DSRestoreMat(nep->ds,DS_MAT_A,&A));
197:   PetscCall(MatDestroy(&H));
198:   PetscCall(DSSolve(nep->ds,nep->eigr,nep->eigi));
199:   PetscCall(NEPDeflationReset(extop));
200:   PetscCall(VecDestroy(&u));
201:   PetscCall(VecDestroy(&r));
202:   PetscCall(BVDestroy(&Vext));
203:   PetscCall(DSDestroy(&ds));
204:   PetscCall(PetscFree(eigr));
205:   PetscFunctionReturn(PETSC_SUCCESS);
206: }

208: static PetscErrorCode NEPNArnoldiSetLagPreconditioner_NArnoldi(NEP nep,PetscInt lag)
209: {
210:   NEP_NARNOLDI *ctx = (NEP_NARNOLDI*)nep->data;

212:   PetscFunctionBegin;
213:   PetscCheck(lag>=0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be non-negative");
214:   ctx->lag = lag;
215:   PetscFunctionReturn(PETSC_SUCCESS);
216: }

218: /*@
219:    NEPNArnoldiSetLagPreconditioner - Determines when the preconditioner is rebuilt in the
220:    nonlinear solve.

222:    Logically Collective

224:    Input Parameters:
225: +  nep - nonlinear eigenvalue solver
226: -  lag - 0 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is
227:           computed within the nonlinear iteration, 2 means every second time
228:           the Jacobian is built, etc.

230:    Options Database Keys:
231: .  -nep_narnoldi_lag_preconditioner <lag> - the lag value

233:    Notes:
234:    The default is 1.
235:    The preconditioner is ALWAYS built in the first iteration of a nonlinear solve.

237:    Level: intermediate

239: .seealso: NEPNArnoldiGetLagPreconditioner()
240: @*/
241: PetscErrorCode NEPNArnoldiSetLagPreconditioner(NEP nep,PetscInt lag)
242: {
243:   PetscFunctionBegin;
246:   PetscTryMethod(nep,"NEPNArnoldiSetLagPreconditioner_C",(NEP,PetscInt),(nep,lag));
247:   PetscFunctionReturn(PETSC_SUCCESS);
248: }

250: static PetscErrorCode NEPNArnoldiGetLagPreconditioner_NArnoldi(NEP nep,PetscInt *lag)
251: {
252:   NEP_NARNOLDI *ctx = (NEP_NARNOLDI*)nep->data;

254:   PetscFunctionBegin;
255:   *lag = ctx->lag;
256:   PetscFunctionReturn(PETSC_SUCCESS);
257: }

259: /*@
260:    NEPNArnoldiGetLagPreconditioner - Indicates how often the preconditioner is rebuilt.

262:    Not Collective

264:    Input Parameter:
265: .  nep - nonlinear eigenvalue solver

267:    Output Parameter:
268: .  lag - the lag parameter

270:    Level: intermediate

272: .seealso: NEPNArnoldiSetLagPreconditioner()
273: @*/
274: PetscErrorCode NEPNArnoldiGetLagPreconditioner(NEP nep,PetscInt *lag)
275: {
276:   PetscFunctionBegin;
278:   PetscAssertPointer(lag,2);
279:   PetscUseMethod(nep,"NEPNArnoldiGetLagPreconditioner_C",(NEP,PetscInt*),(nep,lag));
280:   PetscFunctionReturn(PETSC_SUCCESS);
281: }

283: static PetscErrorCode NEPSetFromOptions_NArnoldi(NEP nep,PetscOptionItems *PetscOptionsObject)
284: {
285:   PetscInt       i;
286:   PetscBool      flg;
287:   NEP_NARNOLDI   *ctx = (NEP_NARNOLDI*)nep->data;

289:   PetscFunctionBegin;
290:   PetscOptionsHeadBegin(PetscOptionsObject,"NEP N-Arnoldi Options");
291:     i = 0;
292:     PetscCall(PetscOptionsInt("-nep_narnoldi_lag_preconditioner","Interval to rebuild preconditioner","NEPNArnoldiSetLagPreconditioner",ctx->lag,&i,&flg));
293:     if (flg) PetscCall(NEPNArnoldiSetLagPreconditioner(nep,i));

295:   PetscOptionsHeadEnd();

297:   if (!ctx->ksp) PetscCall(NEPNArnoldiGetKSP(nep,&ctx->ksp));
298:   PetscCall(KSPSetFromOptions(ctx->ksp));
299:   PetscFunctionReturn(PETSC_SUCCESS);
300: }

302: static PetscErrorCode NEPNArnoldiSetKSP_NArnoldi(NEP nep,KSP ksp)
303: {
304:   NEP_NARNOLDI   *ctx = (NEP_NARNOLDI*)nep->data;

306:   PetscFunctionBegin;
307:   PetscCall(PetscObjectReference((PetscObject)ksp));
308:   PetscCall(KSPDestroy(&ctx->ksp));
309:   ctx->ksp = ksp;
310:   nep->state = NEP_STATE_INITIAL;
311:   PetscFunctionReturn(PETSC_SUCCESS);
312: }

314: /*@
315:    NEPNArnoldiSetKSP - Associate a linear solver object (KSP) to the nonlinear
316:    eigenvalue solver.

318:    Collective

320:    Input Parameters:
321: +  nep - eigenvalue solver
322: -  ksp - the linear solver object

324:    Level: advanced

326: .seealso: NEPNArnoldiGetKSP()
327: @*/
328: PetscErrorCode NEPNArnoldiSetKSP(NEP nep,KSP ksp)
329: {
330:   PetscFunctionBegin;
333:   PetscCheckSameComm(nep,1,ksp,2);
334:   PetscTryMethod(nep,"NEPNArnoldiSetKSP_C",(NEP,KSP),(nep,ksp));
335:   PetscFunctionReturn(PETSC_SUCCESS);
336: }

338: static PetscErrorCode NEPNArnoldiGetKSP_NArnoldi(NEP nep,KSP *ksp)
339: {
340:   NEP_NARNOLDI   *ctx = (NEP_NARNOLDI*)nep->data;

342:   PetscFunctionBegin;
343:   if (!ctx->ksp) {
344:     PetscCall(KSPCreate(PetscObjectComm((PetscObject)nep),&ctx->ksp));
345:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->ksp,(PetscObject)nep,1));
346:     PetscCall(KSPSetOptionsPrefix(ctx->ksp,((PetscObject)nep)->prefix));
347:     PetscCall(KSPAppendOptionsPrefix(ctx->ksp,"nep_narnoldi_"));
348:     PetscCall(PetscObjectSetOptions((PetscObject)ctx->ksp,((PetscObject)nep)->options));
349:     PetscCall(KSPSetErrorIfNotConverged(ctx->ksp,PETSC_TRUE));
350:     PetscCall(KSPSetTolerances(ctx->ksp,SlepcDefaultTol(nep->tol),PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT));
351:   }
352:   *ksp = ctx->ksp;
353:   PetscFunctionReturn(PETSC_SUCCESS);
354: }

356: /*@
357:    NEPNArnoldiGetKSP - Retrieve the linear solver object (KSP) associated with
358:    the nonlinear eigenvalue solver.

360:    Collective

362:    Input Parameter:
363: .  nep - nonlinear eigenvalue solver

365:    Output Parameter:
366: .  ksp - the linear solver object

368:    Level: advanced

370: .seealso: NEPNArnoldiSetKSP()
371: @*/
372: PetscErrorCode NEPNArnoldiGetKSP(NEP nep,KSP *ksp)
373: {
374:   PetscFunctionBegin;
376:   PetscAssertPointer(ksp,2);
377:   PetscUseMethod(nep,"NEPNArnoldiGetKSP_C",(NEP,KSP*),(nep,ksp));
378:   PetscFunctionReturn(PETSC_SUCCESS);
379: }

381: static PetscErrorCode NEPView_NArnoldi(NEP nep,PetscViewer viewer)
382: {
383:   NEP_NARNOLDI   *ctx = (NEP_NARNOLDI*)nep->data;
384:   PetscBool      isascii;

386:   PetscFunctionBegin;
387:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
388:   if (isascii) {
389:     if (ctx->lag) PetscCall(PetscViewerASCIIPrintf(viewer,"  updating the preconditioner every %" PetscInt_FMT " iterations\n",ctx->lag));
390:     if (!ctx->ksp) PetscCall(NEPNArnoldiGetKSP(nep,&ctx->ksp));
391:     PetscCall(PetscViewerASCIIPushTab(viewer));
392:     PetscCall(KSPView(ctx->ksp,viewer));
393:     PetscCall(PetscViewerASCIIPopTab(viewer));
394:   }
395:   PetscFunctionReturn(PETSC_SUCCESS);
396: }

398: static PetscErrorCode NEPReset_NArnoldi(NEP nep)
399: {
400:   NEP_NARNOLDI   *ctx = (NEP_NARNOLDI*)nep->data;

402:   PetscFunctionBegin;
403:   PetscCall(KSPReset(ctx->ksp));
404:   PetscFunctionReturn(PETSC_SUCCESS);
405: }

407: static PetscErrorCode NEPDestroy_NArnoldi(NEP nep)
408: {
409:   NEP_NARNOLDI   *ctx = (NEP_NARNOLDI*)nep->data;

411:   PetscFunctionBegin;
412:   PetscCall(KSPDestroy(&ctx->ksp));
413:   PetscCall(PetscFree(nep->data));
414:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiSetLagPreconditioner_C",NULL));
415:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiGetLagPreconditioner_C",NULL));
416:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiSetKSP_C",NULL));
417:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiGetKSP_C",NULL));
418:   PetscFunctionReturn(PETSC_SUCCESS);
419: }

421: SLEPC_EXTERN PetscErrorCode NEPCreate_NArnoldi(NEP nep)
422: {
423:   NEP_NARNOLDI   *ctx;

425:   PetscFunctionBegin;
426:   PetscCall(PetscNew(&ctx));
427:   nep->data = (void*)ctx;
428:   ctx->lag  = 1;

430:   nep->useds = PETSC_TRUE;

432:   nep->ops->solve          = NEPSolve_NArnoldi;
433:   nep->ops->setup          = NEPSetUp_NArnoldi;
434:   nep->ops->setfromoptions = NEPSetFromOptions_NArnoldi;
435:   nep->ops->reset          = NEPReset_NArnoldi;
436:   nep->ops->destroy        = NEPDestroy_NArnoldi;
437:   nep->ops->view           = NEPView_NArnoldi;
438:   nep->ops->computevectors = NEPComputeVectors_Schur;

440:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiSetLagPreconditioner_C",NEPNArnoldiSetLagPreconditioner_NArnoldi));
441:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiGetLagPreconditioner_C",NEPNArnoldiGetLagPreconditioner_NArnoldi));
442:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiSetKSP_C",NEPNArnoldiSetKSP_NArnoldi));
443:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiGetKSP_C",NEPNArnoldiGetKSP_NArnoldi));
444:   PetscFunctionReturn(PETSC_SUCCESS);
445: }