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 45 : PetscErrorCode EPSDelayedArnoldi(EPS eps,PetscScalar *H,PetscInt ldh,PetscInt k,PetscInt *M,PetscReal *beta,PetscBool *breakdown)
25 : {
26 45 : PetscInt i,j,m=*M;
27 45 : Vec u,t;
28 45 : PetscScalar shh[100],*lhh,dot,dot2;
29 45 : PetscReal norm1=0.0,norm2=1.0;
30 45 : Vec vj,vj1,vj2=NULL;
31 :
32 45 : PetscFunctionBegin;
33 45 : if (m<=100) lhh = shh;
34 0 : else PetscCall(PetscMalloc1(m,&lhh));
35 45 : PetscCall(BVCreateVec(eps->V,&u));
36 45 : PetscCall(BVCreateVec(eps->V,&t));
37 :
38 45 : PetscCall(BVSetActiveColumns(eps->V,0,m));
39 825 : for (j=k;j<m;j++) {
40 780 : PetscCall(BVGetColumn(eps->V,j,&vj));
41 780 : PetscCall(BVGetColumn(eps->V,j+1,&vj1));
42 780 : PetscCall(STApply(eps->st,vj,vj1));
43 780 : PetscCall(BVRestoreColumn(eps->V,j,&vj));
44 780 : PetscCall(BVRestoreColumn(eps->V,j+1,&vj1));
45 :
46 780 : PetscCall(BVDotColumnBegin(eps->V,j+1,H+ldh*j));
47 780 : if (j>k) {
48 735 : PetscCall(BVDotColumnBegin(eps->V,j,lhh));
49 735 : PetscCall(BVGetColumn(eps->V,j,&vj));
50 735 : PetscCall(VecDotBegin(vj,vj,&dot));
51 735 : if (j>k+1) {
52 690 : PetscCall(BVNormVecBegin(eps->V,u,NORM_2,&norm2));
53 690 : PetscCall(BVGetColumn(eps->V,j-2,&vj2));
54 690 : PetscCall(VecDotBegin(u,vj2,&dot2));
55 : }
56 735 : PetscCall(BVDotColumnEnd(eps->V,j+1,H+ldh*j));
57 735 : PetscCall(BVDotColumnEnd(eps->V,j,lhh));
58 735 : PetscCall(VecDotEnd(vj,vj,&dot));
59 735 : PetscCall(BVRestoreColumn(eps->V,j,&vj));
60 735 : if (j>k+1) {
61 690 : PetscCall(BVNormVecEnd(eps->V,u,NORM_2,&norm2));
62 690 : PetscCall(VecDotEnd(u,vj2,&dot2));
63 690 : PetscCall(BVRestoreColumn(eps->V,j-2,&vj2));
64 : }
65 735 : norm1 = PetscSqrtReal(PetscRealPart(dot));
66 8016 : for (i=0;i<j;i++) H[ldh*j+i] = H[ldh*j+i]/norm1;
67 735 : H[ldh*j+j] = H[ldh*j+j]/dot;
68 735 : PetscCall(BVCopyVec(eps->V,j,t));
69 735 : PetscCall(BVScaleColumn(eps->V,j,1.0/norm1));
70 735 : PetscCall(BVScaleColumn(eps->V,j+1,1.0/norm1));
71 45 : } else PetscCall(BVDotColumnEnd(eps->V,j+1,H+ldh*j)); /* j==k */
72 :
73 780 : PetscCall(BVMultColumn(eps->V,-1.0,1.0,j+1,H+ldh*j));
74 :
75 780 : if (j>k) {
76 735 : PetscCall(BVSetActiveColumns(eps->V,0,j));
77 735 : PetscCall(BVMultVec(eps->V,-1.0,1.0,t,lhh));
78 735 : PetscCall(BVSetActiveColumns(eps->V,0,m));
79 8016 : for (i=0;i<j;i++) H[ldh*(j-1)+i] += lhh[i];
80 : }
81 :
82 780 : if (j>k+1) {
83 690 : PetscCall(BVGetColumn(eps->V,j-1,&vj1));
84 690 : PetscCall(VecCopy(u,vj1));
85 690 : PetscCall(BVRestoreColumn(eps->V,j-1,&vj1));
86 690 : PetscCall(BVScaleColumn(eps->V,j-1,1.0/norm2));
87 690 : H[ldh*(j-2)+j-1] = norm2;
88 : }
89 :
90 780 : if (j<m-1) PetscCall(VecCopy(t,u));
91 : }
92 :
93 45 : PetscCall(BVNormVec(eps->V,t,NORM_2,&norm2));
94 45 : PetscCall(VecScale(t,1.0/norm2));
95 45 : PetscCall(BVGetColumn(eps->V,m-1,&vj1));
96 45 : PetscCall(VecCopy(t,vj1));
97 45 : PetscCall(BVRestoreColumn(eps->V,m-1,&vj1));
98 45 : H[ldh*(m-2)+m-1] = norm2;
99 :
100 45 : PetscCall(BVDotColumn(eps->V,m,lhh));
101 :
102 45 : PetscCall(BVMultColumn(eps->V,-1.0,1.0,m,lhh));
103 882 : for (i=0;i<m;i++)
104 837 : H[ldh*(m-1)+i] += lhh[i];
105 :
106 45 : PetscCall(BVNormColumn(eps->V,m,NORM_2,beta));
107 45 : PetscCall(BVScaleColumn(eps->V,m,1.0 / *beta));
108 45 : *breakdown = PETSC_FALSE;
109 :
110 45 : if (m>100) PetscCall(PetscFree(lhh));
111 45 : PetscCall(VecDestroy(&u));
112 45 : PetscCall(VecDestroy(&t));
113 45 : 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 15 : static PetscErrorCode EPSKrylovConvergence_Filter(EPS eps,PetscBool getall,PetscInt kini,PetscInt nits,PetscReal beta,PetscReal gamma,PetscInt *kout)
161 : {
162 15 : PetscInt k,ninside,nconv;
163 15 : PetscScalar re,im;
164 15 : PetscReal resnorm;
165 :
166 15 : PetscFunctionBegin;
167 15 : ninside = 0; /* count how many eigenvalues are located in the interval */
168 117 : for (k=kini;k<kini+nits;k++) {
169 117 : if (PetscRealPart(eps->eigr[k]) < gamma) break;
170 102 : ninside++;
171 : }
172 15 : eps->nev = ninside+kini; /* adjust eigenvalue count */
173 15 : nconv = 0; /* count how many eigenvalues satisfy the convergence criterion */
174 78 : for (k=kini;k<kini+ninside;k++) {
175 : /* eigenvalue */
176 70 : re = eps->eigr[k];
177 70 : im = eps->eigi[k];
178 70 : PetscCall(DSVectors(eps->ds,DS_MAT_X,&k,&resnorm));
179 70 : resnorm *= beta;
180 : /* error estimate */
181 70 : PetscCall((*eps->converged)(eps,re,im,resnorm,&eps->errest[k],eps->convergedctx));
182 70 : if (eps->errest[k] < eps->tol) nconv++;
183 : else break;
184 : }
185 15 : *kout = kini+nconv;
186 15 : 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 15 : 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 7109 : PetscErrorCode EPSKrylovConvergence(EPS eps,PetscBool getall,PetscInt kini,PetscInt nits,PetscReal beta,PetscReal betat,PetscReal corrf,PetscInt *kout)
208 : {
209 7109 : PetscInt k,newk,newk2,marker,ld,inside;
210 7109 : PetscScalar re,im,*Zr,*Zi,*X;
211 7109 : PetscReal resnorm,gamma,lerrest;
212 7109 : PetscBool isshift,isfilter,refined,istrivial;
213 7109 : Vec x=NULL,y=NULL,w[3];
214 :
215 7109 : PetscFunctionBegin;
216 7109 : if (PetscUnlikely(eps->which == EPS_ALL)) {
217 77 : PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilter));
218 77 : if (isfilter) {
219 15 : PetscCall(STFilterGetThreshold(eps->st,&gamma));
220 15 : PetscCall(EPSKrylovConvergence_Filter(eps,getall,kini,nits,beta,gamma,kout));
221 15 : PetscFunctionReturn(PETSC_SUCCESS);
222 : }
223 : }
224 7094 : PetscCall(RGIsTrivial(eps->rg,&istrivial));
225 7094 : if (PetscUnlikely(eps->trueres)) {
226 48 : PetscCall(BVCreateVec(eps->V,&x));
227 48 : PetscCall(BVCreateVec(eps->V,&y));
228 48 : PetscCall(BVCreateVec(eps->V,&w[0]));
229 48 : PetscCall(BVCreateVec(eps->V,&w[2]));
230 : #if !defined(PETSC_USE_COMPLEX)
231 48 : PetscCall(BVCreateVec(eps->V,&w[1]));
232 : #else
233 : w[1] = NULL;
234 : #endif
235 : }
236 7094 : PetscCall(DSGetLeadingDimension(eps->ds,&ld));
237 7094 : PetscCall(DSGetRefined(eps->ds,&refined));
238 7094 : PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STSHIFT,&isshift));
239 7094 : marker = -1;
240 7094 : if (eps->trackall) getall = PETSC_TRUE;
241 15771 : for (k=kini;k<kini+nits;k++) {
242 : /* eigenvalue */
243 15656 : re = eps->eigr[k];
244 15656 : im = eps->eigi[k];
245 15656 : if (!istrivial || eps->trueres || isshift || eps->conv==EPS_CONV_NORM) PetscCall(STBackTransform(eps->st,1,&re,&im));
246 15656 : if (PetscUnlikely(!istrivial)) {
247 109 : PetscCall(RGCheckInside(eps->rg,1,&re,&im,&inside));
248 109 : if (marker==-1 && inside<0) marker = k;
249 109 : if (!(eps->trueres || isshift || eps->conv==EPS_CONV_NORM)) { /* make sure eps->converged below uses the right value */
250 42 : re = eps->eigr[k];
251 42 : im = eps->eigi[k];
252 : }
253 : }
254 15656 : newk = k;
255 15656 : PetscCall(DSVectors(eps->ds,DS_MAT_X,&newk,&resnorm));
256 15656 : if (PetscUnlikely(eps->trueres)) {
257 65 : PetscCall(DSGetArray(eps->ds,DS_MAT_X,&X));
258 65 : Zr = X+k*ld;
259 65 : if (newk==k+1) Zi = X+newk*ld;
260 : else Zi = NULL;
261 65 : PetscCall(EPSComputeRitzVector(eps,Zr,Zi,eps->V,x,y));
262 65 : PetscCall(DSRestoreArray(eps->ds,DS_MAT_X,&X));
263 65 : PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_FALSE,re,im,x,y,w,&resnorm));
264 : }
265 15591 : else if (!refined) resnorm *= beta*corrf;
266 : /* error estimate */
267 15656 : PetscCall((*eps->converged)(eps,re,im,resnorm,&eps->errest[k],eps->convergedctx));
268 15656 : if (marker==-1 && eps->errest[k] >= eps->tol) marker = k;
269 15656 : if (PetscUnlikely(eps->twosided)) {
270 150 : newk2 = k;
271 150 : PetscCall(DSVectors(eps->ds,DS_MAT_Y,&newk2,&resnorm));
272 150 : resnorm *= betat;
273 150 : PetscCall((*eps->converged)(eps,re,im,resnorm,&lerrest,eps->convergedctx));
274 150 : eps->errest[k] = PetscMax(eps->errest[k],lerrest);
275 150 : if (marker==-1 && lerrest >= eps->tol) marker = k;
276 : }
277 15656 : if (PetscUnlikely(newk==k+1)) {
278 402 : eps->errest[k+1] = eps->errest[k];
279 402 : k++;
280 : }
281 15656 : if (marker!=-1 && !getall) break;
282 : }
283 7094 : if (marker!=-1) k = marker;
284 7094 : *kout = k;
285 7094 : if (PetscUnlikely(eps->trueres)) {
286 48 : PetscCall(VecDestroy(&x));
287 48 : PetscCall(VecDestroy(&y));
288 48 : PetscCall(VecDestroy(&w[0]));
289 48 : PetscCall(VecDestroy(&w[2]));
290 : #if !defined(PETSC_USE_COMPLEX)
291 48 : PetscCall(VecDestroy(&w[1]));
292 : #endif
293 : }
294 7094 : PetscFunctionReturn(PETSC_SUCCESS);
295 : }
296 :
297 2440 : PetscErrorCode EPSPseudoLanczos(EPS eps,PetscReal *alpha,PetscReal *beta,PetscReal *omega,PetscInt k,PetscInt *M,PetscBool *breakdown,PetscBool *symmlost,PetscReal *cos,Vec w)
298 : {
299 2440 : PetscInt j,m = *M,i,ld,l;
300 2440 : Vec vj,vj1;
301 2440 : PetscScalar *hwork,lhwork[100];
302 2440 : PetscReal norm,norm1,norm2,t,sym=0.0,fro=0.0;
303 2440 : PetscBLASInt j_,one=1;
304 :
305 2440 : PetscFunctionBegin;
306 2440 : PetscCall(DSGetLeadingDimension(eps->ds,&ld));
307 2440 : PetscCall(DSGetDimensions(eps->ds,NULL,&l,NULL,NULL));
308 2440 : if (cos) *cos = 1.0;
309 2440 : if (m > 100) PetscCall(PetscMalloc1(m,&hwork));
310 2434 : else hwork = lhwork;
311 :
312 2440 : PetscCall(BVSetActiveColumns(eps->V,0,m));
313 42624 : for (j=k;j<m;j++) {
314 40184 : PetscCall(BVGetColumn(eps->V,j,&vj));
315 40184 : PetscCall(BVGetColumn(eps->V,j+1,&vj1));
316 40184 : PetscCall(STApply(eps->st,vj,vj1));
317 40184 : PetscCall(BVRestoreColumn(eps->V,j,&vj));
318 40184 : PetscCall(BVRestoreColumn(eps->V,j+1,&vj1));
319 40184 : PetscCall(BVOrthogonalizeColumn(eps->V,j+1,hwork,&norm,breakdown));
320 40184 : alpha[j] = PetscRealPart(hwork[j]);
321 40184 : beta[j] = PetscAbsReal(norm);
322 40184 : if (j==k) {
323 2440 : PetscReal *f;
324 :
325 2440 : PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&f));
326 3490 : for (i=0;i<l;i++) hwork[i] = 0.0;
327 38284 : for (;i<j-1;i++) hwork[i] -= f[2*ld+i];
328 2440 : PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&f));
329 : }
330 40184 : if (j>0) {
331 40168 : hwork[j-1] -= beta[j-1];
332 40168 : PetscCall(PetscBLASIntCast(j,&j_));
333 40168 : sym = SlepcAbs(BLASnrm2_(&j_,hwork,&one),sym);
334 : }
335 40184 : fro = SlepcAbs(fro,SlepcAbs(alpha[j],beta[j]));
336 40184 : if (j>0) fro = SlepcAbs(fro,beta[j-1]);
337 40259 : if (sym/fro>PetscMax(PETSC_SQRT_MACHINE_EPSILON,10*eps->tol)) { *symmlost = PETSC_TRUE; *M=j; break; }
338 40184 : omega[j+1] = (norm<0.0)? -1.0: 1.0;
339 40184 : PetscCall(BVScaleColumn(eps->V,j+1,1.0/norm));
340 : /* */
341 40184 : 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 2440 : if (m > 100) PetscCall(PetscFree(hwork));
352 2440 : PetscFunctionReturn(PETSC_SUCCESS);
353 : }
|
