Actual source code: epskrylov.c
slepc-3.23.3 2025-09-08
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: if (eps->trackall) getall = PETSC_TRUE;
173: eps->nev = ninside+kini; /* adjust eigenvalue count */
174: nconv = 0; /* count how many eigenvalues satisfy the convergence criterion */
175: for (k=kini;k<kini+(getall?nits:ninside);k++) {
176: /* eigenvalue */
177: re = eps->eigr[k];
178: im = eps->eigi[k];
179: PetscCall(DSVectors(eps->ds,DS_MAT_X,&k,&resnorm));
180: resnorm *= beta;
181: /* error estimate */
182: PetscCall((*eps->converged)(eps,re,im,resnorm,&eps->errest[k],eps->convergedctx));
183: if (eps->errest[k] < eps->tol) nconv++;
184: else if (!getall) break;
185: }
186: *kout = kini+nconv;
187: 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));
188: PetscFunctionReturn(PETSC_SUCCESS);
189: }
191: /*
192: EPSKrylovConvergence - Implements the loop that checks for convergence
193: in Krylov methods.
195: Input Parameters:
196: eps - the eigensolver; some error estimates are updated in eps->errest
197: getall - whether all residuals must be computed
198: kini - initial value of k (the loop variable)
199: nits - number of iterations of the loop
200: V - set of basis vectors (used only if trueresidual is activated)
201: nv - number of vectors to process (dimension of Q, columns of V)
202: beta - norm of f (the residual vector of the Arnoldi/Lanczos factorization)
203: corrf - correction factor for residual estimates (only in harmonic KS)
205: Output Parameters:
206: kout - the first index where the convergence test failed
207: */
208: PetscErrorCode EPSKrylovConvergence(EPS eps,PetscBool getall,PetscInt kini,PetscInt nits,PetscReal beta,PetscReal betat,PetscReal corrf,PetscInt *kout)
209: {
210: PetscInt k,newk,newk2,marker,ld,inside;
211: PetscScalar re,im,*Zr,*Zi,*X;
212: PetscReal resnorm,gamma,lerrest;
213: PetscBool isshift,isfilter,refined,istrivial;
214: Vec x=NULL,y=NULL,w[3];
216: PetscFunctionBegin;
217: if (PetscUnlikely(eps->which == EPS_ALL)) {
218: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilter));
219: if (isfilter) {
220: PetscCall(STFilterGetThreshold(eps->st,&gamma));
221: PetscCall(EPSKrylovConvergence_Filter(eps,getall,kini,nits,beta,gamma,kout));
222: PetscFunctionReturn(PETSC_SUCCESS);
223: }
224: }
225: PetscCall(RGIsTrivial(eps->rg,&istrivial));
226: if (PetscUnlikely(eps->trueres)) {
227: PetscCall(BVCreateVec(eps->V,&x));
228: PetscCall(BVCreateVec(eps->V,&y));
229: PetscCall(BVCreateVec(eps->V,&w[0]));
230: PetscCall(BVCreateVec(eps->V,&w[2]));
231: #if !defined(PETSC_USE_COMPLEX)
232: PetscCall(BVCreateVec(eps->V,&w[1]));
233: #else
234: w[1] = NULL;
235: #endif
236: }
237: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
238: PetscCall(DSGetRefined(eps->ds,&refined));
239: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STSHIFT,&isshift));
240: marker = -1;
241: if (eps->trackall) getall = PETSC_TRUE;
242: for (k=kini;k<kini+nits;k++) {
243: /* eigenvalue */
244: re = eps->eigr[k];
245: im = eps->eigi[k];
246: if (!istrivial || eps->trueres || isshift || eps->conv==EPS_CONV_NORM) PetscCall(STBackTransform(eps->st,1,&re,&im));
247: if (PetscUnlikely(!istrivial)) {
248: PetscCall(RGCheckInside(eps->rg,1,&re,&im,&inside));
249: if (marker==-1 && inside<0) marker = k;
250: if (!(eps->trueres || isshift || eps->conv==EPS_CONV_NORM)) { /* make sure eps->converged below uses the right value */
251: re = eps->eigr[k];
252: im = eps->eigi[k];
253: }
254: }
255: newk = k;
256: PetscCall(DSVectors(eps->ds,DS_MAT_X,&newk,&resnorm));
257: if (PetscUnlikely(eps->trueres)) {
258: PetscCall(DSGetArray(eps->ds,DS_MAT_X,&X));
259: Zr = X+k*ld;
260: if (newk==k+1) Zi = X+newk*ld;
261: else Zi = NULL;
262: PetscCall(EPSComputeRitzVector(eps,Zr,Zi,eps->V,x,y));
263: PetscCall(DSRestoreArray(eps->ds,DS_MAT_X,&X));
264: PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_FALSE,re,im,x,y,w,&resnorm));
265: }
266: else if (!refined) resnorm *= beta*corrf;
267: /* error estimate */
268: PetscCall((*eps->converged)(eps,re,im,resnorm,&eps->errest[k],eps->convergedctx));
269: if (marker==-1 && eps->errest[k] >= eps->tol) marker = k;
270: if (PetscUnlikely(eps->twosided)) {
271: newk2 = k;
272: PetscCall(DSVectors(eps->ds,DS_MAT_Y,&newk2,&resnorm));
273: resnorm *= betat;
274: PetscCall((*eps->converged)(eps,re,im,resnorm,&lerrest,eps->convergedctx));
275: eps->errest[k] = PetscMax(eps->errest[k],lerrest);
276: if (marker==-1 && lerrest >= eps->tol) marker = k;
277: }
278: if (PetscUnlikely(newk==k+1)) {
279: eps->errest[k+1] = eps->errest[k];
280: k++;
281: }
282: if (marker!=-1 && !getall) break;
283: }
284: if (marker!=-1) k = marker;
285: *kout = k;
286: if (PetscUnlikely(eps->trueres)) {
287: PetscCall(VecDestroy(&x));
288: PetscCall(VecDestroy(&y));
289: PetscCall(VecDestroy(&w[0]));
290: PetscCall(VecDestroy(&w[2]));
291: #if !defined(PETSC_USE_COMPLEX)
292: PetscCall(VecDestroy(&w[1]));
293: #endif
294: }
295: PetscFunctionReturn(PETSC_SUCCESS);
296: }
298: PetscErrorCode EPSPseudoLanczos(EPS eps,PetscReal *alpha,PetscReal *beta,PetscReal *omega,PetscInt k,PetscInt *M,PetscBool *breakdown,PetscBool *symmlost,PetscReal *cos,Vec w)
299: {
300: PetscInt j,m = *M,i,ld,l;
301: Vec vj,vj1;
302: PetscScalar *hwork,lhwork[100];
303: PetscReal norm,norm1,norm2,t,sym=0.0,fro=0.0;
304: PetscBLASInt j_,one=1;
306: PetscFunctionBegin;
307: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
308: PetscCall(DSGetDimensions(eps->ds,NULL,&l,NULL,NULL));
309: if (cos) *cos = 1.0;
310: if (m > 100) PetscCall(PetscMalloc1(m,&hwork));
311: else hwork = lhwork;
313: PetscCall(BVSetActiveColumns(eps->V,0,m));
314: for (j=k;j<m;j++) {
315: PetscCall(BVGetColumn(eps->V,j,&vj));
316: PetscCall(BVGetColumn(eps->V,j+1,&vj1));
317: PetscCall(STApply(eps->st,vj,vj1));
318: PetscCall(BVRestoreColumn(eps->V,j,&vj));
319: PetscCall(BVRestoreColumn(eps->V,j+1,&vj1));
320: PetscCall(BVOrthogonalizeColumn(eps->V,j+1,hwork,&norm,breakdown));
321: alpha[j] = PetscRealPart(hwork[j]);
322: beta[j] = PetscAbsReal(norm);
323: if (j==k) {
324: PetscReal *f;
326: PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&f));
327: for (i=0;i<l;i++) hwork[i] = 0.0;
328: for (;i<j-1;i++) hwork[i] -= f[2*ld+i];
329: PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&f));
330: }
331: if (j>0) {
332: hwork[j-1] -= beta[j-1];
333: PetscCall(PetscBLASIntCast(j,&j_));
334: sym = SlepcAbs(BLASnrm2_(&j_,hwork,&one),sym);
335: }
336: fro = SlepcAbs(fro,SlepcAbs(alpha[j],beta[j]));
337: if (j>0) fro = SlepcAbs(fro,beta[j-1]);
338: if (sym/fro>PetscMax(PETSC_SQRT_MACHINE_EPSILON,10*eps->tol)) { *symmlost = PETSC_TRUE; *M=j; break; }
339: omega[j+1] = (norm<0.0)? -1.0: 1.0;
340: PetscCall(BVScaleColumn(eps->V,j+1,1.0/norm));
341: /* */
342: if (cos) {
343: PetscCall(BVGetColumn(eps->V,j+1,&vj1));
344: PetscCall(VecNorm(vj1,NORM_2,&norm1));
345: PetscCall(BVApplyMatrix(eps->V,vj1,w));
346: PetscCall(BVRestoreColumn(eps->V,j+1,&vj1));
347: PetscCall(VecNorm(w,NORM_2,&norm2));
348: t = 1.0/(norm1*norm2);
349: if (*cos>t) *cos = t;
350: }
351: }
352: if (m > 100) PetscCall(PetscFree(hwork));
353: PetscFunctionReturn(PETSC_SUCCESS);
354: }