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