Actual source code: subspace.c
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: MatState fstate; /* state of matrix used to compute the filter */
31: } EPS_SUBSPACE;
33: static PetscErrorCode EPSSetUp_Subspace_Filter(EPS eps)
34: {
35: EPS_SUBSPACE *ctx = (EPS_SUBSPACE*)eps->data;
36: PetscBool estimaterange=PETSC_TRUE,same;
37: PetscReal rleft,rright;
38: Mat A;
40: PetscFunctionBegin;
41: EPSCheckHermitianCondition(eps,PETSC_TRUE," with polynomial filter");
42: EPSCheckStandardCondition(eps,PETSC_TRUE," with polynomial filter");
43: 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");
44: EPSCheckUnsupportedCondition(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION,PETSC_TRUE," with polynomial filter");
45: PetscCall(STFilterSetInterval(eps->st,eps->inta,eps->intb));
46: if (!ctx->estimatedrange) {
47: PetscCall(STFilterGetRange(eps->st,&rleft,&rright));
48: estimaterange = (!rleft && !rright)? PETSC_TRUE: PETSC_FALSE;
49: }
50: PetscCall(STGetMatrix(eps->st,0,&A));
51: PetscCall(MatStateCompareUpdate(A,&ctx->fstate,&same));
52: if (estimaterange && !same) { /* user did not set a range */
53: PetscCall(MatEstimateSpectralRange_EPS(A,&rleft,&rright));
54: PetscCall(PetscInfo(eps,"Setting eigenvalue range to [%g,%g]\n",(double)rleft,(double)rright));
55: PetscCall(STFilterSetRange(eps->st,rleft,rright));
56: ctx->estimatedrange = PETSC_TRUE;
57: }
58: PetscCall(EPSSetStoppingTest(eps,EPS_STOP_THRESHOLD));
59: PetscCall(EPSSetDimensions_Default(eps,&eps->nev,&eps->ncv,&eps->mpd));
60: 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");
61: PetscFunctionReturn(PETSC_SUCCESS);
62: }
64: static PetscErrorCode EPSSetUp_Subspace(EPS eps)
65: {
66: PetscBool isfilt;
68: PetscFunctionBegin;
69: EPSCheckDefinite(eps);
70: EPSCheckNotStructured(eps);
71: if (eps->max_it==PETSC_DETERMINE) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
72: if (!eps->which) PetscCall(EPSSetWhichEigenpairs_Default(eps));
73: if (eps->which==EPS_ALL) {
74: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
75: PetscCheck(isfilt,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Spectrum slicing not supported in this solver");
76: PetscCall(EPSSetUp_Subspace_Filter(eps));
77: } else {
78: 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");
79: PetscCall(EPSSetDimensions_Default(eps,&eps->nev,&eps->ncv,&eps->mpd));
80: }
81: EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_EXTRACTION | EPS_FEATURE_TWOSIDED);
82: PetscCheck(eps->converged==EPSConvergedRelative,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver only supports relative convergence test");
84: PetscCall(EPSAllocateSolution(eps,0));
85: PetscCall(EPS_SetInnerProduct(eps));
86: if (eps->ishermitian) PetscCall(DSSetType(eps->ds,DSHEP));
87: else PetscCall(DSSetType(eps->ds,DSNHEP));
88: PetscCall(DSAllocate(eps->ds,eps->ncv));
89: PetscFunctionReturn(PETSC_SUCCESS);
90: }
92: static PetscErrorCode EPSSetUpSort_Subspace(EPS eps)
93: {
94: SlepcSC sc;
96: PetscFunctionBegin;
97: PetscCall(EPSSetUpSort_Default(eps));
98: if (eps->which==EPS_ALL) {
99: PetscCall(DSGetSlepcSC(eps->ds,&sc));
100: sc->rg = NULL;
101: sc->comparison = SlepcCompareLargestReal;
102: sc->comparisonctx = NULL;
103: sc->map = NULL;
104: sc->mapobj = NULL;
105: }
106: PetscFunctionReturn(PETSC_SUCCESS);
107: }
109: /*
110: EPSSubspaceFindGroup - Find a group of nearly equimodular eigenvalues, provided
111: in arrays wr and wi, according to the tolerance grptol. Also the 2-norms
112: of the residuals must be passed in (rsd). Arrays are processed from index
113: l to index m only. The output information is:
115: ngrp - number of entries of the group
116: ctr - (w(l)+w(l+ngrp-1))/2
117: ae - average of wr(l),...,wr(l+ngrp-1)
118: arsd - average of rsd(l),...,rsd(l+ngrp-1)
119: */
120: static PetscErrorCode EPSSubspaceFindGroup(PetscInt l,PetscInt m,PetscScalar *wr,PetscScalar *wi,PetscReal *rsd,PetscReal grptol,PetscInt *ngrp,PetscReal *ctr,PetscReal *ae,PetscReal *arsd)
121: {
122: PetscInt i;
123: PetscReal rmod,rmod1;
125: PetscFunctionBegin;
126: *ngrp = 0;
127: *ctr = 0;
128: rmod = SlepcAbsEigenvalue(wr[l],wi[l]);
130: for (i=l;i<m;) {
131: rmod1 = SlepcAbsEigenvalue(wr[i],wi[i]);
132: if (PetscAbsReal(rmod-rmod1) > grptol*(rmod+rmod1)) break;
133: *ctr = (rmod+rmod1)/2.0;
134: if (wi[i] == 0.0) {
135: (*ngrp)++;
136: i++;
137: } else {
138: (*ngrp)+=2;
139: i+=2;
140: }
141: }
143: *ae = 0;
144: *arsd = 0;
145: if (*ngrp) {
146: for (i=l;i<l+*ngrp;i++) {
147: (*ae) += PetscRealPart(wr[i]);
148: (*arsd) += rsd[i]*rsd[i];
149: }
150: *ae = *ae / *ngrp;
151: *arsd = PetscSqrtReal(*arsd / *ngrp);
152: }
153: PetscFunctionReturn(PETSC_SUCCESS);
154: }
156: /*
157: EPSSubspaceResidualNorms - Computes the column norms of residual vectors
158: OP*V(1:n,l:m) - V*T(1:m,l:m), where, on entry, OP*V has been computed and
159: stored in R. On exit, rsd(l) to rsd(m) contain the computed norms.
160: */
161: static PetscErrorCode EPSSubspaceResidualNorms(BV R,BV V,Mat T,PetscInt l,PetscInt m,PetscScalar *eigi,PetscReal *rsd)
162: {
163: PetscInt i;
165: PetscFunctionBegin;
166: PetscCall(BVMult(R,-1.0,1.0,V,T));
167: for (i=l;i<m;i++) PetscCall(BVNormColumnBegin(R,i,NORM_2,rsd+i));
168: for (i=l;i<m;i++) PetscCall(BVNormColumnEnd(R,i,NORM_2,rsd+i));
169: #if !PetscDefined(USE_COMPLEX)
170: for (i=l;i<m-1;i++) {
171: if (eigi[i]!=0.0) {
172: rsd[i] = SlepcAbs(rsd[i],rsd[i+1]);
173: rsd[i+1] = rsd[i];
174: i++;
175: }
176: }
177: #endif
178: PetscFunctionReturn(PETSC_SUCCESS);
179: }
181: static PetscErrorCode EPSSolve_Subspace(EPS eps)
182: {
183: Mat H,Q,S,T,B;
184: BV AV,R;
185: PetscBool indef;
186: PetscInt i,k,ngrp,nogrp,*itrsd,*itrsdold,*p_itrsd,*p_itrsdold;
187: PetscInt nxtsrr,idsrr,idort,nxtort,nv,its;
188: PetscReal arsd,oarsd,ctr,octr,ae,oae,*rsd,*orsd,*p_rsd,*p_orsd,tcond=1.0;
189: PetscScalar *oeigr,*oeigi,*p_oeigr,*p_oeigi;
190: /* Parameters */
191: PetscInt init = 5; /* Number of initial iterations */
192: PetscReal stpfac = 1.5; /* Max num of iter before next SRR step */
193: PetscReal alpha = 1.0; /* Used to predict convergence of next residual */
194: PetscReal beta = 1.1; /* Used to predict convergence of next residual */
195: PetscReal grptol = SLEPC_DEFAULT_TOL; /* Tolerance for EPSSubspaceFindGroup */
196: PetscReal cnvtol = 1e-6; /* Convergence criterion for cnv */
197: PetscInt orttol = 2; /* Number of decimal digits whose loss
198: can be tolerated in orthogonalization */
200: PetscFunctionBegin;
201: its = 0;
202: PetscCall(PetscMalloc6(eps->ncv,&rsd,eps->ncv,&orsd,eps->ncv,&oeigr,eps->ncv,&oeigi,eps->ncv,&itrsd,eps->ncv,&itrsdold));
203: PetscCall(BVDuplicate(eps->V,&AV));
204: PetscCall(BVDuplicate(eps->V,&R));
205: PetscCall(STGetOperator(eps->st,&S));
207: for (i=0;i<eps->ncv;i++) {
208: rsd[i] = 0.0;
209: itrsd[i] = -1;
210: }
212: /* Complete the initial basis with random vectors and orthonormalize them */
213: for (k=eps->nini;k<eps->ncv;k++) {
214: PetscCall(BVSetRandomColumn(eps->V,k));
215: PetscCall(BVOrthonormalizeColumn(eps->V,k,PETSC_TRUE,NULL,NULL));
216: }
218: while (eps->reason == EPS_CONVERGED_ITERATING) {
219: eps->its++;
220: nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
221: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,0));
223: for (i=eps->nconv;i<nv;i++) {
224: oeigr[i] = eps->eigr[i];
225: oeigi[i] = eps->eigi[i];
226: orsd[i] = rsd[i];
227: }
229: /* AV(:,idx) = OP * V(:,idx) */
230: PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
231: PetscCall(BVSetActiveColumns(AV,eps->nconv,nv));
232: PetscCall(BVMatMult(eps->V,S,AV));
234: /* T(:,idx) = V' * AV(:,idx) */
235: PetscCall(BVSetActiveColumns(eps->V,0,nv));
236: PetscCall(DSGetMat(eps->ds,DS_MAT_A,&H));
237: PetscCall(BVDot(AV,eps->V,H));
238: PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&H));
239: PetscCall(DSSetState(eps->ds,DS_STATE_RAW));
241: /* Solve projected problem */
242: PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
243: PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
244: PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
246: /* Update vectors V(:,idx) = V * U(:,idx) */
247: PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&Q));
248: PetscCall(BVSetActiveColumns(AV,0,nv));
249: PetscCall(BVSetActiveColumns(R,0,nv));
250: PetscCall(BVMultInPlace(eps->V,Q,eps->nconv,nv));
251: PetscCall(BVMultInPlace(AV,Q,eps->nconv,nv));
252: PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&Q));
253: PetscCall(BVCopy(AV,R));
255: /* Convergence check */
256: PetscCall(DSGetMat(eps->ds,DS_MAT_A,&T));
257: PetscCall(EPSSubspaceResidualNorms(R,eps->V,T,eps->nconv,nv,eps->eigi,rsd));
258: PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&T));
260: for (i=eps->nconv;i<nv;i++) {
261: itrsdold[i] = itrsd[i];
262: itrsd[i] = its;
263: eps->errest[i] = rsd[i];
264: }
266: for (;;) {
267: /* Find clusters of computed eigenvalues */
268: PetscCall(EPSSubspaceFindGroup(eps->nconv,nv,eps->eigr,eps->eigi,eps->errest,grptol,&ngrp,&ctr,&ae,&arsd));
269: PetscCall(EPSSubspaceFindGroup(eps->nconv,nv,oeigr,oeigi,orsd,grptol,&nogrp,&octr,&oae,&oarsd));
270: if (ngrp!=nogrp) break;
271: if (ngrp==0) break;
272: if (PetscAbsReal(ae-oae)>ctr*cnvtol*(itrsd[eps->nconv]-itrsdold[eps->nconv])) break;
273: if (arsd>ctr*eps->tol) break;
274: eps->nconv = eps->nconv + ngrp;
275: if (eps->nconv>=nv) break;
276: }
278: PetscCall(EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,nv));
279: EPSSetCtxThreshold(eps,eps->eigr,eps->eigi,eps->errest,eps->nconv,nv);
280: PetscCall((*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 = 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: PetscCall(DSCond(eps->ds,&tcond));
296: idort = PetscMax(1,(PetscInt)PetscFloorReal(orttol/PetscMax(1,PetscLog10Real(tcond))));
297: nxtort = PetscMin(its+idort,nxtsrr);
298: PetscCall(PetscInfo(eps,"Updated iteration counts: nxtort=%" PetscInt_FMT ", nxtsrr=%" PetscInt_FMT "\n",nxtort,nxtsrr));
300: /* V(:,idx) = AV(:,idx) */
301: PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
302: PetscCall(BVSetActiveColumns(AV,eps->nconv,nv));
303: PetscCall(BVCopy(AV,eps->V));
304: its++;
306: /* Orthogonalization loop */
307: do {
308: PetscCall(BVGetMatrix(eps->V,&B,&indef));
309: PetscCall(BVSetMatrix(eps->V,NULL,PETSC_FALSE));
310: while (its<nxtort) {
311: /* A(:,idx) = OP*V(:,idx) with normalization */
312: PetscCall(BVMatMult(eps->V,S,AV));
313: PetscCall(BVCopy(AV,eps->V));
314: PetscCall(BVNormalize(eps->V,NULL));
315: its++;
316: }
317: PetscCall(BVSetMatrix(eps->V,B,indef));
318: /* Orthonormalize vectors */
319: PetscCall(BVOrthogonalize(eps->V,NULL));
320: nxtort = PetscMin(its+idort,nxtsrr);
321: } while (its<nxtsrr);
323: if (eps->stop==EPS_STOP_THRESHOLD && nv/(PetscReal)eps->nconv<1.3) { /* reallocate */
324: eps->ncv = 1.3*eps->nconv+1;
325: PetscCall(EPSReallocateSolution(eps,eps->ncv));
326: PetscCall(BVResize(AV,eps->ncv,PETSC_TRUE));
327: PetscCall(BVResize(R,eps->ncv,PETSC_TRUE));
328: for (i=nv;i<eps->ncv;i++) eps->perm[i] = i;
329: PetscCall(DSReallocate(eps->ds,eps->ncv));
330: for (k=nv;k<eps->ncv;k++) {
331: PetscCall(BVSetRandomColumn(eps->V,k));
332: PetscCall(BVOrthonormalizeColumn(eps->V,k,PETSC_TRUE,NULL,NULL));
333: }
334: p_rsd = rsd;
335: p_orsd = orsd;
336: p_oeigr = oeigr;
337: p_oeigi = oeigi;
338: p_itrsd = itrsd;
339: p_itrsdold = itrsdold;
340: PetscCall(PetscMalloc6(eps->ncv,&rsd,eps->ncv,&orsd,eps->ncv,&oeigr,eps->ncv,&oeigi,eps->ncv,&itrsd,eps->ncv,&itrsdold));
341: for (i=0;i<nv;i++) {
342: rsd[i] = p_rsd[i];
343: orsd[i] = p_orsd[i];
344: oeigr[i] = p_oeigr[i];
345: oeigi[i] = p_oeigi[i];
346: itrsd[i] = p_itrsd[i];
347: itrsdold[i] = p_itrsdold[i];
348: }
349: PetscCall(PetscFree6(p_rsd,p_orsd,p_oeigr,p_oeigi,p_itrsd,p_itrsdold));
350: }
352: }
354: PetscCall(PetscFree6(rsd,orsd,oeigr,oeigi,itrsd,itrsdold));
355: PetscCall(BVDestroy(&AV));
356: PetscCall(BVDestroy(&R));
357: PetscCall(STRestoreOperator(eps->st,&S));
358: PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
359: PetscFunctionReturn(PETSC_SUCCESS);
360: }
362: static PetscErrorCode EPSDestroy_Subspace(EPS eps)
363: {
364: PetscFunctionBegin;
365: PetscCall(PetscFree(eps->data));
366: PetscFunctionReturn(PETSC_SUCCESS);
367: }
369: /*MC
370: EPSSUBSPACE - EPSSUBSPACE = "subspace" - The simple subspace iteration method.
372: Notes:
373: This solver is very basic and is not recommended in general, since it
374: will not be competitive with respect to other solvers.
376: The implemented method is subspace iteration with Rayleigh-Ritz projection and
377: locking, based on the SRRIT code {cite:p}`Bai97`.
379: Level: beginner
381: .seealso: [](ch:eps), `EPS`, `EPSType`, `EPSSetType()`
382: M*/
383: SLEPC_EXTERN PetscErrorCode EPSCreate_Subspace(EPS eps)
384: {
385: EPS_SUBSPACE *ctx;
387: PetscFunctionBegin;
388: PetscCall(PetscNew(&ctx));
389: eps->data = (void*)ctx;
391: eps->useds = PETSC_TRUE;
392: eps->categ = EPS_CATEGORY_OTHER;
394: eps->ops->solve = EPSSolve_Subspace;
395: eps->ops->setup = EPSSetUp_Subspace;
396: eps->ops->setupsort = EPSSetUpSort_Subspace;
397: eps->ops->destroy = EPSDestroy_Subspace;
398: eps->ops->backtransform = EPSBackTransform_Default;
399: eps->ops->computevectors = EPSComputeVectors_Schur;
400: PetscFunctionReturn(PETSC_SUCCESS);
401: }