Actual source code: lanczos.c
slepc-3.14.2 2021-02-01
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2020, 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;
42: PetscErrorCode ierr;
45: EPSCheckHermitianDefinite(eps);
46: EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd);
47: if (eps->ncv>eps->nev+eps->mpd) SETERRQ(PetscObjectComm((PetscObject)eps),1,"The value of ncv must not be larger than nev+mpd");
48: if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
49: if (!eps->which) { EPSSetWhichEigenpairs_Default(eps); }
50: if (eps->which==EPS_ALL) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support computing all eigenvalues");
51: EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION);
52: 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) {
75: EPSSetWorkVecs(eps,1);
76: }
77: return(0);
78: }
80: /*
81: EPSLocalLanczos - Local reorthogonalization.
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: static PetscErrorCode EPSLocalLanczos(EPS eps,PetscReal *alpha,PetscReal *beta,PetscInt k,PetscInt *M,PetscBool *breakdown)
90: {
92: PetscInt i,j,m = *M;
93: Mat Op;
94: PetscBool *which,lwhich[100];
95: PetscScalar *hwork,lhwork[100];
98: if (m > 100) {
99: PetscMalloc2(m,&which,m,&hwork);
100: } else {
101: which = lwhich;
102: hwork = lhwork;
103: }
104: for (i=0;i<k;i++) which[i] = PETSC_TRUE;
106: BVSetActiveColumns(eps->V,0,m);
107: STGetOperator(eps->st,&Op);
108: for (j=k;j<m;j++) {
109: BVMatMultColumn(eps->V,Op,j);
110: which[j] = PETSC_TRUE;
111: if (j-2>=k) which[j-2] = PETSC_FALSE;
112: BVOrthogonalizeSomeColumn(eps->V,j+1,which,hwork,beta+j,breakdown);
113: alpha[j] = PetscRealPart(hwork[j]);
114: if (*breakdown) {
115: *M = j+1;
116: break;
117: } else {
118: BVScaleColumn(eps->V,j+1,1/beta[j]);
119: }
120: }
121: STRestoreOperator(eps->st,&Op);
122: if (m > 100) {
123: PetscFree2(which,hwork);
124: }
125: return(0);
126: }
128: /*
129: DenseTridiagonal - Solves a real tridiagonal Hermitian Eigenvalue Problem.
131: Input Parameters:
132: + n - dimension of the eigenproblem
133: . D - pointer to the array containing the diagonal elements
134: - E - pointer to the array containing the off-diagonal elements
136: Output Parameters:
137: + w - pointer to the array to store the computed eigenvalues
138: - V - pointer to the array to store the eigenvectors
140: Notes:
141: If V is NULL then the eigenvectors are not computed.
143: This routine use LAPACK routines xSTEVR.
144: */
145: static PetscErrorCode DenseTridiagonal(PetscInt n_,PetscReal *D,PetscReal *E,PetscReal *w,PetscScalar *V)
146: {
148: PetscReal abstol = 0.0,vl,vu,*work;
149: PetscBLASInt il,iu,m,*isuppz,n,lwork,*iwork,liwork,info;
150: const char *jobz;
151: #if defined(PETSC_USE_COMPLEX)
152: PetscInt i,j;
153: PetscReal *VV=NULL;
154: #endif
157: PetscBLASIntCast(n_,&n);
158: PetscBLASIntCast(20*n_,&lwork);
159: PetscBLASIntCast(10*n_,&liwork);
160: if (V) {
161: jobz = "V";
162: #if defined(PETSC_USE_COMPLEX)
163: PetscMalloc1(n*n,&VV);
164: #endif
165: } else jobz = "N";
166: PetscMalloc3(2*n,&isuppz,lwork,&work,liwork,&iwork);
167: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
168: #if defined(PETSC_USE_COMPLEX)
169: PetscStackCallBLAS("LAPACKstevr",LAPACKstevr_(jobz,"A",&n,D,E,&vl,&vu,&il,&iu,&abstol,&m,w,VV,&n,isuppz,work,&lwork,iwork,&liwork,&info));
170: #else
171: PetscStackCallBLAS("LAPACKstevr",LAPACKstevr_(jobz,"A",&n,D,E,&vl,&vu,&il,&iu,&abstol,&m,w,V,&n,isuppz,work,&lwork,iwork,&liwork,&info));
172: #endif
173: PetscFPTrapPop();
174: SlepcCheckLapackInfo("stevr",info);
175: #if defined(PETSC_USE_COMPLEX)
176: if (V) {
177: for (i=0;i<n;i++)
178: for (j=0;j<n;j++)
179: V[i*n+j] = VV[i*n+j];
180: PetscFree(VV);
181: }
182: #endif
183: PetscFree3(isuppz,work,iwork);
184: return(0);
185: }
187: /*
188: EPSSelectiveLanczos - Selective reorthogonalization.
189: */
190: static PetscErrorCode EPSSelectiveLanczos(EPS eps,PetscReal *alpha,PetscReal *beta,PetscInt k,PetscInt *M,PetscBool *breakdown,PetscReal anorm)
191: {
193: EPS_LANCZOS *lanczos = (EPS_LANCZOS*)eps->data;
194: PetscInt i,j,m = *M,n,nritz=0,nritzo;
195: Vec vj1,av;
196: Mat Op;
197: PetscReal *d,*e,*ritz,norm;
198: PetscScalar *Y,*hwork;
199: PetscBool *which;
202: PetscCalloc6(m+1,&d,m,&e,m,&ritz,m*m,&Y,m,&which,m,&hwork);
203: for (i=0;i<k;i++) which[i] = PETSC_TRUE;
204: STGetOperator(eps->st,&Op);
206: for (j=k;j<m;j++) {
207: BVSetActiveColumns(eps->V,0,m);
209: /* Lanczos step */
210: BVMatMultColumn(eps->V,Op,j);
211: which[j] = PETSC_TRUE;
212: if (j-2>=k) which[j-2] = PETSC_FALSE;
213: BVOrthogonalizeSomeColumn(eps->V,j+1,which,hwork,&norm,breakdown);
214: alpha[j] = PetscRealPart(hwork[j]);
215: beta[j] = norm;
216: if (*breakdown) {
217: *M = j+1;
218: break;
219: }
221: /* Compute eigenvalues and eigenvectors Y of the tridiagonal block */
222: n = j-k+1;
223: for (i=0;i<n;i++) {
224: d[i] = alpha[i+k];
225: e[i] = beta[i+k];
226: }
227: DenseTridiagonal(n,d,e,ritz,Y);
229: /* Estimate ||A|| */
230: for (i=0;i<n;i++)
231: if (PetscAbsReal(ritz[i]) > anorm) anorm = PetscAbsReal(ritz[i]);
233: /* Compute nearly converged Ritz vectors */
234: nritzo = 0;
235: for (i=0;i<n;i++) {
236: if (norm*PetscAbsScalar(Y[i*n+n-1]) < PETSC_SQRT_MACHINE_EPSILON*anorm) nritzo++;
237: }
238: if (nritzo>nritz) {
239: nritz = 0;
240: for (i=0;i<n;i++) {
241: if (norm*PetscAbsScalar(Y[i*n+n-1]) < PETSC_SQRT_MACHINE_EPSILON*anorm) {
242: BVSetActiveColumns(eps->V,k,k+n);
243: BVGetColumn(lanczos->AV,nritz,&av);
244: BVMultVec(eps->V,1.0,0.0,av,Y+i*n);
245: BVRestoreColumn(lanczos->AV,nritz,&av);
246: nritz++;
247: }
248: }
249: }
250: if (nritz > 0) {
251: BVGetColumn(eps->V,j+1,&vj1);
252: BVSetActiveColumns(lanczos->AV,0,nritz);
253: BVOrthogonalizeVec(lanczos->AV,vj1,hwork,&norm,breakdown);
254: BVRestoreColumn(eps->V,j+1,&vj1);
255: if (*breakdown) {
256: *M = j+1;
257: break;
258: }
259: }
260: BVScaleColumn(eps->V,j+1,1.0/norm);
261: }
263: STRestoreOperator(eps->st,&Op);
264: PetscFree6(d,e,ritz,Y,which,hwork);
265: return(0);
266: }
268: static void update_omega(PetscReal *omega,PetscReal *omega_old,PetscInt j,PetscReal *alpha,PetscReal *beta,PetscReal eps1,PetscReal anorm)
269: {
270: PetscInt k;
271: PetscReal T,binv;
274: /* Estimate of contribution to roundoff errors from A*v
275: fl(A*v) = A*v + f,
276: where ||f|| \approx eps1*||A||.
277: For a full matrix A, a rule-of-thumb estimate is eps1 = sqrt(n)*eps */
278: T = eps1*anorm;
279: binv = 1.0/beta[j+1];
281: /* Update omega(1) using omega(0)==0 */
282: omega_old[0]= beta[1]*omega[1] + (alpha[0]-alpha[j])*omega[0] - beta[j]*omega_old[0];
283: if (omega_old[0] > 0) omega_old[0] = binv*(omega_old[0] + T);
284: else omega_old[0] = binv*(omega_old[0] - T);
286: /* Update remaining components */
287: for (k=1;k<j-1;k++) {
288: 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];
289: if (omega_old[k] > 0) omega_old[k] = binv*(omega_old[k] + T);
290: else omega_old[k] = binv*(omega_old[k] - T);
291: }
292: omega_old[j-1] = binv*T;
294: /* Swap omega and omega_old */
295: for (k=0;k<j;k++) {
296: omega[k] = omega_old[k];
297: omega_old[k] = omega[k];
298: }
299: omega[j] = eps1;
300: PetscFunctionReturnVoid();
301: }
303: static void compute_int(PetscBool *which,PetscReal *mu,PetscInt j,PetscReal delta,PetscReal eta)
304: {
305: PetscInt i,k,maxpos;
306: PetscReal max;
307: PetscBool found;
310: /* initialize which */
311: found = PETSC_FALSE;
312: maxpos = 0;
313: max = 0.0;
314: for (i=0;i<j;i++) {
315: if (PetscAbsReal(mu[i]) >= delta) {
316: which[i] = PETSC_TRUE;
317: found = PETSC_TRUE;
318: } else which[i] = PETSC_FALSE;
319: if (PetscAbsReal(mu[i]) > max) {
320: maxpos = i;
321: max = PetscAbsReal(mu[i]);
322: }
323: }
324: if (!found) which[maxpos] = PETSC_TRUE;
326: for (i=0;i<j;i++) {
327: if (which[i]) {
328: /* find left interval */
329: for (k=i;k>=0;k--) {
330: if (PetscAbsReal(mu[k])<eta || which[k]) break;
331: else which[k] = PETSC_TRUE;
332: }
333: /* find right interval */
334: for (k=i+1;k<j;k++) {
335: if (PetscAbsReal(mu[k])<eta || which[k]) break;
336: else which[k] = PETSC_TRUE;
337: }
338: }
339: }
340: PetscFunctionReturnVoid();
341: }
343: /*
344: EPSPartialLanczos - Partial reorthogonalization.
345: */
346: static PetscErrorCode EPSPartialLanczos(EPS eps,PetscReal *alpha,PetscReal *beta,PetscInt k,PetscInt *M,PetscBool *breakdown,PetscReal anorm)
347: {
349: EPS_LANCZOS *lanczos = (EPS_LANCZOS*)eps->data;
350: PetscInt i,j,m = *M;
351: Mat Op;
352: PetscReal norm,*omega,lomega[100],*omega_old,lomega_old[100],eps1,delta,eta;
353: PetscBool *which,lwhich[100],*which2,lwhich2[100];
354: PetscBool reorth = PETSC_FALSE,force_reorth = PETSC_FALSE;
355: PetscBool fro = PETSC_FALSE,estimate_anorm = PETSC_FALSE;
356: PetscScalar *hwork,lhwork[100];
359: if (m>100) {
360: PetscMalloc5(m,&omega,m,&omega_old,m,&which,m,&which2,m,&hwork);
361: } else {
362: omega = lomega;
363: omega_old = lomega_old;
364: which = lwhich;
365: which2 = lwhich2;
366: hwork = lhwork;
367: }
369: eps1 = PetscSqrtReal((PetscReal)eps->n)*PETSC_MACHINE_EPSILON/2;
370: delta = PETSC_SQRT_MACHINE_EPSILON/PetscSqrtReal((PetscReal)eps->ncv);
371: eta = PetscPowReal(PETSC_MACHINE_EPSILON,3.0/4.0)/PetscSqrtReal((PetscReal)eps->ncv);
372: if (anorm < 0.0) {
373: anorm = 1.0;
374: estimate_anorm = PETSC_TRUE;
375: }
376: for (i=0;i<PetscMax(100,m);i++) omega[i] = omega_old[i] = 0.0;
377: for (i=0;i<k;i++) which[i] = PETSC_TRUE;
379: BVSetActiveColumns(eps->V,0,m);
380: STGetOperator(eps->st,&Op);
381: for (j=k;j<m;j++) {
382: BVMatMultColumn(eps->V,Op,j);
383: if (fro) {
384: /* Lanczos step with full reorthogonalization */
385: BVOrthogonalizeColumn(eps->V,j+1,hwork,&norm,breakdown);
386: alpha[j] = PetscRealPart(hwork[j]);
387: } else {
388: /* Lanczos step */
389: which[j] = PETSC_TRUE;
390: if (j-2>=k) which[j-2] = PETSC_FALSE;
391: BVOrthogonalizeSomeColumn(eps->V,j+1,which,hwork,&norm,breakdown);
392: alpha[j] = PetscRealPart(hwork[j]);
393: beta[j] = norm;
395: /* Estimate ||A|| if needed */
396: if (estimate_anorm) {
397: if (j>k) anorm = PetscMax(anorm,PetscAbsReal(alpha[j])+norm+beta[j-1]);
398: else anorm = PetscMax(anorm,PetscAbsReal(alpha[j])+norm);
399: }
401: /* Check if reorthogonalization is needed */
402: reorth = PETSC_FALSE;
403: if (j>k) {
404: update_omega(omega,omega_old,j,alpha,beta-1,eps1,anorm);
405: for (i=0;i<j-k;i++) {
406: if (PetscAbsReal(omega[i]) > delta) reorth = PETSC_TRUE;
407: }
408: }
409: if (reorth || force_reorth) {
410: for (i=0;i<k;i++) which2[i] = PETSC_FALSE;
411: for (i=k;i<=j;i++) which2[i] = PETSC_TRUE;
412: if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_PERIODIC) {
413: /* Periodic reorthogonalization */
414: if (force_reorth) force_reorth = PETSC_FALSE;
415: else force_reorth = PETSC_TRUE;
416: for (i=0;i<j-k;i++) omega[i] = eps1;
417: } else {
418: /* Partial reorthogonalization */
419: if (force_reorth) force_reorth = PETSC_FALSE;
420: else {
421: force_reorth = PETSC_TRUE;
422: compute_int(which2+k,omega,j-k,delta,eta);
423: for (i=0;i<j-k;i++) {
424: if (which2[i+k]) omega[i] = eps1;
425: }
426: }
427: }
428: BVOrthogonalizeSomeColumn(eps->V,j+1,which2,hwork,&norm,breakdown);
429: }
430: }
432: if (*breakdown || norm < eps->n*anorm*PETSC_MACHINE_EPSILON) {
433: *M = j+1;
434: break;
435: }
436: if (!fro && norm*delta < anorm*eps1) {
437: fro = PETSC_TRUE;
438: PetscInfo1(eps,"Switching to full reorthogonalization at iteration %D\n",eps->its);
439: }
440: beta[j] = norm;
441: BVScaleColumn(eps->V,j+1,1.0/norm);
442: }
444: STRestoreOperator(eps->st,&Op);
445: if (m>100) {
446: PetscFree5(omega,omega_old,which,which2,hwork);
447: }
448: return(0);
449: }
451: /*
452: EPSBasicLanczos - Computes an m-step Lanczos factorization. The first k
453: columns are assumed to be locked and therefore they are not modified. On
454: exit, the following relation is satisfied:
456: OP * V - V * T = f * e_m^T
458: where the columns of V are the Lanczos vectors, T is a tridiagonal matrix,
459: f is the residual vector and e_m is the m-th vector of the canonical basis.
460: The Lanczos vectors (together with vector f) are B-orthogonal (to working
461: accuracy) if full reorthogonalization is being used, otherwise they are
462: (B-)semi-orthogonal. On exit, beta contains the B-norm of f and the next
463: Lanczos vector can be computed as v_{m+1} = f / beta.
465: This function simply calls another function which depends on the selected
466: reorthogonalization strategy.
467: */
468: static PetscErrorCode EPSBasicLanczos(EPS eps,PetscReal *alpha,PetscReal *beta,PetscInt k,PetscInt *m,PetscBool *breakdown,PetscReal anorm)
469: {
470: PetscErrorCode ierr;
471: EPS_LANCZOS *lanczos = (EPS_LANCZOS*)eps->data;
472: PetscScalar *T;
473: PetscInt i,n=*m;
474: PetscReal betam;
475: BVOrthogRefineType orthog_ref;
476: Mat Op;
479: switch (lanczos->reorthog) {
480: case EPS_LANCZOS_REORTHOG_LOCAL:
481: EPSLocalLanczos(eps,alpha,beta,k,m,breakdown);
482: break;
483: case EPS_LANCZOS_REORTHOG_FULL:
484: STGetOperator(eps->st,&Op);
485: BVMatLanczos(eps->V,Op,alpha,beta,k,m,breakdown);
486: STRestoreOperator(eps->st,&Op);
487: break;
488: case EPS_LANCZOS_REORTHOG_SELECTIVE:
489: EPSSelectiveLanczos(eps,alpha,beta,k,m,breakdown,anorm);
490: break;
491: case EPS_LANCZOS_REORTHOG_PERIODIC:
492: case EPS_LANCZOS_REORTHOG_PARTIAL:
493: EPSPartialLanczos(eps,alpha,beta,k,m,breakdown,anorm);
494: break;
495: case EPS_LANCZOS_REORTHOG_DELAYED:
496: PetscMalloc1(n*n,&T);
497: BVGetOrthogonalization(eps->V,NULL,&orthog_ref,NULL,NULL);
498: if (orthog_ref == BV_ORTHOG_REFINE_NEVER) {
499: EPSDelayedArnoldi1(eps,T,n,k,m,&betam,breakdown);
500: } else {
501: EPSDelayedArnoldi(eps,T,n,k,m,&betam,breakdown);
502: }
503: for (i=k;i<n-1;i++) {
504: alpha[i] = PetscRealPart(T[n*i+i]);
505: beta[i] = PetscRealPart(T[n*i+i+1]);
506: }
507: alpha[n-1] = PetscRealPart(T[n*(n-1)+n-1]);
508: beta[n-1] = betam;
509: PetscFree(T);
510: break;
511: }
512: return(0);
513: }
515: PetscErrorCode EPSSolve_Lanczos(EPS eps)
516: {
517: EPS_LANCZOS *lanczos = (EPS_LANCZOS*)eps->data;
519: PetscInt nconv,i,j,k,l,x,n,*perm,restart,ncv=eps->ncv,r,ld;
520: Vec vi,vj,w;
521: Mat U;
522: PetscScalar *Y,*ritz,stmp;
523: PetscReal *d,*e,*bnd,anorm,beta,norm,rtmp,resnorm;
524: PetscBool breakdown;
525: char *conv,ctmp;
528: DSGetLeadingDimension(eps->ds,&ld);
529: PetscMalloc4(ncv,&ritz,ncv,&bnd,ncv,&perm,ncv,&conv);
531: /* The first Lanczos vector is the normalized initial vector */
532: EPSGetStartVector(eps,0,NULL);
534: anorm = -1.0;
535: nconv = 0;
537: /* Restart loop */
538: while (eps->reason == EPS_CONVERGED_ITERATING) {
539: eps->its++;
541: /* Compute an ncv-step Lanczos factorization */
542: n = PetscMin(nconv+eps->mpd,ncv);
543: DSGetArrayReal(eps->ds,DS_MAT_T,&d);
544: e = d + ld;
545: EPSBasicLanczos(eps,d,e,nconv,&n,&breakdown,anorm);
546: beta = e[n-1];
547: DSRestoreArrayReal(eps->ds,DS_MAT_T,&d);
548: DSSetDimensions(eps->ds,n,0,nconv,0);
549: DSSetState(eps->ds,DS_STATE_INTERMEDIATE);
550: BVSetActiveColumns(eps->V,nconv,n);
552: /* Solve projected problem */
553: DSSolve(eps->ds,ritz,NULL);
554: DSSort(eps->ds,ritz,NULL,NULL,NULL,NULL);
555: DSSynchronize(eps->ds,ritz,NULL);
557: /* Estimate ||A|| */
558: for (i=nconv;i<n;i++)
559: anorm = PetscMax(anorm,PetscAbsReal(PetscRealPart(ritz[i])));
561: /* Compute residual norm estimates as beta*abs(Y(m,:)) + eps*||A|| */
562: DSGetArray(eps->ds,DS_MAT_Q,&Y);
563: for (i=nconv;i<n;i++) {
564: resnorm = beta*PetscAbsScalar(Y[n-1+i*ld]) + PETSC_MACHINE_EPSILON*anorm;
565: (*eps->converged)(eps,ritz[i],eps->eigi[i],resnorm,&bnd[i],eps->convergedctx);
566: if (bnd[i]<eps->tol) conv[i] = 'C';
567: else conv[i] = 'N';
568: }
569: DSRestoreArray(eps->ds,DS_MAT_Q,&Y);
571: /* purge repeated ritz values */
572: if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_LOCAL) {
573: for (i=nconv+1;i<n;i++) {
574: if (conv[i] == 'C' && PetscAbsScalar((ritz[i]-ritz[i-1])/ritz[i]) < eps->tol) conv[i] = 'R';
575: }
576: }
578: /* Compute restart vector */
579: if (breakdown) {
580: PetscInfo2(eps,"Breakdown in Lanczos method (it=%D norm=%g)\n",eps->its,(double)beta);
581: } else {
582: restart = nconv;
583: while (restart<n && conv[restart] != 'N') restart++;
584: if (restart >= n) {
585: breakdown = PETSC_TRUE;
586: } else {
587: for (i=restart+1;i<n;i++) {
588: if (conv[i] == 'N') {
589: SlepcSCCompare(eps->sc,ritz[restart],0.0,ritz[i],0.0,&r);
590: if (r>0) restart = i;
591: }
592: }
593: DSGetArray(eps->ds,DS_MAT_Q,&Y);
594: BVMultColumn(eps->V,1.0,0.0,n,Y+restart*ld+nconv);
595: DSRestoreArray(eps->ds,DS_MAT_Q,&Y);
596: }
597: }
599: /* Count and put converged eigenvalues first */
600: for (i=nconv;i<n;i++) perm[i] = i;
601: for (k=nconv;k<n;k++) {
602: if (conv[perm[k]] != 'C') {
603: j = k + 1;
604: while (j<n && conv[perm[j]] != 'C') j++;
605: if (j>=n) break;
606: l = perm[k]; perm[k] = perm[j]; perm[j] = l;
607: }
608: }
610: /* Sort eigenvectors according to permutation */
611: DSGetArray(eps->ds,DS_MAT_Q,&Y);
612: for (i=nconv;i<k;i++) {
613: x = perm[i];
614: if (x != i) {
615: j = i + 1;
616: while (perm[j] != i) j++;
617: /* swap eigenvalues i and j */
618: stmp = ritz[x]; ritz[x] = ritz[i]; ritz[i] = stmp;
619: rtmp = bnd[x]; bnd[x] = bnd[i]; bnd[i] = rtmp;
620: ctmp = conv[x]; conv[x] = conv[i]; conv[i] = ctmp;
621: perm[j] = x; perm[i] = i;
622: /* swap eigenvectors i and j */
623: for (l=0;l<n;l++) {
624: stmp = Y[l+x*ld]; Y[l+x*ld] = Y[l+i*ld]; Y[l+i*ld] = stmp;
625: }
626: }
627: }
628: DSRestoreArray(eps->ds,DS_MAT_Q,&Y);
630: /* compute converged eigenvectors */
631: DSGetMat(eps->ds,DS_MAT_Q,&U);
632: BVMultInPlace(eps->V,U,nconv,k);
633: MatDestroy(&U);
635: /* purge spurious ritz values */
636: if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_LOCAL) {
637: for (i=nconv;i<k;i++) {
638: BVGetColumn(eps->V,i,&vi);
639: VecNorm(vi,NORM_2,&norm);
640: VecScale(vi,1.0/norm);
641: w = eps->work[0];
642: STApply(eps->st,vi,w);
643: VecAXPY(w,-ritz[i],vi);
644: BVRestoreColumn(eps->V,i,&vi);
645: VecNorm(w,NORM_2,&norm);
646: (*eps->converged)(eps,ritz[i],eps->eigi[i],norm,&bnd[i],eps->convergedctx);
647: if (bnd[i]>=eps->tol) conv[i] = 'S';
648: }
649: for (i=nconv;i<k;i++) {
650: if (conv[i] != 'C') {
651: j = i + 1;
652: while (j<k && conv[j] != 'C') j++;
653: if (j>=k) break;
654: /* swap eigenvalues i and j */
655: stmp = ritz[j]; ritz[j] = ritz[i]; ritz[i] = stmp;
656: rtmp = bnd[j]; bnd[j] = bnd[i]; bnd[i] = rtmp;
657: ctmp = conv[j]; conv[j] = conv[i]; conv[i] = ctmp;
658: /* swap eigenvectors i and j */
659: BVGetColumn(eps->V,i,&vi);
660: BVGetColumn(eps->V,j,&vj);
661: VecSwap(vi,vj);
662: BVRestoreColumn(eps->V,i,&vi);
663: BVRestoreColumn(eps->V,j,&vj);
664: }
665: }
666: k = i;
667: }
669: /* store ritz values and estimated errors */
670: for (i=nconv;i<n;i++) {
671: eps->eigr[i] = ritz[i];
672: eps->errest[i] = bnd[i];
673: }
674: nconv = k;
675: EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,n);
676: (*eps->stopping)(eps,eps->its,eps->max_it,nconv,eps->nev,&eps->reason,eps->stoppingctx);
678: if (eps->reason == EPS_CONVERGED_ITERATING) { /* copy restart vector */
679: BVCopyColumn(eps->V,n,nconv);
680: if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_LOCAL && !breakdown) {
681: /* Reorthonormalize restart vector */
682: BVOrthonormalizeColumn(eps->V,nconv,PETSC_FALSE,NULL,&breakdown);
683: }
684: if (breakdown) {
685: /* Use random vector for restarting */
686: PetscInfo(eps,"Using random vector for restart\n");
687: EPSGetStartVector(eps,nconv,&breakdown);
688: }
689: if (breakdown) { /* give up */
690: eps->reason = EPS_DIVERGED_BREAKDOWN;
691: PetscInfo(eps,"Unable to generate more start vectors\n");
692: }
693: }
694: }
695: eps->nconv = nconv;
697: PetscFree4(ritz,bnd,perm,conv);
698: return(0);
699: }
701: PetscErrorCode EPSSetFromOptions_Lanczos(PetscOptionItems *PetscOptionsObject,EPS eps)
702: {
703: PetscErrorCode ierr;
704: EPS_LANCZOS *lanczos = (EPS_LANCZOS*)eps->data;
705: PetscBool flg;
706: EPSLanczosReorthogType reorthog;
709: PetscOptionsHead(PetscOptionsObject,"EPS Lanczos Options");
711: PetscOptionsEnum("-eps_lanczos_reorthog","Lanczos reorthogonalization","EPSLanczosSetReorthog",EPSLanczosReorthogTypes,(PetscEnum)lanczos->reorthog,(PetscEnum*)&reorthog,&flg);
712: if (flg) { EPSLanczosSetReorthog(eps,reorthog); }
714: PetscOptionsTail();
715: return(0);
716: }
718: static PetscErrorCode EPSLanczosSetReorthog_Lanczos(EPS eps,EPSLanczosReorthogType reorthog)
719: {
720: EPS_LANCZOS *lanczos = (EPS_LANCZOS*)eps->data;
723: switch (reorthog) {
724: case EPS_LANCZOS_REORTHOG_LOCAL:
725: case EPS_LANCZOS_REORTHOG_FULL:
726: case EPS_LANCZOS_REORTHOG_DELAYED:
727: case EPS_LANCZOS_REORTHOG_SELECTIVE:
728: case EPS_LANCZOS_REORTHOG_PERIODIC:
729: case EPS_LANCZOS_REORTHOG_PARTIAL:
730: if (lanczos->reorthog != reorthog) {
731: lanczos->reorthog = reorthog;
732: eps->state = EPS_STATE_INITIAL;
733: }
734: break;
735: default:
736: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid reorthogonalization type");
737: }
738: return(0);
739: }
741: /*@
742: EPSLanczosSetReorthog - Sets the type of reorthogonalization used during the Lanczos
743: iteration.
745: Logically Collective on eps
747: Input Parameters:
748: + eps - the eigenproblem solver context
749: - reorthog - the type of reorthogonalization
751: Options Database Key:
752: . -eps_lanczos_reorthog - Sets the reorthogonalization type (either 'local', 'selective',
753: 'periodic', 'partial', 'full' or 'delayed')
755: Level: advanced
757: .seealso: EPSLanczosGetReorthog(), EPSLanczosReorthogType
758: @*/
759: PetscErrorCode EPSLanczosSetReorthog(EPS eps,EPSLanczosReorthogType reorthog)
760: {
766: PetscTryMethod(eps,"EPSLanczosSetReorthog_C",(EPS,EPSLanczosReorthogType),(eps,reorthog));
767: return(0);
768: }
770: static PetscErrorCode EPSLanczosGetReorthog_Lanczos(EPS eps,EPSLanczosReorthogType *reorthog)
771: {
772: EPS_LANCZOS *lanczos = (EPS_LANCZOS*)eps->data;
775: *reorthog = lanczos->reorthog;
776: return(0);
777: }
779: /*@
780: EPSLanczosGetReorthog - Gets the type of reorthogonalization used during
781: the Lanczos iteration.
783: Not Collective
785: Input Parameter:
786: . eps - the eigenproblem solver context
788: Output Parameter:
789: . reorthog - the type of reorthogonalization
791: Level: advanced
793: .seealso: EPSLanczosSetReorthog(), EPSLanczosReorthogType
794: @*/
795: PetscErrorCode EPSLanczosGetReorthog(EPS eps,EPSLanczosReorthogType *reorthog)
796: {
802: PetscUseMethod(eps,"EPSLanczosGetReorthog_C",(EPS,EPSLanczosReorthogType*),(eps,reorthog));
803: return(0);
804: }
806: PetscErrorCode EPSReset_Lanczos(EPS eps)
807: {
809: EPS_LANCZOS *lanczos = (EPS_LANCZOS*)eps->data;
812: BVDestroy(&lanczos->AV);
813: lanczos->allocsize = 0;
814: return(0);
815: }
817: PetscErrorCode EPSDestroy_Lanczos(EPS eps)
818: {
822: PetscFree(eps->data);
823: PetscObjectComposeFunction((PetscObject)eps,"EPSLanczosSetReorthog_C",NULL);
824: PetscObjectComposeFunction((PetscObject)eps,"EPSLanczosGetReorthog_C",NULL);
825: return(0);
826: }
828: PetscErrorCode EPSView_Lanczos(EPS eps,PetscViewer viewer)
829: {
831: EPS_LANCZOS *lanczos = (EPS_LANCZOS*)eps->data;
832: PetscBool isascii;
835: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
836: if (isascii) {
837: PetscViewerASCIIPrintf(viewer," %s reorthogonalization\n",EPSLanczosReorthogTypes[lanczos->reorthog]);
838: }
839: return(0);
840: }
842: SLEPC_EXTERN PetscErrorCode EPSCreate_Lanczos(EPS eps)
843: {
844: EPS_LANCZOS *ctx;
848: PetscNewLog(eps,&ctx);
849: eps->data = (void*)ctx;
851: eps->useds = PETSC_TRUE;
853: eps->ops->solve = EPSSolve_Lanczos;
854: eps->ops->setup = EPSSetUp_Lanczos;
855: eps->ops->setupsort = EPSSetUpSort_Default;
856: eps->ops->setfromoptions = EPSSetFromOptions_Lanczos;
857: eps->ops->destroy = EPSDestroy_Lanczos;
858: eps->ops->reset = EPSReset_Lanczos;
859: eps->ops->view = EPSView_Lanczos;
860: eps->ops->backtransform = EPSBackTransform_Default;
861: eps->ops->computevectors = EPSComputeVectors_Hermitian;
863: PetscObjectComposeFunction((PetscObject)eps,"EPSLanczosSetReorthog_C",EPSLanczosSetReorthog_Lanczos);
864: PetscObjectComposeFunction((PetscObject)eps,"EPSLanczosGetReorthog_C",EPSLanczosGetReorthog_Lanczos);
865: return(0);
866: }