Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
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 : */
13 :
14 : #include <slepc/private/epsimpl.h>
15 : #include <slepc/private/slepcimpl.h>
16 : #include <slepcblaslapack.h>
17 :
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 41 : PetscErrorCode EPSDelayedArnoldi(EPS eps,PetscScalar *H,PetscInt ldh,PetscInt k,PetscInt *M,PetscReal *beta,PetscBool *breakdown)
25 : {
26 41 : PetscInt i,j,m=*M;
27 41 : Vec u,t;
28 41 : PetscScalar shh[100],*lhh,dot,dot2;
29 41 : PetscReal norm1=0.0,norm2=1.0;
30 41 : Vec vj,vj1,vj2=NULL;
31 :
32 41 : PetscFunctionBegin;
33 41 : if (m<=100) lhh = shh;
34 0 : else PetscCall(PetscMalloc1(m,&lhh));
35 41 : PetscCall(BVCreateVec(eps->V,&u));
36 41 : PetscCall(BVCreateVec(eps->V,&t));
37 :
38 41 : PetscCall(BVSetActiveColumns(eps->V,0,m));
39 752 : for (j=k;j<m;j++) {
40 711 : PetscCall(BVGetColumn(eps->V,j,&vj));
41 711 : PetscCall(BVGetColumn(eps->V,j+1,&vj1));
42 711 : PetscCall(STApply(eps->st,vj,vj1));
43 711 : PetscCall(BVRestoreColumn(eps->V,j,&vj));
44 711 : PetscCall(BVRestoreColumn(eps->V,j+1,&vj1));
45 :
46 711 : PetscCall(BVDotColumnBegin(eps->V,j+1,H+ldh*j));
47 711 : if (j>k) {
48 670 : PetscCall(BVDotColumnBegin(eps->V,j,lhh));
49 670 : PetscCall(BVGetColumn(eps->V,j,&vj));
50 670 : PetscCall(VecDotBegin(vj,vj,&dot));
51 670 : if (j>k+1) {
52 629 : PetscCall(BVNormVecBegin(eps->V,u,NORM_2,&norm2));
53 629 : PetscCall(BVGetColumn(eps->V,j-2,&vj2));
54 629 : PetscCall(VecDotBegin(u,vj2,&dot2));
55 : }
56 670 : PetscCall(BVDotColumnEnd(eps->V,j+1,H+ldh*j));
57 670 : PetscCall(BVDotColumnEnd(eps->V,j,lhh));
58 670 : PetscCall(VecDotEnd(vj,vj,&dot));
59 670 : PetscCall(BVRestoreColumn(eps->V,j,&vj));
60 670 : if (j>k+1) {
61 629 : PetscCall(BVNormVecEnd(eps->V,u,NORM_2,&norm2));
62 629 : PetscCall(VecDotEnd(u,vj2,&dot2));
63 629 : PetscCall(BVRestoreColumn(eps->V,j-2,&vj2));
64 : }
65 670 : norm1 = PetscSqrtReal(PetscRealPart(dot));
66 7278 : for (i=0;i<j;i++) H[ldh*j+i] = H[ldh*j+i]/norm1;
67 670 : H[ldh*j+j] = H[ldh*j+j]/dot;
68 670 : PetscCall(BVCopyVec(eps->V,j,t));
69 670 : PetscCall(BVScaleColumn(eps->V,j,1.0/norm1));
70 670 : PetscCall(BVScaleColumn(eps->V,j+1,1.0/norm1));
71 41 : } else PetscCall(BVDotColumnEnd(eps->V,j+1,H+ldh*j)); /* j==k */
72 :
73 711 : PetscCall(BVMultColumn(eps->V,-1.0,1.0,j+1,H+ldh*j));
74 :
75 711 : if (j>k) {
76 670 : PetscCall(BVSetActiveColumns(eps->V,0,j));
77 670 : PetscCall(BVMultVec(eps->V,-1.0,1.0,t,lhh));
78 670 : PetscCall(BVSetActiveColumns(eps->V,0,m));
79 7278 : for (i=0;i<j;i++) H[ldh*(j-1)+i] += lhh[i];
80 : }
81 :
82 711 : if (j>k+1) {
83 629 : PetscCall(BVGetColumn(eps->V,j-1,&vj1));
84 629 : PetscCall(VecCopy(u,vj1));
85 629 : PetscCall(BVRestoreColumn(eps->V,j-1,&vj1));
86 629 : PetscCall(BVScaleColumn(eps->V,j-1,1.0/norm2));
87 629 : H[ldh*(j-2)+j-1] = norm2;
88 : }
89 :
90 711 : if (j<m-1) PetscCall(VecCopy(t,u));
91 : }
92 :
93 41 : PetscCall(BVNormVec(eps->V,t,NORM_2,&norm2));
94 41 : PetscCall(VecScale(t,1.0/norm2));
95 41 : PetscCall(BVGetColumn(eps->V,m-1,&vj1));
96 41 : PetscCall(VecCopy(t,vj1));
97 41 : PetscCall(BVRestoreColumn(eps->V,m-1,&vj1));
98 41 : H[ldh*(m-2)+m-1] = norm2;
99 :
100 41 : PetscCall(BVDotColumn(eps->V,m,lhh));
101 :
102 41 : PetscCall(BVMultColumn(eps->V,-1.0,1.0,m,lhh));
103 802 : for (i=0;i<m;i++)
104 761 : H[ldh*(m-1)+i] += lhh[i];
105 :
106 41 : PetscCall(BVNormColumn(eps->V,m,NORM_2,beta));
107 41 : PetscCall(BVScaleColumn(eps->V,m,1.0 / *beta));
108 41 : *breakdown = PETSC_FALSE;
109 :
110 41 : if (m>100) PetscCall(PetscFree(lhh));
111 41 : PetscCall(VecDestroy(&u));
112 41 : PetscCall(VecDestroy(&t));
113 41 : PetscFunctionReturn(PETSC_SUCCESS);
114 : }
115 :
116 : /*
117 : EPSDelayedArnoldi1 - This function is similar to EPSDelayedArnoldi,
118 : but without reorthogonalization (only delayed normalization).
119 : */
120 18 : PetscErrorCode EPSDelayedArnoldi1(EPS eps,PetscScalar *H,PetscInt ldh,PetscInt k,PetscInt *M,PetscReal *beta,PetscBool *breakdown)
121 : {
122 18 : PetscInt i,j,m=*M;
123 18 : PetscScalar dot;
124 18 : PetscReal norm=0.0;
125 18 : Vec vj,vj1;
126 :
127 18 : PetscFunctionBegin;
128 18 : PetscCall(BVSetActiveColumns(eps->V,0,m));
129 326 : for (j=k;j<m;j++) {
130 308 : PetscCall(BVGetColumn(eps->V,j,&vj));
131 308 : PetscCall(BVGetColumn(eps->V,j+1,&vj1));
132 308 : PetscCall(STApply(eps->st,vj,vj1));
133 308 : PetscCall(BVRestoreColumn(eps->V,j+1,&vj1));
134 308 : if (j>k) {
135 290 : PetscCall(BVDotColumnBegin(eps->V,j+1,H+ldh*j));
136 290 : PetscCall(VecDotBegin(vj,vj,&dot));
137 290 : PetscCall(BVDotColumnEnd(eps->V,j+1,H+ldh*j));
138 290 : PetscCall(VecDotEnd(vj,vj,&dot));
139 290 : norm = PetscSqrtReal(PetscRealPart(dot));
140 290 : PetscCall(BVScaleColumn(eps->V,j,1.0/norm));
141 290 : H[ldh*(j-1)+j] = norm;
142 3022 : for (i=0;i<j;i++) H[ldh*j+i] = H[ldh*j+i]/norm;
143 290 : H[ldh*j+j] = H[ldh*j+j]/dot;
144 290 : PetscCall(BVScaleColumn(eps->V,j+1,1.0/norm));
145 290 : *beta = norm;
146 : } else { /* j==k */
147 18 : PetscCall(BVDotColumn(eps->V,j+1,H+ldh*j));
148 : }
149 308 : PetscCall(BVRestoreColumn(eps->V,j,&vj));
150 308 : PetscCall(BVMultColumn(eps->V,-1.0,1.0,j+1,H+ldh*j));
151 : }
152 :
153 18 : *breakdown = PETSC_FALSE;
154 18 : PetscFunctionReturn(PETSC_SUCCESS);
155 : }
156 :
157 : /*
158 : EPSKrylovConvergence_Filter - Specialized version for STFILTER.
159 : */
160 13 : static PetscErrorCode EPSKrylovConvergence_Filter(EPS eps,PetscBool getall,PetscInt kini,PetscInt nits,PetscReal beta,PetscReal gamma,PetscInt *kout)
161 : {
162 13 : PetscInt k,ninside,nconv;
163 13 : PetscScalar re,im;
164 13 : PetscReal resnorm;
165 :
166 13 : PetscFunctionBegin;
167 13 : ninside = 0; /* count how many eigenvalues are located in the interval */
168 113 : for (k=kini;k<kini+nits;k++) {
169 113 : if (PetscRealPart(eps->eigr[k]) < gamma) break;
170 100 : ninside++;
171 : }
172 13 : eps->nev = ninside+kini; /* adjust eigenvalue count */
173 13 : nconv = 0; /* count how many eigenvalues satisfy the convergence criterion */
174 76 : for (k=kini;k<kini+ninside;k++) {
175 : /* eigenvalue */
176 68 : re = eps->eigr[k];
177 68 : im = eps->eigi[k];
178 68 : PetscCall(DSVectors(eps->ds,DS_MAT_X,&k,&resnorm));
179 68 : resnorm *= beta;
180 : /* error estimate */
181 68 : PetscCall((*eps->converged)(eps,re,im,resnorm,&eps->errest[k],eps->convergedctx));
182 68 : if (eps->errest[k] < eps->tol) nconv++;
183 : else break;
184 : }
185 13 : *kout = kini+nconv;
186 13 : 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 13 : PetscFunctionReturn(PETSC_SUCCESS);
188 : }
189 :
190 : /*
191 : EPSKrylovConvergence - Implements the loop that checks for convergence
192 : in Krylov methods.
193 :
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)
203 :
204 : Output Parameters:
205 : kout - the first index where the convergence test failed
206 : */
207 3196 : PetscErrorCode EPSKrylovConvergence(EPS eps,PetscBool getall,PetscInt kini,PetscInt nits,PetscReal beta,PetscReal betat,PetscReal corrf,PetscInt *kout)
208 : {
209 3196 : PetscInt k,newk,newk2,marker,ld,inside;
210 3196 : PetscScalar re,im,*Zr,*Zi,*X;
211 3196 : PetscReal resnorm,gamma,lerrest;
212 3196 : PetscBool isshift,isfilter,refined,istrivial;
213 3196 : Vec x=NULL,y=NULL,w[3];
214 :
215 3196 : PetscFunctionBegin;
216 3196 : if (PetscUnlikely(eps->which == EPS_ALL)) {
217 81 : PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilter));
218 81 : if (isfilter) {
219 13 : PetscCall(STFilterGetThreshold(eps->st,&gamma));
220 13 : PetscCall(EPSKrylovConvergence_Filter(eps,getall,kini,nits,beta,gamma,kout));
221 13 : PetscFunctionReturn(PETSC_SUCCESS);
222 : }
223 : }
224 3183 : PetscCall(RGIsTrivial(eps->rg,&istrivial));
225 3183 : if (PetscUnlikely(eps->trueres)) {
226 36 : PetscCall(BVCreateVec(eps->V,&x));
227 36 : PetscCall(BVCreateVec(eps->V,&y));
228 36 : PetscCall(BVCreateVec(eps->V,&w[0]));
229 36 : PetscCall(BVCreateVec(eps->V,&w[2]));
230 : #if !defined(PETSC_USE_COMPLEX)
231 : PetscCall(BVCreateVec(eps->V,&w[1]));
232 : #else
233 36 : w[1] = NULL;
234 : #endif
235 : }
236 3183 : PetscCall(DSGetLeadingDimension(eps->ds,&ld));
237 3183 : PetscCall(DSGetRefined(eps->ds,&refined));
238 3183 : PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STSHIFT,&isshift));
239 3183 : marker = -1;
240 3183 : if (eps->trackall) getall = PETSC_TRUE;
241 11732 : for (k=kini;k<kini+nits;k++) {
242 : /* eigenvalue */
243 11598 : re = eps->eigr[k];
244 11598 : im = eps->eigi[k];
245 11598 : if (!istrivial || eps->trueres || isshift || eps->conv==EPS_CONV_NORM) PetscCall(STBackTransform(eps->st,1,&re,&im));
246 11598 : if (PetscUnlikely(!istrivial)) {
247 121 : PetscCall(RGCheckInside(eps->rg,1,&re,&im,&inside));
248 121 : if (marker==-1 && inside<0) marker = k;
249 121 : if (!(eps->trueres || isshift || eps->conv==EPS_CONV_NORM)) { /* make sure eps->converged below uses the right value */
250 43 : re = eps->eigr[k];
251 43 : im = eps->eigi[k];
252 : }
253 : }
254 11598 : newk = k;
255 11598 : PetscCall(DSVectors(eps->ds,DS_MAT_X,&newk,&resnorm));
256 11598 : if (PetscUnlikely(eps->trueres)) {
257 55 : PetscCall(DSGetArray(eps->ds,DS_MAT_X,&X));
258 55 : Zr = X+k*ld;
259 55 : if (newk==k+1) Zi = X+newk*ld;
260 : else Zi = NULL;
261 55 : PetscCall(EPSComputeRitzVector(eps,Zr,Zi,eps->V,x,y));
262 55 : PetscCall(DSRestoreArray(eps->ds,DS_MAT_X,&X));
263 55 : PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_FALSE,re,im,x,y,w,&resnorm));
264 : }
265 11543 : else if (!refined) resnorm *= beta*corrf;
266 : /* error estimate */
267 11598 : PetscCall((*eps->converged)(eps,re,im,resnorm,&eps->errest[k],eps->convergedctx));
268 11598 : if (marker==-1 && eps->errest[k] >= eps->tol) marker = k;
269 11598 : if (PetscUnlikely(eps->twosided)) {
270 161 : newk2 = k;
271 161 : PetscCall(DSVectors(eps->ds,DS_MAT_Y,&newk2,&resnorm));
272 161 : resnorm *= betat;
273 161 : PetscCall((*eps->converged)(eps,re,im,resnorm,&lerrest,eps->convergedctx));
274 161 : eps->errest[k] = PetscMax(eps->errest[k],lerrest);
275 161 : if (marker==-1 && lerrest >= eps->tol) marker = k;
276 : }
277 11598 : if (PetscUnlikely(newk==k+1)) {
278 9 : eps->errest[k+1] = eps->errest[k];
279 9 : k++;
280 : }
281 11598 : if (marker!=-1 && !getall) break;
282 : }
283 3183 : if (marker!=-1) k = marker;
284 3183 : *kout = k;
285 3183 : if (PetscUnlikely(eps->trueres)) {
286 36 : PetscCall(VecDestroy(&x));
287 36 : PetscCall(VecDestroy(&y));
288 36 : PetscCall(VecDestroy(&w[0]));
289 36 : PetscCall(VecDestroy(&w[2]));
290 : #if !defined(PETSC_USE_COMPLEX)
291 : PetscCall(VecDestroy(&w[1]));
292 : #endif
293 : }
294 3183 : PetscFunctionReturn(PETSC_SUCCESS);
295 : }
296 :
297 61 : PetscErrorCode EPSPseudoLanczos(EPS eps,PetscReal *alpha,PetscReal *beta,PetscReal *omega,PetscInt k,PetscInt *M,PetscBool *breakdown,PetscBool *symmlost,PetscReal *cos,Vec w)
298 : {
299 61 : PetscInt j,m = *M,i,ld,l;
300 61 : Vec vj,vj1;
301 61 : PetscScalar *hwork,lhwork[100];
302 61 : PetscReal norm,norm1,norm2,t,sym=0.0,fro=0.0;
303 61 : PetscBLASInt j_,one=1;
304 :
305 61 : PetscFunctionBegin;
306 61 : PetscCall(DSGetLeadingDimension(eps->ds,&ld));
307 61 : PetscCall(DSGetDimensions(eps->ds,NULL,&l,NULL,NULL));
308 61 : if (cos) *cos = 1.0;
309 61 : if (m > 100) PetscCall(PetscMalloc1(m,&hwork));
310 53 : else hwork = lhwork;
311 :
312 61 : PetscCall(BVSetActiveColumns(eps->V,0,m));
313 2027 : for (j=k;j<m;j++) {
314 1966 : PetscCall(BVGetColumn(eps->V,j,&vj));
315 1966 : PetscCall(BVGetColumn(eps->V,j+1,&vj1));
316 1966 : PetscCall(STApply(eps->st,vj,vj1));
317 1966 : PetscCall(BVRestoreColumn(eps->V,j,&vj));
318 1966 : PetscCall(BVRestoreColumn(eps->V,j+1,&vj1));
319 1966 : PetscCall(BVOrthogonalizeColumn(eps->V,j+1,hwork,&norm,breakdown));
320 1966 : alpha[j] = PetscRealPart(hwork[j]);
321 1966 : beta[j] = PetscAbsReal(norm);
322 1966 : if (j==k) {
323 61 : PetscReal *f;
324 :
325 61 : PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&f));
326 526 : for (i=0;i<l;i++) hwork[i] = 0.0;
327 1028 : for (;i<j-1;i++) hwork[i] -= f[2*ld+i];
328 61 : PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&f));
329 : }
330 1966 : if (j>0) {
331 1953 : hwork[j-1] -= beta[j-1];
332 1953 : PetscCall(PetscBLASIntCast(j,&j_));
333 1953 : sym = SlepcAbs(BLASnrm2_(&j_,hwork,&one),sym);
334 : }
335 1966 : fro = SlepcAbs(fro,SlepcAbs(alpha[j],beta[j]));
336 1966 : if (j>0) fro = SlepcAbs(fro,beta[j-1]);
337 2038 : if (sym/fro>PetscMax(PETSC_SQRT_MACHINE_EPSILON,10*eps->tol)) { *symmlost = PETSC_TRUE; *M=j; break; }
338 1966 : omega[j+1] = (norm<0.0)? -1.0: 1.0;
339 1966 : PetscCall(BVScaleColumn(eps->V,j+1,1.0/norm));
340 : /* */
341 1966 : if (cos) {
342 0 : PetscCall(BVGetColumn(eps->V,j+1,&vj1));
343 0 : PetscCall(VecNorm(vj1,NORM_2,&norm1));
344 0 : PetscCall(BVApplyMatrix(eps->V,vj1,w));
345 0 : PetscCall(BVRestoreColumn(eps->V,j+1,&vj1));
346 0 : PetscCall(VecNorm(w,NORM_2,&norm2));
347 0 : t = 1.0/(norm1*norm2);
348 0 : if (*cos>t) *cos = t;
349 : }
350 : }
351 61 : if (m > 100) PetscCall(PetscFree(hwork));
352 61 : PetscFunctionReturn(PETSC_SUCCESS);
353 : }
|
