Actual source code: subspace.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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 http://slepc.upv.es.
24: */
26: #include <slepc/private/epsimpl.h>
28: PetscErrorCode EPSSetUp_Subspace(EPS eps)
29: {
33: EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd);
34: if (!eps->max_it) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
35: if (!eps->which) { EPSSetWhichEigenpairs_Default(eps); }
36: if (eps->which!=EPS_LARGEST_MAGNITUDE && eps->which!=EPS_TARGET_MAGNITUDE) SETERRQ(PetscObjectComm((PetscObject)eps),1,"Wrong value of eps->which");
37: if (!eps->extraction) {
38: EPSSetExtraction(eps,EPS_RITZ);
39: } else if (eps->extraction!=EPS_RITZ) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
40: if (eps->arbitrary) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Arbitrary selection of eigenpairs not supported in this solver");
42: EPSAllocateSolution(eps,0);
43: EPS_SetInnerProduct(eps);
44: if (eps->ishermitian) {
45: DSSetType(eps->ds,DSHEP);
46: } else {
47: DSSetType(eps->ds,DSNHEP);
48: }
49: DSAllocate(eps->ds,eps->ncv);
50: EPSSetWorkVecs(eps,1);
52: if (eps->isgeneralized && eps->ishermitian && !eps->ispositive) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Requested method does not work for indefinite problems");
53: return(0);
54: }
56: /*
57: EPSSubspaceFindGroup - Find a group of nearly equimodular eigenvalues, provided
58: in arrays wr and wi, according to the tolerance grptol. Also the 2-norms
59: of the residuals must be passed in (rsd). Arrays are processed from index
60: l to index m only. The output information is:
62: ngrp - number of entries of the group
63: ctr - (w(l)+w(l+ngrp-1))/2
64: ae - average of wr(l),...,wr(l+ngrp-1)
65: arsd - average of rsd(l),...,rsd(l+ngrp-1)
66: */
67: static PetscErrorCode EPSSubspaceFindGroup(PetscInt l,PetscInt m,PetscScalar *wr,PetscScalar *wi,PetscReal *rsd,PetscReal grptol,PetscInt *ngrp,PetscReal *ctr,PetscReal *ae,PetscReal *arsd)
68: {
69: PetscInt i;
70: PetscReal rmod,rmod1;
73: *ngrp = 0;
74: *ctr = 0;
75: rmod = SlepcAbsEigenvalue(wr[l],wi[l]);
77: for (i=l;i<m;) {
78: rmod1 = SlepcAbsEigenvalue(wr[i],wi[i]);
79: if (PetscAbsReal(rmod-rmod1) > grptol*(rmod+rmod1)) break;
80: *ctr = (rmod+rmod1)/2.0;
81: if (wi[i] != 0.0) {
82: (*ngrp)+=2;
83: i+=2;
84: } else {
85: (*ngrp)++;
86: i++;
87: }
88: }
90: *ae = 0;
91: *arsd = 0;
92: if (*ngrp) {
93: for (i=l;i<l+*ngrp;i++) {
94: (*ae) += PetscRealPart(wr[i]);
95: (*arsd) += rsd[i]*rsd[i];
96: }
97: *ae = *ae / *ngrp;
98: *arsd = PetscSqrtScalar(*arsd / *ngrp);
99: }
100: return(0);
101: }
103: /*
104: EPSSubspaceResidualNorms - Computes the column norms of residual vectors
105: OP*V(1:n,l:m) - V*T(1:m,l:m), where, on entry, OP*V has been computed and
106: stored in AV. ldt is the leading dimension of T. On exit, rsd(l) to
107: rsd(m) contain the computed norms.
108: */
109: static PetscErrorCode EPSSubspaceResidualNorms(BV V,BV AV,PetscScalar *T,PetscInt l,PetscInt m,PetscInt ldt,Vec w,PetscReal *rsd)
110: {
112: PetscInt i,k;
113: PetscScalar t;
116: for (i=l;i<m;i++) {
117: if (i==m-1 || T[i+1+ldt*i]==0.0) k=i+1;
118: else k=i+2;
119: BVSetActiveColumns(V,0,k);
120: BVCopyVec(AV,i,w);
121: BVMultVec(V,-1.0,1.0,w,T+ldt*i);
122: VecDot(w,w,&t);
123: rsd[i] = PetscRealPart(t);
124: }
125: for (i=l;i<m;i++) {
126: if (i == m-1) {
127: rsd[i] = PetscSqrtReal(rsd[i]);
128: } else if (T[i+1+(ldt*i)]==0.0) {
129: rsd[i] = PetscSqrtReal(rsd[i]);
130: } else {
131: rsd[i] = PetscSqrtReal((rsd[i]+rsd[i+1])/2.0);
132: rsd[i+1] = rsd[i];
133: i++;
134: }
135: }
136: return(0);
137: }
139: PetscErrorCode EPSSolve_Subspace(EPS eps)
140: {
142: Vec w=eps->work[0];
143: Mat H,Q,S;
144: BV AV;
145: PetscInt i,k,ld,ngrp,nogrp,*itrsd,*itrsdold;
146: PetscInt nxtsrr,idsrr,idort,nxtort,nv,ncv = eps->ncv,its;
147: PetscScalar *T;
148: PetscReal arsd,oarsd,ctr,octr,ae,oae,*rsd,norm,tcond=1.0;
149: /* Parameters */
150: PetscInt init = 5; /* Number of initial iterations */
151: PetscReal stpfac = 1.5; /* Max num of iter before next SRR step */
152: PetscReal alpha = 1.0; /* Used to predict convergence of next residual */
153: PetscReal beta = 1.1; /* Used to predict convergence of next residual */
154: PetscReal grptol = 1e-8; /* Tolerance for EPSSubspaceFindGroup */
155: PetscReal cnvtol = 1e-6; /* Convergence criterion for cnv */
156: PetscInt orttol = 2; /* Number of decimal digits whose loss
157: can be tolerated in orthogonalization */
160: its = 0;
161: PetscMalloc3(ncv,&rsd,ncv,&itrsd,ncv,&itrsdold);
162: DSGetLeadingDimension(eps->ds,&ld);
163: BVDuplicate(eps->V,&AV);
164: STGetOperator(eps->st,&S);
166: for (i=0;i<ncv;i++) {
167: rsd[i] = 0.0;
168: itrsd[i] = -1;
169: }
171: /* Complete the initial basis with random vectors and orthonormalize them */
172: for (k=eps->nini;k<ncv;k++) {
173: BVSetRandomColumn(eps->V,k);
174: BVOrthonormalizeColumn(eps->V,k,PETSC_TRUE,NULL,NULL);
175: }
177: while (eps->reason == EPS_CONVERGED_ITERATING) {
178: eps->its++;
179: nv = PetscMin(eps->nconv+eps->mpd,ncv);
180: DSSetDimensions(eps->ds,nv,0,eps->nconv,0);
182: /* Find group in previously computed eigenvalues */
183: EPSSubspaceFindGroup(eps->nconv,nv,eps->eigr,eps->eigi,rsd,grptol,&nogrp,&octr,&oae,&oarsd);
185: /* AV(:,idx) = OP * V(:,idx) */
186: BVSetActiveColumns(eps->V,eps->nconv,nv);
187: BVSetActiveColumns(AV,eps->nconv,nv);
188: BVMatMult(eps->V,S,AV);
190: /* T(:,idx) = V' * AV(:,idx) */
191: BVSetActiveColumns(eps->V,0,nv);
192: DSGetMat(eps->ds,DS_MAT_A,&H);
193: BVDot(AV,eps->V,H);
194: DSRestoreMat(eps->ds,DS_MAT_A,&H);
195: DSSetState(eps->ds,DS_STATE_RAW);
197: /* Solve projected problem */
198: DSSolve(eps->ds,eps->eigr,eps->eigi);
199: DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL);
200: DSSynchronize(eps->ds,eps->eigr,eps->eigi);
202: /* Update vectors V(:,idx) = V * U(:,idx) */
203: DSGetMat(eps->ds,DS_MAT_Q,&Q);
204: BVSetActiveColumns(AV,0,nv);
205: BVMultInPlace(eps->V,Q,eps->nconv,nv);
206: BVMultInPlace(AV,Q,eps->nconv,nv);
207: MatDestroy(&Q);
209: /* Convergence check */
210: DSGetArray(eps->ds,DS_MAT_A,&T);
211: EPSSubspaceResidualNorms(eps->V,AV,T,eps->nconv,nv,ld,w,rsd);
212: DSRestoreArray(eps->ds,DS_MAT_A,&T);
214: for (i=eps->nconv;i<nv;i++) {
215: itrsdold[i] = itrsd[i];
216: itrsd[i] = its;
217: eps->errest[i] = rsd[i];
218: }
220: for (;;) {
221: /* Find group in currently computed eigenvalues */
222: EPSSubspaceFindGroup(eps->nconv,nv,eps->eigr,eps->eigi,eps->errest,grptol,&ngrp,&ctr,&ae,&arsd);
223: if (ngrp!=nogrp) break;
224: if (ngrp==0) break;
225: if (PetscAbsReal(ae-oae)>ctr*cnvtol*(itrsd[eps->nconv]-itrsdold[eps->nconv])) break;
226: if (arsd>ctr*eps->tol) break;
227: eps->nconv = eps->nconv + ngrp;
228: if (eps->nconv>=nv) break;
229: }
231: EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,nv);
232: (*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx);
233: if (eps->reason != EPS_CONVERGED_ITERATING) break;
235: /* Compute nxtsrr (iteration of next projection step) */
236: nxtsrr = PetscMin(eps->max_it,PetscMax((PetscInt)PetscFloorReal(stpfac*its),init));
238: if (ngrp!=nogrp || ngrp==0 || arsd>=oarsd) {
239: idsrr = nxtsrr - its;
240: } else {
241: idsrr = (PetscInt)PetscFloorReal(alpha+beta*(itrsdold[eps->nconv]-itrsd[eps->nconv])*PetscLogReal(arsd/eps->tol)/PetscLogReal(arsd/oarsd));
242: idsrr = PetscMax(1,idsrr);
243: }
244: nxtsrr = PetscMin(nxtsrr,its+idsrr);
246: /* Compute nxtort (iteration of next orthogonalization step) */
247: DSCond(eps->ds,&tcond);
248: idort = PetscMax(1,(PetscInt)PetscFloorReal(orttol/PetscMax(1,PetscLog10Real(tcond))));
249: nxtort = PetscMin(its+idort,nxtsrr);
251: /* V(:,idx) = AV(:,idx) */
252: BVSetActiveColumns(eps->V,eps->nconv,nv);
253: BVSetActiveColumns(AV,eps->nconv,nv);
254: BVCopy(AV,eps->V);
255: its++;
257: /* Orthogonalization loop */
258: do {
259: while (its<nxtort) {
260: /* A(:,idx) = OP*V(:,idx) with normalization */
261: BVMatMult(eps->V,S,AV);
262: BVCopy(AV,eps->V);
263: for (i=eps->nconv;i<nv;i++) {
264: BVNormColumn(eps->V,i,NORM_INFINITY,&norm);
265: BVScaleColumn(eps->V,i,1/norm);
266: }
267: its++;
268: }
269: /* Orthonormalize vectors */
270: BVOrthogonalize(eps->V,NULL);
271: nxtort = PetscMin(its+idort,nxtsrr);
272: } while (its<nxtsrr);
273: }
275: PetscFree3(rsd,itrsd,itrsdold);
276: BVDestroy(&AV);
277: MatDestroy(&S);
278: /* truncate Schur decomposition and change the state to raw so that
279: DSVectors() computes eigenvectors from scratch */
280: DSSetDimensions(eps->ds,eps->nconv,0,0,0);
281: DSSetState(eps->ds,DS_STATE_RAW);
282: return(0);
283: }
285: PetscErrorCode EPSDestroy_Subspace(EPS eps)
286: {
290: PetscFree(eps->data);
291: return(0);
292: }
294: SLEPC_EXTERN PetscErrorCode EPSCreate_Subspace(EPS eps)
295: {
297: eps->useds = PETSC_TRUE;
298: eps->categ = EPS_CATEGORY_OTHER;
300: eps->ops->solve = EPSSolve_Subspace;
301: eps->ops->setup = EPSSetUp_Subspace;
302: eps->ops->destroy = EPSDestroy_Subspace;
303: eps->ops->backtransform = EPSBackTransform_Default;
304: eps->ops->computevectors = EPSComputeVectors_Schur;
305: return(0);
306: }