Actual source code: epskrylov.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: 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 - Implements the loop that checks for convergence
159: in Krylov methods.
161: Input Parameters:
162: eps - the eigensolver; some error estimates are updated in eps->errest
163: getall - whether all residuals must be computed
164: kini - initial value of k (the loop variable)
165: nits - number of iterations of the loop
166: beta - norm of f (the residual vector of the Arnoldi/Lanczos factorization)
167: betat - norm of the residual of the second Krylov relation in two-sided methods
168: corrf - correction factor for residual estimates (only in harmonic KS)
170: Output Parameter:
171: kout - the first index where the convergence test failed
172: */
173: PetscErrorCode EPSKrylovConvergence(EPS eps,PetscBool getall,PetscInt kini,PetscInt nits,PetscReal beta,PetscReal betat,PetscReal corrf,PetscInt *kout)
174: {
175: PetscInt k,newk,newk2,marker,ld,inside;
176: PetscScalar re,im,*Zr,*Zi,*X;
177: PetscReal resnorm,lerrest;
178: PetscBool isshift,refined,istrivial;
179: Vec x=NULL,y=NULL,w[3];
181: PetscFunctionBegin;
182: PetscCall(RGIsTrivial(eps->rg,&istrivial));
183: if (PetscUnlikely(eps->trueres)) {
184: PetscCall(BVCreateVec(eps->V,&x));
185: PetscCall(BVCreateVec(eps->V,&y));
186: PetscCall(BVCreateVec(eps->V,&w[0]));
187: PetscCall(BVCreateVec(eps->V,&w[2]));
188: #if !defined(PETSC_USE_COMPLEX)
189: PetscCall(BVCreateVec(eps->V,&w[1]));
190: #else
191: w[1] = NULL;
192: #endif
193: }
194: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
195: PetscCall(DSGetRefined(eps->ds,&refined));
196: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STSHIFT,&isshift));
197: marker = -1;
198: if (eps->trackall) getall = PETSC_TRUE;
199: for (k=kini;k<kini+nits;k++) {
200: /* eigenvalue */
201: re = eps->eigr[k];
202: im = eps->eigi[k];
203: if (!istrivial || eps->trueres || isshift || eps->conv==EPS_CONV_NORM) PetscCall(STBackTransform(eps->st,1,&re,&im));
204: if (PetscUnlikely(!istrivial)) {
205: PetscCall(RGCheckInside(eps->rg,1,&re,&im,&inside));
206: if (marker==-1 && inside<0) marker = k;
207: if (!(eps->trueres || isshift || eps->conv==EPS_CONV_NORM)) { /* make sure eps->converged below uses the right value */
208: re = eps->eigr[k];
209: im = eps->eigi[k];
210: }
211: }
212: newk = k;
213: PetscCall(DSVectors(eps->ds,DS_MAT_X,&newk,&resnorm));
214: if (PetscUnlikely(eps->trueres)) {
215: PetscCall(DSGetArray(eps->ds,DS_MAT_X,&X));
216: Zr = X+k*ld;
217: if (newk==k+1) Zi = X+newk*ld;
218: else Zi = NULL;
219: PetscCall(EPSComputeRitzVector(eps,Zr,Zi,eps->V,x,y));
220: PetscCall(DSRestoreArray(eps->ds,DS_MAT_X,&X));
221: PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_FALSE,re,im,x,y,w,&resnorm));
222: }
223: else if (!refined) resnorm *= beta*corrf;
224: /* error estimate */
225: PetscCall((*eps->converged)(eps,re,im,resnorm,&eps->errest[k],eps->convergedctx));
226: if (marker==-1 && eps->errest[k] >= eps->tol) marker = k;
227: if (PetscUnlikely(eps->twosided)) {
228: newk2 = k;
229: PetscCall(DSVectors(eps->ds,DS_MAT_Y,&newk2,&resnorm));
230: resnorm *= betat;
231: PetscCall((*eps->converged)(eps,re,im,resnorm,&lerrest,eps->convergedctx));
232: eps->errest[k] = PetscMax(eps->errest[k],lerrest);
233: if (marker==-1 && lerrest >= eps->tol) marker = k;
234: }
235: if (PetscUnlikely(newk==k+1)) {
236: eps->errest[k+1] = eps->errest[k];
237: k++;
238: }
239: if (marker!=-1 && !getall) break;
240: }
241: if (marker!=-1) k = marker;
242: *kout = k;
243: if (PetscUnlikely(eps->trueres)) {
244: PetscCall(VecDestroy(&x));
245: PetscCall(VecDestroy(&y));
246: PetscCall(VecDestroy(&w[0]));
247: PetscCall(VecDestroy(&w[2]));
248: #if !defined(PETSC_USE_COMPLEX)
249: PetscCall(VecDestroy(&w[1]));
250: #endif
251: }
252: PetscFunctionReturn(PETSC_SUCCESS);
253: }
255: PetscErrorCode EPSPseudoLanczos(EPS eps,PetscReal *alpha,PetscReal *beta,PetscReal *omega,PetscInt k,PetscInt *M,PetscBool *breakdown,PetscBool *symmlost,PetscReal *cos,Vec w)
256: {
257: PetscInt j,m = *M,i,ld,l;
258: Vec vj,vj1;
259: PetscScalar *hwork,lhwork[100];
260: PetscReal norm,norm1,norm2,t,sym=0.0,fro=0.0;
261: PetscBLASInt j_,one=1;
263: PetscFunctionBegin;
264: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
265: PetscCall(DSGetDimensions(eps->ds,NULL,&l,NULL,NULL));
266: if (cos) *cos = 1.0;
267: if (m > 100) PetscCall(PetscMalloc1(m,&hwork));
268: else hwork = lhwork;
270: PetscCall(BVSetActiveColumns(eps->V,0,m));
271: for (j=k;j<m;j++) {
272: PetscCall(BVGetColumn(eps->V,j,&vj));
273: PetscCall(BVGetColumn(eps->V,j+1,&vj1));
274: PetscCall(STApply(eps->st,vj,vj1));
275: PetscCall(BVRestoreColumn(eps->V,j,&vj));
276: PetscCall(BVRestoreColumn(eps->V,j+1,&vj1));
277: PetscCall(BVOrthogonalizeColumn(eps->V,j+1,hwork,&norm,breakdown));
278: alpha[j] = PetscRealPart(hwork[j]);
279: beta[j] = PetscAbsReal(norm);
280: if (j==k) {
281: PetscReal *f;
283: PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&f));
284: for (i=0;i<l;i++) hwork[i] = 0.0;
285: for (;i<j-1;i++) hwork[i] -= f[2*ld+i];
286: PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&f));
287: }
288: if (j>0) {
289: hwork[j-1] -= beta[j-1];
290: PetscCall(PetscBLASIntCast(j,&j_));
291: sym = SlepcAbs(BLASnrm2_(&j_,hwork,&one),sym);
292: }
293: fro = SlepcAbs(fro,SlepcAbs(alpha[j],beta[j]));
294: if (j>0) fro = SlepcAbs(fro,beta[j-1]);
295: if (sym/fro>PetscMax(PETSC_SQRT_MACHINE_EPSILON,10*eps->tol)) { *symmlost = PETSC_TRUE; *M=j; break; }
296: omega[j+1] = (norm<0.0)? -1.0: 1.0;
297: PetscCall(BVScaleColumn(eps->V,j+1,1.0/norm));
298: /* */
299: if (cos) {
300: PetscCall(BVGetColumn(eps->V,j+1,&vj1));
301: PetscCall(VecNorm(vj1,NORM_2,&norm1));
302: PetscCall(BVApplyMatrix(eps->V,vj1,w));
303: PetscCall(BVRestoreColumn(eps->V,j+1,&vj1));
304: PetscCall(VecNorm(w,NORM_2,&norm2));
305: t = 1.0/(norm1*norm2);
306: if (*cos>t) *cos = t;
307: }
308: }
309: if (m > 100) PetscCall(PetscFree(hwork));
310: PetscFunctionReturn(PETSC_SUCCESS);
311: }