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