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 514 : static PetscErrorCode EPSLocalLanczos(EPS eps,PetscReal *alpha,PetscReal *beta,PetscInt k,PetscInt *M,PetscBool *breakdown)
89 : {
90 514 : PetscInt i,j,m = *M;
91 514 : Mat Op;
92 514 : PetscBool *which,lwhich[100];
93 514 : PetscScalar *hwork,lhwork[100];
94 :
95 514 : PetscFunctionBegin;
96 514 : if (m > 100) PetscCall(PetscMalloc2(m,&which,m,&hwork));
97 : else {
98 514 : which = lwhich;
99 514 : hwork = lhwork;
100 : }
101 1784 : for (i=0;i<k;i++) which[i] = PETSC_TRUE;
102 :
103 514 : PetscCall(BVSetActiveColumns(eps->V,0,m));
104 514 : PetscCall(STGetOperator(eps->st,&Op));
105 7480 : for (j=k;j<m;j++) {
106 6966 : PetscCall(BVMatMultColumn(eps->V,Op,j));
107 6966 : which[j] = PETSC_TRUE;
108 6966 : if (j-2>=k) which[j-2] = PETSC_FALSE;
109 6966 : PetscCall(BVOrthogonalizeSomeColumn(eps->V,j+1,which,hwork,beta+j,breakdown));
110 6966 : alpha[j] = PetscRealPart(hwork[j]);
111 6966 : if (PetscUnlikely(*breakdown)) {
112 0 : *M = j+1;
113 0 : break;
114 6966 : } else PetscCall(BVScaleColumn(eps->V,j+1,1/beta[j]));
115 : }
116 514 : PetscCall(STRestoreOperator(eps->st,&Op));
117 514 : if (m > 100) PetscCall(PetscFree2(which,hwork));
118 514 : 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 676 : static PetscErrorCode DenseTridiagonal(PetscInt n_,PetscReal *D,PetscReal *E,PetscReal *w,PetscScalar *V)
139 : {
140 676 : PetscReal abstol = 0.0,vl,vu,*work;
141 676 : PetscBLASInt il,iu,m,*isuppz,n,lwork,*iwork,liwork,info;
142 676 : const char *jobz;
143 : #if defined(PETSC_USE_COMPLEX)
144 676 : PetscInt i,j;
145 676 : PetscReal *VV=NULL;
146 : #endif
147 :
148 676 : PetscFunctionBegin;
149 676 : PetscCall(PetscBLASIntCast(n_,&n));
150 676 : PetscCall(PetscBLASIntCast(20*n_,&lwork));
151 676 : PetscCall(PetscBLASIntCast(10*n_,&liwork));
152 676 : if (V) {
153 676 : jobz = "V";
154 : #if defined(PETSC_USE_COMPLEX)
155 676 : PetscCall(PetscMalloc1(n*n,&VV));
156 : #endif
157 : } else jobz = "N";
158 676 : PetscCall(PetscMalloc3(2*n,&isuppz,lwork,&work,liwork,&iwork));
159 676 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
160 : #if defined(PETSC_USE_COMPLEX)
161 676 : 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 : 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 676 : PetscCall(PetscFPTrapPop());
166 676 : SlepcCheckLapackInfo("stevr",info);
167 : #if defined(PETSC_USE_COMPLEX)
168 676 : if (V) {
169 7053 : for (i=0;i<n;i++)
170 84762 : for (j=0;j<n;j++)
171 78385 : V[i*n+j] = VV[i*n+j];
172 676 : PetscCall(PetscFree(VV));
173 : }
174 : #endif
175 676 : PetscCall(PetscFree3(isuppz,work,iwork));
176 676 : 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 84 : for (i=0;i<k;i++) which[i] = PETSC_TRUE;
195 38 : PetscCall(STGetOperator(eps->st,&Op));
196 :
197 714 : for (j=k;j<m;j++) {
198 676 : PetscCall(BVSetActiveColumns(eps->V,0,m));
199 :
200 : /* Lanczos step */
201 676 : PetscCall(BVMatMultColumn(eps->V,Op,j));
202 676 : which[j] = PETSC_TRUE;
203 676 : if (j-2>=k) which[j-2] = PETSC_FALSE;
204 676 : PetscCall(BVOrthogonalizeSomeColumn(eps->V,j+1,which,hwork,&norm,breakdown));
205 676 : alpha[j] = PetscRealPart(hwork[j]);
206 676 : beta[j] = norm;
207 676 : if (PetscUnlikely(*breakdown)) {
208 0 : *M = j+1;
209 0 : break;
210 : }
211 :
212 : /* Compute eigenvalues and eigenvectors Y of the tridiagonal block */
213 676 : n = j-k+1;
214 7053 : for (i=0;i<n;i++) {
215 6377 : d[i] = alpha[i+k];
216 6377 : e[i] = beta[i+k];
217 : }
218 676 : PetscCall(DenseTridiagonal(n,d,e,ritz,Y));
219 :
220 : /* Estimate ||A|| */
221 7053 : for (i=0;i<n;i++)
222 6377 : if (PetscAbsReal(ritz[i]) > anorm) anorm = PetscAbsReal(ritz[i]);
223 :
224 : /* Compute nearly converged Ritz vectors */
225 : nritzo = 0;
226 7053 : for (i=0;i<n;i++) {
227 6377 : if (norm*PetscAbsScalar(Y[i*n+n-1]) < PETSC_SQRT_MACHINE_EPSILON*anorm) nritzo++;
228 : }
229 676 : if (nritzo>nritz) {
230 : nritz = 0;
231 160 : for (i=0;i<n;i++) {
232 138 : if (norm*PetscAbsScalar(Y[i*n+n-1]) < PETSC_SQRT_MACHINE_EPSILON*anorm) {
233 32 : PetscCall(BVSetActiveColumns(eps->V,k,k+n));
234 32 : PetscCall(BVGetColumn(lanczos->AV,nritz,&av));
235 32 : PetscCall(BVMultVec(eps->V,1.0,0.0,av,Y+i*n));
236 32 : PetscCall(BVRestoreColumn(lanczos->AV,nritz,&av));
237 32 : nritz++;
238 : }
239 : }
240 : }
241 676 : if (nritz > 0) {
242 255 : PetscCall(BVGetColumn(eps->V,j+1,&vj1));
243 255 : PetscCall(BVSetActiveColumns(lanczos->AV,0,nritz));
244 255 : PetscCall(BVOrthogonalizeVec(lanczos->AV,vj1,hwork,&norm,breakdown));
245 255 : PetscCall(BVRestoreColumn(eps->V,j+1,&vj1));
246 255 : if (PetscUnlikely(*breakdown)) {
247 0 : *M = j+1;
248 0 : break;
249 : }
250 : }
251 676 : 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 842 : static void update_omega(PetscReal *omega,PetscReal *omega_old,PetscInt j,PetscReal *alpha,PetscReal *beta,PetscReal eps1,PetscReal anorm)
260 : {
261 842 : PetscInt k;
262 842 : PetscReal T,binv;
263 :
264 842 : 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 842 : T = eps1*anorm;
270 842 : binv = 1.0/beta[j+1];
271 :
272 : /* Update omega(1) using omega(0)==0 */
273 842 : omega_old[0]= beta[1]*omega[1] + (alpha[0]-alpha[j])*omega[0] - beta[j]*omega_old[0];
274 842 : if (omega_old[0] > 0) omega_old[0] = binv*(omega_old[0] + T);
275 506 : else omega_old[0] = binv*(omega_old[0] - T);
276 :
277 : /* Update remaining components */
278 7632 : for (k=1;k<j-1;k++) {
279 6790 : 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 6790 : if (omega_old[k] > 0) omega_old[k] = binv*(omega_old[k] + T);
281 3594 : else omega_old[k] = binv*(omega_old[k] - T);
282 : }
283 842 : omega_old[j-1] = binv*T;
284 :
285 : /* Swap omega and omega_old */
286 9296 : for (k=0;k<j;k++) {
287 8454 : omega[k] = omega_old[k];
288 8454 : omega_old[k] = omega[k];
289 : }
290 842 : omega[j] = eps1;
291 842 : PetscFunctionReturnVoid();
292 : }
293 :
294 53 : static void compute_int(PetscBool *which,PetscReal *mu,PetscInt j,PetscReal delta,PetscReal eta)
295 : {
296 53 : PetscInt i,k,maxpos;
297 53 : PetscReal max;
298 53 : PetscBool found;
299 :
300 53 : PetscFunctionBegin;
301 : /* initialize which */
302 53 : found = PETSC_FALSE;
303 53 : maxpos = 0;
304 53 : max = 0.0;
305 843 : for (i=0;i<j;i++) {
306 790 : if (PetscAbsReal(mu[i]) >= delta) {
307 189 : which[i] = PETSC_TRUE;
308 189 : found = PETSC_TRUE;
309 601 : } else which[i] = PETSC_FALSE;
310 790 : if (PetscAbsReal(mu[i]) > max) {
311 88 : maxpos = i;
312 88 : max = PetscAbsReal(mu[i]);
313 : }
314 : }
315 53 : if (!found) which[maxpos] = PETSC_TRUE;
316 :
317 843 : for (i=0;i<j;i++) {
318 790 : if (which[i]) {
319 : /* find left interval */
320 689 : for (k=i;k>=0;k--) {
321 689 : if (PetscAbsReal(mu[k])<eta || which[k]) break;
322 0 : else which[k] = PETSC_TRUE;
323 : }
324 : /* find right interval */
325 1189 : for (k=i+1;k<j;k++) {
326 1159 : if (PetscAbsReal(mu[k])<eta || which[k]) break;
327 500 : else which[k] = PETSC_TRUE;
328 : }
329 : }
330 : }
331 53 : 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 168 : 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 1428 : for (j=k;j<m;j++) {
371 1352 : PetscCall(BVMatMultColumn(eps->V,Op,j));
372 1352 : if (fro) {
373 : /* Lanczos step with full reorthogonalization */
374 434 : PetscCall(BVOrthogonalizeColumn(eps->V,j+1,hwork,&norm,breakdown));
375 434 : alpha[j] = PetscRealPart(hwork[j]);
376 : } else {
377 : /* Lanczos step */
378 918 : which[j] = PETSC_TRUE;
379 918 : if (j-2>=k) which[j-2] = PETSC_FALSE;
380 918 : PetscCall(BVOrthogonalizeSomeColumn(eps->V,j+1,which,hwork,&norm,breakdown));
381 918 : alpha[j] = PetscRealPart(hwork[j]);
382 918 : beta[j] = norm;
383 :
384 : /* Estimate ||A|| if needed */
385 918 : 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 918 : reorth = PETSC_FALSE;
392 918 : if (j>k) {
393 842 : update_omega(omega,omega_old,j,alpha,beta-1,eps1,anorm);
394 9228 : for (i=0;i<j-k;i++) {
395 7544 : if (PetscAbsReal(omega[i]) > delta) reorth = PETSC_TRUE;
396 : }
397 : }
398 918 : if (reorth || force_reorth) {
399 298 : for (i=0;i<k;i++) which2[i] = PETSC_FALSE;
400 2905 : for (i=k;i<=j;i++) which2[i] = PETSC_TRUE;
401 173 : if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_PERIODIC) {
402 : /* Periodic reorthogonalization */
403 85 : if (force_reorth) force_reorth = PETSC_FALSE;
404 50 : else force_reorth = PETSC_TRUE;
405 1338 : for (i=0;i<j-k;i++) omega[i] = eps1;
406 : } else {
407 : /* Partial reorthogonalization */
408 88 : if (force_reorth) force_reorth = PETSC_FALSE;
409 : else {
410 53 : force_reorth = PETSC_TRUE;
411 53 : compute_int(which2+k,omega,j-k,delta,eta);
412 896 : for (i=0;i<j-k;i++) {
413 790 : if (which2[i+k]) omega[i] = eps1;
414 : }
415 : }
416 : }
417 173 : PetscCall(BVOrthogonalizeSomeColumn(eps->V,j+1,which2,hwork,&norm,breakdown));
418 : }
419 : }
420 :
421 1352 : if (PetscUnlikely(*breakdown || norm < eps->n*anorm*PETSC_MACHINE_EPSILON)) {
422 0 : *M = j+1;
423 0 : break;
424 : }
425 1352 : if (!fro && norm*delta < anorm*eps1) {
426 26 : fro = PETSC_TRUE;
427 26 : PetscCall(PetscInfo(eps,"Switching to full reorthogonalization at iteration %" PetscInt_FMT "\n",eps->its));
428 : }
429 1352 : beta[j] = norm;
430 1352 : 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 938 : static PetscErrorCode EPSBasicLanczos(EPS eps,PetscInt k,PetscInt *m,PetscReal *betam,PetscBool *breakdown,PetscReal anorm)
456 : {
457 938 : EPS_LANCZOS *lanczos = (EPS_LANCZOS*)eps->data;
458 938 : PetscScalar *T;
459 938 : PetscInt i,n=*m,ld;
460 938 : PetscReal *alpha,*beta;
461 938 : BVOrthogRefineType orthog_ref;
462 938 : Mat Op,M;
463 :
464 938 : PetscFunctionBegin;
465 938 : PetscCall(DSGetLeadingDimension(eps->ds,&ld));
466 938 : switch (lanczos->reorthog) {
467 514 : case EPS_LANCZOS_REORTHOG_LOCAL:
468 514 : PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&alpha));
469 514 : beta = alpha + ld;
470 514 : PetscCall(EPSLocalLanczos(eps,alpha,beta,k,m,breakdown));
471 514 : *betam = beta[*m-1];
472 514 : PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&alpha));
473 : break;
474 287 : case EPS_LANCZOS_REORTHOG_FULL:
475 287 : PetscCall(STGetOperator(eps->st,&Op));
476 287 : PetscCall(DSGetMat(eps->ds,DS_MAT_T,&M));
477 287 : PetscCall(BVMatLanczos(eps->V,Op,M,k,m,betam,breakdown));
478 287 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_T,&M));
479 287 : 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 23 : case EPS_LANCZOS_REORTHOG_DELAYED:
497 23 : PetscCall(PetscMalloc1(n*n,&T));
498 23 : PetscCall(BVGetOrthogonalization(eps->V,NULL,&orthog_ref,NULL,NULL));
499 23 : if (orthog_ref == BV_ORTHOG_REFINE_NEVER) PetscCall(EPSDelayedArnoldi1(eps,T,n,k,m,betam,breakdown));
500 23 : else PetscCall(EPSDelayedArnoldi(eps,T,n,k,m,betam,breakdown));
501 23 : n = *m;
502 23 : PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&alpha));
503 23 : beta = alpha + ld;
504 403 : for (i=k;i<n-1;i++) {
505 380 : alpha[i] = PetscRealPart(T[n*i+i]);
506 380 : beta[i] = PetscRealPart(T[n*i+i+1]);
507 : }
508 23 : alpha[n-1] = PetscRealPart(T[n*(n-1)+n-1]);
509 23 : beta[n-1] = *betam;
510 23 : PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&alpha));
511 23 : PetscCall(PetscFree(T));
512 : break;
513 : }
514 938 : 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 974 : while (eps->reason == EPS_CONVERGED_ITERATING) {
540 938 : eps->its++;
541 :
542 : /* Compute an ncv-step Lanczos factorization */
543 938 : n = PetscMin(nconv+eps->mpd,ncv);
544 938 : PetscCall(DSSetDimensions(eps->ds,n,nconv,PETSC_DEFAULT));
545 938 : PetscCall(EPSBasicLanczos(eps,nconv,&n,&beta,&breakdown,anorm));
546 938 : PetscCall(DSSetDimensions(eps->ds,n,nconv,0));
547 938 : PetscCall(DSSetState(eps->ds,DS_STATE_INTERMEDIATE));
548 938 : PetscCall(BVSetActiveColumns(eps->V,nconv,n));
549 :
550 : /* Solve projected problem */
551 938 : PetscCall(DSSolve(eps->ds,ritz,NULL));
552 938 : PetscCall(DSSort(eps->ds,ritz,NULL,NULL,NULL,NULL));
553 938 : PetscCall(DSSynchronize(eps->ds,ritz,NULL));
554 :
555 : /* Estimate ||A|| */
556 13337 : for (i=nconv;i<n;i++)
557 12830 : anorm = PetscMax(anorm,PetscAbsReal(PetscRealPart(ritz[i])));
558 :
559 : /* Compute residual norm estimates as beta*abs(Y(m,:)) + eps*||A|| */
560 938 : PetscCall(DSGetArray(eps->ds,DS_MAT_Q,&Y));
561 13337 : for (i=nconv;i<n;i++) {
562 12399 : resnorm = beta*PetscAbsScalar(Y[n-1+i*ld]) + PETSC_MACHINE_EPSILON*anorm;
563 12399 : PetscCall((*eps->converged)(eps,ritz[i],eps->eigi[i],resnorm,&bnd[i],eps->convergedctx));
564 12399 : if (bnd[i]<eps->tol) conv[i] = 'C';
565 12230 : else conv[i] = 'N';
566 : }
567 938 : PetscCall(DSRestoreArray(eps->ds,DS_MAT_Q,&Y));
568 :
569 : /* purge repeated ritz values */
570 938 : if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_LOCAL) {
571 6966 : for (i=nconv+1;i<n;i++) {
572 6452 : 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 938 : 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 1084 : while (restart<n && conv[restart] != 'N') restart++;
581 937 : if (restart >= n) {
582 0 : breakdown = PETSC_TRUE;
583 : } else {
584 12234 : for (i=restart+1;i<n;i++) {
585 11297 : if (conv[i] == 'N') {
586 11293 : PetscCall(SlepcSCCompare(eps->sc,ritz[restart],0.0,ritz[i],0.0,&r));
587 11293 : if (r>0) restart = i;
588 : }
589 : }
590 937 : PetscCall(DSGetArray(eps->ds,DS_MAT_Q,&Y));
591 937 : PetscCall(BVMultColumn(eps->V,1.0,0.0,n,Y+restart*ld+nconv));
592 937 : PetscCall(DSRestoreArray(eps->ds,DS_MAT_Q,&Y));
593 : }
594 : }
595 :
596 : /* Count and put converged eigenvalues first */
597 13337 : for (i=nconv;i<n;i++) perm[i] = i;
598 1104 : for (k=nconv;k<n;k++) {
599 1103 : if (conv[perm[k]] != 'C') {
600 940 : j = k + 1;
601 12240 : while (j<n && conv[perm[j]] != 'C') j++;
602 940 : 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 938 : PetscCall(DSGetArray(eps->ds,DS_MAT_Q,&Y));
609 1104 : 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 938 : PetscCall(DSRestoreArray(eps->ds,DS_MAT_Q,&Y));
626 :
627 : /* compute converged eigenvectors */
628 938 : PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&U));
629 938 : PetscCall(BVMultInPlace(eps->V,U,nconv,k));
630 938 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&U));
631 :
632 : /* purge spurious ritz values */
633 938 : if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_LOCAL) {
634 589 : 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 588 : 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 13337 : for (i=nconv;i<n;i++) {
668 12399 : eps->eigr[i] = ritz[i];
669 12399 : eps->errest[i] = bnd[i];
670 : }
671 938 : nconv = k;
672 938 : PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,n));
673 938 : PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,nconv,eps->nev,&eps->reason,eps->stoppingctx));
674 :
675 938 : if (eps->reason == EPS_CONVERGED_ITERATING) { /* copy restart vector */
676 902 : PetscCall(BVCopyColumn(eps->V,n,nconv));
677 902 : if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_LOCAL && !breakdown) {
678 : /* Reorthonormalize restart vector */
679 496 : PetscCall(BVOrthonormalizeColumn(eps->V,nconv,PETSC_FALSE,NULL,&breakdown));
680 : }
681 902 : 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 902 : if (PetscUnlikely(breakdown)) { /* give up */
687 0 : eps->reason = EPS_DIVERGED_BREAKDOWN;
688 974 : 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 : }
|