Actual source code: epskrylov.c
slepc-3.21.1 2024-04-26
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: Common subroutines for all Krylov-type solvers
12: */
14: #include <slepc/private/epsimpl.h>
15: #include <slepc/private/slepcimpl.h>
16: #include <slepcblaslapack.h>
18: /*
19: EPSDelayedArnoldi - This function is equivalent to BVMatArnoldi but
20: performs the computation in a different way. The main idea is that
21: reorthogonalization is delayed to the next Arnoldi step. This version is
22: more scalable but in some cases convergence may stagnate.
23: */
24: PetscErrorCode EPSDelayedArnoldi(EPS eps,PetscScalar *H,PetscInt ldh,PetscInt k,PetscInt *M,PetscReal *beta,PetscBool *breakdown)
25: {
26: PetscInt i,j,m=*M;
27: Vec u,t;
28: PetscScalar shh[100],*lhh,dot,dot2;
29: PetscReal norm1=0.0,norm2=1.0;
30: Vec vj,vj1,vj2=NULL;
32: PetscFunctionBegin;
33: if (m<=100) lhh = shh;
34: else PetscCall(PetscMalloc1(m,&lhh));
35: PetscCall(BVCreateVec(eps->V,&u));
36: PetscCall(BVCreateVec(eps->V,&t));
38: PetscCall(BVSetActiveColumns(eps->V,0,m));
39: for (j=k;j<m;j++) {
40: PetscCall(BVGetColumn(eps->V,j,&vj));
41: PetscCall(BVGetColumn(eps->V,j+1,&vj1));
42: PetscCall(STApply(eps->st,vj,vj1));
43: PetscCall(BVRestoreColumn(eps->V,j,&vj));
44: PetscCall(BVRestoreColumn(eps->V,j+1,&vj1));
46: PetscCall(BVDotColumnBegin(eps->V,j+1,H+ldh*j));
47: if (j>k) {
48: PetscCall(BVDotColumnBegin(eps->V,j,lhh));
49: PetscCall(BVGetColumn(eps->V,j,&vj));
50: PetscCall(VecDotBegin(vj,vj,&dot));
51: if (j>k+1) {
52: PetscCall(BVNormVecBegin(eps->V,u,NORM_2,&norm2));
53: PetscCall(BVGetColumn(eps->V,j-2,&vj2));
54: PetscCall(VecDotBegin(u,vj2,&dot2));
55: }
56: PetscCall(BVDotColumnEnd(eps->V,j+1,H+ldh*j));
57: PetscCall(BVDotColumnEnd(eps->V,j,lhh));
58: PetscCall(VecDotEnd(vj,vj,&dot));
59: PetscCall(BVRestoreColumn(eps->V,j,&vj));
60: if (j>k+1) {
61: PetscCall(BVNormVecEnd(eps->V,u,NORM_2,&norm2));
62: PetscCall(VecDotEnd(u,vj2,&dot2));
63: PetscCall(BVRestoreColumn(eps->V,j-2,&vj2));
64: }
65: norm1 = PetscSqrtReal(PetscRealPart(dot));
66: for (i=0;i<j;i++) H[ldh*j+i] = H[ldh*j+i]/norm1;
67: H[ldh*j+j] = H[ldh*j+j]/dot;
68: PetscCall(BVCopyVec(eps->V,j,t));
69: PetscCall(BVScaleColumn(eps->V,j,1.0/norm1));
70: PetscCall(BVScaleColumn(eps->V,j+1,1.0/norm1));
71: } else PetscCall(BVDotColumnEnd(eps->V,j+1,H+ldh*j)); /* j==k */
73: PetscCall(BVMultColumn(eps->V,-1.0,1.0,j+1,H+ldh*j));
75: if (j>k) {
76: PetscCall(BVSetActiveColumns(eps->V,0,j));
77: PetscCall(BVMultVec(eps->V,-1.0,1.0,t,lhh));
78: PetscCall(BVSetActiveColumns(eps->V,0,m));
79: for (i=0;i<j;i++) H[ldh*(j-1)+i] += lhh[i];
80: }
82: if (j>k+1) {
83: PetscCall(BVGetColumn(eps->V,j-1,&vj1));
84: PetscCall(VecCopy(u,vj1));
85: PetscCall(BVRestoreColumn(eps->V,j-1,&vj1));
86: PetscCall(BVScaleColumn(eps->V,j-1,1.0/norm2));
87: H[ldh*(j-2)+j-1] = norm2;
88: }
90: if (j<m-1) PetscCall(VecCopy(t,u));
91: }
93: PetscCall(BVNormVec(eps->V,t,NORM_2,&norm2));
94: PetscCall(VecScale(t,1.0/norm2));
95: PetscCall(BVGetColumn(eps->V,m-1,&vj1));
96: PetscCall(VecCopy(t,vj1));
97: PetscCall(BVRestoreColumn(eps->V,m-1,&vj1));
98: H[ldh*(m-2)+m-1] = norm2;
100: PetscCall(BVDotColumn(eps->V,m,lhh));
102: PetscCall(BVMultColumn(eps->V,-1.0,1.0,m,lhh));
103: for (i=0;i<m;i++)
104: H[ldh*(m-1)+i] += lhh[i];
106: PetscCall(BVNormColumn(eps->V,m,NORM_2,beta));
107: PetscCall(BVScaleColumn(eps->V,m,1.0 / *beta));
108: *breakdown = PETSC_FALSE;
110: if (m>100) PetscCall(PetscFree(lhh));
111: PetscCall(VecDestroy(&u));
112: PetscCall(VecDestroy(&t));
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: /*
117: EPSDelayedArnoldi1 - This function is similar to EPSDelayedArnoldi,
118: but without reorthogonalization (only delayed normalization).
119: */
120: PetscErrorCode EPSDelayedArnoldi1(EPS eps,PetscScalar *H,PetscInt ldh,PetscInt k,PetscInt *M,PetscReal *beta,PetscBool *breakdown)
121: {
122: PetscInt i,j,m=*M;
123: PetscScalar dot;
124: PetscReal norm=0.0;
125: Vec vj,vj1;
127: PetscFunctionBegin;
128: PetscCall(BVSetActiveColumns(eps->V,0,m));
129: for (j=k;j<m;j++) {
130: PetscCall(BVGetColumn(eps->V,j,&vj));
131: PetscCall(BVGetColumn(eps->V,j+1,&vj1));
132: PetscCall(STApply(eps->st,vj,vj1));
133: PetscCall(BVRestoreColumn(eps->V,j+1,&vj1));
134: if (j>k) {
135: PetscCall(BVDotColumnBegin(eps->V,j+1,H+ldh*j));
136: PetscCall(VecDotBegin(vj,vj,&dot));
137: PetscCall(BVDotColumnEnd(eps->V,j+1,H+ldh*j));
138: PetscCall(VecDotEnd(vj,vj,&dot));
139: norm = PetscSqrtReal(PetscRealPart(dot));
140: PetscCall(BVScaleColumn(eps->V,j,1.0/norm));
141: H[ldh*(j-1)+j] = norm;
142: for (i=0;i<j;i++) H[ldh*j+i] = H[ldh*j+i]/norm;
143: H[ldh*j+j] = H[ldh*j+j]/dot;
144: PetscCall(BVScaleColumn(eps->V,j+1,1.0/norm));
145: *beta = norm;
146: } else { /* j==k */
147: PetscCall(BVDotColumn(eps->V,j+1,H+ldh*j));
148: }
149: PetscCall(BVRestoreColumn(eps->V,j,&vj));
150: PetscCall(BVMultColumn(eps->V,-1.0,1.0,j+1,H+ldh*j));
151: }
153: *breakdown = PETSC_FALSE;
154: PetscFunctionReturn(PETSC_SUCCESS);
155: }
157: /*
158: EPSKrylovConvergence_Filter - Specialized version for STFILTER.
159: */
160: static PetscErrorCode EPSKrylovConvergence_Filter(EPS eps,PetscBool getall,PetscInt kini,PetscInt nits,PetscReal beta,PetscReal gamma,PetscInt *kout)
161: {
162: PetscInt k,ninside,nconv;
163: PetscScalar re,im;
164: PetscReal resnorm;
166: PetscFunctionBegin;
167: ninside = 0; /* count how many eigenvalues are located in the interval */
168: for (k=kini;k<kini+nits;k++) {
169: if (PetscRealPart(eps->eigr[k]) < gamma) break;
170: ninside++;
171: }
172: eps->nev = ninside+kini; /* adjust eigenvalue count */
173: nconv = 0; /* count how many eigenvalues satisfy the convergence criterion */
174: for (k=kini;k<kini+ninside;k++) {
175: /* eigenvalue */
176: re = eps->eigr[k];
177: im = eps->eigi[k];
178: PetscCall(DSVectors(eps->ds,DS_MAT_X,&k,&resnorm));
179: resnorm *= beta;
180: /* error estimate */
181: PetscCall((*eps->converged)(eps,re,im,resnorm,&eps->errest[k],eps->convergedctx));
182: if (eps->errest[k] < eps->tol) nconv++;
183: else break;
184: }
185: *kout = kini+nconv;
186: PetscCall(PetscInfo(eps,"Found %" PetscInt_FMT " eigenvalue approximations inside the interval (gamma=%g), k=%" PetscInt_FMT " nconv=%" PetscInt_FMT "\n",ninside,(double)gamma,k,nconv));
187: PetscFunctionReturn(PETSC_SUCCESS);
188: }
190: /*
191: EPSKrylovConvergence - Implements the loop that checks for convergence
192: in Krylov methods.
194: Input Parameters:
195: eps - the eigensolver; some error estimates are updated in eps->errest
196: getall - whether all residuals must be computed
197: kini - initial value of k (the loop variable)
198: nits - number of iterations of the loop
199: V - set of basis vectors (used only if trueresidual is activated)
200: nv - number of vectors to process (dimension of Q, columns of V)
201: beta - norm of f (the residual vector of the Arnoldi/Lanczos factorization)
202: corrf - correction factor for residual estimates (only in harmonic KS)
204: Output Parameters:
205: kout - the first index where the convergence test failed
206: */
207: PetscErrorCode EPSKrylovConvergence(EPS eps,PetscBool getall,PetscInt kini,PetscInt nits,PetscReal beta,PetscReal betat,PetscReal corrf,PetscInt *kout)
208: {
209: PetscInt k,newk,newk2,marker,ld,inside;
210: PetscScalar re,im,*Zr,*Zi,*X;
211: PetscReal resnorm,gamma,lerrest;
212: PetscBool isshift,isfilter,refined,istrivial;
213: Vec x=NULL,y=NULL,w[3];
215: PetscFunctionBegin;
216: if (PetscUnlikely(eps->which == EPS_ALL)) {
217: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilter));
218: if (isfilter) {
219: PetscCall(STFilterGetThreshold(eps->st,&gamma));
220: PetscCall(EPSKrylovConvergence_Filter(eps,getall,kini,nits,beta,gamma,kout));
221: PetscFunctionReturn(PETSC_SUCCESS);
222: }
223: }
224: PetscCall(RGIsTrivial(eps->rg,&istrivial));
225: if (PetscUnlikely(eps->trueres)) {
226: PetscCall(BVCreateVec(eps->V,&x));
227: PetscCall(BVCreateVec(eps->V,&y));
228: PetscCall(BVCreateVec(eps->V,&w[0]));
229: PetscCall(BVCreateVec(eps->V,&w[2]));
230: #if !defined(PETSC_USE_COMPLEX)
231: PetscCall(BVCreateVec(eps->V,&w[1]));
232: #else
233: w[1] = NULL;
234: #endif
235: }
236: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
237: PetscCall(DSGetRefined(eps->ds,&refined));
238: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STSHIFT,&isshift));
239: marker = -1;
240: if (eps->trackall) getall = PETSC_TRUE;
241: for (k=kini;k<kini+nits;k++) {
242: /* eigenvalue */
243: re = eps->eigr[k];
244: im = eps->eigi[k];
245: if (!istrivial || eps->trueres || isshift || eps->conv==EPS_CONV_NORM) PetscCall(STBackTransform(eps->st,1,&re,&im));
246: if (PetscUnlikely(!istrivial)) {
247: PetscCall(RGCheckInside(eps->rg,1,&re,&im,&inside));
248: if (marker==-1 && inside<0) marker = k;
249: if (!(eps->trueres || isshift || eps->conv==EPS_CONV_NORM)) { /* make sure eps->converged below uses the right value */
250: re = eps->eigr[k];
251: im = eps->eigi[k];
252: }
253: }
254: newk = k;
255: PetscCall(DSVectors(eps->ds,DS_MAT_X,&newk,&resnorm));
256: if (PetscUnlikely(eps->trueres)) {
257: PetscCall(DSGetArray(eps->ds,DS_MAT_X,&X));
258: Zr = X+k*ld;
259: if (newk==k+1) Zi = X+newk*ld;
260: else Zi = NULL;
261: PetscCall(EPSComputeRitzVector(eps,Zr,Zi,eps->V,x,y));
262: PetscCall(DSRestoreArray(eps->ds,DS_MAT_X,&X));
263: PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_FALSE,re,im,x,y,w,&resnorm));
264: }
265: else if (!refined) resnorm *= beta*corrf;
266: /* error estimate */
267: PetscCall((*eps->converged)(eps,re,im,resnorm,&eps->errest[k],eps->convergedctx));
268: if (marker==-1 && eps->errest[k] >= eps->tol) marker = k;
269: if (PetscUnlikely(eps->twosided)) {
270: newk2 = k;
271: PetscCall(DSVectors(eps->ds,DS_MAT_Y,&newk2,&resnorm));
272: resnorm *= betat;
273: PetscCall((*eps->converged)(eps,re,im,resnorm,&lerrest,eps->convergedctx));
274: eps->errest[k] = PetscMax(eps->errest[k],lerrest);
275: if (marker==-1 && lerrest >= eps->tol) marker = k;
276: }
277: if (PetscUnlikely(newk==k+1)) {
278: eps->errest[k+1] = eps->errest[k];
279: k++;
280: }
281: if (marker!=-1 && !getall) break;
282: }
283: if (marker!=-1) k = marker;
284: *kout = k;
285: if (PetscUnlikely(eps->trueres)) {
286: PetscCall(VecDestroy(&x));
287: PetscCall(VecDestroy(&y));
288: PetscCall(VecDestroy(&w[0]));
289: PetscCall(VecDestroy(&w[2]));
290: #if !defined(PETSC_USE_COMPLEX)
291: PetscCall(VecDestroy(&w[1]));
292: #endif
293: }
294: PetscFunctionReturn(PETSC_SUCCESS);
295: }
297: PetscErrorCode EPSPseudoLanczos(EPS eps,PetscReal *alpha,PetscReal *beta,PetscReal *omega,PetscInt k,PetscInt *M,PetscBool *breakdown,PetscBool *symmlost,PetscReal *cos,Vec w)
298: {
299: PetscInt j,m = *M,i,ld,l;
300: Vec vj,vj1;
301: PetscScalar *hwork,lhwork[100];
302: PetscReal norm,norm1,norm2,t,sym=0.0,fro=0.0;
303: PetscBLASInt j_,one=1;
305: PetscFunctionBegin;
306: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
307: PetscCall(DSGetDimensions(eps->ds,NULL,&l,NULL,NULL));
308: if (cos) *cos = 1.0;
309: if (m > 100) PetscCall(PetscMalloc1(m,&hwork));
310: else hwork = lhwork;
312: PetscCall(BVSetActiveColumns(eps->V,0,m));
313: for (j=k;j<m;j++) {
314: PetscCall(BVGetColumn(eps->V,j,&vj));
315: PetscCall(BVGetColumn(eps->V,j+1,&vj1));
316: PetscCall(STApply(eps->st,vj,vj1));
317: PetscCall(BVRestoreColumn(eps->V,j,&vj));
318: PetscCall(BVRestoreColumn(eps->V,j+1,&vj1));
319: PetscCall(BVOrthogonalizeColumn(eps->V,j+1,hwork,&norm,breakdown));
320: alpha[j] = PetscRealPart(hwork[j]);
321: beta[j] = PetscAbsReal(norm);
322: if (j==k) {
323: PetscReal *f;
325: PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&f));
326: for (i=0;i<l;i++) hwork[i] = 0.0;
327: for (;i<j-1;i++) hwork[i] -= f[2*ld+i];
328: PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&f));
329: }
330: if (j>0) {
331: hwork[j-1] -= beta[j-1];
332: PetscCall(PetscBLASIntCast(j,&j_));
333: sym = SlepcAbs(BLASnrm2_(&j_,hwork,&one),sym);
334: }
335: fro = SlepcAbs(fro,SlepcAbs(alpha[j],beta[j]));
336: if (j>0) fro = SlepcAbs(fro,beta[j-1]);
337: if (sym/fro>PetscMax(PETSC_SQRT_MACHINE_EPSILON,10*eps->tol)) { *symmlost = PETSC_TRUE; *M=j; break; }
338: omega[j+1] = (norm<0.0)? -1.0: 1.0;
339: PetscCall(BVScaleColumn(eps->V,j+1,1.0/norm));
340: /* */
341: if (cos) {
342: PetscCall(BVGetColumn(eps->V,j+1,&vj1));
343: PetscCall(VecNorm(vj1,NORM_2,&norm1));
344: PetscCall(BVApplyMatrix(eps->V,vj1,w));
345: PetscCall(BVRestoreColumn(eps->V,j+1,&vj1));
346: PetscCall(VecNorm(w,NORM_2,&norm2));
347: t = 1.0/(norm1*norm2);
348: if (*cos>t) *cos = t;
349: }
350: }
351: if (m > 100) PetscCall(PetscFree(hwork));
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }