Actual source code: ks-twosided.c
slepc-3.18.2 2023-01-26
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: "krylovschur"
13: Method: Two-sided Arnoldi with Krylov-Schur restart (for left eigenvectors)
15: References:
17: [1] I.N. Zwaan and M.E. Hochstenbach, "Krylov-Schur-type restarts
18: for the two-sided Arnoldi method", SIAM J. Matrix Anal. Appl.
19: 38(2):297-321, 2017.
21: */
23: #include <slepc/private/epsimpl.h>
24: #include "krylovschur.h"
25: #include <slepcblaslapack.h>
27: static PetscErrorCode EPSTwoSidedRQUpdate1(EPS eps,Mat M,PetscInt nv,PetscReal beta,PetscReal betat)
28: {
29: PetscScalar *T,*S,*A,*w;
30: const PetscScalar *pM;
31: Vec u;
32: PetscInt ld,ncv=eps->ncv,i,l,nnv;
33: PetscBLASInt info,n_,ncv_,*p,one=1;
35: DSGetLeadingDimension(eps->ds,&ld);
36: PetscMalloc3(nv,&p,ncv*ncv,&A,ncv,&w);
37: BVGetActiveColumns(eps->V,&l,&nnv);
38: BVSetActiveColumns(eps->V,0,nv);
39: BVSetActiveColumns(eps->W,0,nv);
40: BVGetColumn(eps->V,nv,&u);
41: BVDotVec(eps->W,u,w);
42: BVRestoreColumn(eps->V,nv,&u);
43: MatDenseGetArrayRead(M,&pM);
44: PetscArraycpy(A,pM,ncv*ncv);
45: MatDenseRestoreArrayRead(M,&pM);
46: PetscBLASIntCast(nv,&n_);
47: PetscBLASIntCast(ncv,&ncv_);
48: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
49: PetscCallBLAS("LAPACKgetrf",LAPACKgetrf_(&n_,&n_,A,&ncv_,p,&info));
50: SlepcCheckLapackInfo("getrf",info);
51: PetscLogFlops(2.0*n_*n_*n_/3.0);
52: PetscCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&n_,&one,A,&ncv_,p,w,&ncv_,&info));
53: SlepcCheckLapackInfo("getrs",info);
54: PetscLogFlops(2.0*n_*n_-n_);
55: BVMultColumn(eps->V,-1.0,1.0,nv,w);
56: DSGetArray(eps->ds,DS_MAT_A,&S);
57: for (i=0;i<nv;i++) S[(nv-1)*ld+i] += beta*w[i];
58: DSRestoreArray(eps->ds,DS_MAT_A,&S);
59: BVGetColumn(eps->W,nv,&u);
60: BVDotVec(eps->V,u,w);
61: BVRestoreColumn(eps->W,nv,&u);
62: PetscCallBLAS("LAPACKgetrs",LAPACKgetrs_("C",&n_,&one,A,&ncv_,p,w,&ncv_,&info));
63: PetscFPTrapPop();
64: BVMultColumn(eps->W,-1.0,1.0,nv,w);
65: DSGetArray(eps->ds,DS_MAT_B,&T);
66: for (i=0;i<nv;i++) T[(nv-1)*ld+i] += betat*w[i];
67: DSRestoreArray(eps->ds,DS_MAT_B,&T);
68: PetscFree3(p,A,w);
69: BVSetActiveColumns(eps->V,l,nnv);
70: BVSetActiveColumns(eps->W,l,nnv);
71: return 0;
72: }
74: static PetscErrorCode EPSTwoSidedRQUpdate2(EPS eps,Mat M,PetscInt k)
75: {
76: PetscScalar *Q,*pM,*w,zero=0.0,sone=1.0,*c,*A;
77: PetscBLASInt n_,ncv_,ld_;
78: PetscReal norm;
79: PetscInt l,nv,ncv=eps->ncv,ld,i,j;
81: DSGetLeadingDimension(eps->ds,&ld);
82: BVGetActiveColumns(eps->V,&l,&nv);
83: BVSetActiveColumns(eps->V,0,nv);
84: BVSetActiveColumns(eps->W,0,nv);
85: PetscMalloc2(ncv*ncv,&w,ncv,&c);
86: /* u = u - V*V'*u */
87: BVOrthogonalizeColumn(eps->V,k,c,&norm,NULL);
88: BVScaleColumn(eps->V,k,1.0/norm);
89: DSGetArray(eps->ds,DS_MAT_A,&A);
90: /* H = H + V'*u*b' */
91: for (j=l;j<k;j++) {
92: for (i=0;i<k;i++) A[i+j*ld] += c[i]*A[k+j*ld];
93: A[k+j*ld] *= norm;
94: }
95: DSRestoreArray(eps->ds,DS_MAT_A,&A);
96: BVOrthogonalizeColumn(eps->W,k,c,&norm,NULL);
97: BVScaleColumn(eps->W,k,1.0/norm);
98: DSGetArray(eps->ds,DS_MAT_B,&A);
99: /* H = H + V'*u*b' */
100: for (j=l;j<k;j++) {
101: for (i=0;i<k;i++) A[i+j*ld] += c[i]*A[k+j*ld];
102: A[k+j*ld] *= norm;
103: }
104: DSRestoreArray(eps->ds,DS_MAT_B,&A);
106: /* M = Q'*M*Q */
107: MatDenseGetArray(M,&pM);
108: PetscBLASIntCast(ncv,&ncv_);
109: PetscBLASIntCast(nv,&n_);
110: PetscBLASIntCast(ld,&ld_);
111: DSGetArray(eps->ds,DS_MAT_Q,&Q);
112: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,pM,&ncv_,Q,&ld_,&zero,w,&ncv_));
113: DSRestoreArray(eps->ds,DS_MAT_Q,&Q);
114: DSGetArray(eps->ds,DS_MAT_Z,&Q);
115: PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&n_,&n_,&n_,&sone,Q,&ld_,w,&ncv_,&zero,pM,&ncv_));
116: DSRestoreArray(eps->ds,DS_MAT_Z,&Q);
117: MatDenseRestoreArray(M,&pM);
118: PetscFree2(w,c);
119: BVSetActiveColumns(eps->V,l,nv);
120: BVSetActiveColumns(eps->W,l,nv);
121: return 0;
122: }
124: PetscErrorCode EPSSolve_KrylovSchur_TwoSided(EPS eps)
125: {
126: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
127: Mat M,U,Op,OpHT,S,T;
128: PetscReal norm,norm2,beta,betat;
129: PetscInt ld,l,nv,nvt,k,nconv,dsn,dsk;
130: PetscBool breakdownt,breakdown,breakdownl;
132: DSGetLeadingDimension(eps->ds,&ld);
133: EPSGetStartVector(eps,0,NULL);
134: EPSGetLeftStartVector(eps,0,NULL);
135: l = 0;
136: MatCreateSeqDense(PETSC_COMM_SELF,eps->ncv,eps->ncv,NULL,&M);
138: STGetOperator(eps->st,&Op);
139: MatCreateHermitianTranspose(Op,&OpHT);
141: /* Restart loop */
142: while (eps->reason == EPS_CONVERGED_ITERATING) {
143: eps->its++;
145: /* Compute an nv-step Arnoldi factorization for Op */
146: nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
147: DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l);
148: DSGetMat(eps->ds,DS_MAT_A,&S);
149: BVMatArnoldi(eps->V,Op,S,eps->nconv+l,&nv,&beta,&breakdown);
150: DSRestoreMat(eps->ds,DS_MAT_A,&S);
152: /* Compute an nv-step Arnoldi factorization for Op' */
153: nvt = nv;
154: DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l);
155: DSGetMat(eps->ds,DS_MAT_B,&T);
156: BVMatArnoldi(eps->W,OpHT,T,eps->nconv+l,&nvt,&betat,&breakdownt);
157: DSRestoreMat(eps->ds,DS_MAT_B,&T);
159: /* Make sure both factorizations have the same length */
160: nv = PetscMin(nv,nvt);
161: DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l);
162: if (l==0) DSSetState(eps->ds,DS_STATE_INTERMEDIATE);
163: else DSSetState(eps->ds,DS_STATE_RAW);
164: breakdown = (breakdown || breakdownt)? PETSC_TRUE: PETSC_FALSE;
166: /* Update M, modify Rayleigh quotients S and T */
167: BVSetActiveColumns(eps->V,eps->nconv+l,nv);
168: BVSetActiveColumns(eps->W,eps->nconv+l,nv);
169: BVMatProject(eps->V,NULL,eps->W,M);
171: EPSTwoSidedRQUpdate1(eps,M,nv,beta,betat);
173: /* Solve projected problem */
174: DSSolve(eps->ds,eps->eigr,eps->eigi);
175: DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL);
176: DSSynchronize(eps->ds,eps->eigr,eps->eigi);
177: DSUpdateExtraRow(eps->ds);
179: /* Check convergence */
180: BVNormColumn(eps->V,nv,NORM_2,&norm);
181: BVNormColumn(eps->W,nv,NORM_2,&norm2);
182: EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta*norm,betat*norm2,1.0,&k);
183: (*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx);
184: nconv = k;
186: /* Update l */
187: if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
188: else {
189: l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
190: DSGetTruncateSize(eps->ds,k,nv,&l);
191: }
192: if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
193: if (l) PetscInfo(eps,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l);
195: /* Update the corresponding vectors V(:,idx) = V*Q(:,idx) */
196: BVSetActiveColumns(eps->V,eps->nconv,nv);
197: BVSetActiveColumns(eps->W,eps->nconv,nv);
198: DSGetMat(eps->ds,DS_MAT_Q,&U);
199: BVMultInPlace(eps->V,U,eps->nconv,k+l);
200: DSRestoreMat(eps->ds,DS_MAT_Q,&U);
201: DSGetMat(eps->ds,DS_MAT_Z,&U);
202: BVMultInPlace(eps->W,U,eps->nconv,k+l);
203: DSRestoreMat(eps->ds,DS_MAT_Z,&U);
204: if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
205: BVCopyColumn(eps->V,nv,k+l);
206: BVCopyColumn(eps->W,nv,k+l);
207: }
209: if (eps->reason == EPS_CONVERGED_ITERATING) {
210: if (breakdown || k==nv) {
211: /* Start a new Arnoldi factorization */
212: PetscInfo(eps,"Breakdown in Krylov-Schur method (it=%" PetscInt_FMT " norm=%g)\n",eps->its,(double)beta);
213: if (k<eps->nev) {
214: EPSGetStartVector(eps,k,&breakdown);
215: EPSGetLeftStartVector(eps,k,&breakdownl);
216: if (breakdown || breakdownl) {
217: eps->reason = EPS_DIVERGED_BREAKDOWN;
218: PetscInfo(eps,"Unable to generate more start vectors\n");
219: }
220: }
221: } else {
222: DSGetDimensions(eps->ds,&dsn,NULL,&dsk,NULL);
223: DSSetDimensions(eps->ds,dsn,k,dsk);
224: DSTruncate(eps->ds,k+l,PETSC_FALSE);
225: }
226: EPSTwoSidedRQUpdate2(eps,M,k+l);
227: }
228: eps->nconv = k;
229: EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv);
230: }
232: STRestoreOperator(eps->st,&Op);
233: MatDestroy(&OpHT);
235: DSTruncate(eps->ds,eps->nconv,PETSC_TRUE);
236: MatDestroy(&M);
237: return 0;
238: }