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 1857 : 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 7578 : for (j=k;j<m;j++) {
107 7058 : PetscCall(BVMatMultColumn(eps->V,Op,j));
108 7058 : which[j] = PETSC_TRUE;
109 7058 : if (j-2>=k) which[j-2] = PETSC_FALSE;
110 7058 : PetscCall(BVOrthogonalizeSomeColumn(eps->V,j+1,which,hwork,beta+j,breakdown));
111 7058 : alpha[j] = PetscRealPart(hwork[j]);
112 7058 : if (PetscUnlikely(*breakdown)) {
113 0 : *M = j+1;
114 0 : break;
115 7058 : } 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 673 : static PetscErrorCode DenseTridiagonal(PetscInt n_,PetscReal *D,PetscReal *E,PetscReal *w,PetscScalar *V)
140 : {
141 673 : PetscReal abstol = 0.0,vl,vu,*work;
142 673 : PetscBLASInt il,iu,m,*isuppz,n,lwork,*iwork,liwork,info;
143 673 : const char *jobz;
144 : #if defined(PETSC_USE_COMPLEX)
145 : PetscInt i,j;
146 : PetscReal *VV=NULL;
147 : #endif
148 :
149 673 : PetscFunctionBegin;
150 673 : PetscCall(PetscBLASIntCast(n_,&n));
151 673 : PetscCall(PetscBLASIntCast(20*n_,&lwork));
152 673 : PetscCall(PetscBLASIntCast(10*n_,&liwork));
153 673 : if (V) {
154 : jobz = "V";
155 : #if defined(PETSC_USE_COMPLEX)
156 : PetscCall(PetscMalloc1(n*n,&VV));
157 : #endif
158 0 : } else jobz = "N";
159 673 : PetscCall(PetscMalloc3(2*n,&isuppz,lwork,&work,liwork,&iwork));
160 673 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
161 : #if defined(PETSC_USE_COMPLEX)
162 : 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 673 : 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 673 : PetscCall(PetscFPTrapPop());
167 673 : SlepcCheckLapackInfo("stevr",info);
168 : #if defined(PETSC_USE_COMPLEX)
169 : if (V) {
170 : for (i=0;i<n;i++)
171 : for (j=0;j<n;j++)
172 : V[i*n+j] = VV[i*n+j];
173 : PetscCall(PetscFree(VV));
174 : }
175 : #endif
176 673 : PetscCall(PetscFree3(isuppz,work,iwork));
177 673 : 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 87 : for (i=0;i<k;i++) which[i] = PETSC_TRUE;
196 38 : PetscCall(STGetOperator(eps->st,&Op));
197 :
198 711 : for (j=k;j<m;j++) {
199 673 : PetscCall(BVSetActiveColumns(eps->V,0,m));
200 :
201 : /* Lanczos step */
202 673 : PetscCall(BVMatMultColumn(eps->V,Op,j));
203 673 : which[j] = PETSC_TRUE;
204 673 : if (j-2>=k) which[j-2] = PETSC_FALSE;
205 673 : PetscCall(BVOrthogonalizeSomeColumn(eps->V,j+1,which,hwork,&norm,breakdown));
206 673 : alpha[j] = PetscRealPart(hwork[j]);
207 673 : beta[j] = norm;
208 673 : if (PetscUnlikely(*breakdown)) {
209 0 : *M = j+1;
210 0 : break;
211 : }
212 :
213 : /* Compute eigenvalues and eigenvectors Y of the tridiagonal block */
214 673 : n = j-k+1;
215 6994 : for (i=0;i<n;i++) {
216 6321 : d[i] = alpha[i+k];
217 6321 : e[i] = beta[i+k];
218 : }
219 673 : PetscCall(DenseTridiagonal(n,d,e,ritz,Y));
220 :
221 : /* Estimate ||A|| */
222 6994 : for (i=0;i<n;i++)
223 6321 : if (PetscAbsReal(ritz[i]) > anorm) anorm = PetscAbsReal(ritz[i]);
224 :
225 : /* Compute nearly converged Ritz vectors */
226 : nritzo = 0;
227 6994 : for (i=0;i<n;i++) {
228 6321 : if (norm*PetscAbsScalar(Y[i*n+n-1]) < PETSC_SQRT_MACHINE_EPSILON*anorm) nritzo++;
229 : }
230 673 : if (nritzo>nritz) {
231 : nritz = 0;
232 196 : for (i=0;i<n;i++) {
233 173 : if (norm*PetscAbsScalar(Y[i*n+n-1]) < PETSC_SQRT_MACHINE_EPSILON*anorm) {
234 37 : PetscCall(BVSetActiveColumns(eps->V,k,k+n));
235 37 : PetscCall(BVGetColumn(lanczos->AV,nritz,&av));
236 37 : PetscCall(BVMultVec(eps->V,1.0,0.0,av,Y+i*n));
237 37 : PetscCall(BVRestoreColumn(lanczos->AV,nritz,&av));
238 37 : nritz++;
239 : }
240 : }
241 : }
242 673 : if (nritz > 0) {
243 247 : PetscCall(BVGetColumn(eps->V,j+1,&vj1));
244 247 : PetscCall(BVSetActiveColumns(lanczos->AV,0,nritz));
245 247 : PetscCall(BVOrthogonalizeVec(lanczos->AV,vj1,hwork,&norm,breakdown));
246 247 : PetscCall(BVRestoreColumn(eps->V,j+1,&vj1));
247 247 : if (PetscUnlikely(*breakdown)) {
248 0 : *M = j+1;
249 0 : break;
250 : }
251 : }
252 673 : 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 768 : static void update_omega(PetscReal *omega,PetscReal *omega_old,PetscInt j,PetscReal *alpha,PetscReal *beta,PetscReal eps1,PetscReal anorm)
261 : {
262 768 : PetscInt k;
263 768 : PetscReal T,binv;
264 :
265 768 : 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 768 : T = eps1*anorm;
271 768 : binv = 1.0/beta[j+1];
272 :
273 : /* Update omega(1) using omega(0)==0 */
274 768 : omega_old[0]= beta[1]*omega[1] + (alpha[0]-alpha[j])*omega[0] - beta[j]*omega_old[0];
275 768 : if (omega_old[0] > 0) omega_old[0] = binv*(omega_old[0] + T);
276 461 : else omega_old[0] = binv*(omega_old[0] - T);
277 :
278 : /* Update remaining components */
279 7014 : for (k=1;k<j-1;k++) {
280 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];
281 6246 : if (omega_old[k] > 0) omega_old[k] = binv*(omega_old[k] + T);
282 3192 : else omega_old[k] = binv*(omega_old[k] - T);
283 : }
284 768 : omega_old[j-1] = binv*T;
285 :
286 : /* Swap omega and omega_old */
287 8534 : for (k=0;k<j;k++) {
288 7766 : omega[k] = omega_old[k];
289 7766 : omega_old[k] = omega[k];
290 : }
291 768 : omega[j] = eps1;
292 768 : PetscFunctionReturnVoid();
293 : }
294 :
295 48 : static void compute_int(PetscBool *which,PetscReal *mu,PetscInt j,PetscReal delta,PetscReal eta)
296 : {
297 48 : PetscInt i,k,maxpos;
298 48 : PetscReal max;
299 48 : PetscBool found;
300 :
301 48 : PetscFunctionBegin;
302 : /* initialize which */
303 48 : found = PETSC_FALSE;
304 48 : maxpos = 0;
305 48 : max = 0.0;
306 755 : for (i=0;i<j;i++) {
307 707 : if (PetscAbsReal(mu[i]) >= delta) {
308 175 : which[i] = PETSC_TRUE;
309 175 : found = PETSC_TRUE;
310 532 : } else which[i] = PETSC_FALSE;
311 707 : if (PetscAbsReal(mu[i]) > max) {
312 78 : maxpos = i;
313 78 : max = PetscAbsReal(mu[i]);
314 : }
315 : }
316 48 : if (!found) which[maxpos] = PETSC_TRUE;
317 :
318 755 : for (i=0;i<j;i++) {
319 707 : if (which[i]) {
320 : /* find left interval */
321 623 : for (k=i;k>=0;k--) {
322 623 : if (PetscAbsReal(mu[k])<eta || which[k]) break;
323 0 : else which[k] = PETSC_TRUE;
324 : }
325 : /* find right interval */
326 1071 : for (k=i+1;k<j;k++) {
327 1037 : if (PetscAbsReal(mu[k])<eta || which[k]) break;
328 448 : else which[k] = PETSC_TRUE;
329 : }
330 : }
331 : }
332 48 : 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 174 : 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 1422 : for (j=k;j<m;j++) {
372 1346 : PetscCall(BVMatMultColumn(eps->V,Op,j));
373 1346 : if (fro) {
374 : /* Lanczos step with full reorthogonalization */
375 502 : PetscCall(BVOrthogonalizeColumn(eps->V,j+1,hwork,&norm,breakdown));
376 502 : alpha[j] = PetscRealPart(hwork[j]);
377 : } else {
378 : /* Lanczos step */
379 844 : which[j] = PETSC_TRUE;
380 844 : if (j-2>=k) which[j-2] = PETSC_FALSE;
381 844 : PetscCall(BVOrthogonalizeSomeColumn(eps->V,j+1,which,hwork,&norm,breakdown));
382 844 : alpha[j] = PetscRealPart(hwork[j]);
383 844 : beta[j] = norm;
384 :
385 : /* Estimate ||A|| if needed */
386 844 : 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 844 : reorth = PETSC_FALSE;
393 844 : if (j>k) {
394 768 : update_omega(omega,omega_old,j,alpha,beta-1,eps1,anorm);
395 8362 : for (i=0;i<j-k;i++) {
396 6826 : if (PetscAbsReal(omega[i]) > delta) reorth = PETSC_TRUE;
397 : }
398 : }
399 844 : if (reorth || force_reorth) {
400 290 : for (i=0;i<k;i++) which2[i] = PETSC_FALSE;
401 2653 : for (i=k;i<=j;i++) which2[i] = PETSC_TRUE;
402 159 : if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_PERIODIC) {
403 : /* Periodic reorthogonalization */
404 79 : if (force_reorth) force_reorth = PETSC_FALSE;
405 47 : else force_reorth = PETSC_TRUE;
406 1238 : for (i=0;i<j-k;i++) omega[i] = eps1;
407 : } else {
408 : /* Partial reorthogonalization */
409 80 : if (force_reorth) force_reorth = PETSC_FALSE;
410 : else {
411 48 : force_reorth = PETSC_TRUE;
412 48 : compute_int(which2+k,omega,j-k,delta,eta);
413 803 : for (i=0;i<j-k;i++) {
414 707 : if (which2[i+k]) omega[i] = eps1;
415 : }
416 : }
417 : }
418 159 : PetscCall(BVOrthogonalizeSomeColumn(eps->V,j+1,which2,hwork,&norm,breakdown));
419 : }
420 : }
421 :
422 1346 : if (PetscUnlikely(*breakdown || norm < eps->n*anorm*PETSC_MACHINE_EPSILON)) {
423 0 : *M = j+1;
424 0 : break;
425 : }
426 1346 : if (!fro && norm*delta < anorm*eps1) {
427 30 : fro = PETSC_TRUE;
428 30 : PetscCall(PetscInfo(eps,"Switching to full reorthogonalization at iteration %" PetscInt_FMT "\n",eps->its));
429 : }
430 1346 : beta[j] = norm;
431 1346 : 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 932 : static PetscErrorCode EPSBasicLanczos(EPS eps,PetscInt k,PetscInt *m,PetscReal *betam,PetscBool *breakdown,PetscReal anorm)
457 : {
458 932 : EPS_LANCZOS *lanczos = (EPS_LANCZOS*)eps->data;
459 932 : PetscScalar *T;
460 932 : PetscInt i,n=*m,ld;
461 932 : PetscReal *alpha,*beta;
462 932 : BVOrthogRefineType orthog_ref;
463 932 : Mat Op,M;
464 :
465 932 : PetscFunctionBegin;
466 932 : PetscCall(DSGetLeadingDimension(eps->ds,&ld));
467 932 : 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 271 : case EPS_LANCZOS_REORTHOG_FULL:
476 271 : PetscCall(STGetOperator(eps->st,&Op));
477 271 : PetscCall(DSGetMat(eps->ds,DS_MAT_T,&M));
478 271 : PetscCall(BVMatLanczos(eps->V,Op,M,k,m,betam,breakdown));
479 271 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_T,&M));
480 271 : 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 27 : case EPS_LANCZOS_REORTHOG_DELAYED:
498 27 : PetscCall(PetscMalloc1(n*n,&T));
499 27 : PetscCall(BVGetOrthogonalization(eps->V,NULL,&orthog_ref,NULL,NULL));
500 27 : if (orthog_ref == BV_ORTHOG_REFINE_NEVER) PetscCall(EPSDelayedArnoldi1(eps,T,n,k,m,betam,breakdown));
501 27 : else PetscCall(EPSDelayedArnoldi(eps,T,n,k,m,betam,breakdown));
502 27 : n = *m;
503 27 : PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&alpha));
504 27 : beta = alpha + ld;
505 472 : for (i=k;i<n-1;i++) {
506 445 : alpha[i] = PetscRealPart(T[n*i+i]);
507 445 : beta[i] = PetscRealPart(T[n*i+i+1]);
508 : }
509 27 : alpha[n-1] = PetscRealPart(T[n*(n-1)+n-1]);
510 27 : beta[n-1] = *betam;
511 27 : PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&alpha));
512 27 : PetscCall(PetscFree(T));
513 : break;
514 : }
515 932 : 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 968 : while (eps->reason == EPS_CONVERGED_ITERATING) {
541 932 : eps->its++;
542 :
543 : /* Compute an ncv-step Lanczos factorization */
544 932 : n = PetscMin(nconv+eps->mpd,ncv);
545 932 : PetscCall(DSSetDimensions(eps->ds,n,nconv,PETSC_DETERMINE));
546 932 : PetscCall(EPSBasicLanczos(eps,nconv,&n,&beta,&breakdown,anorm));
547 932 : PetscCall(DSSetDimensions(eps->ds,n,nconv,0));
548 932 : PetscCall(DSSetState(eps->ds,DS_STATE_INTERMEDIATE));
549 932 : PetscCall(BVSetActiveColumns(eps->V,nconv,n));
550 :
551 : /* Solve projected problem */
552 932 : PetscCall(DSSolve(eps->ds,ritz,NULL));
553 932 : PetscCall(DSSort(eps->ds,ritz,NULL,NULL,NULL,NULL));
554 932 : PetscCall(DSSynchronize(eps->ds,ritz,NULL));
555 :
556 : /* Estimate ||A|| */
557 13272 : for (i=nconv;i<n;i++)
558 12761 : anorm = PetscMax(anorm,PetscAbsReal(PetscRealPart(ritz[i])));
559 :
560 : /* Compute residual norm estimates as beta*abs(Y(m,:)) + eps*||A|| */
561 932 : PetscCall(DSGetArray(eps->ds,DS_MAT_Q,&Y));
562 13272 : for (i=nconv;i<n;i++) {
563 12340 : resnorm = beta*PetscAbsScalar(Y[n-1+i*ld]) + PETSC_MACHINE_EPSILON*anorm;
564 12340 : PetscCall((*eps->converged)(eps,ritz[i],eps->eigi[i],resnorm,&bnd[i],eps->convergedctx));
565 12340 : if (bnd[i]<eps->tol) conv[i] = 'C';
566 12171 : else conv[i] = 'N';
567 : }
568 932 : PetscCall(DSRestoreArray(eps->ds,DS_MAT_Q,&Y));
569 :
570 : /* purge repeated ritz values */
571 932 : if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_LOCAL) {
572 7058 : for (i=nconv+1;i<n;i++) {
573 6538 : 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 932 : 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 1078 : while (restart<n && conv[restart] != 'N') restart++;
582 931 : if (restart >= n) {
583 0 : breakdown = PETSC_TRUE;
584 : } else {
585 12175 : for (i=restart+1;i<n;i++) {
586 11244 : if (conv[i] == 'N') {
587 11240 : PetscCall(SlepcSCCompare(eps->sc,ritz[restart],0.0,ritz[i],0.0,&r));
588 11240 : if (r>0) restart = i;
589 : }
590 : }
591 931 : PetscCall(DSGetArray(eps->ds,DS_MAT_Q,&Y));
592 931 : PetscCall(BVMultColumn(eps->V,1.0,0.0,n,Y+restart*ld+nconv));
593 931 : PetscCall(DSRestoreArray(eps->ds,DS_MAT_Q,&Y));
594 : }
595 : }
596 :
597 : /* Count and put converged eigenvalues first */
598 13272 : for (i=nconv;i<n;i++) perm[i] = i;
599 1098 : for (k=nconv;k<n;k++) {
600 1097 : if (conv[perm[k]] != 'C') {
601 934 : j = k + 1;
602 12181 : while (j<n && conv[perm[j]] != 'C') j++;
603 934 : if (j>=n) break;
604 3 : l = perm[k]; perm[k] = perm[j]; perm[j] = l;
605 : }
606 : }
607 :
608 : /* Sort eigenvectors according to permutation */
609 932 : PetscCall(DSGetArray(eps->ds,DS_MAT_Q,&Y));
610 1098 : for (i=nconv;i<k;i++) {
611 166 : x = perm[i];
612 166 : if (x != i) {
613 3 : j = i + 1;
614 7 : while (perm[j] != i) j++;
615 : /* swap eigenvalues i and j */
616 3 : stmp = ritz[x]; ritz[x] = ritz[i]; ritz[i] = stmp;
617 3 : rtmp = bnd[x]; bnd[x] = bnd[i]; bnd[i] = rtmp;
618 3 : ctmp = conv[x]; conv[x] = conv[i]; conv[i] = ctmp;
619 3 : perm[j] = x; perm[i] = i;
620 : /* swap eigenvectors i and j */
621 60 : for (l=0;l<n;l++) {
622 57 : stmp = Y[l+x*ld]; Y[l+x*ld] = Y[l+i*ld]; Y[l+i*ld] = stmp;
623 : }
624 : }
625 : }
626 932 : PetscCall(DSRestoreArray(eps->ds,DS_MAT_Q,&Y));
627 :
628 : /* compute converged eigenvectors */
629 932 : PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&U));
630 932 : PetscCall(BVMultInPlace(eps->V,U,nconv,k));
631 932 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&U));
632 :
633 : /* purge spurious ritz values */
634 932 : 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 13272 : for (i=nconv;i<n;i++) {
669 12340 : eps->eigr[i] = ritz[i];
670 12340 : eps->errest[i] = bnd[i];
671 : }
672 932 : nconv = k;
673 932 : PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,n));
674 932 : PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,nconv,eps->nev,&eps->reason,eps->stoppingctx));
675 :
676 932 : if (eps->reason == EPS_CONVERGED_ITERATING) { /* copy restart vector */
677 896 : PetscCall(BVCopyColumn(eps->V,n,nconv));
678 896 : 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 896 : 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 896 : if (PetscUnlikely(breakdown)) { /* give up */
688 0 : eps->reason = EPS_DIVERGED_BREAKDOWN;
689 968 : 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 : }
|
