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 : SLEPc eigensolver: "lanczos"
12 :
13 : Method: Explicitly Restarted Symmetric/Hermitian Lanczos
14 :
15 : Algorithm:
16 :
17 : Lanczos method for symmetric (Hermitian) problems, with explicit
18 : restart and deflation. Several reorthogonalization strategies can
19 : be selected.
20 :
21 : References:
22 :
23 : [1] "Lanczos Methods in SLEPc", SLEPc Technical Report STR-5,
24 : available at https://slepc.upv.es.
25 : */
26 :
27 : #include <slepc/private/epsimpl.h> /*I "slepceps.h" I*/
28 : #include <slepcblaslapack.h>
29 :
30 : typedef struct {
31 : EPSLanczosReorthogType reorthog; /* user-provided reorthogonalization parameter */
32 : PetscInt allocsize; /* number of columns of work BV's allocated at setup */
33 : BV AV; /* work BV used in selective reorthogonalization */
34 : } EPS_LANCZOS;
35 :
36 36 : static PetscErrorCode EPSSetUp_Lanczos(EPS eps)
37 : {
38 36 : EPS_LANCZOS *lanczos = (EPS_LANCZOS*)eps->data;
39 36 : BVOrthogRefineType refine;
40 36 : BVOrthogBlockType btype;
41 36 : PetscReal eta;
42 :
43 36 : PetscFunctionBegin;
44 36 : EPSCheckHermitianDefinite(eps);
45 36 : PetscCall(EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd));
46 36 : PetscCheck(eps->ncv<=eps->nev+eps->mpd,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must not be larger than nev+mpd");
47 36 : if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
48 36 : if (!eps->which) PetscCall(EPSSetWhichEigenpairs_Default(eps));
49 36 : PetscCheck(eps->which!=EPS_ALL,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support computing all eigenvalues");
50 36 : EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION);
51 36 : EPSCheckIgnored(eps,EPS_FEATURE_BALANCE);
52 :
53 36 : PetscCheck(lanczos->reorthog!=(EPSLanczosReorthogType)-1,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"You should explicitly provide the reorthogonalization type, e.g., -eps_lanczos_reorthog local. Note that the EPSLANCZOS solver is *NOT RECOMMENDED* for general use, because it uses explicit restart which typically has slow convergence. The recommended solver is EPSKRYLOVSCHUR (the default), which implements Lanczos with thick restart in the case of symmetric/Hermitian problems");
54 :
55 36 : PetscCall(EPSAllocateSolution(eps,1));
56 36 : PetscCall(EPS_SetInnerProduct(eps));
57 36 : if (lanczos->reorthog != EPS_LANCZOS_REORTHOG_FULL) {
58 29 : PetscCall(BVGetOrthogonalization(eps->V,NULL,&refine,&eta,&btype));
59 29 : PetscCall(BVSetOrthogonalization(eps->V,BV_ORTHOG_MGS,refine,eta,btype));
60 29 : PetscCall(PetscInfo(eps,"Switching to MGS orthogonalization\n"));
61 : }
62 36 : if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_SELECTIVE) {
63 3 : if (!lanczos->allocsize) {
64 1 : PetscCall(BVDuplicate(eps->V,&lanczos->AV));
65 1 : PetscCall(BVGetSizes(lanczos->AV,NULL,NULL,&lanczos->allocsize));
66 : } else { /* make sure V and AV have the same size */
67 2 : PetscCall(BVGetSizes(eps->V,NULL,NULL,&lanczos->allocsize));
68 2 : PetscCall(BVResize(lanczos->AV,lanczos->allocsize,PETSC_FALSE));
69 : }
70 : }
71 :
72 36 : PetscCall(DSSetType(eps->ds,DSHEP));
73 36 : PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
74 36 : PetscCall(DSAllocate(eps->ds,eps->ncv+1));
75 36 : if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_LOCAL) PetscCall(EPSSetWorkVecs(eps,1));
76 36 : PetscFunctionReturn(PETSC_SUCCESS);
77 : }
78 :
79 : /*
80 : EPSLocalLanczos - Local reorthogonalization.
81 :
82 : This is the simplest variant. At each Lanczos step, the corresponding Lanczos vector
83 : is orthogonalized with respect to the two previous Lanczos vectors, according to
84 : the three term Lanczos recurrence. WARNING: This variant does not track the loss of
85 : orthogonality that occurs in finite-precision arithmetic and, therefore, the
86 : generated vectors are not guaranteed to be (semi-)orthogonal.
87 : */
88 520 : static PetscErrorCode EPSLocalLanczos(EPS eps,PetscReal *alpha,PetscReal *beta,PetscInt k,PetscInt *M,PetscBool *breakdown)
89 : {
90 520 : PetscInt i,j,m = *M;
91 520 : Mat Op;
92 520 : PetscBool *which,lwhich[100];
93 520 : PetscScalar *hwork,lhwork[100];
94 :
95 520 : PetscFunctionBegin;
96 520 : if (m > 100) PetscCall(PetscMalloc2(m,&which,m,&hwork));
97 : else {
98 520 : which = lwhich;
99 520 : hwork = lhwork;
100 : }
101 1857 : for (i=0;i<k;i++) which[i] = PETSC_TRUE;
102 :
103 520 : PetscCall(BVSetActiveColumns(eps->V,0,m));
104 520 : PetscCall(STGetOperator(eps->st,&Op));
105 7578 : for (j=k;j<m;j++) {
106 7058 : PetscCall(BVMatMultColumn(eps->V,Op,j));
107 7058 : which[j] = PETSC_TRUE;
108 7058 : if (j-2>=k) which[j-2] = PETSC_FALSE;
109 7058 : PetscCall(BVOrthogonalizeSomeColumn(eps->V,j+1,which,hwork,beta+j,breakdown));
110 7058 : alpha[j] = PetscRealPart(hwork[j]);
111 7058 : if (PetscUnlikely(*breakdown)) {
112 0 : *M = j+1;
113 0 : break;
114 7058 : } else PetscCall(BVScaleColumn(eps->V,j+1,1/beta[j]));
115 : }
116 520 : PetscCall(STRestoreOperator(eps->st,&Op));
117 520 : if (m > 100) PetscCall(PetscFree2(which,hwork));
118 520 : PetscFunctionReturn(PETSC_SUCCESS);
119 : }
120 :
121 : /*
122 : DenseTridiagonal - Solves a real tridiagonal Hermitian Eigenvalue Problem.
123 :
124 : Input Parameters:
125 : + n - dimension of the eigenproblem
126 : . D - pointer to the array containing the diagonal elements
127 : - E - pointer to the array containing the off-diagonal elements
128 :
129 : Output Parameters:
130 : + w - pointer to the array to store the computed eigenvalues
131 : - V - pointer to the array to store the eigenvectors
132 :
133 : Notes:
134 : If V is NULL then the eigenvectors are not computed.
135 :
136 : This routine use LAPACK routines xSTEVR.
137 : */
138 673 : static PetscErrorCode DenseTridiagonal(PetscInt n_,PetscReal *D,PetscReal *E,PetscReal *w,PetscScalar *V)
139 : {
140 673 : PetscReal abstol = 0.0,vl,vu,*work;
141 673 : PetscBLASInt il,iu,m,*isuppz,n,lwork,*iwork,liwork,info;
142 673 : const char *jobz;
143 : #if defined(PETSC_USE_COMPLEX)
144 : PetscInt i,j;
145 : PetscReal *VV=NULL;
146 : #endif
147 :
148 673 : PetscFunctionBegin;
149 673 : PetscCall(PetscBLASIntCast(n_,&n));
150 673 : PetscCall(PetscBLASIntCast(20*n_,&lwork));
151 673 : PetscCall(PetscBLASIntCast(10*n_,&liwork));
152 673 : if (V) {
153 : jobz = "V";
154 : #if defined(PETSC_USE_COMPLEX)
155 : PetscCall(PetscMalloc1(n*n,&VV));
156 : #endif
157 0 : } else jobz = "N";
158 673 : PetscCall(PetscMalloc3(2*n,&isuppz,lwork,&work,liwork,&iwork));
159 673 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
160 : #if defined(PETSC_USE_COMPLEX)
161 : PetscCallBLAS("LAPACKstevr",LAPACKstevr_(jobz,"A",&n,D,E,&vl,&vu,&il,&iu,&abstol,&m,w,VV,&n,isuppz,work,&lwork,iwork,&liwork,&info));
162 : #else
163 673 : PetscCallBLAS("LAPACKstevr",LAPACKstevr_(jobz,"A",&n,D,E,&vl,&vu,&il,&iu,&abstol,&m,w,V,&n,isuppz,work,&lwork,iwork,&liwork,&info));
164 : #endif
165 673 : PetscCall(PetscFPTrapPop());
166 673 : SlepcCheckLapackInfo("stevr",info);
167 : #if defined(PETSC_USE_COMPLEX)
168 : if (V) {
169 : for (i=0;i<n;i++)
170 : for (j=0;j<n;j++)
171 : V[i*n+j] = VV[i*n+j];
172 : PetscCall(PetscFree(VV));
173 : }
174 : #endif
175 673 : PetscCall(PetscFree3(isuppz,work,iwork));
176 673 : PetscFunctionReturn(PETSC_SUCCESS);
177 : }
178 :
179 : /*
180 : EPSSelectiveLanczos - Selective reorthogonalization.
181 : */
182 38 : static PetscErrorCode EPSSelectiveLanczos(EPS eps,PetscReal *alpha,PetscReal *beta,PetscInt k,PetscInt *M,PetscBool *breakdown,PetscReal anorm)
183 : {
184 38 : EPS_LANCZOS *lanczos = (EPS_LANCZOS*)eps->data;
185 38 : PetscInt i,j,m = *M,n,nritz=0,nritzo;
186 38 : Vec vj1,av;
187 38 : Mat Op;
188 38 : PetscReal *d,*e,*ritz,norm;
189 38 : PetscScalar *Y,*hwork;
190 38 : PetscBool *which;
191 :
192 38 : PetscFunctionBegin;
193 38 : PetscCall(PetscCalloc6(m+1,&d,m,&e,m,&ritz,m*m,&Y,m,&which,m,&hwork));
194 87 : for (i=0;i<k;i++) which[i] = PETSC_TRUE;
195 38 : PetscCall(STGetOperator(eps->st,&Op));
196 :
197 711 : for (j=k;j<m;j++) {
198 673 : PetscCall(BVSetActiveColumns(eps->V,0,m));
199 :
200 : /* Lanczos step */
201 673 : PetscCall(BVMatMultColumn(eps->V,Op,j));
202 673 : which[j] = PETSC_TRUE;
203 673 : if (j-2>=k) which[j-2] = PETSC_FALSE;
204 673 : PetscCall(BVOrthogonalizeSomeColumn(eps->V,j+1,which,hwork,&norm,breakdown));
205 673 : alpha[j] = PetscRealPart(hwork[j]);
206 673 : beta[j] = norm;
207 673 : if (PetscUnlikely(*breakdown)) {
208 0 : *M = j+1;
209 0 : break;
210 : }
211 :
212 : /* Compute eigenvalues and eigenvectors Y of the tridiagonal block */
213 673 : n = j-k+1;
214 6994 : for (i=0;i<n;i++) {
215 6321 : d[i] = alpha[i+k];
216 6321 : e[i] = beta[i+k];
217 : }
218 673 : PetscCall(DenseTridiagonal(n,d,e,ritz,Y));
219 :
220 : /* Estimate ||A|| */
221 6994 : for (i=0;i<n;i++)
222 6321 : if (PetscAbsReal(ritz[i]) > anorm) anorm = PetscAbsReal(ritz[i]);
223 :
224 : /* Compute nearly converged Ritz vectors */
225 : nritzo = 0;
226 6994 : for (i=0;i<n;i++) {
227 6321 : if (norm*PetscAbsScalar(Y[i*n+n-1]) < PETSC_SQRT_MACHINE_EPSILON*anorm) nritzo++;
228 : }
229 673 : if (nritzo>nritz) {
230 : nritz = 0;
231 196 : for (i=0;i<n;i++) {
232 173 : if (norm*PetscAbsScalar(Y[i*n+n-1]) < PETSC_SQRT_MACHINE_EPSILON*anorm) {
233 37 : PetscCall(BVSetActiveColumns(eps->V,k,k+n));
234 37 : PetscCall(BVGetColumn(lanczos->AV,nritz,&av));
235 37 : PetscCall(BVMultVec(eps->V,1.0,0.0,av,Y+i*n));
236 37 : PetscCall(BVRestoreColumn(lanczos->AV,nritz,&av));
237 37 : nritz++;
238 : }
239 : }
240 : }
241 673 : if (nritz > 0) {
242 247 : PetscCall(BVGetColumn(eps->V,j+1,&vj1));
243 247 : PetscCall(BVSetActiveColumns(lanczos->AV,0,nritz));
244 247 : PetscCall(BVOrthogonalizeVec(lanczos->AV,vj1,hwork,&norm,breakdown));
245 247 : PetscCall(BVRestoreColumn(eps->V,j+1,&vj1));
246 247 : if (PetscUnlikely(*breakdown)) {
247 0 : *M = j+1;
248 0 : break;
249 : }
250 : }
251 673 : PetscCall(BVScaleColumn(eps->V,j+1,1.0/norm));
252 : }
253 :
254 38 : PetscCall(STRestoreOperator(eps->st,&Op));
255 38 : PetscCall(PetscFree6(d,e,ritz,Y,which,hwork));
256 38 : PetscFunctionReturn(PETSC_SUCCESS);
257 : }
258 :
259 768 : static void update_omega(PetscReal *omega,PetscReal *omega_old,PetscInt j,PetscReal *alpha,PetscReal *beta,PetscReal eps1,PetscReal anorm)
260 : {
261 768 : PetscInt k;
262 768 : PetscReal T,binv;
263 :
264 768 : PetscFunctionBegin;
265 : /* Estimate of contribution to roundoff errors from A*v
266 : fl(A*v) = A*v + f,
267 : where ||f|| \approx eps1*||A||.
268 : For a full matrix A, a rule-of-thumb estimate is eps1 = sqrt(n)*eps */
269 768 : T = eps1*anorm;
270 768 : binv = 1.0/beta[j+1];
271 :
272 : /* Update omega(1) using omega(0)==0 */
273 768 : omega_old[0]= beta[1]*omega[1] + (alpha[0]-alpha[j])*omega[0] - beta[j]*omega_old[0];
274 768 : if (omega_old[0] > 0) omega_old[0] = binv*(omega_old[0] + T);
275 461 : else omega_old[0] = binv*(omega_old[0] - T);
276 :
277 : /* Update remaining components */
278 7014 : for (k=1;k<j-1;k++) {
279 6246 : omega_old[k] = beta[k+1]*omega[k+1] + (alpha[k]-alpha[j])*omega[k] + beta[k]*omega[k-1] - beta[j]*omega_old[k];
280 6246 : if (omega_old[k] > 0) omega_old[k] = binv*(omega_old[k] + T);
281 3192 : else omega_old[k] = binv*(omega_old[k] - T);
282 : }
283 768 : omega_old[j-1] = binv*T;
284 :
285 : /* Swap omega and omega_old */
286 8534 : for (k=0;k<j;k++) {
287 7766 : omega[k] = omega_old[k];
288 7766 : omega_old[k] = omega[k];
289 : }
290 768 : omega[j] = eps1;
291 768 : PetscFunctionReturnVoid();
292 : }
293 :
294 48 : static void compute_int(PetscBool *which,PetscReal *mu,PetscInt j,PetscReal delta,PetscReal eta)
295 : {
296 48 : PetscInt i,k,maxpos;
297 48 : PetscReal max;
298 48 : PetscBool found;
299 :
300 48 : PetscFunctionBegin;
301 : /* initialize which */
302 48 : found = PETSC_FALSE;
303 48 : maxpos = 0;
304 48 : max = 0.0;
305 755 : for (i=0;i<j;i++) {
306 707 : if (PetscAbsReal(mu[i]) >= delta) {
307 175 : which[i] = PETSC_TRUE;
308 175 : found = PETSC_TRUE;
309 532 : } else which[i] = PETSC_FALSE;
310 707 : if (PetscAbsReal(mu[i]) > max) {
311 78 : maxpos = i;
312 78 : max = PetscAbsReal(mu[i]);
313 : }
314 : }
315 48 : if (!found) which[maxpos] = PETSC_TRUE;
316 :
317 755 : for (i=0;i<j;i++) {
318 707 : if (which[i]) {
319 : /* find left interval */
320 623 : for (k=i;k>=0;k--) {
321 623 : if (PetscAbsReal(mu[k])<eta || which[k]) break;
322 0 : else which[k] = PETSC_TRUE;
323 : }
324 : /* find right interval */
325 1071 : for (k=i+1;k<j;k++) {
326 1037 : if (PetscAbsReal(mu[k])<eta || which[k]) break;
327 448 : else which[k] = PETSC_TRUE;
328 : }
329 : }
330 : }
331 48 : PetscFunctionReturnVoid();
332 : }
333 :
334 : /*
335 : EPSPartialLanczos - Partial reorthogonalization.
336 : */
337 76 : static PetscErrorCode EPSPartialLanczos(EPS eps,PetscReal *alpha,PetscReal *beta,PetscInt k,PetscInt *M,PetscBool *breakdown,PetscReal anorm)
338 : {
339 76 : EPS_LANCZOS *lanczos = (EPS_LANCZOS*)eps->data;
340 76 : PetscInt i,j,m = *M;
341 76 : Mat Op;
342 76 : PetscReal norm,*omega,lomega[100],*omega_old,lomega_old[100],eps1,delta,eta;
343 76 : PetscBool *which,lwhich[100],*which2,lwhich2[100];
344 76 : PetscBool reorth = PETSC_FALSE,force_reorth = PETSC_FALSE;
345 76 : PetscBool fro = PETSC_FALSE,estimate_anorm = PETSC_FALSE;
346 76 : PetscScalar *hwork,lhwork[100];
347 :
348 76 : PetscFunctionBegin;
349 76 : if (m>100) PetscCall(PetscMalloc5(m,&omega,m,&omega_old,m,&which,m,&which2,m,&hwork));
350 : else {
351 76 : omega = lomega;
352 76 : omega_old = lomega_old;
353 76 : which = lwhich;
354 76 : which2 = lwhich2;
355 76 : hwork = lhwork;
356 : }
357 :
358 76 : eps1 = PetscSqrtReal((PetscReal)eps->n)*PETSC_MACHINE_EPSILON/2;
359 76 : delta = PETSC_SQRT_MACHINE_EPSILON/PetscSqrtReal((PetscReal)eps->ncv);
360 76 : eta = PetscPowReal(PETSC_MACHINE_EPSILON,3.0/4.0)/PetscSqrtReal((PetscReal)eps->ncv);
361 76 : if (anorm < 0.0) {
362 6 : anorm = 1.0;
363 6 : estimate_anorm = PETSC_TRUE;
364 : }
365 7676 : for (i=0;i<PetscMax(100,m);i++) omega[i] = omega_old[i] = 0.0;
366 174 : for (i=0;i<k;i++) which[i] = PETSC_TRUE;
367 :
368 76 : PetscCall(BVSetActiveColumns(eps->V,0,m));
369 76 : PetscCall(STGetOperator(eps->st,&Op));
370 1422 : for (j=k;j<m;j++) {
371 1346 : PetscCall(BVMatMultColumn(eps->V,Op,j));
372 1346 : if (fro) {
373 : /* Lanczos step with full reorthogonalization */
374 502 : PetscCall(BVOrthogonalizeColumn(eps->V,j+1,hwork,&norm,breakdown));
375 502 : alpha[j] = PetscRealPart(hwork[j]);
376 : } else {
377 : /* Lanczos step */
378 844 : which[j] = PETSC_TRUE;
379 844 : if (j-2>=k) which[j-2] = PETSC_FALSE;
380 844 : PetscCall(BVOrthogonalizeSomeColumn(eps->V,j+1,which,hwork,&norm,breakdown));
381 844 : alpha[j] = PetscRealPart(hwork[j]);
382 844 : beta[j] = norm;
383 :
384 : /* Estimate ||A|| if needed */
385 844 : if (estimate_anorm) {
386 114 : if (j>k) anorm = PetscMax(anorm,PetscAbsReal(alpha[j])+norm+beta[j-1]);
387 6 : else anorm = PetscMax(anorm,PetscAbsReal(alpha[j])+norm);
388 : }
389 :
390 : /* Check if reorthogonalization is needed */
391 844 : reorth = PETSC_FALSE;
392 844 : if (j>k) {
393 768 : update_omega(omega,omega_old,j,alpha,beta-1,eps1,anorm);
394 8362 : for (i=0;i<j-k;i++) {
395 6826 : if (PetscAbsReal(omega[i]) > delta) reorth = PETSC_TRUE;
396 : }
397 : }
398 844 : if (reorth || force_reorth) {
399 290 : for (i=0;i<k;i++) which2[i] = PETSC_FALSE;
400 2653 : for (i=k;i<=j;i++) which2[i] = PETSC_TRUE;
401 159 : if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_PERIODIC) {
402 : /* Periodic reorthogonalization */
403 79 : if (force_reorth) force_reorth = PETSC_FALSE;
404 47 : else force_reorth = PETSC_TRUE;
405 1238 : for (i=0;i<j-k;i++) omega[i] = eps1;
406 : } else {
407 : /* Partial reorthogonalization */
408 80 : if (force_reorth) force_reorth = PETSC_FALSE;
409 : else {
410 48 : force_reorth = PETSC_TRUE;
411 48 : compute_int(which2+k,omega,j-k,delta,eta);
412 803 : for (i=0;i<j-k;i++) {
413 707 : if (which2[i+k]) omega[i] = eps1;
414 : }
415 : }
416 : }
417 159 : PetscCall(BVOrthogonalizeSomeColumn(eps->V,j+1,which2,hwork,&norm,breakdown));
418 : }
419 : }
420 :
421 1346 : if (PetscUnlikely(*breakdown || norm < eps->n*anorm*PETSC_MACHINE_EPSILON)) {
422 0 : *M = j+1;
423 0 : break;
424 : }
425 1346 : if (!fro && norm*delta < anorm*eps1) {
426 30 : fro = PETSC_TRUE;
427 30 : PetscCall(PetscInfo(eps,"Switching to full reorthogonalization at iteration %" PetscInt_FMT "\n",eps->its));
428 : }
429 1346 : beta[j] = norm;
430 1346 : PetscCall(BVScaleColumn(eps->V,j+1,1.0/norm));
431 : }
432 :
433 76 : PetscCall(STRestoreOperator(eps->st,&Op));
434 76 : if (m>100) PetscCall(PetscFree5(omega,omega_old,which,which2,hwork));
435 76 : PetscFunctionReturn(PETSC_SUCCESS);
436 : }
437 :
438 : /*
439 : EPSBasicLanczos - Computes an m-step Lanczos factorization. The first k
440 : columns are assumed to be locked and therefore they are not modified. On
441 : exit, the following relation is satisfied:
442 :
443 : OP * V - V * T = f * e_m^T
444 :
445 : where the columns of V are the Lanczos vectors, T is a tridiagonal matrix,
446 : f is the residual vector and e_m is the m-th vector of the canonical basis.
447 : The Lanczos vectors (together with vector f) are B-orthogonal (to working
448 : accuracy) if full reorthogonalization is being used, otherwise they are
449 : (B-)semi-orthogonal. On exit, beta contains the B-norm of f and the next
450 : Lanczos vector can be computed as v_{m+1} = f / beta.
451 :
452 : This function simply calls another function which depends on the selected
453 : reorthogonalization strategy.
454 : */
455 932 : static PetscErrorCode EPSBasicLanczos(EPS eps,PetscInt k,PetscInt *m,PetscReal *betam,PetscBool *breakdown,PetscReal anorm)
456 : {
457 932 : EPS_LANCZOS *lanczos = (EPS_LANCZOS*)eps->data;
458 932 : PetscScalar *T;
459 932 : PetscInt i,n=*m,ld;
460 932 : PetscReal *alpha,*beta;
461 932 : BVOrthogRefineType orthog_ref;
462 932 : Mat Op,M;
463 :
464 932 : PetscFunctionBegin;
465 932 : PetscCall(DSGetLeadingDimension(eps->ds,&ld));
466 932 : switch (lanczos->reorthog) {
467 520 : case EPS_LANCZOS_REORTHOG_LOCAL:
468 520 : PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&alpha));
469 520 : beta = alpha + ld;
470 520 : PetscCall(EPSLocalLanczos(eps,alpha,beta,k,m,breakdown));
471 520 : *betam = beta[*m-1];
472 520 : PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&alpha));
473 : break;
474 271 : case EPS_LANCZOS_REORTHOG_FULL:
475 271 : PetscCall(STGetOperator(eps->st,&Op));
476 271 : PetscCall(DSGetMat(eps->ds,DS_MAT_T,&M));
477 271 : PetscCall(BVMatLanczos(eps->V,Op,M,k,m,betam,breakdown));
478 271 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_T,&M));
479 271 : PetscCall(STRestoreOperator(eps->st,&Op));
480 : break;
481 38 : case EPS_LANCZOS_REORTHOG_SELECTIVE:
482 38 : PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&alpha));
483 38 : beta = alpha + ld;
484 38 : PetscCall(EPSSelectiveLanczos(eps,alpha,beta,k,m,breakdown,anorm));
485 38 : *betam = beta[*m-1];
486 38 : PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&alpha));
487 : break;
488 76 : case EPS_LANCZOS_REORTHOG_PERIODIC:
489 : case EPS_LANCZOS_REORTHOG_PARTIAL:
490 76 : PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&alpha));
491 76 : beta = alpha + ld;
492 76 : PetscCall(EPSPartialLanczos(eps,alpha,beta,k,m,breakdown,anorm));
493 76 : *betam = beta[*m-1];
494 76 : PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&alpha));
495 : break;
496 27 : case EPS_LANCZOS_REORTHOG_DELAYED:
497 27 : PetscCall(PetscMalloc1(n*n,&T));
498 27 : PetscCall(BVGetOrthogonalization(eps->V,NULL,&orthog_ref,NULL,NULL));
499 27 : if (orthog_ref == BV_ORTHOG_REFINE_NEVER) PetscCall(EPSDelayedArnoldi1(eps,T,n,k,m,betam,breakdown));
500 27 : else PetscCall(EPSDelayedArnoldi(eps,T,n,k,m,betam,breakdown));
501 27 : n = *m;
502 27 : PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&alpha));
503 27 : beta = alpha + ld;
504 472 : for (i=k;i<n-1;i++) {
505 445 : alpha[i] = PetscRealPart(T[n*i+i]);
506 445 : beta[i] = PetscRealPart(T[n*i+i+1]);
507 : }
508 27 : alpha[n-1] = PetscRealPart(T[n*(n-1)+n-1]);
509 27 : beta[n-1] = *betam;
510 27 : PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&alpha));
511 27 : PetscCall(PetscFree(T));
512 : break;
513 : }
514 932 : PetscFunctionReturn(PETSC_SUCCESS);
515 : }
516 :
517 36 : static PetscErrorCode EPSSolve_Lanczos(EPS eps)
518 : {
519 36 : EPS_LANCZOS *lanczos = (EPS_LANCZOS*)eps->data;
520 36 : PetscInt nconv,i,j,k,l,x,n,*perm,restart,ncv=eps->ncv,r,ld;
521 36 : Vec vi,vj,w;
522 36 : Mat U;
523 36 : PetscScalar *Y,*ritz,stmp;
524 36 : PetscReal *bnd,anorm,beta,norm,rtmp,resnorm;
525 36 : PetscBool breakdown;
526 36 : char *conv,ctmp;
527 :
528 36 : PetscFunctionBegin;
529 36 : PetscCall(DSGetLeadingDimension(eps->ds,&ld));
530 36 : PetscCall(PetscMalloc4(ncv,&ritz,ncv,&bnd,ncv,&perm,ncv,&conv));
531 :
532 : /* The first Lanczos vector is the normalized initial vector */
533 36 : PetscCall(EPSGetStartVector(eps,0,NULL));
534 :
535 : anorm = -1.0;
536 : nconv = 0;
537 :
538 : /* Restart loop */
539 968 : while (eps->reason == EPS_CONVERGED_ITERATING) {
540 932 : eps->its++;
541 :
542 : /* Compute an ncv-step Lanczos factorization */
543 932 : n = PetscMin(nconv+eps->mpd,ncv);
544 932 : PetscCall(DSSetDimensions(eps->ds,n,nconv,PETSC_DEFAULT));
545 932 : PetscCall(EPSBasicLanczos(eps,nconv,&n,&beta,&breakdown,anorm));
546 932 : PetscCall(DSSetDimensions(eps->ds,n,nconv,0));
547 932 : PetscCall(DSSetState(eps->ds,DS_STATE_INTERMEDIATE));
548 932 : PetscCall(BVSetActiveColumns(eps->V,nconv,n));
549 :
550 : /* Solve projected problem */
551 932 : PetscCall(DSSolve(eps->ds,ritz,NULL));
552 932 : PetscCall(DSSort(eps->ds,ritz,NULL,NULL,NULL,NULL));
553 932 : PetscCall(DSSynchronize(eps->ds,ritz,NULL));
554 :
555 : /* Estimate ||A|| */
556 13272 : for (i=nconv;i<n;i++)
557 12761 : anorm = PetscMax(anorm,PetscAbsReal(PetscRealPart(ritz[i])));
558 :
559 : /* Compute residual norm estimates as beta*abs(Y(m,:)) + eps*||A|| */
560 932 : PetscCall(DSGetArray(eps->ds,DS_MAT_Q,&Y));
561 13272 : for (i=nconv;i<n;i++) {
562 12340 : resnorm = beta*PetscAbsScalar(Y[n-1+i*ld]) + PETSC_MACHINE_EPSILON*anorm;
563 12340 : PetscCall((*eps->converged)(eps,ritz[i],eps->eigi[i],resnorm,&bnd[i],eps->convergedctx));
564 12340 : if (bnd[i]<eps->tol) conv[i] = 'C';
565 12171 : else conv[i] = 'N';
566 : }
567 932 : PetscCall(DSRestoreArray(eps->ds,DS_MAT_Q,&Y));
568 :
569 : /* purge repeated ritz values */
570 932 : if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_LOCAL) {
571 7058 : for (i=nconv+1;i<n;i++) {
572 6538 : if (conv[i] == 'C' && PetscAbsScalar((ritz[i]-ritz[i-1])/ritz[i]) < eps->tol) conv[i] = 'R';
573 : }
574 : }
575 :
576 : /* Compute restart vector */
577 932 : if (breakdown) PetscCall(PetscInfo(eps,"Breakdown in Lanczos method (it=%" PetscInt_FMT " norm=%g)\n",eps->its,(double)beta));
578 : else {
579 : restart = nconv;
580 1078 : while (restart<n && conv[restart] != 'N') restart++;
581 931 : if (restart >= n) {
582 0 : breakdown = PETSC_TRUE;
583 : } else {
584 12175 : for (i=restart+1;i<n;i++) {
585 11244 : if (conv[i] == 'N') {
586 11240 : PetscCall(SlepcSCCompare(eps->sc,ritz[restart],0.0,ritz[i],0.0,&r));
587 11240 : if (r>0) restart = i;
588 : }
589 : }
590 931 : PetscCall(DSGetArray(eps->ds,DS_MAT_Q,&Y));
591 931 : PetscCall(BVMultColumn(eps->V,1.0,0.0,n,Y+restart*ld+nconv));
592 931 : PetscCall(DSRestoreArray(eps->ds,DS_MAT_Q,&Y));
593 : }
594 : }
595 :
596 : /* Count and put converged eigenvalues first */
597 13272 : for (i=nconv;i<n;i++) perm[i] = i;
598 1098 : for (k=nconv;k<n;k++) {
599 1097 : if (conv[perm[k]] != 'C') {
600 934 : j = k + 1;
601 12181 : while (j<n && conv[perm[j]] != 'C') j++;
602 934 : if (j>=n) break;
603 3 : l = perm[k]; perm[k] = perm[j]; perm[j] = l;
604 : }
605 : }
606 :
607 : /* Sort eigenvectors according to permutation */
608 932 : PetscCall(DSGetArray(eps->ds,DS_MAT_Q,&Y));
609 1098 : for (i=nconv;i<k;i++) {
610 166 : x = perm[i];
611 166 : if (x != i) {
612 3 : j = i + 1;
613 7 : while (perm[j] != i) j++;
614 : /* swap eigenvalues i and j */
615 3 : stmp = ritz[x]; ritz[x] = ritz[i]; ritz[i] = stmp;
616 3 : rtmp = bnd[x]; bnd[x] = bnd[i]; bnd[i] = rtmp;
617 3 : ctmp = conv[x]; conv[x] = conv[i]; conv[i] = ctmp;
618 3 : perm[j] = x; perm[i] = i;
619 : /* swap eigenvectors i and j */
620 60 : for (l=0;l<n;l++) {
621 57 : stmp = Y[l+x*ld]; Y[l+x*ld] = Y[l+i*ld]; Y[l+i*ld] = stmp;
622 : }
623 : }
624 : }
625 932 : PetscCall(DSRestoreArray(eps->ds,DS_MAT_Q,&Y));
626 :
627 : /* compute converged eigenvectors */
628 932 : PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&U));
629 932 : PetscCall(BVMultInPlace(eps->V,U,nconv,k));
630 932 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&U));
631 :
632 : /* purge spurious ritz values */
633 932 : if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_LOCAL) {
634 595 : for (i=nconv;i<k;i++) {
635 75 : PetscCall(BVGetColumn(eps->V,i,&vi));
636 75 : PetscCall(VecNorm(vi,NORM_2,&norm));
637 75 : PetscCall(VecScale(vi,1.0/norm));
638 75 : w = eps->work[0];
639 75 : PetscCall(STApply(eps->st,vi,w));
640 75 : PetscCall(VecAXPY(w,-ritz[i],vi));
641 75 : PetscCall(BVRestoreColumn(eps->V,i,&vi));
642 75 : PetscCall(VecNorm(w,NORM_2,&norm));
643 75 : PetscCall((*eps->converged)(eps,ritz[i],eps->eigi[i],norm,&bnd[i],eps->convergedctx));
644 75 : if (bnd[i]>=eps->tol) conv[i] = 'S';
645 : }
646 594 : for (i=nconv;i<k;i++) {
647 75 : if (conv[i] != 'C') {
648 1 : j = i + 1;
649 1 : while (j<k && conv[j] != 'C') j++;
650 1 : if (j>=k) break;
651 : /* swap eigenvalues i and j */
652 0 : stmp = ritz[j]; ritz[j] = ritz[i]; ritz[i] = stmp;
653 0 : rtmp = bnd[j]; bnd[j] = bnd[i]; bnd[i] = rtmp;
654 0 : ctmp = conv[j]; conv[j] = conv[i]; conv[i] = ctmp;
655 : /* swap eigenvectors i and j */
656 0 : PetscCall(BVGetColumn(eps->V,i,&vi));
657 0 : PetscCall(BVGetColumn(eps->V,j,&vj));
658 0 : PetscCall(VecSwap(vi,vj));
659 0 : PetscCall(BVRestoreColumn(eps->V,i,&vi));
660 74 : PetscCall(BVRestoreColumn(eps->V,j,&vj));
661 : }
662 : }
663 : k = i;
664 : }
665 :
666 : /* store ritz values and estimated errors */
667 13272 : for (i=nconv;i<n;i++) {
668 12340 : eps->eigr[i] = ritz[i];
669 12340 : eps->errest[i] = bnd[i];
670 : }
671 932 : nconv = k;
672 932 : PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,n));
673 932 : PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,nconv,eps->nev,&eps->reason,eps->stoppingctx));
674 :
675 932 : if (eps->reason == EPS_CONVERGED_ITERATING) { /* copy restart vector */
676 896 : PetscCall(BVCopyColumn(eps->V,n,nconv));
677 896 : if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_LOCAL && !breakdown) {
678 : /* Reorthonormalize restart vector */
679 502 : PetscCall(BVOrthonormalizeColumn(eps->V,nconv,PETSC_FALSE,NULL,&breakdown));
680 : }
681 896 : if (breakdown) {
682 : /* Use random vector for restarting */
683 0 : PetscCall(PetscInfo(eps,"Using random vector for restart\n"));
684 0 : PetscCall(EPSGetStartVector(eps,nconv,&breakdown));
685 : }
686 896 : if (PetscUnlikely(breakdown)) { /* give up */
687 0 : eps->reason = EPS_DIVERGED_BREAKDOWN;
688 968 : PetscCall(PetscInfo(eps,"Unable to generate more start vectors\n"));
689 : }
690 : }
691 : }
692 36 : eps->nconv = nconv;
693 :
694 36 : PetscCall(PetscFree4(ritz,bnd,perm,conv));
695 36 : PetscFunctionReturn(PETSC_SUCCESS);
696 : }
697 :
698 18 : static PetscErrorCode EPSSetFromOptions_Lanczos(EPS eps,PetscOptionItems *PetscOptionsObject)
699 : {
700 18 : EPS_LANCZOS *lanczos = (EPS_LANCZOS*)eps->data;
701 18 : PetscBool flg;
702 18 : EPSLanczosReorthogType reorthog=EPS_LANCZOS_REORTHOG_LOCAL,curval;
703 :
704 18 : PetscFunctionBegin;
705 18 : PetscOptionsHeadBegin(PetscOptionsObject,"EPS Lanczos Options");
706 :
707 18 : curval = (lanczos->reorthog==(EPSLanczosReorthogType)-1)? EPS_LANCZOS_REORTHOG_LOCAL: lanczos->reorthog;
708 18 : PetscCall(PetscOptionsEnum("-eps_lanczos_reorthog","Lanczos reorthogonalization","EPSLanczosSetReorthog",EPSLanczosReorthogTypes,(PetscEnum)curval,(PetscEnum*)&reorthog,&flg));
709 18 : if (flg) PetscCall(EPSLanczosSetReorthog(eps,reorthog));
710 :
711 18 : PetscOptionsHeadEnd();
712 18 : PetscFunctionReturn(PETSC_SUCCESS);
713 : }
714 :
715 19 : static PetscErrorCode EPSLanczosSetReorthog_Lanczos(EPS eps,EPSLanczosReorthogType reorthog)
716 : {
717 19 : EPS_LANCZOS *lanczos = (EPS_LANCZOS*)eps->data;
718 :
719 19 : PetscFunctionBegin;
720 19 : switch (reorthog) {
721 19 : case EPS_LANCZOS_REORTHOG_LOCAL:
722 : case EPS_LANCZOS_REORTHOG_FULL:
723 : case EPS_LANCZOS_REORTHOG_DELAYED:
724 : case EPS_LANCZOS_REORTHOG_SELECTIVE:
725 : case EPS_LANCZOS_REORTHOG_PERIODIC:
726 : case EPS_LANCZOS_REORTHOG_PARTIAL:
727 19 : if (lanczos->reorthog != reorthog) {
728 19 : lanczos->reorthog = reorthog;
729 19 : eps->state = EPS_STATE_INITIAL;
730 : }
731 19 : break;
732 0 : default:
733 0 : SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid reorthogonalization type");
734 : }
735 19 : PetscFunctionReturn(PETSC_SUCCESS);
736 : }
737 :
738 : /*@
739 : EPSLanczosSetReorthog - Sets the type of reorthogonalization used during the Lanczos
740 : iteration.
741 :
742 : Logically Collective
743 :
744 : Input Parameters:
745 : + eps - the eigenproblem solver context
746 : - reorthog - the type of reorthogonalization
747 :
748 : Options Database Key:
749 : . -eps_lanczos_reorthog - Sets the reorthogonalization type (either 'local', 'selective',
750 : 'periodic', 'partial', 'full' or 'delayed')
751 :
752 : Level: advanced
753 :
754 : .seealso: EPSLanczosGetReorthog(), EPSLanczosReorthogType
755 : @*/
756 19 : PetscErrorCode EPSLanczosSetReorthog(EPS eps,EPSLanczosReorthogType reorthog)
757 : {
758 19 : PetscFunctionBegin;
759 19 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
760 76 : PetscValidLogicalCollectiveEnum(eps,reorthog,2);
761 19 : PetscTryMethod(eps,"EPSLanczosSetReorthog_C",(EPS,EPSLanczosReorthogType),(eps,reorthog));
762 19 : PetscFunctionReturn(PETSC_SUCCESS);
763 : }
764 :
765 5 : static PetscErrorCode EPSLanczosGetReorthog_Lanczos(EPS eps,EPSLanczosReorthogType *reorthog)
766 : {
767 5 : EPS_LANCZOS *lanczos = (EPS_LANCZOS*)eps->data;
768 :
769 5 : PetscFunctionBegin;
770 5 : *reorthog = lanczos->reorthog;
771 5 : PetscFunctionReturn(PETSC_SUCCESS);
772 : }
773 :
774 : /*@
775 : EPSLanczosGetReorthog - Gets the type of reorthogonalization used during
776 : the Lanczos iteration.
777 :
778 : Not Collective
779 :
780 : Input Parameter:
781 : . eps - the eigenproblem solver context
782 :
783 : Output Parameter:
784 : . reorthog - the type of reorthogonalization
785 :
786 : Level: advanced
787 :
788 : .seealso: EPSLanczosSetReorthog(), EPSLanczosReorthogType
789 : @*/
790 5 : PetscErrorCode EPSLanczosGetReorthog(EPS eps,EPSLanczosReorthogType *reorthog)
791 : {
792 5 : PetscFunctionBegin;
793 5 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
794 5 : PetscAssertPointer(reorthog,2);
795 5 : PetscUseMethod(eps,"EPSLanczosGetReorthog_C",(EPS,EPSLanczosReorthogType*),(eps,reorthog));
796 5 : PetscFunctionReturn(PETSC_SUCCESS);
797 : }
798 :
799 22 : static PetscErrorCode EPSReset_Lanczos(EPS eps)
800 : {
801 22 : EPS_LANCZOS *lanczos = (EPS_LANCZOS*)eps->data;
802 :
803 22 : PetscFunctionBegin;
804 22 : PetscCall(BVDestroy(&lanczos->AV));
805 22 : lanczos->allocsize = 0;
806 22 : PetscFunctionReturn(PETSC_SUCCESS);
807 : }
808 :
809 19 : static PetscErrorCode EPSDestroy_Lanczos(EPS eps)
810 : {
811 19 : PetscFunctionBegin;
812 19 : PetscCall(PetscFree(eps->data));
813 19 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLanczosSetReorthog_C",NULL));
814 19 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLanczosGetReorthog_C",NULL));
815 19 : PetscFunctionReturn(PETSC_SUCCESS);
816 : }
817 :
818 0 : static PetscErrorCode EPSView_Lanczos(EPS eps,PetscViewer viewer)
819 : {
820 0 : EPS_LANCZOS *lanczos = (EPS_LANCZOS*)eps->data;
821 0 : PetscBool isascii;
822 :
823 0 : PetscFunctionBegin;
824 0 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
825 0 : if (isascii) {
826 0 : if (lanczos->reorthog != (EPSLanczosReorthogType)-1) PetscCall(PetscViewerASCIIPrintf(viewer," %s reorthogonalization\n",EPSLanczosReorthogTypes[lanczos->reorthog]));
827 : }
828 0 : PetscFunctionReturn(PETSC_SUCCESS);
829 : }
830 :
831 19 : SLEPC_EXTERN PetscErrorCode EPSCreate_Lanczos(EPS eps)
832 : {
833 19 : EPS_LANCZOS *ctx;
834 :
835 19 : PetscFunctionBegin;
836 19 : PetscCall(PetscNew(&ctx));
837 19 : eps->data = (void*)ctx;
838 19 : ctx->reorthog = (EPSLanczosReorthogType)-1;
839 :
840 19 : eps->useds = PETSC_TRUE;
841 :
842 19 : eps->ops->solve = EPSSolve_Lanczos;
843 19 : eps->ops->setup = EPSSetUp_Lanczos;
844 19 : eps->ops->setupsort = EPSSetUpSort_Default;
845 19 : eps->ops->setfromoptions = EPSSetFromOptions_Lanczos;
846 19 : eps->ops->destroy = EPSDestroy_Lanczos;
847 19 : eps->ops->reset = EPSReset_Lanczos;
848 19 : eps->ops->view = EPSView_Lanczos;
849 19 : eps->ops->backtransform = EPSBackTransform_Default;
850 19 : eps->ops->computevectors = EPSComputeVectors_Hermitian;
851 :
852 19 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLanczosSetReorthog_C",EPSLanczosSetReorthog_Lanczos));
853 19 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLanczosGetReorthog_C",EPSLanczosGetReorthog_Lanczos));
854 19 : PetscFunctionReturn(PETSC_SUCCESS);
855 : }
|