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