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