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