Actual source code: ks-lrep.c
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: thick-restarted Lanczos for Linear Response eigenvalue problems
15: References:
17: [1] Z. Teng, R.-C. Li, "Convergence analysis of Lanczos-type methods for the
18: linear response eigenvalue problem", J. Comput. Appl. Math. 247, 2013.
20: [2] F. Alvarruiz, B. Mellado-Pinto, J. E. Roman, "Restarted Lanczos methods
21: for the linear response eigenvalue problem", in preparation, 2026.
23: */
24: #include <slepc/private/epsimpl.h>
25: #include "krylovschur.h"
27: static PetscErrorCode Orthog_Teng(Vec x,BV U,BV V,PetscInt j,PetscScalar *h,PetscScalar *c)
28: {
29: PetscInt i;
31: PetscFunctionBegin;
32: PetscCall(BVSetActiveColumns(U,0,j));
33: PetscCall(BVSetActiveColumns(V,0,j));
34: /* c = U^* x */
35: PetscCall(BVDotVec(U,x,c));
36: /* x = x-V*c */
37: PetscCall(BVMultVec(V,-1.0,1.0,x,c));
38: /* accumulate orthog coeffs into h */
39: for (i=0;i<2*j;i++) h[i] += c[i];
40: PetscFunctionReturn(PETSC_SUCCESS);
41: }
43: /* Orthogonalize vector x against first j vectors in U and V
44: v is column j-1 of V */
45: static PetscErrorCode OrthogonalizeVector_Teng(Vec x,BV U,BV V,PetscInt j,Vec u,PetscReal *beta,PetscInt k,PetscScalar *h)
46: {
47: PetscReal alpha;
48: PetscInt i,l;
50: PetscFunctionBegin;
51: PetscCall(PetscArrayzero(h,2*j));
53: /* Local orthogonalization */
54: l = j==k+1?0:j-2; /* 1st column to orthogonalize against */
55: PetscCall(VecDotRealPart(x,u,&alpha));
56: for (i=l;i<j-1;i++) h[i] = beta[i];
57: h[j-1] = alpha;
58: /* x = x-V(:,l:j-1)*h(l:j-1) */
59: PetscCall(BVSetActiveColumns(V,l,j));
60: PetscCall(BVMultVec(V,-1.0,1.0,x,h+l));
62: /* Full orthogonalization */
63: PetscCall(Orthog_Teng(x,U,V,j,h,h+2*j));
64: PetscFunctionReturn(PETSC_SUCCESS);
65: }
67: static PetscErrorCode EPSLREPLanczos_Teng(EPS eps,Mat K,Mat M,BV U,BV V,PetscReal *alpha,PetscReal *beta,PetscInt k,PetscInt *min,PetscBool *breakdown)
68: {
69: PetscInt j,m = *min;
70: Vec u,v,uh,vh;
71: PetscReal beta0;
72: PetscScalar *hwork,lhwork[100],gamma;
73: PetscBool alloc=PETSC_FALSE;
75: PetscFunctionBegin;
76: if (4*m > 100) {
77: PetscCall(PetscMalloc1(4*m,&hwork));
78: alloc = PETSC_TRUE;
79: } else hwork = lhwork;
81: /* Normalize initial vector */
82: if (k==0) {
83: if (eps->nini==0) PetscCall(BVSetRandomColumn(V,0));
84: PetscCall(BVGetColumn(U,0,&u));
85: PetscCall(BVGetColumn(V,0,&v));
86: PetscCall(MatMult(M,v,u));
87: PetscCall(VecDot(u,v,&gamma));
88: beta0 = PetscSqrtReal(PetscRealPart(gamma));
89: if (beta0==0.0) {
90: if (breakdown) *breakdown = PETSC_TRUE;
91: *min = 1; m = 0;
92: } else {
93: PetscCall(VecScale(u,1.0/beta0));
94: PetscCall(VecScale(v,1.0/beta0));
95: }
96: PetscCall(BVRestoreColumn(U,0,&u));
97: PetscCall(BVRestoreColumn(V,0,&v));
98: }
100: for (j=k;j<m;j++) {
101: /* j+1 columns (indices 0 to j) have been computed */
102: PetscCall(BVGetColumn(U,j+1,&uh));
103: PetscCall(BVGetColumn(V,j+1,&vh));
104: PetscCall(BVGetColumn(U,j,&u));
105: PetscCall(MatMult(K,u,vh));
106: PetscCall(OrthogonalizeVector_Teng(vh,U,V,j+1,u,beta,k,hwork));
107: alpha[j] = PetscRealPart(hwork[j]);
108: PetscCall(MatMult(M,vh,uh));
109: PetscCall(VecDot(uh,vh,&gamma));
110: beta[j] = PetscSqrtReal(PetscRealPart(gamma));
111: if (beta[j]==0.0) {
112: if (breakdown) *breakdown = PETSC_TRUE;
113: *min = j+1; m = j;
114: } else {
115: PetscCall(VecScale(uh,1.0/beta[j]));
116: PetscCall(VecScale(vh,1.0/beta[j]));
117: }
118: PetscCall(BVRestoreColumn(U,j+1,&uh));
119: PetscCall(BVRestoreColumn(V,j+1,&vh));
120: PetscCall(BVRestoreColumn(U,j,&u));
121: }
122: if (alloc) PetscCall(PetscFree(hwork));
123: PetscFunctionReturn(PETSC_SUCCESS);
124: }
126: static PetscErrorCode EPSComputeVectors_LREP_Teng(EPS eps)
127: {
128: Mat H;
129: Vec u,v,w;
130: BV U,V;
131: IS is[2];
132: PetscInt k;
133: PetscScalar lambda;
134: PetscBool reduced;
136: PetscFunctionBegin;
137: PetscCall(STGetMatrix(eps->st,0,&H));
138: PetscCall(MatNestGetISs(H,is,NULL));
139: PetscCall(SlepcCheckMatLREPReduced(H,&reduced));
140: PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&V,&U));
141: for (k=0;k<eps->nconv;k++) {
142: PetscCall(BVGetColumn(V,k,&v));
143: /* approx eigenvector is [eigr[k]*v; u] */
144: lambda = eps->eigr[k];
145: PetscCall(STBackTransform(eps->st,1,&lambda,&eps->eigi[k]));
146: PetscCall(VecScale(v,lambda));
147: PetscCall(BVRestoreColumn(V,k,&v));
148: }
149: if (!reduced) {
150: /* the eigenvector [v;u] = J*[y;x] where [y;x] is the reduced eigenvector
151: and J = 1/sqrt(2)[I I; I -I], i.e, v=1/sqrt(2)*(y+x) u=1/sqrt(2)*(y-x) */
152: PetscCall(BVCreateVec(V,&w));
153: for (k=0;k<eps->nconv;k++) {
154: PetscCall(BVGetColumn(U,k,&u));
155: PetscCall(BVGetColumn(V,k,&v));
156: PetscCall(VecCopy(u,w));
157: PetscCall(VecCopy(v,u));
158: PetscCall(VecAXPY(u,-1.0,w));
159: PetscCall(VecAXPY(v,1.0,w));
160: PetscCall(BVRestoreColumn(U,k,&u));
161: PetscCall(BVRestoreColumn(V,k,&v));
162: }
163: PetscCall(VecDestroy(&w));
164: }
165: PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&V,&U));
166: /* Normalize eigenvectors */
167: PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
168: PetscCall(BVNormalize(eps->V,NULL));
169: PetscFunctionReturn(PETSC_SUCCESS);
170: }
172: PetscErrorCode EPSSetUp_KrylovSchur_LREP(EPS eps)
173: {
174: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
175: PetscBool flg;
177: PetscFunctionBegin;
178: PetscCheck((eps->problem_type==EPS_LREP),PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Problem type should be LREP");
179: EPSCheckUnsupportedCondition(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION | EPS_FEATURE_BALANCE,PETSC_TRUE," with LREP structure");
180: PetscCall(EPSSetDimensions_Default(eps,&eps->nev,&eps->ncv,&eps->mpd));
181: PetscCheck(eps->ncv<=eps->nev+eps->mpd,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must not be larger than nev+mpd");
182: if (eps->max_it==PETSC_DETERMINE) eps->max_it = PetscMax(100,2*eps->n/eps->ncv)*((eps->stop==EPS_STOP_THRESHOLD)?10:1);
184: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STSHIFT,&flg));
185: PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Krylov-Schur LREP only supports shift ST");
186: if (!eps->which) eps->which = EPS_SMALLEST_MAGNITUDE;
188: if (!ctx->keep) ctx->keep = 0.5;
189: PetscCall(STSetStructured(eps->st,PETSC_FALSE));
191: PetscCall(EPSAllocateSolution(eps,1));
192: /* Teng */
193: eps->ops->solve = EPSSolve_KrylovSchur_LREP_Teng;
194: eps->ops->computevectors = EPSComputeVectors_LREP_Teng;
195: PetscCall(DSSetType(eps->ds,DSHEP));
196: PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
197: PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
198: PetscCall(DSAllocate(eps->ds,eps->ncv+1));
199: PetscFunctionReturn(PETSC_SUCCESS);
200: }
202: PetscErrorCode EPSSolve_KrylovSchur_LREP_Teng(EPS eps)
203: {
204: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
205: PetscInt i,k,l,ld,nv,nconv=0,nevsave,ma,na,Ma,Na;
206: Mat H,Q,K,M,A,B;
207: BV U,V;
208: IS is[2];
209: PetscReal *a,*b,beta;
210: PetscBool reduced,breakdown=PETSC_FALSE;
211: const PetscScalar scal[] = { 1.0, -1.0 };
213: PetscFunctionBegin;
214: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
216: /* Extract matrix blocks */
217: PetscCall(STGetMatrix(eps->st,0,&H));
218: PetscCall(MatNestGetISs(H,is,NULL));
219: PetscCall(SlepcCheckMatLREPReduced(H,&reduced));
220: if (reduced) {
221: PetscCall(MatNestGetSubMat(H,0,1,&K));
222: PetscCall(MatNestGetSubMat(H,1,0,&M));
223: } else {
224: PetscCall(MatNestGetSubMat(H,0,0,&A));
225: PetscCall(MatNestGetSubMat(H,0,1,&B));
226: PetscCall(MatGetSize(A,&Ma,&Na));
227: PetscCall(MatGetLocalSize(A,&ma,&na));
228: /* K = A-B */
229: PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&K));
230: PetscCall(MatSetSizes(K,ma,na,Ma,Na));
231: PetscCall(MatSetType(K,MATCOMPOSITE));
232: PetscCall(MatCompositeAddMat(K,A));
233: PetscCall(MatCompositeAddMat(K,B));
234: PetscCall(MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY));
235: PetscCall(MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY));
236: PetscCall(MatCompositeSetScalings(K,scal));
237: /* M = A+B */
238: PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&M));
239: PetscCall(MatSetSizes(M,ma,na,Ma,Na));
240: PetscCall(MatSetType(M,MATCOMPOSITE));
241: PetscCall(MatCompositeAddMat(M,A));
242: PetscCall(MatCompositeAddMat(M,B));
243: PetscCall(MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY));
244: PetscCall(MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY));
245: }
247: /* Get the split bases */
248: PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&V,&U));
250: nevsave = eps->nev;
251: eps->nev = (eps->nev+1)/2;
252: l = 0;
254: /* Restart loop */
255: while (eps->reason == EPS_CONVERGED_ITERATING) {
256: eps->its++;
258: /* Compute an nv-step Lanczos factorization */
259: nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
260: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
261: PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&a));
262: b = a + ld;
263: PetscCall(EPSLREPLanczos_Teng(eps,K,M,U,V,a,b,eps->nconv+l,&nv,&breakdown));
264: beta = b[nv-1];
265: PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&a));
266: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
267: PetscCall(DSSetState(eps->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
268: PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
270: /* Solve projected problem */
271: PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
272: PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
273: PetscCall(DSUpdateExtraRow(eps->ds));
274: PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
276: /* Check convergence */
277: for (i=0;i<nv;i++) eps->eigr[i] = PetscSqrtReal(PetscRealPart(eps->eigr[i]));
278: PetscCall(EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta,0.0,1.0,&k));
279: EPSSetCtxThreshold(eps,eps->eigr,eps->eigi,eps->errest,k,nv);
280: PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
281: nconv = k;
283: /* Update l */
284: if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
285: else l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
286: if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
287: if (l) PetscCall(PetscInfo(eps,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
289: if (eps->reason == EPS_CONVERGED_ITERATING) {
290: PetscCheck(!breakdown,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Breakdown in LREP Krylov-Schur (beta=%g)",(double)beta);
291: /* Prepare the Rayleigh quotient for restart */
292: PetscCall(DSTruncate(eps->ds,k+l,PETSC_FALSE));
293: }
294: /* Update the corresponding vectors
295: U(:,idx) = U*Q(:,idx), V(:,idx) = V*Q(:,idx) */
296: PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&Q));
297: PetscCall(BVMultInPlace(U,Q,eps->nconv,k+l));
298: PetscCall(BVMultInPlace(V,Q,eps->nconv,k+l));
299: PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&Q));
301: if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
302: PetscCall(BVCopyColumn(eps->V,nv,k+l));
303: if (eps->stop==EPS_STOP_THRESHOLD && nv-k<5) { /* reallocate */
304: eps->ncv = eps->mpd+k;
305: PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&V,&U));
306: PetscCall(EPSReallocateSolution(eps,eps->ncv+1));
307: PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&V,&U));
308: for (i=nv;i<eps->ncv;i++) eps->perm[i] = i;
309: PetscCall(DSReallocate(eps->ds,eps->ncv+1));
310: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
311: }
312: }
313: eps->nconv = k;
314: PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv));
315: }
317: eps->nev = nevsave;
319: PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
320: PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&V,&U));
321: if (!reduced) {
322: PetscCall(MatDestroy(&K));
323: PetscCall(MatDestroy(&M));
324: }
325: PetscFunctionReturn(PETSC_SUCCESS);
326: }