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: "subspace"
12 :
13 : Method: Subspace Iteration
14 :
15 : Algorithm:
16 :
17 : Subspace iteration with Rayleigh-Ritz projection and locking,
18 : based on the SRRIT implementation.
19 :
20 : References:
21 :
22 : [1] "Subspace Iteration in SLEPc", SLEPc Technical Report STR-3,
23 : available at https://slepc.upv.es.
24 : */
25 :
26 : #include <slepc/private/epsimpl.h>
27 :
28 : typedef struct {
29 : PetscBool estimatedrange; /* the filter range was not set by the user */
30 : } EPS_SUBSPACE;
31 :
32 1 : static PetscErrorCode EPSSetUp_Subspace_Filter(EPS eps)
33 : {
34 1 : EPS_SUBSPACE *ctx = (EPS_SUBSPACE*)eps->data;
35 1 : PetscBool estimaterange=PETSC_TRUE;
36 1 : PetscReal rleft,rright;
37 1 : Mat A;
38 :
39 1 : PetscFunctionBegin;
40 1 : EPSCheckHermitianCondition(eps,PETSC_TRUE," with polynomial filter");
41 1 : EPSCheckStandardCondition(eps,PETSC_TRUE," with polynomial filter");
42 1 : PetscCheck(eps->intb<PETSC_MAX_REAL || eps->inta>PETSC_MIN_REAL,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The defined computational interval should have at least one of their sides bounded");
43 1 : EPSCheckUnsupportedCondition(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION,PETSC_TRUE," with polynomial filter");
44 1 : PetscCall(STFilterSetInterval(eps->st,eps->inta,eps->intb));
45 1 : if (!ctx->estimatedrange) {
46 1 : PetscCall(STFilterGetRange(eps->st,&rleft,&rright));
47 1 : estimaterange = (!rleft && !rright)? PETSC_TRUE: PETSC_FALSE;
48 : }
49 1 : if (estimaterange) { /* user did not set a range */
50 1 : PetscCall(STGetMatrix(eps->st,0,&A));
51 1 : PetscCall(MatEstimateSpectralRange_EPS(A,&rleft,&rright));
52 1 : PetscCall(PetscInfo(eps,"Setting eigenvalue range to [%g,%g]\n",(double)rleft,(double)rright));
53 1 : PetscCall(STFilterSetRange(eps->st,rleft,rright));
54 1 : ctx->estimatedrange = PETSC_TRUE;
55 : }
56 1 : if (eps->ncv==PETSC_DEFAULT && eps->nev==1) eps->nev = 40; /* user did not provide nev estimation */
57 1 : PetscCall(EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd));
58 1 : 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");
59 1 : PetscFunctionReturn(PETSC_SUCCESS);
60 : }
61 :
62 18 : static PetscErrorCode EPSSetUp_Subspace(EPS eps)
63 : {
64 18 : PetscBool isfilt;
65 :
66 18 : PetscFunctionBegin;
67 18 : EPSCheckDefinite(eps);
68 18 : if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
69 18 : if (!eps->which) PetscCall(EPSSetWhichEigenpairs_Default(eps));
70 18 : if (eps->which==EPS_ALL) {
71 1 : PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
72 1 : PetscCheck(isfilt,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Spectrum slicing not supported in this solver");
73 1 : PetscCall(EPSSetUp_Subspace_Filter(eps));
74 : } else {
75 17 : PetscCheck(eps->which==EPS_LARGEST_MAGNITUDE || eps->which==EPS_TARGET_MAGNITUDE,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver supports only largest magnitude or target magnitude eigenvalues");
76 17 : PetscCall(EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd));
77 : }
78 18 : EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_EXTRACTION | EPS_FEATURE_TWOSIDED);
79 18 : PetscCheck(eps->converged==EPSConvergedRelative,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver only supports relative convergence test");
80 :
81 18 : PetscCall(EPSAllocateSolution(eps,0));
82 18 : PetscCall(EPS_SetInnerProduct(eps));
83 18 : if (eps->ishermitian) PetscCall(DSSetType(eps->ds,DSHEP));
84 6 : else PetscCall(DSSetType(eps->ds,DSNHEP));
85 18 : PetscCall(DSAllocate(eps->ds,eps->ncv));
86 18 : PetscFunctionReturn(PETSC_SUCCESS);
87 : }
88 :
89 18 : static PetscErrorCode EPSSetUpSort_Subspace(EPS eps)
90 : {
91 18 : SlepcSC sc;
92 :
93 18 : PetscFunctionBegin;
94 18 : PetscCall(EPSSetUpSort_Default(eps));
95 18 : if (eps->which==EPS_ALL) {
96 1 : PetscCall(DSGetSlepcSC(eps->ds,&sc));
97 1 : sc->rg = NULL;
98 1 : sc->comparison = SlepcCompareLargestReal;
99 1 : sc->comparisonctx = NULL;
100 1 : sc->map = NULL;
101 1 : sc->mapobj = NULL;
102 : }
103 18 : PetscFunctionReturn(PETSC_SUCCESS);
104 : }
105 :
106 : /*
107 : EPSSubspaceFindGroup - Find a group of nearly equimodular eigenvalues, provided
108 : in arrays wr and wi, according to the tolerance grptol. Also the 2-norms
109 : of the residuals must be passed in (rsd). Arrays are processed from index
110 : l to index m only. The output information is:
111 :
112 : ngrp - number of entries of the group
113 : ctr - (w(l)+w(l+ngrp-1))/2
114 : ae - average of wr(l),...,wr(l+ngrp-1)
115 : arsd - average of rsd(l),...,rsd(l+ngrp-1)
116 : */
117 544 : static PetscErrorCode EPSSubspaceFindGroup(PetscInt l,PetscInt m,PetscScalar *wr,PetscScalar *wi,PetscReal *rsd,PetscReal grptol,PetscInt *ngrp,PetscReal *ctr,PetscReal *ae,PetscReal *arsd)
118 : {
119 544 : PetscInt i;
120 544 : PetscReal rmod,rmod1;
121 :
122 544 : PetscFunctionBegin;
123 544 : *ngrp = 0;
124 544 : *ctr = 0;
125 544 : rmod = SlepcAbsEigenvalue(wr[l],wi[l]);
126 :
127 1921 : for (i=l;i<m;) {
128 1357 : rmod1 = SlepcAbsEigenvalue(wr[i],wi[i]);
129 1357 : if (PetscAbsReal(rmod-rmod1) > grptol*(rmod+rmod1)) break;
130 833 : *ctr = (rmod+rmod1)/2.0;
131 833 : if (wi[i] == 0.0) {
132 833 : (*ngrp)++;
133 833 : i++;
134 : } else {
135 0 : (*ngrp)+=2;
136 0 : i+=2;
137 : }
138 : }
139 :
140 544 : *ae = 0;
141 544 : *arsd = 0;
142 544 : if (*ngrp) {
143 1377 : for (i=l;i<l+*ngrp;i++) {
144 833 : (*ae) += PetscRealPart(wr[i]);
145 833 : (*arsd) += rsd[i]*rsd[i];
146 : }
147 544 : *ae = *ae / *ngrp;
148 544 : *arsd = PetscSqrtReal(*arsd / *ngrp);
149 : }
150 544 : PetscFunctionReturn(PETSC_SUCCESS);
151 : }
152 :
153 : /*
154 : EPSSubspaceResidualNorms - Computes the column norms of residual vectors
155 : OP*V(1:n,l:m) - V*T(1:m,l:m), where, on entry, OP*V has been computed and
156 : stored in R. On exit, rsd(l) to rsd(m) contain the computed norms.
157 : */
158 158 : static PetscErrorCode EPSSubspaceResidualNorms(BV R,BV V,Mat T,PetscInt l,PetscInt m,PetscScalar *eigi,PetscReal *rsd)
159 : {
160 158 : PetscInt i;
161 :
162 158 : PetscFunctionBegin;
163 158 : PetscCall(BVMult(R,-1.0,1.0,V,T));
164 2650 : for (i=l;i<m;i++) PetscCall(BVNormColumnBegin(R,i,NORM_2,rsd+i));
165 2650 : for (i=l;i<m;i++) PetscCall(BVNormColumnEnd(R,i,NORM_2,rsd+i));
166 : #if !defined(PETSC_USE_COMPLEX)
167 2492 : for (i=l;i<m-1;i++) {
168 2334 : if (eigi[i]!=0.0) {
169 0 : rsd[i] = SlepcAbs(rsd[i],rsd[i+1]);
170 0 : rsd[i+1] = rsd[i];
171 0 : i++;
172 : }
173 : }
174 : #endif
175 158 : PetscFunctionReturn(PETSC_SUCCESS);
176 : }
177 :
178 18 : static PetscErrorCode EPSSolve_Subspace(EPS eps)
179 : {
180 18 : Mat H,Q,S,T,B;
181 18 : BV AV,R;
182 18 : PetscBool indef;
183 18 : PetscInt i,k,ld,ngrp,nogrp,*itrsd,*itrsdold;
184 18 : PetscInt nxtsrr,idsrr,idort,nxtort,nv,ncv = eps->ncv,its,ninside;
185 18 : PetscReal arsd,oarsd,ctr,octr,ae,oae,*rsd,*orsd,tcond=1.0,gamma;
186 18 : PetscScalar *oeigr,*oeigi;
187 : /* Parameters */
188 18 : PetscInt init = 5; /* Number of initial iterations */
189 18 : PetscReal stpfac = 1.5; /* Max num of iter before next SRR step */
190 18 : PetscReal alpha = 1.0; /* Used to predict convergence of next residual */
191 18 : PetscReal beta = 1.1; /* Used to predict convergence of next residual */
192 18 : PetscReal grptol = SLEPC_DEFAULT_TOL; /* Tolerance for EPSSubspaceFindGroup */
193 18 : PetscReal cnvtol = 1e-6; /* Convergence criterion for cnv */
194 18 : PetscInt orttol = 2; /* Number of decimal digits whose loss
195 : can be tolerated in orthogonalization */
196 :
197 18 : PetscFunctionBegin;
198 18 : its = 0;
199 18 : PetscCall(PetscMalloc6(ncv,&rsd,ncv,&orsd,ncv,&oeigr,ncv,&oeigi,ncv,&itrsd,ncv,&itrsdold));
200 18 : PetscCall(DSGetLeadingDimension(eps->ds,&ld));
201 18 : PetscCall(BVDuplicate(eps->V,&AV));
202 18 : PetscCall(BVDuplicate(eps->V,&R));
203 18 : PetscCall(STGetOperator(eps->st,&S));
204 :
205 331 : for (i=0;i<ncv;i++) {
206 313 : rsd[i] = 0.0;
207 313 : itrsd[i] = -1;
208 : }
209 :
210 : /* Complete the initial basis with random vectors and orthonormalize them */
211 329 : for (k=eps->nini;k<ncv;k++) {
212 311 : PetscCall(BVSetRandomColumn(eps->V,k));
213 311 : PetscCall(BVOrthonormalizeColumn(eps->V,k,PETSC_TRUE,NULL,NULL));
214 : }
215 :
216 158 : while (eps->reason == EPS_CONVERGED_ITERATING) {
217 158 : eps->its++;
218 158 : nv = PetscMin(eps->nconv+eps->mpd,ncv);
219 158 : PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,0));
220 :
221 2650 : for (i=eps->nconv;i<nv;i++) {
222 2492 : oeigr[i] = eps->eigr[i];
223 2492 : oeigi[i] = eps->eigi[i];
224 2492 : orsd[i] = rsd[i];
225 : }
226 :
227 : /* AV(:,idx) = OP * V(:,idx) */
228 158 : PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
229 158 : PetscCall(BVSetActiveColumns(AV,eps->nconv,nv));
230 158 : PetscCall(BVMatMult(eps->V,S,AV));
231 :
232 : /* T(:,idx) = V' * AV(:,idx) */
233 158 : PetscCall(BVSetActiveColumns(eps->V,0,nv));
234 158 : PetscCall(DSGetMat(eps->ds,DS_MAT_A,&H));
235 158 : PetscCall(BVDot(AV,eps->V,H));
236 158 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&H));
237 158 : PetscCall(DSSetState(eps->ds,DS_STATE_RAW));
238 :
239 : /* Solve projected problem */
240 158 : PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
241 158 : PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
242 158 : PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
243 :
244 : /* Update vectors V(:,idx) = V * U(:,idx) */
245 158 : PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&Q));
246 158 : PetscCall(BVSetActiveColumns(AV,0,nv));
247 158 : PetscCall(BVSetActiveColumns(R,0,nv));
248 158 : PetscCall(BVMultInPlace(eps->V,Q,eps->nconv,nv));
249 158 : PetscCall(BVMultInPlace(AV,Q,eps->nconv,nv));
250 158 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&Q));
251 158 : PetscCall(BVCopy(AV,R));
252 :
253 : /* Convergence check */
254 158 : PetscCall(DSGetMat(eps->ds,DS_MAT_A,&T));
255 158 : PetscCall(EPSSubspaceResidualNorms(R,eps->V,T,eps->nconv,nv,eps->eigi,rsd));
256 158 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&T));
257 :
258 158 : if (eps->which==EPS_ALL && eps->its>1) { /* adjust eigenvalue count */
259 10 : ninside = 0;
260 10 : PetscCall(STFilterGetThreshold(eps->st,&gamma));
261 91 : for (i=eps->nconv;i<nv;i++) {
262 91 : if (PetscRealPart(eps->eigr[i]) < gamma) break;
263 81 : ninside++;
264 : }
265 10 : eps->nev = eps->nconv+ninside;
266 : }
267 2650 : for (i=eps->nconv;i<nv;i++) {
268 2492 : itrsdold[i] = itrsd[i];
269 2492 : itrsd[i] = its;
270 2492 : eps->errest[i] = rsd[i];
271 : }
272 :
273 272 : for (;;) {
274 : /* Find clusters of computed eigenvalues */
275 272 : PetscCall(EPSSubspaceFindGroup(eps->nconv,nv,eps->eigr,eps->eigi,eps->errest,grptol,&ngrp,&ctr,&ae,&arsd));
276 272 : PetscCall(EPSSubspaceFindGroup(eps->nconv,nv,oeigr,oeigi,orsd,grptol,&nogrp,&octr,&oae,&oarsd));
277 272 : if (ngrp!=nogrp) break;
278 254 : if (ngrp==0) break;
279 254 : if (PetscAbsReal(ae-oae)>ctr*cnvtol*(itrsd[eps->nconv]-itrsdold[eps->nconv])) break;
280 184 : if (arsd>ctr*eps->tol) break;
281 115 : eps->nconv = eps->nconv + ngrp;
282 115 : if (eps->nconv>=nv) break;
283 : }
284 :
285 158 : PetscCall(EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,nv));
286 158 : PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx));
287 158 : if (eps->reason != EPS_CONVERGED_ITERATING) break;
288 :
289 : /* Compute nxtsrr (iteration of next projection step) */
290 140 : nxtsrr = PetscMin(eps->max_it,PetscMax((PetscInt)PetscFloorReal(stpfac*its),init));
291 :
292 140 : if (ngrp!=nogrp || ngrp==0 || arsd>=oarsd) {
293 21 : idsrr = nxtsrr - its;
294 : } else {
295 119 : idsrr = (PetscInt)PetscFloorReal(alpha+beta*(itrsdold[eps->nconv]-itrsd[eps->nconv])*PetscLogReal(arsd/eps->tol)/PetscLogReal(arsd/oarsd));
296 119 : idsrr = PetscMax(1,idsrr);
297 : }
298 140 : nxtsrr = PetscMin(nxtsrr,its+idsrr);
299 :
300 : /* Compute nxtort (iteration of next orthogonalization step) */
301 140 : PetscCall(DSCond(eps->ds,&tcond));
302 140 : idort = PetscMax(1,(PetscInt)PetscFloorReal(orttol/PetscMax(1,PetscLog10Real(tcond))));
303 140 : nxtort = PetscMin(its+idort,nxtsrr);
304 140 : PetscCall(PetscInfo(eps,"Updated iteration counts: nxtort=%" PetscInt_FMT ", nxtsrr=%" PetscInt_FMT "\n",nxtort,nxtsrr));
305 :
306 : /* V(:,idx) = AV(:,idx) */
307 140 : PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
308 140 : PetscCall(BVSetActiveColumns(AV,eps->nconv,nv));
309 140 : PetscCall(BVCopy(AV,eps->V));
310 140 : its++;
311 :
312 : /* Orthogonalization loop */
313 1742 : do {
314 1742 : PetscCall(BVGetMatrix(eps->V,&B,&indef));
315 1742 : PetscCall(BVSetMatrix(eps->V,NULL,PETSC_FALSE));
316 4905 : while (its<nxtort) {
317 : /* A(:,idx) = OP*V(:,idx) with normalization */
318 3163 : PetscCall(BVMatMult(eps->V,S,AV));
319 3163 : PetscCall(BVCopy(AV,eps->V));
320 3163 : PetscCall(BVNormalize(eps->V,NULL));
321 3163 : its++;
322 : }
323 1742 : PetscCall(BVSetMatrix(eps->V,B,indef));
324 : /* Orthonormalize vectors */
325 1742 : PetscCall(BVOrthogonalize(eps->V,NULL));
326 1742 : nxtort = PetscMin(its+idort,nxtsrr);
327 1742 : } while (its<nxtsrr);
328 : }
329 :
330 18 : PetscCall(PetscFree6(rsd,orsd,oeigr,oeigi,itrsd,itrsdold));
331 18 : PetscCall(BVDestroy(&AV));
332 18 : PetscCall(BVDestroy(&R));
333 18 : PetscCall(STRestoreOperator(eps->st,&S));
334 18 : PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
335 18 : PetscFunctionReturn(PETSC_SUCCESS);
336 : }
337 :
338 16 : static PetscErrorCode EPSDestroy_Subspace(EPS eps)
339 : {
340 16 : PetscFunctionBegin;
341 16 : PetscCall(PetscFree(eps->data));
342 16 : PetscFunctionReturn(PETSC_SUCCESS);
343 : }
344 :
345 16 : SLEPC_EXTERN PetscErrorCode EPSCreate_Subspace(EPS eps)
346 : {
347 16 : EPS_SUBSPACE *ctx;
348 :
349 16 : PetscFunctionBegin;
350 16 : PetscCall(PetscNew(&ctx));
351 16 : eps->data = (void*)ctx;
352 :
353 16 : eps->useds = PETSC_TRUE;
354 16 : eps->categ = EPS_CATEGORY_OTHER;
355 :
356 16 : eps->ops->solve = EPSSolve_Subspace;
357 16 : eps->ops->setup = EPSSetUp_Subspace;
358 16 : eps->ops->setupsort = EPSSetUpSort_Subspace;
359 16 : eps->ops->destroy = EPSDestroy_Subspace;
360 16 : eps->ops->backtransform = EPSBackTransform_Default;
361 16 : eps->ops->computevectors = EPSComputeVectors_Schur;
362 16 : PetscFunctionReturn(PETSC_SUCCESS);
363 : }
|