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