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