Actual source code: narnoldi.c
slepc-3.11.2 2019-07-30
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: "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> /*I "slepcnep.h" I*/
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: {
36: PetscBool istrivial;
39: NEPSetDimensions_Default(nep,nep->nev,&nep->ncv,&nep->mpd);
40: if (nep->ncv>nep->nev+nep->mpd) SETERRQ(PetscObjectComm((PetscObject)nep),1,"The value of ncv must not be larger than nev+mpd");
41: if (!nep->max_it) nep->max_it = nep->nev*nep->ncv;
42: if (nep->which && nep->which!=NEP_TARGET_MAGNITUDE) SETERRQ(PetscObjectComm((PetscObject)nep),1,"Wrong value of which");
43: if (nep->fui!=NEP_USER_INTERFACE_SPLIT) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"NARNOLDI only available for split operator");
45: RGIsTrivial(nep->rg,&istrivial);
46: if (!istrivial) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"This solver does not support region filtering");
48: NEPAllocateSolution(nep,0);
49: NEPSetWorkVecs(nep,3);
50: return(0);
51: }
53: PetscErrorCode NEPSolve_NArnoldi(NEP nep)
54: {
55: PetscErrorCode ierr;
56: NEP_NARNOLDI *ctx = (NEP_NARNOLDI*)nep->data;
57: Mat T,H;
58: Vec f,r,u,uu;
59: PetscScalar *X,lambda=0.0,lambda2=0.0,*eigr,*Hp,*Ap,sigma;
60: PetscReal beta,resnorm=0.0,nrm,perr=0.0;
61: PetscInt n,i,j,ldds,ldh;
62: PetscBool breakdown,skip=PETSC_FALSE;
63: BV Vext;
64: DS ds;
65: NEP_EXT_OP extop=NULL;
66: SlepcSC sc;
67: KSPConvergedReason kspreason;
70: /* get initial space and shift */
71: NEPGetDefaultShift(nep,&sigma);
72: if (!nep->nini) {
73: BVSetRandomColumn(nep->V,0);
74: BVNormColumn(nep->V,0,NORM_2,&nrm);
75: BVScaleColumn(nep->V,0,1.0/nrm);
76: n = 1;
77: } else n = nep->nini;
79: if (!ctx->ksp) { NEPNArnoldiGetKSP(nep,&ctx->ksp); }
80: NEPDeflationInitialize(nep,nep->V,ctx->ksp,PETSC_FALSE,nep->nev,&extop);
81: NEPDeflationCreateBV(extop,nep->ncv,&Vext);
83: /* prepare linear solver */
84: NEPDeflationSolveSetUp(extop,sigma);
86: BVGetColumn(Vext,0,&f);
87: VecDuplicate(f,&r);
88: VecDuplicate(f,&u);
89: BVGetColumn(nep->V,0,&uu);
90: NEPDeflationCopyToExtendedVec(extop,uu,NULL,f,PETSC_FALSE);
91: BVRestoreColumn(nep->V,0,&uu);
92: VecCopy(f,r);
93: NEPDeflationFunctionSolve(extop,r,f);
94: VecNorm(f,NORM_2,&nrm);
95: VecScale(f,1.0/nrm);
96: BVRestoreColumn(Vext,0,&f);
98: DSCreate(PetscObjectComm((PetscObject)nep),&ds);
99: PetscLogObjectParent((PetscObject)nep,(PetscObject)ds);
100: DSSetType(ds,DSNEP);
101: DSNEPSetFN(ds,nep->nt,nep->f);
102: DSAllocate(ds,nep->ncv);
103: DSGetSlepcSC(ds,&sc);
104: sc->comparison = nep->sc->comparison;
105: sc->comparisonctx = nep->sc->comparisonctx;
106: DSSetFromOptions(ds);
108: /* build projected matrices for initial space */
109: DSSetDimensions(ds,n,0,0,0);
110: NEPDeflationProjectOperator(extop,Vext,ds,0,n);
112: PetscMalloc1(nep->ncv,&eigr);
114: /* Restart loop */
115: while (nep->reason == NEP_CONVERGED_ITERATING) {
116: nep->its++;
118: /* solve projected problem */
119: DSSetDimensions(ds,n,0,0,0);
120: DSSetState(ds,DS_STATE_RAW);
121: DSSolve(ds,eigr,NULL);
122: DSSynchronize(ds,eigr,NULL);
123: if (nep->its>1) lambda2 = lambda;
124: lambda = eigr[0];
125: nep->eigr[nep->nconv] = lambda;
127: /* compute Ritz vector, x = V*s */
128: DSGetArray(ds,DS_MAT_X,&X);
129: BVSetActiveColumns(Vext,0,n);
130: BVMultVec(Vext,1.0,0.0,u,X);
131: DSRestoreArray(ds,DS_MAT_X,&X);
133: /* compute the residual, r = T(lambda)*x */
134: NEPDeflationComputeFunction(extop,lambda,&T);
135: MatMult(T,u,r);
137: /* convergence test */
138: VecNorm(r,NORM_2,&resnorm);
139: if (nep->its>1) perr=nep->errest[nep->nconv];
140: (*nep->converged)(nep,lambda,0,resnorm,&nep->errest[nep->nconv],nep->convergedctx);
141: if (nep->errest[nep->nconv]<=nep->tol) {
142: nep->nconv = nep->nconv + 1;
143: NEPDeflationLocking(extop,u,lambda);
144: skip = PETSC_TRUE;
145: }
146: (*nep->stopping)(nep,nep->its,nep->max_it,nep->nconv,nep->nev,&nep->reason,nep->stoppingctx);
147: if (!skip || nep->reason>0) {
148: NEPMonitor(nep,nep->its,nep->nconv,nep->eigr,nep->eigi,nep->errest,(nep->reason>0)?nep->nconv:nep->nconv+1);
149: }
151: if (nep->reason == NEP_CONVERGED_ITERATING) {
152: if (!skip) {
153: if (n>=nep->ncv) {
154: nep->reason = NEP_DIVERGED_SUBSPACE_EXHAUSTED;
155: break;
156: }
157: if (ctx->lag && !(nep->its%ctx->lag) && nep->its>=2*ctx->lag && perr && nep->errest[nep->nconv]>.5*perr) {
158: NEPDeflationSolveSetUp(extop,lambda2);
159: }
161: /* continuation vector: f = T(sigma)\r */
162: BVGetColumn(Vext,n,&f);
163: NEPDeflationFunctionSolve(extop,r,f);
164: BVRestoreColumn(Vext,n,&f);
165: KSPGetConvergedReason(ctx->ksp,&kspreason);
166: if (kspreason<0) {
167: PetscInfo1(nep,"iter=%D, linear solve failed, stopping solve\n",nep->its);
168: nep->reason = NEP_DIVERGED_LINEAR_SOLVE;
169: break;
170: }
172: /* orthonormalize */
173: BVOrthonormalizeColumn(Vext,n,PETSC_FALSE,&beta,&breakdown);
174: if (breakdown || beta==0.0) {
175: PetscInfo1(nep,"iter=%D, orthogonalization failed, stopping solve\n",nep->its);
176: nep->reason = NEP_DIVERGED_BREAKDOWN;
177: break;
178: }
180: /* update projected matrices */
181: DSSetDimensions(ds,n+1,0,0,0);
182: NEPDeflationProjectOperator(extop,Vext,ds,n,n+1);
183: n++;
184: } else {
185: nep->its--; /* do not count this as a full iteration */
186: BVGetColumn(Vext,0,&f);
187: NEPDeflationSetRandomVec(extop,f);
188: NEPDeflationSolveSetUp(extop,sigma);
189: VecCopy(f,r);
190: NEPDeflationFunctionSolve(extop,r,f);
191: VecNorm(f,NORM_2,&nrm);
192: VecScale(f,1.0/nrm);
193: BVRestoreColumn(Vext,0,&f);
194: n = 1;
195: DSSetDimensions(ds,n,0,0,0);
196: NEPDeflationProjectOperator(extop,Vext,ds,n-1,n);
197: skip = PETSC_FALSE;
198: }
199: }
200: }
202: NEPDeflationGetInvariantPair(extop,NULL,&H);
203: MatGetSize(H,NULL,&ldh);
204: DSSetType(nep->ds,DSNHEP);
205: DSAllocate(nep->ds,PetscMax(nep->nconv,1));
206: DSGetLeadingDimension(nep->ds,&ldds);
207: MatDenseGetArray(H,&Hp);
208: DSGetArray(nep->ds,DS_MAT_A,&Ap);
209: for (j=0;j<nep->nconv;j++)
210: for (i=0;i<nep->nconv;i++) Ap[j*ldds+i] = Hp[j*ldh+i];
211: DSRestoreArray(nep->ds,DS_MAT_A,&Ap);
212: MatDenseRestoreArray(H,&Hp);
213: MatDestroy(&H);
214: DSSetDimensions(nep->ds,nep->nconv,0,0,nep->nconv);
215: DSSolve(nep->ds,nep->eigr,nep->eigi);
216: NEPDeflationReset(extop);
217: VecDestroy(&u);
218: VecDestroy(&r);
219: BVDestroy(&Vext);
220: DSDestroy(&ds);
221: PetscFree(eigr);
222: return(0);
223: }
225: static PetscErrorCode NEPNArnoldiSetLagPreconditioner_NArnoldi(NEP nep,PetscInt lag)
226: {
227: NEP_NARNOLDI *ctx = (NEP_NARNOLDI*)nep->data;
230: if (lag<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be non-negative");
231: ctx->lag = lag;
232: return(0);
233: }
235: /*@
236: NEPNArnoldiSetLagPreconditioner - Determines when the preconditioner is rebuilt in the
237: nonlinear solve.
239: Logically Collective on NEP
241: Input Parameters:
242: + nep - nonlinear eigenvalue solver
243: - lag - 0 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is
244: computed within the nonlinear iteration, 2 means every second time
245: the Jacobian is built, etc.
247: Options Database Keys:
248: . -nep_narnoldi_lag_preconditioner <lag>
250: Notes:
251: The default is 1.
252: The preconditioner is ALWAYS built in the first iteration of a nonlinear solve.
254: Level: intermediate
256: .seealso: NEPNArnoldiGetLagPreconditioner()
257: @*/
258: PetscErrorCode NEPNArnoldiSetLagPreconditioner(NEP nep,PetscInt lag)
259: {
265: PetscTryMethod(nep,"NEPNArnoldiSetLagPreconditioner_C",(NEP,PetscInt),(nep,lag));
266: return(0);
267: }
269: static PetscErrorCode NEPNArnoldiGetLagPreconditioner_NArnoldi(NEP nep,PetscInt *lag)
270: {
271: NEP_NARNOLDI *ctx = (NEP_NARNOLDI*)nep->data;
274: *lag = ctx->lag;
275: return(0);
276: }
278: /*@
279: NEPNArnoldiGetLagPreconditioner - Indicates how often the preconditioner is rebuilt.
281: Not Collective
283: Input Parameter:
284: . nep - nonlinear eigenvalue solver
286: Output Parameter:
287: . lag - the lag parameter
289: Level: intermediate
291: .seealso: NEPNArnoldiSetLagPreconditioner()
292: @*/
293: PetscErrorCode NEPNArnoldiGetLagPreconditioner(NEP nep,PetscInt *lag)
294: {
300: PetscUseMethod(nep,"NEPNArnoldiGetLagPreconditioner_C",(NEP,PetscInt*),(nep,lag));
301: return(0);
302: }
304: PetscErrorCode NEPSetFromOptions_NArnoldi(PetscOptionItems *PetscOptionsObject,NEP nep)
305: {
307: PetscInt i;
308: PetscBool flg;
309: NEP_NARNOLDI *ctx = (NEP_NARNOLDI*)nep->data;
312: PetscOptionsHead(PetscOptionsObject,"NEP N-Arnoldi Options");
313: i = 0;
314: PetscOptionsInt("-nep_narnoldi_lag_preconditioner","Interval to rebuild preconditioner","NEPNArnoldiSetLagPreconditioner",ctx->lag,&i,&flg);
315: if (flg) { NEPNArnoldiSetLagPreconditioner(nep,i); }
317: PetscOptionsTail();
319: if (!ctx->ksp) { NEPNArnoldiGetKSP(nep,&ctx->ksp); }
320: KSPSetFromOptions(ctx->ksp);
321: return(0);
322: }
324: static PetscErrorCode NEPNArnoldiSetKSP_NArnoldi(NEP nep,KSP ksp)
325: {
327: NEP_NARNOLDI *ctx = (NEP_NARNOLDI*)nep->data;
330: PetscObjectReference((PetscObject)ksp);
331: KSPDestroy(&ctx->ksp);
332: ctx->ksp = ksp;
333: PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->ksp);
334: nep->state = NEP_STATE_INITIAL;
335: return(0);
336: }
338: /*@
339: NEPNArnoldiSetKSP - Associate a linear solver object (KSP) to the nonlinear
340: eigenvalue solver.
342: Collective on NEP
344: Input Parameters:
345: + nep - eigenvalue solver
346: - ksp - the linear solver object
348: Level: advanced
350: .seealso: NEPNArnoldiGetKSP()
351: @*/
352: PetscErrorCode NEPNArnoldiSetKSP(NEP nep,KSP ksp)
353: {
360: PetscTryMethod(nep,"NEPNArnoldiSetKSP_C",(NEP,KSP),(nep,ksp));
361: return(0);
362: }
364: static PetscErrorCode NEPNArnoldiGetKSP_NArnoldi(NEP nep,KSP *ksp)
365: {
367: NEP_NARNOLDI *ctx = (NEP_NARNOLDI*)nep->data;
370: if (!ctx->ksp) {
371: KSPCreate(PetscObjectComm((PetscObject)nep),&ctx->ksp);
372: PetscObjectIncrementTabLevel((PetscObject)ctx->ksp,(PetscObject)nep,1);
373: KSPSetOptionsPrefix(ctx->ksp,((PetscObject)nep)->prefix);
374: KSPAppendOptionsPrefix(ctx->ksp,"nep_narnoldi_");
375: PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->ksp);
376: PetscObjectSetOptions((PetscObject)ctx->ksp,((PetscObject)nep)->options);
377: KSPSetErrorIfNotConverged(ctx->ksp,PETSC_TRUE);
378: KSPSetTolerances(ctx->ksp,SLEPC_DEFAULT_TOL,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
379: }
380: *ksp = ctx->ksp;
381: return(0);
382: }
384: /*@
385: NEPNArnoldiGetKSP - Retrieve the linear solver object (KSP) associated with
386: the nonlinear eigenvalue solver.
388: Not Collective
390: Input Parameter:
391: . nep - nonlinear eigenvalue solver
393: Output Parameter:
394: . ksp - the linear solver object
396: Level: advanced
398: .seealso: NEPNArnoldiSetKSP()
399: @*/
400: PetscErrorCode NEPNArnoldiGetKSP(NEP nep,KSP *ksp)
401: {
407: PetscUseMethod(nep,"NEPNArnoldiGetKSP_C",(NEP,KSP*),(nep,ksp));
408: return(0);
409: }
411: PetscErrorCode NEPView_NArnoldi(NEP nep,PetscViewer viewer)
412: {
414: NEP_NARNOLDI *ctx = (NEP_NARNOLDI*)nep->data;
415: PetscBool isascii;
418: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
419: if (isascii) {
420: if (ctx->lag) {
421: PetscViewerASCIIPrintf(viewer," updating the preconditioner every %D iterations\n",ctx->lag);
422: }
423: if (!ctx->ksp) { NEPNArnoldiGetKSP(nep,&ctx->ksp); }
424: PetscViewerASCIIPushTab(viewer);
425: KSPView(ctx->ksp,viewer);
426: PetscViewerASCIIPopTab(viewer);
427: }
428: return(0);
429: }
431: PetscErrorCode NEPReset_NArnoldi(NEP nep)
432: {
434: NEP_NARNOLDI *ctx = (NEP_NARNOLDI*)nep->data;
437: KSPReset(ctx->ksp);
438: return(0);
439: }
441: PetscErrorCode NEPDestroy_NArnoldi(NEP nep)
442: {
444: NEP_NARNOLDI *ctx = (NEP_NARNOLDI*)nep->data;
447: KSPDestroy(&ctx->ksp);
448: PetscFree(nep->data);
449: PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiSetLagPreconditioner_C",NULL);
450: PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiGetLagPreconditioner_C",NULL);
451: PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiSetKSP_C",NULL);
452: PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiGetKSP_C",NULL);
453: return(0);
454: }
456: SLEPC_EXTERN PetscErrorCode NEPCreate_NArnoldi(NEP nep)
457: {
459: NEP_NARNOLDI *ctx;
462: PetscNewLog(nep,&ctx);
463: nep->data = (void*)ctx;
464: ctx->lag = 1;
466: nep->useds = PETSC_TRUE;
468: nep->ops->solve = NEPSolve_NArnoldi;
469: nep->ops->setup = NEPSetUp_NArnoldi;
470: nep->ops->setfromoptions = NEPSetFromOptions_NArnoldi;
471: nep->ops->reset = NEPReset_NArnoldi;
472: nep->ops->destroy = NEPDestroy_NArnoldi;
473: nep->ops->view = NEPView_NArnoldi;
474: nep->ops->computevectors = NEPComputeVectors_Schur;
476: PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiSetLagPreconditioner_C",NEPNArnoldiSetLagPreconditioner_NArnoldi);
477: PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiGetLagPreconditioner_C",NEPNArnoldiGetLagPreconditioner_NArnoldi);
478: PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiSetKSP_C",NEPNArnoldiSetKSP_NArnoldi);
479: PetscObjectComposeFunction((PetscObject)nep,"NEPNArnoldiGetKSP_C",NEPNArnoldiGetKSP_NArnoldi);
480: return(0);
481: }