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