Actual source code: narnoldi.c
slepc-3.17.0 2022-03-31
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: PetscErrorCode NEPSetUp_NArnoldi(NEP nep)
34: {
35: NEPSetDimensions_Default(nep,nep->nev,&nep->ncv,&nep->mpd);
37: if (nep->max_it==PETSC_DEFAULT) nep->max_it = nep->nev*nep->ncv;
38: if (!nep->which) nep->which = NEP_TARGET_MAGNITUDE;
40: NEPCheckUnsupported(nep,NEP_FEATURE_CALLBACK | NEP_FEATURE_REGION | NEP_FEATURE_TWOSIDED);
41: NEPAllocateSolution(nep,0);
42: NEPSetWorkVecs(nep,3);
43: PetscFunctionReturn(0);
44: }
46: PetscErrorCode NEPSolve_NArnoldi(NEP nep)
47: {
48: NEP_NARNOLDI *ctx = (NEP_NARNOLDI*)nep->data;
49: Mat T,H;
50: Vec f,r,u,uu;
51: PetscScalar *X,lambda=0.0,lambda2=0.0,*eigr,*Ap,sigma;
52: const PetscScalar *Hp;
53: PetscReal beta,resnorm=0.0,nrm,perr=0.0;
54: PetscInt n,i,j,ldds,ldh;
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: /* get initial space and shift */
63: NEPGetDefaultShift(nep,&sigma);
64: if (!nep->nini) {
65: BVSetRandomColumn(nep->V,0);
66: BVNormColumn(nep->V,0,NORM_2,&nrm);
67: BVScaleColumn(nep->V,0,1.0/nrm);
68: n = 1;
69: } else n = nep->nini;
71: if (!ctx->ksp) NEPNArnoldiGetKSP(nep,&ctx->ksp);
72: NEPDeflationInitialize(nep,nep->V,ctx->ksp,PETSC_FALSE,nep->nev,&extop);
73: NEPDeflationCreateBV(extop,nep->ncv,&Vext);
75: /* prepare linear solver */
76: NEPDeflationSolveSetUp(extop,sigma);
78: BVGetColumn(Vext,0,&f);
79: VecDuplicate(f,&r);
80: VecDuplicate(f,&u);
81: BVGetColumn(nep->V,0,&uu);
82: NEPDeflationCopyToExtendedVec(extop,uu,NULL,f,PETSC_FALSE);
83: BVRestoreColumn(nep->V,0,&uu);
84: VecCopy(f,r);
85: NEPDeflationFunctionSolve(extop,r,f);
86: VecNorm(f,NORM_2,&nrm);
87: VecScale(f,1.0/nrm);
88: BVRestoreColumn(Vext,0,&f);
90: DSCreate(PetscObjectComm((PetscObject)nep),&ds);
91: PetscLogObjectParent((PetscObject)nep,(PetscObject)ds);
92: DSSetType(ds,DSNEP);
93: DSNEPSetFN(ds,nep->nt,nep->f);
94: DSAllocate(ds,nep->ncv);
95: DSGetSlepcSC(ds,&sc);
96: sc->comparison = nep->sc->comparison;
97: sc->comparisonctx = nep->sc->comparisonctx;
98: DSSetFromOptions(ds);
100: /* build projected matrices for initial space */
101: DSSetDimensions(ds,n,0,0);
102: NEPDeflationProjectOperator(extop,Vext,ds,0,n);
104: PetscMalloc1(nep->ncv,&eigr);
106: /* Restart loop */
107: while (nep->reason == NEP_CONVERGED_ITERATING) {
108: nep->its++;
110: /* solve projected problem */
111: DSSetDimensions(ds,n,0,0);
112: DSSetState(ds,DS_STATE_RAW);
113: DSSolve(ds,eigr,NULL);
114: 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: DSGetArray(ds,DS_MAT_X,&X);
121: BVSetActiveColumns(Vext,0,n);
122: BVMultVec(Vext,1.0,0.0,u,X);
123: DSRestoreArray(ds,DS_MAT_X,&X);
125: /* compute the residual, r = T(lambda)*x */
126: NEPDeflationComputeFunction(extop,lambda,&T);
127: MatMult(T,u,r);
129: /* convergence test */
130: VecNorm(r,NORM_2,&resnorm);
131: if (nep->its>1) perr=nep->errest[nep->nconv];
132: (*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: NEPDeflationLocking(extop,u,lambda);
136: skip = PETSC_TRUE;
137: }
138: (*nep->stopping)(nep,nep->its,nep->max_it,nep->nconv,nep->nev,&nep->reason,nep->stoppingctx);
139: if (!skip || nep->reason>0) 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) NEPDeflationSolveSetUp(extop,lambda2);
149: /* continuation vector: f = T(sigma)\r */
150: BVGetColumn(Vext,n,&f);
151: NEPDeflationFunctionSolve(extop,r,f);
152: BVRestoreColumn(Vext,n,&f);
153: KSPGetConvergedReason(ctx->ksp,&kspreason);
154: if (kspreason<0) {
155: 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: BVOrthonormalizeColumn(Vext,n,PETSC_FALSE,&beta,&breakdown);
162: if (breakdown || beta==0.0) {
163: 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: DSSetDimensions(ds,n+1,0,0);
170: NEPDeflationProjectOperator(extop,Vext,ds,n,n+1);
171: n++;
172: } else {
173: nep->its--; /* do not count this as a full iteration */
174: BVGetColumn(Vext,0,&f);
175: NEPDeflationSetRandomVec(extop,f);
176: NEPDeflationSolveSetUp(extop,sigma);
177: VecCopy(f,r);
178: NEPDeflationFunctionSolve(extop,r,f);
179: VecNorm(f,NORM_2,&nrm);
180: VecScale(f,1.0/nrm);
181: BVRestoreColumn(Vext,0,&f);
182: n = 1;
183: DSSetDimensions(ds,n,0,0);
184: NEPDeflationProjectOperator(extop,Vext,ds,n-1,n);
185: skip = PETSC_FALSE;
186: }
187: }
188: }
190: NEPDeflationGetInvariantPair(extop,NULL,&H);
191: MatGetSize(H,NULL,&ldh);
192: DSSetType(nep->ds,DSNHEP);
193: DSAllocate(nep->ds,PetscMax(nep->nconv,1));
194: DSGetLeadingDimension(nep->ds,&ldds);
195: MatDenseGetArrayRead(H,&Hp);
196: DSGetArray(nep->ds,DS_MAT_A,&Ap);
197: for (j=0;j<nep->nconv;j++)
198: for (i=0;i<nep->nconv;i++) Ap[j*ldds+i] = Hp[j*ldh+i];
199: DSRestoreArray(nep->ds,DS_MAT_A,&Ap);
200: MatDenseRestoreArrayRead(H,&Hp);
201: MatDestroy(&H);
202: DSSetDimensions(nep->ds,nep->nconv,0,nep->nconv);
203: DSSolve(nep->ds,nep->eigr,nep->eigi);
204: NEPDeflationReset(extop);
205: VecDestroy(&u);
206: VecDestroy(&r);
207: BVDestroy(&Vext);
208: DSDestroy(&ds);
209: PetscFree(eigr);
210: PetscFunctionReturn(0);
211: }
213: static PetscErrorCode NEPNArnoldiSetLagPreconditioner_NArnoldi(NEP nep,PetscInt lag)
214: {
215: NEP_NARNOLDI *ctx = (NEP_NARNOLDI*)nep->data;
218: ctx->lag = lag;
219: PetscFunctionReturn(0);
220: }
222: /*@
223: NEPNArnoldiSetLagPreconditioner - Determines when the preconditioner is rebuilt in the
224: nonlinear solve.
226: Logically Collective on nep
228: Input Parameters:
229: + nep - nonlinear eigenvalue solver
230: - lag - 0 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is
231: computed within the nonlinear iteration, 2 means every second time
232: the Jacobian is built, etc.
234: Options Database Keys:
235: . -nep_narnoldi_lag_preconditioner <lag> - the lag value
237: Notes:
238: The default is 1.
239: The preconditioner is ALWAYS built in the first iteration of a nonlinear solve.
241: Level: intermediate
243: .seealso: NEPNArnoldiGetLagPreconditioner()
244: @*/
245: PetscErrorCode NEPNArnoldiSetLagPreconditioner(NEP nep,PetscInt lag)
246: {
249: PetscTryMethod(nep,"NEPNArnoldiSetLagPreconditioner_C",(NEP,PetscInt),(nep,lag));
250: PetscFunctionReturn(0);
251: }
253: static PetscErrorCode NEPNArnoldiGetLagPreconditioner_NArnoldi(NEP nep,PetscInt *lag)
254: {
255: NEP_NARNOLDI *ctx = (NEP_NARNOLDI*)nep->data;
257: *lag = ctx->lag;
258: PetscFunctionReturn(0);
259: }
261: /*@
262: NEPNArnoldiGetLagPreconditioner - Indicates how often the preconditioner is rebuilt.
264: Not Collective
266: Input Parameter:
267: . nep - nonlinear eigenvalue solver
269: Output Parameter:
270: . lag - the lag parameter
272: Level: intermediate
274: .seealso: NEPNArnoldiSetLagPreconditioner()
275: @*/
276: PetscErrorCode NEPNArnoldiGetLagPreconditioner(NEP nep,PetscInt *lag)
277: {
280: PetscUseMethod(nep,"NEPNArnoldiGetLagPreconditioner_C",(NEP,PetscInt*),(nep,lag));
281: PetscFunctionReturn(0);
282: }
284: PetscErrorCode NEPSetFromOptions_NArnoldi(PetscOptionItems *PetscOptionsObject,NEP nep)
285: {
286: PetscInt i;
287: PetscBool flg;
288: NEP_NARNOLDI *ctx = (NEP_NARNOLDI*)nep->data;
290: PetscOptionsHead(PetscOptionsObject,"NEP N-Arnoldi Options");
291: i = 0;
292: PetscOptionsInt("-nep_narnoldi_lag_preconditioner","Interval to rebuild preconditioner","NEPNArnoldiSetLagPreconditioner",ctx->lag,&i,&flg);
293: if (flg) NEPNArnoldiSetLagPreconditioner(nep,i);
295: PetscOptionsTail();
297: if (!ctx->ksp) NEPNArnoldiGetKSP(nep,&ctx->ksp);
298: KSPSetFromOptions(ctx->ksp);
299: PetscFunctionReturn(0);
300: }
302: static PetscErrorCode NEPNArnoldiSetKSP_NArnoldi(NEP nep,KSP ksp)
303: {
304: NEP_NARNOLDI *ctx = (NEP_NARNOLDI*)nep->data;
306: PetscObjectReference((PetscObject)ksp);
307: KSPDestroy(&ctx->ksp);
308: ctx->ksp = ksp;
309: PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->ksp);
310: nep->state = NEP_STATE_INITIAL;
311: PetscFunctionReturn(0);
312: }
314: /*@
315: NEPNArnoldiSetKSP - Associate a linear solver object (KSP) to the nonlinear
316: eigenvalue solver.
318: Collective on nep
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: {
333: PetscTryMethod(nep,"NEPNArnoldiSetKSP_C",(NEP,KSP),(nep,ksp));
334: PetscFunctionReturn(0);
335: }
337: static PetscErrorCode NEPNArnoldiGetKSP_NArnoldi(NEP nep,KSP *ksp)
338: {
339: NEP_NARNOLDI *ctx = (NEP_NARNOLDI*)nep->data;
341: if (!ctx->ksp) {
342: KSPCreate(PetscObjectComm((PetscObject)nep),&ctx->ksp);
343: PetscObjectIncrementTabLevel((PetscObject)ctx->ksp,(PetscObject)nep,1);
344: KSPSetOptionsPrefix(ctx->ksp,((PetscObject)nep)->prefix);
345: KSPAppendOptionsPrefix(ctx->ksp,"nep_narnoldi_");
346: PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->ksp);
347: PetscObjectSetOptions((PetscObject)ctx->ksp,((PetscObject)nep)->options);
348: KSPSetErrorIfNotConverged(ctx->ksp,PETSC_TRUE);
349: KSPSetTolerances(ctx->ksp,SlepcDefaultTol(nep->tol),PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
350: }
351: *ksp = ctx->ksp;
352: PetscFunctionReturn(0);
353: }
355: /*@
356: NEPNArnoldiGetKSP - Retrieve the linear solver object (KSP) associated with
357: the nonlinear eigenvalue solver.
359: Not Collective
361: Input Parameter:
362: . nep - nonlinear eigenvalue solver
364: Output Parameter:
365: . ksp - the linear solver object
367: Level: advanced
369: .seealso: NEPNArnoldiSetKSP()
370: @*/
371: PetscErrorCode NEPNArnoldiGetKSP(NEP nep,KSP *ksp)
372: {
375: PetscUseMethod(nep,"NEPNArnoldiGetKSP_C",(NEP,KSP*),(nep,ksp));
376: PetscFunctionReturn(0);
377: }
379: PetscErrorCode NEPView_NArnoldi(NEP nep,PetscViewer viewer)
380: {
381: NEP_NARNOLDI *ctx = (NEP_NARNOLDI*)nep->data;
382: PetscBool isascii;
384: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
385: if (isascii) {
386: if (ctx->lag) PetscViewerASCIIPrintf(viewer," updating the preconditioner every %" PetscInt_FMT " iterations\n",ctx->lag);
387: if (!ctx->ksp) NEPNArnoldiGetKSP(nep,&ctx->ksp);
388: PetscViewerASCIIPushTab(viewer);
389: KSPView(ctx->ksp,viewer);
390: PetscViewerASCIIPopTab(viewer);
391: }
392: PetscFunctionReturn(0);
393: }
395: PetscErrorCode NEPReset_NArnoldi(NEP nep)
396: {
397: NEP_NARNOLDI *ctx = (NEP_NARNOLDI*)nep->data;
399: KSPReset(ctx->ksp);
400: PetscFunctionReturn(0);
401: }
403: PetscErrorCode NEPDestroy_NArnoldi(NEP nep)
404: {
405: NEP_NARNOLDI *ctx = (NEP_NARNOLDI*)nep->data;
407: KSPDestroy(&ctx->ksp);
408: PetscFree(nep->data);
409: PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiSetLagPreconditioner_C",NULL);
410: PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiGetLagPreconditioner_C",NULL);
411: PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiSetKSP_C",NULL);
412: PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiGetKSP_C",NULL);
413: PetscFunctionReturn(0);
414: }
416: SLEPC_EXTERN PetscErrorCode NEPCreate_NArnoldi(NEP nep)
417: {
418: NEP_NARNOLDI *ctx;
420: PetscNewLog(nep,&ctx);
421: nep->data = (void*)ctx;
422: ctx->lag = 1;
424: nep->useds = PETSC_TRUE;
426: nep->ops->solve = NEPSolve_NArnoldi;
427: nep->ops->setup = NEPSetUp_NArnoldi;
428: nep->ops->setfromoptions = NEPSetFromOptions_NArnoldi;
429: nep->ops->reset = NEPReset_NArnoldi;
430: nep->ops->destroy = NEPDestroy_NArnoldi;
431: nep->ops->view = NEPView_NArnoldi;
432: nep->ops->computevectors = NEPComputeVectors_Schur;
434: PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiSetLagPreconditioner_C",NEPNArnoldiSetLagPreconditioner_NArnoldi);
435: PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiGetLagPreconditioner_C",NEPNArnoldiGetLagPreconditioner_NArnoldi);
436: PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiSetKSP_C",NEPNArnoldiSetKSP_NArnoldi);
437: PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiGetKSP_C",NEPNArnoldiGetKSP_NArnoldi);
438: PetscFunctionReturn(0);
439: }