Actual source code: narnoldi.c
slepc-3.21.1 2024-04-26
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_DEFAULT) 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_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT));
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: }