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 : 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 19 : static PetscErrorCode EPSSetUp_Subspace(EPS eps)
63 : {
64 19 : PetscBool isfilt;
65 :
66 19 : PetscFunctionBegin;
67 19 : EPSCheckDefinite(eps);
68 19 : if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
69 19 : if (!eps->which) PetscCall(EPSSetWhichEigenpairs_Default(eps));
70 19 : 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 18 : 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 18 : PetscCall(EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd));
77 : }
78 19 : EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_EXTRACTION | EPS_FEATURE_TWOSIDED);
79 19 : PetscCheck(eps->converged==EPSConvergedRelative,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver only supports relative convergence test");
80 :
81 19 : PetscCall(EPSAllocateSolution(eps,0));
82 19 : PetscCall(EPS_SetInnerProduct(eps));
83 19 : if (eps->ishermitian) PetscCall(DSSetType(eps->ds,DSHEP));
84 6 : else PetscCall(DSSetType(eps->ds,DSNHEP));
85 19 : PetscCall(DSAllocate(eps->ds,eps->ncv));
86 19 : PetscFunctionReturn(PETSC_SUCCESS);
87 : }
88 :
89 19 : static PetscErrorCode EPSSetUpSort_Subspace(EPS eps)
90 : {
91 19 : SlepcSC sc;
92 :
93 19 : PetscFunctionBegin;
94 19 : PetscCall(EPSSetUpSort_Default(eps));
95 19 : 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 19 : 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 588 : 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 588 : PetscInt i;
120 588 : PetscReal rmod,rmod1;
121 :
122 588 : PetscFunctionBegin;
123 588 : *ngrp = 0;
124 588 : *ctr = 0;
125 588 : rmod = SlepcAbsEigenvalue(wr[l],wi[l]);
126 :
127 1483 : for (i=l;i<m;) {
128 1462 : rmod1 = SlepcAbsEigenvalue(wr[i],wi[i]);
129 1462 : if (PetscAbsReal(rmod-rmod1) > grptol*(rmod+rmod1)) break;
130 895 : *ctr = (rmod+rmod1)/2.0;
131 895 : if (wi[i] == 0.0) {
132 895 : (*ngrp)++;
133 895 : i++;
134 : } else {
135 0 : (*ngrp)+=2;
136 0 : i+=2;
137 : }
138 : }
139 :
140 588 : *ae = 0;
141 588 : *arsd = 0;
142 588 : if (*ngrp) {
143 1483 : for (i=l;i<l+*ngrp;i++) {
144 895 : (*ae) += PetscRealPart(wr[i]);
145 895 : (*arsd) += rsd[i]*rsd[i];
146 : }
147 588 : *ae = *ae / *ngrp;
148 588 : *arsd = PetscSqrtReal(*arsd / *ngrp);
149 : }
150 588 : 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 172 : static PetscErrorCode EPSSubspaceResidualNorms(BV R,BV V,Mat T,PetscInt l,PetscInt m,PetscScalar *eigi,PetscReal *rsd)
159 : {
160 172 : PetscInt i;
161 :
162 172 : PetscFunctionBegin;
163 172 : PetscCall(BVMult(R,-1.0,1.0,V,T));
164 2917 : for (i=l;i<m;i++) PetscCall(BVNormColumnBegin(R,i,NORM_2,rsd+i));
165 2917 : for (i=l;i<m;i++) PetscCall(BVNormColumnEnd(R,i,NORM_2,rsd+i));
166 : #if !defined(PETSC_USE_COMPLEX)
167 : for (i=l;i<m-1;i++) {
168 : if (eigi[i]!=0.0) {
169 : rsd[i] = SlepcAbs(rsd[i],rsd[i+1]);
170 : rsd[i+1] = rsd[i];
171 : i++;
172 : }
173 : }
174 : #endif
175 172 : PetscFunctionReturn(PETSC_SUCCESS);
176 : }
177 :
178 19 : static PetscErrorCode EPSSolve_Subspace(EPS eps)
179 : {
180 19 : Mat H,Q,S,T,B;
181 19 : BV AV,R;
182 19 : PetscBool indef;
183 19 : PetscInt i,k,ld,ngrp,nogrp,*itrsd,*itrsdold;
184 19 : PetscInt nxtsrr,idsrr,idort,nxtort,nv,ncv = eps->ncv,its,ninside;
185 19 : PetscReal arsd,oarsd,ctr,octr,ae,oae,*rsd,*orsd,tcond=1.0,gamma;
186 19 : PetscScalar *oeigr,*oeigi;
187 : /* Parameters */
188 19 : PetscInt init = 5; /* Number of initial iterations */
189 19 : PetscReal stpfac = 1.5; /* Max num of iter before next SRR step */
190 19 : PetscReal alpha = 1.0; /* Used to predict convergence of next residual */
191 19 : PetscReal beta = 1.1; /* Used to predict convergence of next residual */
192 19 : PetscReal grptol = SLEPC_DEFAULT_TOL; /* Tolerance for EPSSubspaceFindGroup */
193 19 : PetscReal cnvtol = 1e-6; /* Convergence criterion for cnv */
194 19 : PetscInt orttol = 2; /* Number of decimal digits whose loss
195 : can be tolerated in orthogonalization */
196 :
197 19 : PetscFunctionBegin;
198 19 : its = 0;
199 19 : PetscCall(PetscMalloc6(ncv,&rsd,ncv,&orsd,ncv,&oeigr,ncv,&oeigi,ncv,&itrsd,ncv,&itrsdold));
200 19 : PetscCall(DSGetLeadingDimension(eps->ds,&ld));
201 19 : PetscCall(BVDuplicate(eps->V,&AV));
202 19 : PetscCall(BVDuplicate(eps->V,&R));
203 19 : PetscCall(STGetOperator(eps->st,&S));
204 :
205 351 : for (i=0;i<ncv;i++) {
206 332 : rsd[i] = 0.0;
207 332 : itrsd[i] = -1;
208 : }
209 :
210 : /* Complete the initial basis with random vectors and orthonormalize them */
211 349 : for (k=eps->nini;k<ncv;k++) {
212 330 : PetscCall(BVSetRandomColumn(eps->V,k));
213 330 : PetscCall(BVOrthonormalizeColumn(eps->V,k,PETSC_TRUE,NULL,NULL));
214 : }
215 :
216 172 : while (eps->reason == EPS_CONVERGED_ITERATING) {
217 172 : eps->its++;
218 172 : nv = PetscMin(eps->nconv+eps->mpd,ncv);
219 172 : PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,0));
220 :
221 2917 : for (i=eps->nconv;i<nv;i++) {
222 2745 : oeigr[i] = eps->eigr[i];
223 2745 : oeigi[i] = eps->eigi[i];
224 2745 : orsd[i] = rsd[i];
225 : }
226 :
227 : /* AV(:,idx) = OP * V(:,idx) */
228 172 : PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
229 172 : PetscCall(BVSetActiveColumns(AV,eps->nconv,nv));
230 172 : PetscCall(BVMatMult(eps->V,S,AV));
231 :
232 : /* T(:,idx) = V' * AV(:,idx) */
233 172 : PetscCall(BVSetActiveColumns(eps->V,0,nv));
234 172 : PetscCall(DSGetMat(eps->ds,DS_MAT_A,&H));
235 172 : PetscCall(BVDot(AV,eps->V,H));
236 172 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&H));
237 172 : PetscCall(DSSetState(eps->ds,DS_STATE_RAW));
238 :
239 : /* Solve projected problem */
240 172 : PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
241 172 : PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
242 172 : PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
243 :
244 : /* Update vectors V(:,idx) = V * U(:,idx) */
245 172 : PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&Q));
246 172 : PetscCall(BVSetActiveColumns(AV,0,nv));
247 172 : PetscCall(BVSetActiveColumns(R,0,nv));
248 172 : PetscCall(BVMultInPlace(eps->V,Q,eps->nconv,nv));
249 172 : PetscCall(BVMultInPlace(AV,Q,eps->nconv,nv));
250 172 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&Q));
251 172 : PetscCall(BVCopy(AV,R));
252 :
253 : /* Convergence check */
254 172 : PetscCall(DSGetMat(eps->ds,DS_MAT_A,&T));
255 172 : PetscCall(EPSSubspaceResidualNorms(R,eps->V,T,eps->nconv,nv,eps->eigi,rsd));
256 172 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&T));
257 :
258 172 : if (eps->which==EPS_ALL && eps->its>1) { /* adjust eigenvalue count */
259 9 : ninside = 0;
260 9 : PetscCall(STFilterGetThreshold(eps->st,&gamma));
261 83 : for (i=eps->nconv;i<nv;i++) {
262 83 : if (PetscRealPart(eps->eigr[i]) < gamma) break;
263 74 : ninside++;
264 : }
265 9 : eps->nev = eps->nconv+ninside;
266 : }
267 2917 : for (i=eps->nconv;i<nv;i++) {
268 2745 : itrsdold[i] = itrsd[i];
269 2745 : itrsd[i] = its;
270 2745 : eps->errest[i] = rsd[i];
271 : }
272 :
273 294 : for (;;) {
274 : /* Find clusters of computed eigenvalues */
275 294 : PetscCall(EPSSubspaceFindGroup(eps->nconv,nv,eps->eigr,eps->eigi,eps->errest,grptol,&ngrp,&ctr,&ae,&arsd));
276 294 : PetscCall(EPSSubspaceFindGroup(eps->nconv,nv,oeigr,oeigi,orsd,grptol,&nogrp,&octr,&oae,&oarsd));
277 294 : if (ngrp!=nogrp) break;
278 275 : if (ngrp==0) break;
279 275 : if (PetscAbsReal(ae-oae)>ctr*cnvtol*(itrsd[eps->nconv]-itrsdold[eps->nconv])) break;
280 199 : if (arsd>ctr*eps->tol) break;
281 123 : eps->nconv = eps->nconv + ngrp;
282 123 : if (eps->nconv>=nv) break;
283 : }
284 :
285 172 : PetscCall(EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,nv));
286 172 : PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx));
287 172 : if (eps->reason != EPS_CONVERGED_ITERATING) break;
288 :
289 : /* Compute nxtsrr (iteration of next projection step) */
290 153 : nxtsrr = PetscMin(eps->max_it,PetscMax((PetscInt)PetscFloorReal(stpfac*its),init));
291 :
292 153 : if (ngrp!=nogrp || ngrp==0 || arsd>=oarsd) {
293 22 : idsrr = nxtsrr - its;
294 : } else {
295 131 : idsrr = (PetscInt)PetscFloorReal(alpha+beta*(itrsdold[eps->nconv]-itrsd[eps->nconv])*PetscLogReal(arsd/eps->tol)/PetscLogReal(arsd/oarsd));
296 131 : idsrr = PetscMax(1,idsrr);
297 : }
298 153 : nxtsrr = PetscMin(nxtsrr,its+idsrr);
299 :
300 : /* Compute nxtort (iteration of next orthogonalization step) */
301 153 : PetscCall(DSCond(eps->ds,&tcond));
302 153 : idort = PetscMax(1,(PetscInt)PetscFloorReal(orttol/PetscMax(1,PetscLog10Real(tcond))));
303 153 : nxtort = PetscMin(its+idort,nxtsrr);
304 153 : PetscCall(PetscInfo(eps,"Updated iteration counts: nxtort=%" PetscInt_FMT ", nxtsrr=%" PetscInt_FMT "\n",nxtort,nxtsrr));
305 :
306 : /* V(:,idx) = AV(:,idx) */
307 153 : PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
308 153 : PetscCall(BVSetActiveColumns(AV,eps->nconv,nv));
309 153 : PetscCall(BVCopy(AV,eps->V));
310 153 : its++;
311 :
312 : /* Orthogonalization loop */
313 1939 : do {
314 1939 : PetscCall(BVGetMatrix(eps->V,&B,&indef));
315 1939 : PetscCall(BVSetMatrix(eps->V,NULL,PETSC_FALSE));
316 5489 : while (its<nxtort) {
317 : /* A(:,idx) = OP*V(:,idx) with normalization */
318 3550 : PetscCall(BVMatMult(eps->V,S,AV));
319 3550 : PetscCall(BVCopy(AV,eps->V));
320 3550 : PetscCall(BVNormalize(eps->V,NULL));
321 3550 : its++;
322 : }
323 1939 : PetscCall(BVSetMatrix(eps->V,B,indef));
324 : /* Orthonormalize vectors */
325 1939 : PetscCall(BVOrthogonalize(eps->V,NULL));
326 1939 : nxtort = PetscMin(its+idort,nxtsrr);
327 1939 : } while (its<nxtsrr);
328 : }
329 :
330 19 : PetscCall(PetscFree6(rsd,orsd,oeigr,oeigi,itrsd,itrsdold));
331 19 : PetscCall(BVDestroy(&AV));
332 19 : PetscCall(BVDestroy(&R));
333 19 : PetscCall(STRestoreOperator(eps->st,&S));
334 19 : PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
335 19 : PetscFunctionReturn(PETSC_SUCCESS);
336 : }
337 :
338 17 : static PetscErrorCode EPSDestroy_Subspace(EPS eps)
339 : {
340 17 : PetscFunctionBegin;
341 17 : PetscCall(PetscFree(eps->data));
342 17 : PetscFunctionReturn(PETSC_SUCCESS);
343 : }
344 :
345 17 : SLEPC_EXTERN PetscErrorCode EPSCreate_Subspace(EPS eps)
346 : {
347 17 : EPS_SUBSPACE *ctx;
348 :
349 17 : PetscFunctionBegin;
350 17 : PetscCall(PetscNew(&ctx));
351 17 : eps->data = (void*)ctx;
352 :
353 17 : eps->useds = PETSC_TRUE;
354 17 : eps->categ = EPS_CATEGORY_OTHER;
355 :
356 17 : eps->ops->solve = EPSSolve_Subspace;
357 17 : eps->ops->setup = EPSSetUp_Subspace;
358 17 : eps->ops->setupsort = EPSSetUpSort_Subspace;
359 17 : eps->ops->destroy = EPSDestroy_Subspace;
360 17 : eps->ops->backtransform = EPSBackTransform_Default;
361 17 : eps->ops->computevectors = EPSComputeVectors_Schur;
362 17 : PetscFunctionReturn(PETSC_SUCCESS);
363 : }
|