Actual source code: ks-bse.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 Bethe-Salpeter pseudo-Hermitan matrices
15: References:
17: [1] M. Shao et al, "A structure preserving Lanczos algorithm for computing
18: the optical absorption spectrum", SIAM J. Matrix Anal. App. 39(2), 2018.
20: [2] F. Alvarruiz, B. Mellado-Pinto, J. E. Roman, "Variants of thick-restart
21: Lanczos for the Bethe-Salpeter eigenvalue problem", arXiv:2503.20920,
22: 2025.
24: */
25: #include <slepc/private/epsimpl.h>
26: #include "krylovschur.h"
28: static PetscBool cited = PETSC_FALSE;
29: static const char citation[] =
30: "@Misc{slepc-bse,\n"
31: " author = \"F. Alvarruiz and B. Mellado-Pinto and J. E. Roman\",\n"
32: " title = \"Variants of thick-restart {Lanczos} for the {Bethe--Salpeter} eigenvalue problem\",\n"
33: " eprint = \"2503.20920\",\n"
34: " archivePrefix = \"arXiv\",\n"
35: " primaryClass = \"mathematics.numerical analysis\",\n"
36: " year = \"2025,\"\n"
37: " doi = \"https://doi.org/10.48550/arXiv.2503.20920\"\n"
38: "}\n";
40: static PetscErrorCode Orthog_Shao(Vec x,BV U,BV V,PetscInt j,PetscScalar *h,PetscScalar *c)
41: {
42: PetscInt i;
44: PetscFunctionBegin;
45: PetscCall(BVSetActiveColumns(U,0,j));
46: PetscCall(BVSetActiveColumns(V,0,j));
47: /* c = real(V^* x) ; c2 = imag(U^* x)*1i */
48: #if defined(PETSC_USE_COMPLEX)
49: PetscCall(BVDotVecBegin(V,x,c));
50: PetscCall(BVDotVecBegin(U,x,c+j));
51: PetscCall(BVDotVecEnd(V,x,c));
52: PetscCall(BVDotVecEnd(U,x,c+j));
53: #else
54: PetscCall(BVDotVec(V,x,c));
55: #endif
56: #if defined(PETSC_USE_COMPLEX)
57: for (i=0; i<j; i++) {
58: c[i] = PetscRealPart(c[i]);
59: c[j+i] = PetscCMPLX(0.0,PetscImaginaryPart(c[j+i]));
60: }
61: #endif
62: /* x = x-U*c-V*c2 */
63: PetscCall(BVMultVec(U,-1.0,1.0,x,c));
64: #if defined(PETSC_USE_COMPLEX)
65: PetscCall(BVMultVec(V,-1.0,1.0,x,c+j));
66: #endif
67: /* accumulate orthog coeffs into h */
68: for (i=0; i<2*j; i++) h[i] += c[i];
69: PetscFunctionReturn(PETSC_SUCCESS);
70: }
72: /* Orthogonalize vector x against first j vectors in U and V
73: v is column j-1 of V */
74: static PetscErrorCode OrthogonalizeVector_Shao(Vec x,BV U,BV V,PetscInt j,Vec v,PetscReal *beta,PetscInt k,PetscScalar *h)
75: {
76: PetscReal alpha;
77: PetscInt i,l;
79: PetscFunctionBegin;
80: PetscCall(PetscArrayzero(h,2*j));
82: /* Local orthogonalization */
83: l = j==k+1?0:j-2; /* 1st column to orthogonalize against */
84: PetscCall(VecDotRealPart(x,v,&alpha));
85: for (i=l; i<j-1; i++) h[i] = beta[i];
86: h[j-1] = alpha;
87: /* x = x-U(:,l:j-1)*h(l:j-1) */
88: PetscCall(BVSetActiveColumns(U,l,j));
89: PetscCall(BVMultVec(U,-1.0,1.0,x,h+l));
91: /* Full orthogonalization */
92: PetscCall(Orthog_Shao(x,U,V,j,h,h+2*j));
93: PetscFunctionReturn(PETSC_SUCCESS);
94: }
96: static PetscErrorCode EPSBSELanczos_Shao(EPS eps,BV U,BV V,PetscReal *alpha,PetscReal *beta,PetscInt k,PetscInt *M,PetscBool *breakdown)
97: {
98: PetscInt j,m = *M;
99: Vec v,x,y,w,f,g,vecs[2];
100: Mat H;
101: IS is[2];
102: PetscReal nrm;
103: PetscScalar *hwork,lhwork[100],gamma;
104: PetscBool alloc=PETSC_FALSE;
105: PetscContainer container;
106: SlepcMatStruct mctx;
108: PetscFunctionBegin;
109: if (4*m > 100) {
110: PetscCall(PetscMalloc1(4*m,&hwork));
111: alloc = PETSC_TRUE;
112: } else hwork = lhwork;
113: PetscCall(STGetMatrix(eps->st,0,&H));
114: PetscCall(MatNestGetISs(H,is,NULL));
115: PetscCall(PetscObjectQuery((PetscObject)H,"SlepcMatStruct",(PetscObject*)&container));
116: PetscCall(PetscContainerGetPointer(container,(void**)&mctx));
118: /* create work vectors */
119: PetscCall(BVGetColumn(V,0,&v));
120: PetscCall(VecDuplicate(v,&w));
121: vecs[0] = v;
122: vecs[1] = w;
123: PetscCall(VecCreateNest(PetscObjectComm((PetscObject)eps),2,is,vecs,&f));
124: PetscCall(VecCreateNest(PetscObjectComm((PetscObject)eps),2,is,vecs,&g));
125: PetscCall(BVRestoreColumn(V,0,&v));
127: /* Normalize initial vector */
128: if (k==0) {
129: if (eps->nini==0) PetscCall(BVSetRandomColumn(eps->V,0));
130: PetscCall(BVGetColumn(U,0,&x));
131: PetscCall(BVGetColumn(V,0,&y));
132: PetscCall(VecNestSetSubVec(f,0,x));
133: PetscCall(VecNestSetSubVec(g,0,y));
134: mctx->s = 1.0;
135: PetscCall(STApply(eps->st,f,g));
136: PetscCall(VecDot(y,x,&gamma));
137: nrm = PetscSqrtReal(PetscRealPart(gamma));
138: if (nrm==0.0) {
139: if (breakdown) *breakdown = PETSC_TRUE;
140: *M = 1; m = 0;
141: } else {
142: PetscCall(VecScale(x,1.0/nrm));
143: PetscCall(VecScale(y,1.0/nrm));
144: }
145: PetscCall(BVRestoreColumn(U,0,&x));
146: PetscCall(BVRestoreColumn(V,0,&y));
147: }
149: for (j=k;j<m;j++) {
150: /* j+1 columns (indexes 0 to j) have been computed */
151: PetscCall(BVGetColumn(V,j,&v));
152: PetscCall(BVGetColumn(U,j+1,&x));
153: PetscCall(BVGetColumn(V,j+1,&y));
154: PetscCall(VecNestSetSubVec(f,0,v));
155: PetscCall(VecNestSetSubVec(g,0,x));
156: mctx->s = -1.0;
157: PetscCall(STApply(eps->st,f,g));
158: PetscCall(OrthogonalizeVector_Shao(x,U,V,j+1,v,beta,k,hwork));
159: alpha[j] = PetscRealPart(hwork[j]);
160: PetscCall(VecNestSetSubVec(f,0,x));
161: PetscCall(VecNestSetSubVec(g,0,y));
162: mctx->s = 1.0;
163: PetscCall(STApply(eps->st,f,g));
164: PetscCall(VecDot(x,y,&gamma));
165: beta[j] = PetscSqrtReal(PetscRealPart(gamma));
166: if (beta[j]==0.0) {
167: if (breakdown) *breakdown = PETSC_TRUE;
168: *M = j+1; m = j;
169: } else {
170: PetscCall(VecScale(x,1.0/beta[j]));
171: PetscCall(VecScale(y,1.0/beta[j]));
172: }
173: PetscCall(BVRestoreColumn(V,j,&v));
174: PetscCall(BVRestoreColumn(U,j+1,&x));
175: PetscCall(BVRestoreColumn(V,j+1,&y));
176: }
177: if (alloc) PetscCall(PetscFree(hwork));
178: PetscCall(VecDestroy(&w));
179: PetscCall(VecDestroy(&f));
180: PetscCall(VecDestroy(&g));
181: PetscFunctionReturn(PETSC_SUCCESS);
182: }
184: static PetscErrorCode EPSComputeVectors_BSE_Shao(EPS eps)
185: {
186: Mat H;
187: Vec u1,v1;
188: BV U,V;
189: IS is[2];
190: PetscInt k;
191: PetscScalar lambda;
193: PetscFunctionBegin;
194: PetscCall(STGetMatrix(eps->st,0,&H));
195: PetscCall(MatNestGetISs(H,is,NULL));
196: PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
197: for (k=0; k<eps->nconv; k++) {
198: PetscCall(BVGetColumn(U,k,&u1));
199: PetscCall(BVGetColumn(V,k,&v1));
200: /* approx eigenvector is [ (eigr[k]*u1+v1)]
201: [conj(eigr[k]*u1-v1)] */
202: lambda = eps->eigr[k];
203: PetscCall(STBackTransform(eps->st,1,&lambda,&eps->eigi[k]));
204: PetscCall(VecAYPX(u1,lambda,v1));
205: PetscCall(VecAYPX(v1,-2.0,u1));
206: PetscCall(VecConjugate(v1));
207: PetscCall(BVRestoreColumn(U,k,&u1));
208: PetscCall(BVRestoreColumn(V,k,&v1));
209: }
210: PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
211: /* Normalize eigenvectors */
212: PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
213: PetscCall(BVNormalize(eps->V,NULL));
214: PetscFunctionReturn(PETSC_SUCCESS);
215: }
217: static PetscErrorCode Orthog_Gruning(Vec x,BV U,BV V,BV HU,BV HV,PetscInt j,PetscScalar *h,PetscScalar *c)
218: {
219: PetscInt i;
221: PetscFunctionBegin;
222: PetscCall(BVSetActiveColumns(U,0,j));
223: PetscCall(BVSetActiveColumns(HU,0,j));
224: PetscCall(BVSetActiveColumns(V,0,j));
225: PetscCall(BVSetActiveColumns(HV,0,j));
226: #if defined(PETSC_USE_COMPLEX)
227: PetscCall(BVDotVecBegin(HU,x,c));
228: PetscCall(BVDotVecBegin(HV,x,c+j));
229: PetscCall(BVDotVecEnd(HU,x,c));
230: PetscCall(BVDotVecEnd(HV,x,c+j));
231: #else
232: PetscCall(BVDotVec(HU,x,c));
233: #endif
234: for (i=0; i<j; i++) {
235: /* c1 = 2*real(HU^* x) ; c2 = 2*imag(HV^* x)*1i */
236: #if defined(PETSC_USE_COMPLEX)
237: c[i] = PetscRealPart(c[i]);
238: c[j+i] = PetscCMPLX(0.0,PetscImaginaryPart(c[j+i]));
239: #else
240: c[j+i] = 0.0;
241: #endif
242: }
243: /* x = x-U*c1-V*c2 */
244: #if defined(PETSC_USE_COMPLEX)
245: PetscCall(BVMultVec(U,-2.0,1.0,x,c));
246: PetscCall(BVMultVec(V,-2.0,1.0,x,c+j));
247: #else
248: PetscCall(BVMultVec(U,-2.0,1.0,x,c));
249: #endif
250: /* accumulate orthog coeffs into h */
251: for (i=0; i<2*j; i++) h[i] += 2*c[i];
252: PetscFunctionReturn(PETSC_SUCCESS);
253: }
255: /* Orthogonalize vector x against first j vectors in U and V */
256: static PetscErrorCode OrthogonalizeVector_Gruning(Vec x,BV U,BV V,BV HU,BV HV,PetscInt j,PetscReal *beta,PetscInt k,PetscScalar *h,PetscBool s)
257: {
258: PetscInt l,i;
259: Vec u;
261: PetscFunctionBegin;
262: PetscCall(PetscArrayzero(h,4*j));
264: if (s) {
265: /* Local orthogonalization */
266: PetscCall(BVGetColumn(U,j-1,&u));
267: PetscCall(VecAXPY(x,-*beta,u));
268: PetscCall(BVRestoreColumn(U,j-1,&u));
269: h[j-1] = *beta;
270: /* Full orthogonalization */
271: PetscCall(Orthog_Gruning(x,U,V,HU,HV,j,h,h+2*j));
272: } else {
273: /* Local orthogonalization */
274: l = j==k+1?0:j-2; /* 1st column to orthogonalize against */
275: for (i=l; i<j-1; i++) h[j+i] = beta[i];
276: /* x = x-V(:,l:j-2)*h(l:j-2) */
277: PetscCall(BVSetActiveColumns(V,l,j-1));
278: PetscCall(BVMultVec(V,-1.0,1.0,x,h+j+l));
279: }
280: PetscFunctionReturn(PETSC_SUCCESS);
281: }
283: static PetscErrorCode EPSBSELanczos_Gruning(EPS eps,BV U,BV V,BV HU,BV HV,PetscReal *beta1,PetscReal *beta2,PetscInt k,PetscInt *M,PetscBool *breakdown)
284: {
285: PetscInt j,m = *M;
286: Vec v,x,y,w,f,g,vecs[2];
287: Mat H;
288: IS is[2];
289: PetscReal nrm;
290: PetscScalar *hwork,lhwork[100],dot;
291: PetscBool alloc=PETSC_FALSE;
292: PetscContainer container;
293: SlepcMatStruct mctx;
295: PetscFunctionBegin;
296: if (4*m > 100) {
297: PetscCall(PetscMalloc1(4*m,&hwork));
298: alloc = PETSC_TRUE;
299: } else hwork = lhwork;
300: PetscCall(STGetMatrix(eps->st,0,&H));
301: PetscCall(MatNestGetISs(H,is,NULL));
302: PetscCall(PetscObjectQuery((PetscObject)H,"SlepcMatStruct",(PetscObject*)&container));
303: PetscCall(PetscContainerGetPointer(container,(void**)&mctx));
305: /* create work vectors */
306: PetscCall(BVGetColumn(V,0,&v));
307: PetscCall(VecDuplicate(v,&w));
308: vecs[0] = v;
309: vecs[1] = w;
310: PetscCall(VecCreateNest(PetscObjectComm((PetscObject)eps),2,is,vecs,&f));
311: PetscCall(VecCreateNest(PetscObjectComm((PetscObject)eps),2,is,vecs,&g));
312: PetscCall(BVRestoreColumn(V,0,&v));
314: /* Normalize initial vector */
315: if (k==0) {
316: if (eps->nini==0) PetscCall(BVSetRandomColumn(eps->V,0));
317: /* y = Hmult(v1,1) */
318: PetscCall(BVGetColumn(U,k,&x));
319: PetscCall(BVGetColumn(HU,k,&y));
320: PetscCall(VecNestSetSubVec(f,0,x));
321: PetscCall(VecNestSetSubVec(g,0,y));
322: mctx->s = 1.0;
323: PetscCall(STApply(eps->st,f,g));
324: /* nrm = sqrt(2*real(u1'*y)); */
325: PetscCall(VecDot(x,y,&dot));
326: nrm = PetscSqrtReal(PetscRealPart(2*dot));
327: if (nrm==0.0) {
328: if (breakdown) *breakdown = PETSC_TRUE;
329: *M = 1; m = 0;
330: } else {
331: /* U(:,j) = u1/nrm; */
332: /* HU(:,j) = y/nrm; */
333: PetscCall(VecScale(x,1.0/nrm));
334: PetscCall(VecScale(y,1.0/nrm));
335: }
336: PetscCall(BVRestoreColumn(U,k,&x));
337: PetscCall(BVRestoreColumn(HU,k,&y));
338: }
340: for (j=k;j<m;j++) {
341: /* j+1 columns (indexes 0 to j) have been computed */
342: PetscCall(BVGetColumn(HU,j,&x));
343: PetscCall(BVGetColumn(V,j,&v));
344: PetscCall(BVGetColumn(HV,j,&y));
345: PetscCall(VecCopy(x,v));
346: PetscCall(BVRestoreColumn(HU,j,&x));
347: /* v = Orthogonalize HU(:,j) */
348: PetscCall(OrthogonalizeVector_Gruning(v,U,V,HU,HV,j+1,beta2,k,hwork,PETSC_FALSE));
349: /* y = Hmult(v,-1) */
350: PetscCall(VecNestSetSubVec(f,0,v));
351: PetscCall(VecNestSetSubVec(g,0,y));
352: mctx->s = -1.0;
353: PetscCall(STApply(eps->st,f,g));
354: /* beta = sqrt(2*real(v'*y)); */
355: PetscCall(VecDot(v,y,&dot));
356: beta1[j] = PetscSqrtReal(PetscRealPart(2*dot));
357: if (beta1[j]==0.0) {
358: if (breakdown) *breakdown = PETSC_TRUE;
359: *M = j+1; m = j;
360: } else {
361: /* V(:,j) = v/beta1; */
362: /* HV(:,j) = y/beta1; */
363: PetscCall(VecScale(v,1.0/beta1[j]));
364: PetscCall(VecScale(y,1.0/beta1[j]));
365: }
366: PetscCall(BVRestoreColumn(V,j,&v));
367: PetscCall(BVRestoreColumn(HV,j,&y));
368: if (breakdown && *breakdown) continue;
370: PetscCall(BVGetColumn(HV,j,&x));
371: PetscCall(BVGetColumn(U,j+1,&v));
372: PetscCall(BVGetColumn(HU,j+1,&y));
373: PetscCall(VecCopy(x,v));
374: PetscCall(BVRestoreColumn(HV,j,&x));
375: /* v = Orthogonalize HV(:,j) */
376: PetscCall(OrthogonalizeVector_Gruning(v,U,V,HU,HV,j+1,&beta1[j],k,hwork,PETSC_TRUE));
377: /* y = Hmult(v,1) */
378: PetscCall(VecNestSetSubVec(f,0,v));
379: PetscCall(VecNestSetSubVec(g,0,y));
380: mctx->s = 1.0;
381: PetscCall(STApply(eps->st,f,g));
382: /* beta = sqrt(2*real(v'*y)); */
383: PetscCall(VecDot(v,y,&dot));
384: beta2[j] = PetscSqrtReal(PetscRealPart(2*dot));
385: if (beta2[j]==0.0) {
386: if (breakdown) *breakdown = PETSC_TRUE;
387: *M = j+1; m = j;
388: } else {
389: /* U(:,j) = v/beta2; */
390: /* HU(:,j) = y/beta2; */
391: PetscCall(VecScale(v,1.0/beta2[j]));
392: PetscCall(VecScale(y,1.0/beta2[j]));
393: }
394: PetscCall(BVRestoreColumn(U,j+1,&v));
395: PetscCall(BVRestoreColumn(HU,j+1,&y));
396: }
397: if (alloc) PetscCall(PetscFree(hwork));
398: PetscCall(VecDestroy(&w));
399: PetscCall(VecDestroy(&f));
400: PetscCall(VecDestroy(&g));
401: PetscFunctionReturn(PETSC_SUCCESS);
402: }
404: static PetscErrorCode EPSComputeVectors_BSE_Gruning(EPS eps)
405: {
406: Mat H;
407: Vec u1,v1;
408: BV U,V;
409: IS is[2];
410: PetscInt k;
412: PetscFunctionBegin;
413: PetscCall(STGetMatrix(eps->st,0,&H));
414: PetscCall(MatNestGetISs(H,is,NULL));
415: PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
416: /* approx eigenvector [x1] is [ u1+v1 ]
417: [x2] [conj(u1)-conj(v1)] */
418: for (k=0; k<eps->nconv; k++) {
419: PetscCall(BVGetColumn(U,k,&u1));
420: PetscCall(BVGetColumn(V,k,&v1));
421: /* x1 = u1 + v1 */
422: PetscCall(VecAXPY(u1,1.0,v1));
423: /* x2 = conj(u1) - conj(v1) = conj(u1 - v1) = conj((u1 + v1) - 2*v1) */
424: PetscCall(VecAYPX(v1,-2.0,u1));
425: PetscCall(VecConjugate(v1));
426: PetscCall(BVRestoreColumn(U,k,&u1));
427: PetscCall(BVRestoreColumn(V,k,&v1));
428: }
429: PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
430: /* Normalize eigenvectors */
431: PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
432: PetscCall(BVNormalize(eps->V,NULL));
433: PetscFunctionReturn(PETSC_SUCCESS);
434: }
436: /* Full orthogonalization of vector [hx, conj(hx)] against first j vectors in X and Y */
437: static PetscErrorCode Orthog_ProjectedBSE(Vec hx,BV X,BV Y,PetscInt j,PetscScalar *h,PetscScalar *c)
438: {
439: PetscInt i;
441: PetscFunctionBegin;
442: PetscCall(BVSetActiveColumns(X,0,j));
443: PetscCall(BVSetActiveColumns(Y,0,j));
444: /* c1 = X^* hx */
445: PetscCall(BVDotVec(X,hx,c));
446: /* c2 = Y^* conj(hx) */
447: PetscCall(VecConjugate(hx));
448: PetscCall(BVDotVec(Y,hx,c+j));
449: /* c = c1 - c2 */
450: for (i=0;i<j;i++) c[i] -= c[i+j];
451: /* hx = hx - conj(Y*c) */
452: PetscCall(BVMultVec(Y,-1.0,1.0,hx,c));
453: PetscCall(VecConjugate(hx));
454: /* hx = hx - X*c */
455: PetscCall(BVMultVec(X,-1.0,1.0,hx,c));
456: /* accumulate orthog coeffs into h */
457: for (i=0;i<j;i++) h[i] += c[i];
458: PetscFunctionReturn(PETSC_SUCCESS);
459: }
461: /* Orthogonalize vector [hx; hy] against first j vectors in X and Y
462: The result is a vector [u; conj(u)]. Vector hx is overwritten with u. */
463: static PetscErrorCode OrthogonalizeVector_ProjectedBSE(Vec hx,Vec hy,BV X,BV Y,PetscInt j,PetscReal *beta,PetscInt k,PetscScalar *h)
464: {
465: PetscInt l,i;
466: PetscScalar alpha,alpha1,alpha2;
467: Vec x,y;
469: PetscFunctionBegin;
470: PetscCall(PetscArrayzero(h,2*j));
472: /* Local orthogonalization */
473: l = j==k+1?0:j-2; /* 1st column to orthogonalize against */
474: /* alpha = X(:,j-1)'*hx-Y(:,j-1)'*hy */
475: PetscCall(BVGetColumn(X,j-1,&x));
476: PetscCall(BVGetColumn(Y,j-1,&y));
477: PetscCall(VecDotBegin(hx,x,&alpha1));
478: PetscCall(VecDotBegin(hy,y,&alpha2));
479: PetscCall(VecDotEnd(hx,x,&alpha1));
480: PetscCall(VecDotEnd(hy,y,&alpha2));
481: PetscCall(BVRestoreColumn(X,j-1,&x));
482: PetscCall(BVRestoreColumn(Y,j-1,&y));
483: alpha = alpha1-alpha2;
484: /* Store coeffs into h */
485: for (i=l; i<j-1; i++) h[i] = h[j+i] = beta[i]/2.0;
486: h[j-1] = alpha;
488: /* Orthogonalize: hx = hx - X(:,l:j-1)*h1 - conj(Y(:,l:j-1))*h2 */
489: /* hx = hx - X(:,l:j-1)*h1 */
490: PetscCall(BVSetActiveColumns(X,l,j));
491: PetscCall(BVSetActiveColumns(Y,l,j));
492: PetscCall(BVMultVec(X,-1.0,1.0,hx,h+l));
493: /* hx = conj(hx) */
494: PetscCall(VecConjugate(hx));
495: /* hx = hx - Y(:,l:j-1)*conj(h2) */
496: h[2*j-1] = PetscConj(alpha-1.0);
497: PetscCall(BVMultVec(Y,-1.0,1.0,hx,h+j+l));
498: h[2*j-1] = alpha-1.0;
499: /* hx = conj(hx) */
500: PetscCall(VecConjugate(hx));
502: /* Full orthogonalization */
503: PetscCall(Orthog_ProjectedBSE(hx,X,Y,j,h,h+2*j));
504: PetscFunctionReturn(PETSC_SUCCESS);
505: }
507: static PetscErrorCode EPSBSELanczos_ProjectedBSE(EPS eps,BV X,BV Y,Vec v,PetscReal *alpha,PetscReal *beta,PetscInt k,PetscInt *M,PetscBool *breakdown)
508: {
509: PetscInt j,m = *M;
510: Vec u,x,y,w,f,g,vecs[2];
511: Mat H;
512: IS is[2];
513: PetscReal nrm;
514: PetscScalar *hwork,lhwork[100],gamma;
515: PetscBool alloc=PETSC_FALSE;
516: PetscContainer container;
517: SlepcMatStruct mctx;
519: PetscFunctionBegin;
520: if (4*m > 100) {
521: PetscCall(PetscMalloc1(4*m,&hwork));
522: alloc = PETSC_TRUE;
523: } else hwork = lhwork;
524: PetscCall(STGetMatrix(eps->st,0,&H));
525: PetscCall(MatNestGetISs(H,is,NULL));
526: PetscCall(PetscObjectQuery((PetscObject)H,"SlepcMatStruct",(PetscObject*)&container));
527: PetscCall(PetscContainerGetPointer(container,(void**)&mctx));
529: /* create work vectors */
530: PetscCall(BVGetColumn(Y,0,&u));
531: PetscCall(VecDuplicate(u,&w));
532: vecs[0] = u;
533: vecs[1] = w;
534: PetscCall(VecCreateNest(PetscObjectComm((PetscObject)eps),2,is,vecs,&f));
535: PetscCall(VecCreateNest(PetscObjectComm((PetscObject)eps),2,is,vecs,&g));
536: PetscCall(BVRestoreColumn(Y,0,&u));
538: /* Normalize initial vector */
539: if (k==0) {
540: if (eps->nini==0) PetscCall(BVSetRandomColumn(eps->V,0));
541: PetscCall(BVGetColumn(X,0,&x));
542: /* v = Hmult(u,1) */
543: PetscCall(BVGetColumn(Y,0,&y));
544: PetscCall(VecNestSetSubVec(f,0,x));
545: PetscCall(VecNestSetSubVec(g,0,y));
546: mctx->s = 1.0;
547: PetscCall(STApply(eps->st,f,g));
548: /* nrm = sqrt(real(u'*v)) */
549: PetscCall(VecDot(y,x,&gamma));
550: nrm = PetscSqrtReal(PetscRealPart(gamma));
551: if (nrm==0.0) {
552: if (breakdown) *breakdown = PETSC_TRUE;
553: *M = 1; m = 0;
554: } else {
555: /* u = u /(nrm*2) */
556: PetscCall(VecScale(x,1.0/(2.0*nrm)));
557: /* v = v /(nrm*2) */
558: PetscCall(VecScale(y,1.0/(2.0*nrm)));
559: PetscCall(VecCopy(y,v));
560: /* X(:,1) = (u+v) */
561: PetscCall(VecAXPY(x,1,y));
562: /* Y(:,1) = conj(u-v) */
563: PetscCall(VecAYPX(y,-2,x));
564: PetscCall(VecConjugate(y));
565: }
566: PetscCall(BVRestoreColumn(X,0,&x));
567: PetscCall(BVRestoreColumn(Y,0,&y));
568: }
570: for (j=k;j<m;j++) {
571: /* j+1 columns (indexes 0 to j) have been computed */
572: PetscCall(BVGetColumn(X,j+1,&x));
573: PetscCall(BVGetColumn(Y,j+1,&y));
574: /* u = Hmult(v,-1)*/
575: PetscCall(VecNestSetSubVec(f,0,v));
576: PetscCall(VecNestSetSubVec(g,0,x));
577: mctx->s = -1.0;
578: PetscCall(STApply(eps->st,f,g));
579: /* hx = (u+v) */
580: PetscCall(VecCopy(x,y));
581: PetscCall(VecAXPY(x,1,v));
582: /* hy = conj(u-v) */
583: PetscCall(VecAXPY(y,-1,v));
584: PetscCall(VecConjugate(y));
585: /* [u,cd] = orthog(hx,hy,X(:,1:j),Y(:,1:j),opt)*/
586: PetscCall(OrthogonalizeVector_ProjectedBSE(x,y,X,Y,j+1,beta,k,hwork));
587: /* alpha(j) = 2*(real(cd(j))-1/2) */
588: alpha[j] = 2*(PetscRealPart(hwork[j]) - 0.5);
589: /* v = Hmult(u,1) */
590: PetscCall(VecNestSetSubVec(f,0,x));
591: PetscCall(VecNestSetSubVec(g,0,y));
592: mctx->s = 1.0;
593: PetscCall(STApply(eps->st,f,g));
594: /* nrm = sqrt(real(u'*v)) */
595: /* beta(j) = 2*nrm */
596: PetscCall(VecDot(x,y,&gamma));
597: beta[j] = 2.0*PetscSqrtReal(PetscRealPart(gamma));
598: if (beta[j]==0.0) {
599: if (breakdown) *breakdown = PETSC_TRUE;
600: *M = j+1; m = j;
601: } else {
602: /* u = u/(nrm*2) */
603: PetscCall(VecScale(x,1.0/beta[j]));
604: /* v = v/(nrm*2) */
605: PetscCall(VecScale(y,1.0/beta[j]));
606: PetscCall(VecCopy(y,v));
607: /* X(:,j+1) = (u+v) */
608: PetscCall(VecAXPY(x,1,y));
609: /* Y(:,j+1) = conj(u-v) */
610: PetscCall(VecAYPX(y,-2,x));
611: PetscCall(VecConjugate(y));
612: }
613: PetscCall(BVRestoreColumn(X,j+1,&x));
614: PetscCall(BVRestoreColumn(Y,j+1,&y));
615: }
616: if (alloc) PetscCall(PetscFree(hwork));
617: PetscCall(VecDestroy(&w));
618: PetscCall(VecDestroy(&f));
619: PetscCall(VecDestroy(&g));
620: PetscFunctionReturn(PETSC_SUCCESS);
621: }
623: static PetscErrorCode EPSComputeVectors_BSE_ProjectedBSE(EPS eps)
624: {
625: Mat H;
626: Vec x1,y1,cx1,cy1;
627: BV X,Y;
628: IS is[2];
629: PetscInt k;
630: PetscScalar delta1,delta2,lambda;
632: PetscFunctionBegin;
633: PetscCall(STGetMatrix(eps->st,0,&H));
634: PetscCall(MatNestGetISs(H,is,NULL));
635: PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&X,&Y));
636: PetscCall(BVCreateVec(X,&cx1));
637: PetscCall(BVCreateVec(Y,&cy1));
638: for (k=0; k<eps->nconv; k++) {
639: /* approx eigenvector is [ delta1*x1 + delta2*conj(y1) ]
640: [ delta1*y1 + delta2*conj(x1) ] */
641: lambda = eps->eigr[k];
642: PetscCall(STBackTransform(eps->st,1,&lambda,&eps->eigi[k]));
643: delta1 = lambda+1.0;
644: delta2 = lambda-1.0;
645: PetscCall(BVGetColumn(X,k,&x1));
646: PetscCall(BVGetColumn(Y,k,&y1));
647: PetscCall(VecCopy(x1,cx1));
648: PetscCall(VecCopy(y1,cy1));
649: PetscCall(VecConjugate(cx1));
650: PetscCall(VecConjugate(cy1));
651: PetscCall(VecScale(x1,delta1));
652: PetscCall(VecScale(y1,delta1));
653: PetscCall(VecAXPY(x1,delta2,cy1));
654: PetscCall(VecAXPY(y1,delta2,cx1));
655: PetscCall(BVRestoreColumn(X,k,&x1));
656: PetscCall(BVRestoreColumn(Y,k,&y1));
657: }
658: PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&X,&Y));
659: PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
660: /* Normalize eigenvector */
661: PetscCall(BVNormalize(eps->V,NULL));
662: PetscCall(VecDestroy(&cx1));
663: PetscCall(VecDestroy(&cy1));
664: PetscFunctionReturn(PETSC_SUCCESS);
665: }
667: PetscErrorCode EPSSetUp_KrylovSchur_BSE(EPS eps)
668: {
669: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
670: PetscBool flg,sinvert;
672: PetscFunctionBegin;
673: PetscCheck((eps->problem_type==EPS_BSE),PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Problem type should be BSE");
674: EPSCheckUnsupportedCondition(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION | EPS_FEATURE_BALANCE,PETSC_TRUE," with BSE structure");
675: PetscCall(PetscCitationsRegister(citation,&cited));
676: PetscCall(EPSSetDimensions_Default(eps,&eps->nev,&eps->ncv,&eps->mpd));
677: 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");
678: if (eps->max_it==PETSC_DETERMINE) eps->max_it = PetscMax(100,2*eps->n/eps->ncv)*((eps->stop==EPS_STOP_THRESHOLD)?10:1);
680: PetscCall(PetscObjectTypeCompareAny((PetscObject)eps->st,&flg,STSINVERT,STSHIFT,""));
681: PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Krylov-Schur BSE only supports shift and shift-and-invert ST");
682: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STSINVERT,&sinvert));
683: if (!eps->which) {
684: if (sinvert) eps->which = EPS_TARGET_MAGNITUDE;
685: else eps->which = EPS_SMALLEST_MAGNITUDE;
686: }
688: if (!ctx->keep) ctx->keep = 0.5;
689: PetscCall(STSetStructured(eps->st,PETSC_TRUE));
691: PetscCall(EPSAllocateSolution(eps,1));
692: switch (ctx->bse) {
693: case EPS_KRYLOVSCHUR_BSE_SHAO:
694: eps->ops->solve = EPSSolve_KrylovSchur_BSE_Shao;
695: eps->ops->computevectors = EPSComputeVectors_BSE_Shao;
696: PetscCall(DSSetType(eps->ds,DSHEP));
697: PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
698: PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
699: PetscCall(DSAllocate(eps->ds,eps->ncv+1));
700: break;
701: case EPS_KRYLOVSCHUR_BSE_GRUNING:
702: eps->ops->solve = EPSSolve_KrylovSchur_BSE_Gruning;
703: eps->ops->computevectors = EPSComputeVectors_BSE_Gruning;
704: PetscCall(DSSetType(eps->ds,DSSVD));
705: PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
706: PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
707: PetscCall(DSAllocate(eps->ds,eps->ncv+1));
708: break;
709: case EPS_KRYLOVSCHUR_BSE_PROJECTEDBSE:
710: eps->ops->solve = EPSSolve_KrylovSchur_BSE_ProjectedBSE;
711: eps->ops->computevectors = EPSComputeVectors_BSE_ProjectedBSE;
712: PetscCall(DSSetType(eps->ds,DSHEP));
713: PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
714: PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
715: PetscCall(DSAllocate(eps->ds,eps->ncv+1));
716: break;
717: default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Unexpected error");
718: }
719: PetscFunctionReturn(PETSC_SUCCESS);
720: }
722: PetscErrorCode EPSSolve_KrylovSchur_BSE_Shao(EPS eps)
723: {
724: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
725: PetscInt i,k,l,ld,nv,nconv=0,nevsave;
726: Mat H,Q;
727: BV U,V;
728: IS is[2];
729: PetscReal *a,*b,beta;
730: PetscBool breakdown=PETSC_FALSE;
732: PetscFunctionBegin;
733: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
735: /* Extract matrix blocks */
736: PetscCall(STGetMatrix(eps->st,0,&H));
737: PetscCall(MatNestGetISs(H,is,NULL));
739: /* Get the split bases */
740: PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
742: nevsave = eps->nev;
743: eps->nev = (eps->nev+1)/2;
744: l = 0;
746: /* Restart loop */
747: while (eps->reason == EPS_CONVERGED_ITERATING) {
748: eps->its++;
750: /* Compute an nv-step Lanczos factorization */
751: nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
752: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
753: PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&a));
754: b = a + ld;
755: PetscCall(EPSBSELanczos_Shao(eps,U,V,a,b,eps->nconv+l,&nv,&breakdown));
756: beta = b[nv-1];
757: PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&a));
758: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
759: PetscCall(DSSetState(eps->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
760: PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
762: /* Solve projected problem */
763: PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
764: PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
765: PetscCall(DSUpdateExtraRow(eps->ds));
766: PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
768: /* Check convergence */
769: for (i=0;i<nv;i++) eps->eigr[i] = PetscSqrtReal(PetscRealPart(eps->eigr[i]));
770: PetscCall(EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta,0.0,1.0,&k));
771: EPSSetCtxThreshold(eps,eps->eigr,eps->eigi,eps->errest,k,nv);
772: PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
773: nconv = k;
775: /* Update l */
776: if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
777: else l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
778: if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
779: if (l) PetscCall(PetscInfo(eps,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
781: if (eps->reason == EPS_CONVERGED_ITERATING) {
782: PetscCheck(!breakdown,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Breakdown in BSE Krylov-Schur (beta=%g)",(double)beta);
783: /* Prepare the Rayleigh quotient for restart */
784: PetscCall(DSTruncate(eps->ds,k+l,PETSC_FALSE));
785: }
786: /* Update the corresponding vectors
787: U(:,idx) = U*Q(:,idx), V(:,idx) = V*Q(:,idx) */
788: PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&Q));
789: PetscCall(BVMultInPlace(U,Q,eps->nconv,k+l));
790: PetscCall(BVMultInPlace(V,Q,eps->nconv,k+l));
791: PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&Q));
793: if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
794: PetscCall(BVCopyColumn(eps->V,nv,k+l));
795: if (eps->stop==EPS_STOP_THRESHOLD && nv-k<5) { /* reallocate */
796: eps->ncv = eps->mpd+k;
797: PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
798: PetscCall(EPSReallocateSolution(eps,eps->ncv+1));
799: PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
800: for (i=nv;i<eps->ncv;i++) eps->perm[i] = i;
801: PetscCall(DSReallocate(eps->ds,eps->ncv+1));
802: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
803: }
804: }
805: eps->nconv = k;
806: PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv));
807: }
809: eps->nev = nevsave;
811: PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
812: PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
813: PetscFunctionReturn(PETSC_SUCCESS);
814: }
816: /*
817: EPSConvergence_Gruning - convergence check based on SVDKrylovConvergence().
818: */
819: static PetscErrorCode EPSConvergence_Gruning(EPS eps,PetscBool getall,PetscInt kini,PetscInt nits,PetscInt *kout)
820: {
821: PetscInt k,marker,ld;
822: PetscReal *alpha,*beta,resnorm;
823: PetscBool extra;
825: PetscFunctionBegin;
826: *kout = 0;
827: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
828: PetscCall(DSGetExtraRow(eps->ds,&extra));
829: PetscCheck(extra,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Only implemented for DS with extra row");
830: marker = -1;
831: if (eps->trackall) getall = PETSC_TRUE;
832: PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&alpha));
833: beta = alpha + ld;
834: for (k=kini;k<kini+nits;k++) {
835: resnorm = PetscAbsReal(beta[k]);
836: PetscCall((*eps->converged)(eps,eps->eigr[k],eps->eigi[k],resnorm,&eps->errest[k],eps->convergedctx));
837: if (marker==-1 && eps->errest[k] >= eps->tol) marker = k;
838: if (marker!=-1 && !getall) break;
839: }
840: PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&alpha));
841: if (marker!=-1) k = marker;
842: *kout = k;
843: PetscFunctionReturn(PETSC_SUCCESS);
844: }
846: PetscErrorCode EPSSolve_KrylovSchur_BSE_Gruning(EPS eps)
847: {
848: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
849: PetscInt i,k,l,ld,nv,nconv=0,nevsave;
850: Mat H,Q,Z;
851: BV U,V,HU,HV;
852: IS is[2];
853: PetscReal *d1,*d2,beta;
854: PetscBool breakdown=PETSC_FALSE;
856: PetscFunctionBegin;
857: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
859: /* Extract matrix blocks */
860: PetscCall(STGetMatrix(eps->st,0,&H));
861: PetscCall(MatNestGetISs(H,is,NULL));
863: /* Get the split bases */
864: PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
866: /* Create HU and HV */
867: PetscCall(BVDuplicate(U,&HU));
868: PetscCall(BVDuplicate(V,&HV));
870: nevsave = eps->nev;
871: eps->nev = (eps->nev+1)/2;
872: l = 0;
874: /* Restart loop */
875: while (eps->reason == EPS_CONVERGED_ITERATING) {
876: eps->its++;
878: /* Compute an nv-step Lanczos factorization */
879: nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
880: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
881: PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&d1));
882: d2 = d1 + ld;
883: PetscCall(EPSBSELanczos_Gruning(eps,U,V,HU,HV,d1,d2,eps->nconv+l,&nv,&breakdown));
884: beta = d1[nv-1];
885: PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&d1));
887: /* Compute SVD */
888: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
889: PetscCall(DSSVDSetDimensions(eps->ds,nv));
890: PetscCall(DSSetState(eps->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
891: PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
893: PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
894: PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
895: PetscCall(DSUpdateExtraRow(eps->ds));
896: PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
898: /* Check convergence */
899: PetscCall(EPSConvergence_Gruning(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,&k));
900: EPSSetCtxThreshold(eps,eps->eigr,eps->eigi,eps->errest,k,nv);
901: PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
902: nconv = k;
904: /* Update l */
905: if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
906: else l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
907: if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
908: if (l) PetscCall(PetscInfo(eps,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
910: if (eps->reason == EPS_CONVERGED_ITERATING) {
911: PetscCheck(!breakdown,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Breakdown in BSE Krylov-Schur (beta=%g)",(double)beta);
912: /* Prepare the Rayleigh quotient for restart */
913: PetscCall(DSTruncate(eps->ds,k+l,PETSC_FALSE));
914: }
915: /* Update the corresponding vectors
916: U(:,idx) = U*Q(:,idx), V(:,idx) = V*Z(:,idx) */
917: PetscCall(DSGetMat(eps->ds,DS_MAT_U,&Q));
918: PetscCall(DSGetMat(eps->ds,DS_MAT_V,&Z));
919: PetscCall(BVMultInPlace(U,Z,eps->nconv,k+l));
920: PetscCall(BVMultInPlace(V,Q,eps->nconv,k+l));
921: PetscCall(BVMultInPlace(HU,Z,eps->nconv,k+l));
922: PetscCall(BVMultInPlace(HV,Q,eps->nconv,k+l));
923: PetscCall(DSRestoreMat(eps->ds,DS_MAT_U,&Q));
924: PetscCall(DSRestoreMat(eps->ds,DS_MAT_V,&Z));
926: if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
927: PetscCall(BVCopyColumn(U,nv,k+l));
928: PetscCall(BVCopyColumn(HU,nv,k+l));
929: if (eps->stop==EPS_STOP_THRESHOLD && nv-k<5) { /* reallocate */
930: eps->ncv = eps->mpd+k;
931: PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
932: PetscCall(EPSReallocateSolution(eps,eps->ncv+1));
933: PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
934: PetscCall(BVResize(HU,eps->ncv+1,PETSC_TRUE));
935: PetscCall(BVResize(HV,eps->ncv+1,PETSC_TRUE));
936: for (i=nv;i<eps->ncv;i++) { eps->perm[i] = i; eps->eigi[i] = 0.0; }
937: PetscCall(DSReallocate(eps->ds,eps->ncv+1));
938: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
939: }
940: }
941: eps->nconv = k;
942: PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv));
943: }
945: eps->nev = nevsave;
947: PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
948: PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
950: PetscCall(BVDestroy(&HU));
951: PetscCall(BVDestroy(&HV));
952: PetscFunctionReturn(PETSC_SUCCESS);
953: }
955: PetscErrorCode EPSSolve_KrylovSchur_BSE_ProjectedBSE(EPS eps)
956: {
957: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
958: PetscInt i,k,l,ld,nv,nconv=0,nevsave;
959: Mat H,Q;
960: Vec v;
961: BV U,V;
962: IS is[2];
963: PetscReal *a,*b,beta;
964: PetscBool breakdown=PETSC_FALSE;
966: PetscFunctionBegin;
967: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
969: /* Extract matrix blocks */
970: PetscCall(STGetMatrix(eps->st,0,&H));
971: PetscCall(MatNestGetISs(H,is,NULL));
973: /* Get the split bases */
974: PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
976: /* Vector v */
977: PetscCall(BVCreateVec(V,&v));
979: nevsave = eps->nev;
980: eps->nev = (eps->nev+1)/2;
981: l = 0;
983: /* Restart loop */
984: while (eps->reason == EPS_CONVERGED_ITERATING) {
985: eps->its++;
987: /* Compute an nv-step Lanczos factorization */
988: nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
989: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
990: PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&a));
991: b = a + ld;
992: PetscCall(EPSBSELanczos_ProjectedBSE(eps,U,V,v,a,b,eps->nconv+l,&nv,&breakdown));
993: beta = b[nv-1];
994: PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&a));
995: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
996: PetscCall(DSSetState(eps->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
997: PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
999: /* Solve projected problem */
1000: PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
1001: PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
1002: PetscCall(DSUpdateExtraRow(eps->ds));
1003: PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
1005: /* Check convergence */
1006: for (i=0;i<nv;i++) eps->eigr[i] = PetscSqrtReal(PetscRealPart(eps->eigr[i]));
1007: PetscCall(EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta,0.0,1.0,&k));
1008: EPSSetCtxThreshold(eps,eps->eigr,eps->eigi,eps->errest,k,nv);
1009: PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
1010: nconv = k;
1012: /* Update l */
1013: if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
1014: else l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
1015: if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
1016: if (l) PetscCall(PetscInfo(eps,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
1018: if (eps->reason == EPS_CONVERGED_ITERATING) {
1019: PetscCheck(!breakdown,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Breakdown in BSE Krylov-Schur (beta=%g)",(double)beta);
1020: /* Prepare the Rayleigh quotient for restart */
1021: PetscCall(DSTruncate(eps->ds,k+l,PETSC_FALSE));
1022: }
1023: /* Update the corresponding vectors
1024: U(:,idx) = U*Q(:,idx), V(:,idx) = V*Q(:,idx) */
1025: PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&Q));
1026: PetscCall(BVMultInPlace(U,Q,eps->nconv,k+l));
1027: PetscCall(BVMultInPlace(V,Q,eps->nconv,k+l));
1028: PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&Q));
1030: if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
1031: PetscCall(BVCopyColumn(eps->V,nv,k+l));
1032: if (eps->stop==EPS_STOP_THRESHOLD && nv-k<5) { /* reallocate */
1033: eps->ncv = eps->mpd+k;
1034: PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
1035: PetscCall(EPSReallocateSolution(eps,eps->ncv+1));
1036: PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
1037: for (i=nv;i<eps->ncv;i++) eps->perm[i] = i;
1038: PetscCall(DSReallocate(eps->ds,eps->ncv+1));
1039: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
1040: }
1041: }
1042: eps->nconv = k;
1043: PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv));
1044: }
1046: eps->nev = nevsave;
1048: PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
1049: PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
1051: PetscCall(VecDestroy(&v));
1052: PetscFunctionReturn(PETSC_SUCCESS);
1053: }