Actual source code: arnoldi.c
slepc-3.21.0 2024-03-30
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: PetscCall(EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd));
36: 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");
37: if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
38: if (!eps->which) PetscCall(EPSSetWhichEigenpairs_Default(eps));
39: PetscCheck(eps->which!=EPS_ALL,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support computing all eigenvalues");
40: EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_TWOSIDED);
42: PetscCall(EPSAllocateSolution(eps,1));
43: PetscCall(EPS_SetInnerProduct(eps));
44: PetscCall(DSSetType(eps->ds,DSNHEP));
45: if (eps->extraction==EPS_REFINED || eps->extraction==EPS_REFINED_HARMONIC) PetscCall(DSSetRefined(eps->ds,PETSC_TRUE));
46: PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
47: PetscCall(DSAllocate(eps->ds,eps->ncv+1));
48: PetscFunctionReturn(PETSC_SUCCESS);
49: }
51: static PetscErrorCode EPSSolve_Arnoldi(EPS eps)
52: {
53: PetscInt k,nv,ld;
54: Mat U,Op,H;
55: PetscScalar *Harray;
56: PetscReal beta,gamma=1.0;
57: PetscBool breakdown,harmonic,refined;
58: BVOrthogRefineType orthog_ref;
59: EPS_ARNOLDI *arnoldi = (EPS_ARNOLDI*)eps->data;
61: PetscFunctionBegin;
62: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
63: PetscCall(DSGetRefined(eps->ds,&refined));
64: harmonic = (eps->extraction==EPS_HARMONIC || eps->extraction==EPS_REFINED_HARMONIC)?PETSC_TRUE:PETSC_FALSE;
65: PetscCall(BVGetOrthogonalization(eps->V,NULL,&orthog_ref,NULL,NULL));
67: /* Get the starting Arnoldi vector */
68: PetscCall(EPSGetStartVector(eps,0,NULL));
70: /* Restart loop */
71: while (eps->reason == EPS_CONVERGED_ITERATING) {
72: eps->its++;
74: /* Compute an nv-step Arnoldi factorization */
75: nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
76: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,0));
77: if (!arnoldi->delayed) {
78: PetscCall(STGetOperator(eps->st,&Op));
79: PetscCall(DSGetMat(eps->ds,DS_MAT_A,&H));
80: PetscCall(BVMatArnoldi(eps->V,Op,H,eps->nconv,&nv,&beta,&breakdown));
81: PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&H));
82: PetscCall(STRestoreOperator(eps->st,&Op));
83: } else if (orthog_ref == BV_ORTHOG_REFINE_NEVER) {
84: PetscCall(DSGetArray(eps->ds,DS_MAT_A,&Harray));
85: PetscCall(EPSDelayedArnoldi1(eps,Harray,ld,eps->nconv,&nv,&beta,&breakdown));
86: PetscCall(DSRestoreArray(eps->ds,DS_MAT_A,&Harray));
87: } else {
88: PetscCall(DSGetArray(eps->ds,DS_MAT_A,&Harray));
89: PetscCall(EPSDelayedArnoldi(eps,Harray,ld,eps->nconv,&nv,&beta,&breakdown));
90: PetscCall(DSRestoreArray(eps->ds,DS_MAT_A,&Harray));
91: }
92: PetscCall(DSSetState(eps->ds,DS_STATE_INTERMEDIATE));
93: PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
95: /* Compute translation of Krylov decomposition if harmonic extraction used */
96: if (harmonic) PetscCall(DSTranslateHarmonic(eps->ds,eps->target,beta,PETSC_FALSE,NULL,&gamma));
98: /* Solve projected problem */
99: PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
100: PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
101: PetscCall(DSUpdateExtraRow(eps->ds));
102: PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
104: /* Check convergence */
105: PetscCall(EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta,0.0,gamma,&k));
106: if (refined) {
107: PetscCall(DSGetMat(eps->ds,DS_MAT_X,&U));
108: PetscCall(BVMultInPlace(eps->V,U,eps->nconv,k+1));
109: PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&U));
110: PetscCall(BVOrthonormalizeColumn(eps->V,k,PETSC_FALSE,NULL,NULL));
111: } else {
112: PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&U));
113: PetscCall(BVMultInPlace(eps->V,U,eps->nconv,PetscMin(k+1,nv)));
114: PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&U));
115: }
116: PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
117: if (eps->reason == EPS_CONVERGED_ITERATING && breakdown) {
118: PetscCall(PetscInfo(eps,"Breakdown in Arnoldi method (it=%" PetscInt_FMT " norm=%g)\n",eps->its,(double)beta));
119: PetscCall(EPSGetStartVector(eps,k,&breakdown));
120: if (breakdown) {
121: eps->reason = EPS_DIVERGED_BREAKDOWN;
122: PetscCall(PetscInfo(eps,"Unable to generate more start vectors\n"));
123: }
124: }
125: eps->nconv = k;
126: PetscCall(EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,nv));
127: }
128: PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
129: PetscFunctionReturn(PETSC_SUCCESS);
130: }
132: static PetscErrorCode EPSSetFromOptions_Arnoldi(EPS eps,PetscOptionItems *PetscOptionsObject)
133: {
134: PetscBool set,val;
135: EPS_ARNOLDI *arnoldi = (EPS_ARNOLDI*)eps->data;
137: PetscFunctionBegin;
138: PetscOptionsHeadBegin(PetscOptionsObject,"EPS Arnoldi Options");
140: PetscCall(PetscOptionsBool("-eps_arnoldi_delayed","Use delayed reorthogonalization","EPSArnoldiSetDelayed",arnoldi->delayed,&val,&set));
141: if (set) PetscCall(EPSArnoldiSetDelayed(eps,val));
143: PetscOptionsHeadEnd();
144: PetscFunctionReturn(PETSC_SUCCESS);
145: }
147: static PetscErrorCode EPSArnoldiSetDelayed_Arnoldi(EPS eps,PetscBool delayed)
148: {
149: EPS_ARNOLDI *arnoldi = (EPS_ARNOLDI*)eps->data;
151: PetscFunctionBegin;
152: arnoldi->delayed = delayed;
153: PetscFunctionReturn(PETSC_SUCCESS);
154: }
156: /*@
157: EPSArnoldiSetDelayed - Activates or deactivates delayed reorthogonalization
158: in the Arnoldi iteration.
160: Logically Collective
162: Input Parameters:
163: + eps - the eigenproblem solver context
164: - delayed - boolean flag
166: Options Database Key:
167: . -eps_arnoldi_delayed - Activates delayed reorthogonalization in Arnoldi
169: Note:
170: Delayed reorthogonalization is an aggressive optimization for the Arnoldi
171: eigensolver than may provide better scalability, but sometimes makes the
172: solver converge less than the default algorithm.
174: Level: advanced
176: .seealso: EPSArnoldiGetDelayed()
177: @*/
178: PetscErrorCode EPSArnoldiSetDelayed(EPS eps,PetscBool delayed)
179: {
180: PetscFunctionBegin;
183: PetscTryMethod(eps,"EPSArnoldiSetDelayed_C",(EPS,PetscBool),(eps,delayed));
184: PetscFunctionReturn(PETSC_SUCCESS);
185: }
187: static PetscErrorCode EPSArnoldiGetDelayed_Arnoldi(EPS eps,PetscBool *delayed)
188: {
189: EPS_ARNOLDI *arnoldi = (EPS_ARNOLDI*)eps->data;
191: PetscFunctionBegin;
192: *delayed = arnoldi->delayed;
193: PetscFunctionReturn(PETSC_SUCCESS);
194: }
196: /*@
197: EPSArnoldiGetDelayed - Gets the type of reorthogonalization used during the Arnoldi
198: iteration.
200: Not Collective
202: Input Parameter:
203: . eps - the eigenproblem solver context
205: Output Parameter:
206: . delayed - boolean flag indicating if delayed reorthogonalization has been enabled
208: Level: advanced
210: .seealso: EPSArnoldiSetDelayed()
211: @*/
212: PetscErrorCode EPSArnoldiGetDelayed(EPS eps,PetscBool *delayed)
213: {
214: PetscFunctionBegin;
216: PetscAssertPointer(delayed,2);
217: PetscUseMethod(eps,"EPSArnoldiGetDelayed_C",(EPS,PetscBool*),(eps,delayed));
218: PetscFunctionReturn(PETSC_SUCCESS);
219: }
221: static PetscErrorCode EPSDestroy_Arnoldi(EPS eps)
222: {
223: PetscFunctionBegin;
224: PetscCall(PetscFree(eps->data));
225: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSArnoldiSetDelayed_C",NULL));
226: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSArnoldiGetDelayed_C",NULL));
227: PetscFunctionReturn(PETSC_SUCCESS);
228: }
230: static PetscErrorCode EPSView_Arnoldi(EPS eps,PetscViewer viewer)
231: {
232: PetscBool isascii;
233: EPS_ARNOLDI *arnoldi = (EPS_ARNOLDI*)eps->data;
235: PetscFunctionBegin;
236: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
237: if (isascii && arnoldi->delayed) PetscCall(PetscViewerASCIIPrintf(viewer," using delayed reorthogonalization\n"));
238: PetscFunctionReturn(PETSC_SUCCESS);
239: }
241: SLEPC_EXTERN PetscErrorCode EPSCreate_Arnoldi(EPS eps)
242: {
243: EPS_ARNOLDI *ctx;
245: PetscFunctionBegin;
246: PetscCall(PetscNew(&ctx));
247: eps->data = (void*)ctx;
249: eps->useds = PETSC_TRUE;
251: eps->ops->solve = EPSSolve_Arnoldi;
252: eps->ops->setup = EPSSetUp_Arnoldi;
253: eps->ops->setupsort = EPSSetUpSort_Default;
254: eps->ops->setfromoptions = EPSSetFromOptions_Arnoldi;
255: eps->ops->destroy = EPSDestroy_Arnoldi;
256: eps->ops->view = EPSView_Arnoldi;
257: eps->ops->backtransform = EPSBackTransform_Default;
258: eps->ops->computevectors = EPSComputeVectors_Schur;
260: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSArnoldiSetDelayed_C",EPSArnoldiSetDelayed_Arnoldi));
261: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSArnoldiGetDelayed_C",EPSArnoldiGetDelayed_Arnoldi));
262: PetscFunctionReturn(PETSC_SUCCESS);
263: }