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