Actual source code: arnoldi.c
slepc-3.22.1 2024-10-28
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 eigensolver: "arnoldi"
13: Method: Explicitly Restarted Arnoldi
15: Algorithm:
17: Arnoldi method with explicit restart and deflation.
19: References:
21: [1] "Arnoldi Methods in SLEPc", SLEPc Technical Report STR-4,
22: available at https://slepc.upv.es.
23: */
25: #include <slepc/private/epsimpl.h>
27: typedef struct {
28: PetscBool delayed;
29: } EPS_ARNOLDI;
31: static PetscErrorCode EPSSetUp_Arnoldi(EPS eps)
32: {
33: PetscFunctionBegin;
34: EPSCheckDefinite(eps);
35: EPSCheckNotStructured(eps);
36: PetscCall(EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd));
37: PetscCheck(eps->ncv<=eps->nev+eps->mpd,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must not be larger than nev+mpd");
38: if (eps->max_it==PETSC_DETERMINE) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
39: if (!eps->which) PetscCall(EPSSetWhichEigenpairs_Default(eps));
40: PetscCheck(eps->which!=EPS_ALL,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support computing all eigenvalues");
41: EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_TWOSIDED);
43: PetscCall(EPSAllocateSolution(eps,1));
44: PetscCall(EPS_SetInnerProduct(eps));
45: PetscCall(DSSetType(eps->ds,DSNHEP));
46: if (eps->extraction==EPS_REFINED || eps->extraction==EPS_REFINED_HARMONIC) PetscCall(DSSetRefined(eps->ds,PETSC_TRUE));
47: PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
48: PetscCall(DSAllocate(eps->ds,eps->ncv+1));
49: PetscFunctionReturn(PETSC_SUCCESS);
50: }
52: static PetscErrorCode EPSSolve_Arnoldi(EPS eps)
53: {
54: PetscInt k,nv,ld;
55: Mat U,Op,H;
56: PetscScalar *Harray;
57: PetscReal beta,gamma=1.0;
58: PetscBool breakdown,harmonic,refined;
59: BVOrthogRefineType orthog_ref;
60: EPS_ARNOLDI *arnoldi = (EPS_ARNOLDI*)eps->data;
62: PetscFunctionBegin;
63: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
64: PetscCall(DSGetRefined(eps->ds,&refined));
65: harmonic = (eps->extraction==EPS_HARMONIC || eps->extraction==EPS_REFINED_HARMONIC)?PETSC_TRUE:PETSC_FALSE;
66: PetscCall(BVGetOrthogonalization(eps->V,NULL,&orthog_ref,NULL,NULL));
68: /* Get the starting Arnoldi vector */
69: PetscCall(EPSGetStartVector(eps,0,NULL));
71: /* Restart loop */
72: while (eps->reason == EPS_CONVERGED_ITERATING) {
73: eps->its++;
75: /* Compute an nv-step Arnoldi factorization */
76: nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
77: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,0));
78: if (!arnoldi->delayed) {
79: PetscCall(STGetOperator(eps->st,&Op));
80: PetscCall(DSGetMat(eps->ds,DS_MAT_A,&H));
81: PetscCall(BVMatArnoldi(eps->V,Op,H,eps->nconv,&nv,&beta,&breakdown));
82: PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&H));
83: PetscCall(STRestoreOperator(eps->st,&Op));
84: } else if (orthog_ref == BV_ORTHOG_REFINE_NEVER) {
85: PetscCall(DSGetArray(eps->ds,DS_MAT_A,&Harray));
86: PetscCall(EPSDelayedArnoldi1(eps,Harray,ld,eps->nconv,&nv,&beta,&breakdown));
87: PetscCall(DSRestoreArray(eps->ds,DS_MAT_A,&Harray));
88: } else {
89: PetscCall(DSGetArray(eps->ds,DS_MAT_A,&Harray));
90: PetscCall(EPSDelayedArnoldi(eps,Harray,ld,eps->nconv,&nv,&beta,&breakdown));
91: PetscCall(DSRestoreArray(eps->ds,DS_MAT_A,&Harray));
92: }
93: PetscCall(DSSetState(eps->ds,DS_STATE_INTERMEDIATE));
94: PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
96: /* Compute translation of Krylov decomposition if harmonic extraction used */
97: if (harmonic) PetscCall(DSTranslateHarmonic(eps->ds,eps->target,beta,PETSC_FALSE,NULL,&gamma));
99: /* Solve projected problem */
100: PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
101: PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
102: PetscCall(DSUpdateExtraRow(eps->ds));
103: PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
105: /* Check convergence */
106: PetscCall(EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta,0.0,gamma,&k));
107: if (refined) {
108: PetscCall(DSGetMat(eps->ds,DS_MAT_X,&U));
109: PetscCall(BVMultInPlace(eps->V,U,eps->nconv,k+1));
110: PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&U));
111: PetscCall(BVOrthonormalizeColumn(eps->V,k,PETSC_FALSE,NULL,NULL));
112: } else {
113: PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&U));
114: PetscCall(BVMultInPlace(eps->V,U,eps->nconv,PetscMin(k+1,nv)));
115: PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&U));
116: }
117: PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
118: if (eps->reason == EPS_CONVERGED_ITERATING && breakdown) {
119: PetscCall(PetscInfo(eps,"Breakdown in Arnoldi method (it=%" PetscInt_FMT " norm=%g)\n",eps->its,(double)beta));
120: PetscCall(EPSGetStartVector(eps,k,&breakdown));
121: if (breakdown) {
122: eps->reason = EPS_DIVERGED_BREAKDOWN;
123: PetscCall(PetscInfo(eps,"Unable to generate more start vectors\n"));
124: }
125: }
126: eps->nconv = k;
127: PetscCall(EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,nv));
128: }
129: PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
130: PetscFunctionReturn(PETSC_SUCCESS);
131: }
133: static PetscErrorCode EPSSetFromOptions_Arnoldi(EPS eps,PetscOptionItems *PetscOptionsObject)
134: {
135: PetscBool set,val;
136: EPS_ARNOLDI *arnoldi = (EPS_ARNOLDI*)eps->data;
138: PetscFunctionBegin;
139: PetscOptionsHeadBegin(PetscOptionsObject,"EPS Arnoldi Options");
141: PetscCall(PetscOptionsBool("-eps_arnoldi_delayed","Use delayed reorthogonalization","EPSArnoldiSetDelayed",arnoldi->delayed,&val,&set));
142: if (set) PetscCall(EPSArnoldiSetDelayed(eps,val));
144: PetscOptionsHeadEnd();
145: PetscFunctionReturn(PETSC_SUCCESS);
146: }
148: static PetscErrorCode EPSArnoldiSetDelayed_Arnoldi(EPS eps,PetscBool delayed)
149: {
150: EPS_ARNOLDI *arnoldi = (EPS_ARNOLDI*)eps->data;
152: PetscFunctionBegin;
153: arnoldi->delayed = delayed;
154: PetscFunctionReturn(PETSC_SUCCESS);
155: }
157: /*@
158: EPSArnoldiSetDelayed - Activates or deactivates delayed reorthogonalization
159: in the Arnoldi iteration.
161: Logically Collective
163: Input Parameters:
164: + eps - the eigenproblem solver context
165: - delayed - boolean flag
167: Options Database Key:
168: . -eps_arnoldi_delayed - Activates delayed reorthogonalization in Arnoldi
170: Note:
171: Delayed reorthogonalization is an aggressive optimization for the Arnoldi
172: eigensolver than may provide better scalability, but sometimes makes the
173: solver converge less than the default algorithm.
175: Level: advanced
177: .seealso: EPSArnoldiGetDelayed()
178: @*/
179: PetscErrorCode EPSArnoldiSetDelayed(EPS eps,PetscBool delayed)
180: {
181: PetscFunctionBegin;
184: PetscTryMethod(eps,"EPSArnoldiSetDelayed_C",(EPS,PetscBool),(eps,delayed));
185: PetscFunctionReturn(PETSC_SUCCESS);
186: }
188: static PetscErrorCode EPSArnoldiGetDelayed_Arnoldi(EPS eps,PetscBool *delayed)
189: {
190: EPS_ARNOLDI *arnoldi = (EPS_ARNOLDI*)eps->data;
192: PetscFunctionBegin;
193: *delayed = arnoldi->delayed;
194: PetscFunctionReturn(PETSC_SUCCESS);
195: }
197: /*@
198: EPSArnoldiGetDelayed - Gets the type of reorthogonalization used during the Arnoldi
199: iteration.
201: Not Collective
203: Input Parameter:
204: . eps - the eigenproblem solver context
206: Output Parameter:
207: . delayed - boolean flag indicating if delayed reorthogonalization has been enabled
209: Level: advanced
211: .seealso: EPSArnoldiSetDelayed()
212: @*/
213: PetscErrorCode EPSArnoldiGetDelayed(EPS eps,PetscBool *delayed)
214: {
215: PetscFunctionBegin;
217: PetscAssertPointer(delayed,2);
218: PetscUseMethod(eps,"EPSArnoldiGetDelayed_C",(EPS,PetscBool*),(eps,delayed));
219: PetscFunctionReturn(PETSC_SUCCESS);
220: }
222: static PetscErrorCode EPSDestroy_Arnoldi(EPS eps)
223: {
224: PetscFunctionBegin;
225: PetscCall(PetscFree(eps->data));
226: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSArnoldiSetDelayed_C",NULL));
227: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSArnoldiGetDelayed_C",NULL));
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }
231: static PetscErrorCode EPSView_Arnoldi(EPS eps,PetscViewer viewer)
232: {
233: PetscBool isascii;
234: EPS_ARNOLDI *arnoldi = (EPS_ARNOLDI*)eps->data;
236: PetscFunctionBegin;
237: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
238: if (isascii && arnoldi->delayed) PetscCall(PetscViewerASCIIPrintf(viewer," using delayed reorthogonalization\n"));
239: PetscFunctionReturn(PETSC_SUCCESS);
240: }
242: SLEPC_EXTERN PetscErrorCode EPSCreate_Arnoldi(EPS eps)
243: {
244: EPS_ARNOLDI *ctx;
246: PetscFunctionBegin;
247: PetscCall(PetscNew(&ctx));
248: eps->data = (void*)ctx;
250: eps->useds = PETSC_TRUE;
252: eps->ops->solve = EPSSolve_Arnoldi;
253: eps->ops->setup = EPSSetUp_Arnoldi;
254: eps->ops->setupsort = EPSSetUpSort_Default;
255: eps->ops->setfromoptions = EPSSetFromOptions_Arnoldi;
256: eps->ops->destroy = EPSDestroy_Arnoldi;
257: eps->ops->view = EPSView_Arnoldi;
258: eps->ops->backtransform = EPSBackTransform_Default;
259: eps->ops->computevectors = EPSComputeVectors_Schur;
261: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSArnoldiSetDelayed_C",EPSArnoldiSetDelayed_Arnoldi));
262: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSArnoldiGetDelayed_C",EPSArnoldiGetDelayed_Arnoldi));
263: PetscFunctionReturn(PETSC_SUCCESS);
264: }