Actual source code: ks-bse.c
slepc-3.22.1 2024-10-28
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: */
21: #include <slepc/private/epsimpl.h>
22: #include "krylovschur.h"
24: static PetscErrorCode Orthog_Shao(Vec x,BV U,BV V,PetscInt j,PetscScalar *h,PetscScalar *c,PetscBool *breakdown)
25: {
26: PetscInt i;
28: PetscFunctionBegin;
29: PetscCall(BVSetActiveColumns(U,0,j));
30: PetscCall(BVSetActiveColumns(V,0,j));
31: /* c = real(V^* x) ; c2 = imag(U^* x)*1i */
32: #if defined(PETSC_USE_COMPLEX)
33: PetscCall(BVDotVecBegin(V,x,c));
34: PetscCall(BVDotVecBegin(U,x,c+j));
35: PetscCall(BVDotVecEnd(V,x,c));
36: PetscCall(BVDotVecEnd(U,x,c+j));
37: #else
38: PetscCall(BVDotVec(V,x,c));
39: #endif
40: #if defined(PETSC_USE_COMPLEX)
41: for (i=0; i<j; i++) {
42: c[i] = PetscRealPart(c[i]);
43: c[j+i] = PetscCMPLX(0.0,PetscImaginaryPart(c[j+i]));
44: }
45: #endif
46: /* x = x-U*c-V*c2 */
47: PetscCall(BVMultVec(U,-1.0,1.0,x,c));
48: #if defined(PETSC_USE_COMPLEX)
49: PetscCall(BVMultVec(V,-1.0,1.0,x,c+j));
50: #endif
51: /* accumulate orthog coeffs into h */
52: for (i=0; i<2*j; i++) h[i] += c[i];
53: PetscFunctionReturn(PETSC_SUCCESS);
54: }
56: /* Orthogonalize vector x against first j vectors in U and V
57: v is column j-1 of V */
58: static PetscErrorCode OrthogonalizeVector_Shao(Vec x,BV U,BV V,PetscInt j,Vec v,PetscReal *beta,PetscInt k,PetscScalar *h,PetscBool *breakdown)
59: {
60: PetscReal alpha;
61: PetscInt i,l;
63: PetscFunctionBegin;
64: PetscCall(PetscArrayzero(h,2*j));
66: /* Local orthogonalization */
67: l = j==k+1?0:j-2; /* 1st column to orthogonalize against */
68: PetscCall(VecDotRealPart(x,v,&alpha));
69: for (i=l; i<j-1; i++) h[i] = beta[i];
70: h[j-1] = alpha;
71: /* x = x-U(:,l:j-1)*h(l:j-1) */
72: PetscCall(BVSetActiveColumns(U,l,j));
73: PetscCall(BVMultVec(U,-1.0,1.0,x,h+l));
75: /* Full orthogonalization */
76: PetscCall(Orthog_Shao(x,U,V,j,h,h+2*j,breakdown));
77: PetscFunctionReturn(PETSC_SUCCESS);
78: }
80: static PetscErrorCode EPSBSELanczos_Shao(EPS eps,BV U,BV V,PetscReal *alpha,PetscReal *beta,PetscInt k,PetscInt *M,PetscBool *breakdown)
81: {
82: PetscInt j,m = *M;
83: Vec v,x,y,w,f,g,vecs[2];
84: Mat H;
85: IS is[2];
86: PetscReal nrm;
87: PetscScalar *hwork,lhwork[100],gamma;
89: PetscFunctionBegin;
90: if (4*m > 100) PetscCall(PetscMalloc1(4*m,&hwork));
91: else hwork = lhwork;
92: PetscCall(STGetMatrix(eps->st,0,&H));
93: PetscCall(MatNestGetISs(H,is,NULL));
95: /* create work vectors */
96: PetscCall(BVGetColumn(V,0,&v));
97: PetscCall(VecDuplicate(v,&w));
98: vecs[0] = v;
99: vecs[1] = w;
100: PetscCall(VecCreateNest(PetscObjectComm((PetscObject)eps),2,is,vecs,&f));
101: PetscCall(VecCreateNest(PetscObjectComm((PetscObject)eps),2,is,vecs,&g));
102: PetscCall(BVRestoreColumn(V,0,&v));
104: /* Normalize initial vector */
105: if (k==0) {
106: if (eps->nini==0) PetscCall(BVSetRandomColumn(eps->V,0));
107: PetscCall(BVGetColumn(U,0,&x));
108: PetscCall(BVGetColumn(V,0,&y));
109: PetscCall(VecCopy(x,w));
110: PetscCall(VecConjugate(w));
111: PetscCall(VecNestSetSubVec(f,0,x));
112: PetscCall(VecNestSetSubVec(g,0,y));
113: PetscCall(STApply(eps->st,f,g));
114: PetscCall(VecDot(y,x,&gamma));
115: nrm = PetscSqrtReal(PetscRealPart(gamma));
116: PetscCall(VecScale(x,1.0/nrm));
117: PetscCall(VecScale(y,1.0/nrm));
118: PetscCall(BVRestoreColumn(U,0,&x));
119: PetscCall(BVRestoreColumn(V,0,&y));
120: }
122: for (j=k;j<m;j++) {
123: /* j+1 columns (indexes 0 to j) have been computed */
124: PetscCall(BVGetColumn(V,j,&v));
125: PetscCall(BVGetColumn(U,j+1,&x));
126: PetscCall(BVGetColumn(V,j+1,&y));
127: PetscCall(VecCopy(v,w));
128: PetscCall(VecConjugate(w));
129: PetscCall(VecScale(w,-1.0));
130: PetscCall(VecNestSetSubVec(f,0,v));
131: PetscCall(VecNestSetSubVec(g,0,x));
132: PetscCall(STApply(eps->st,f,g));
133: PetscCall(OrthogonalizeVector_Shao(x,U,V,j+1,v,beta,k,hwork,breakdown));
134: alpha[j] = PetscRealPart(hwork[j]);
135: PetscCall(VecCopy(x,w));
136: PetscCall(VecConjugate(w));
137: PetscCall(VecNestSetSubVec(f,0,x));
138: PetscCall(VecNestSetSubVec(g,0,y));
139: PetscCall(STApply(eps->st,f,g));
140: PetscCall(VecDot(x,y,&gamma));
141: beta[j] = PetscSqrtReal(PetscRealPart(gamma));
142: PetscCall(VecScale(x,1.0/beta[j]));
143: PetscCall(VecScale(y,1.0/beta[j]));
144: PetscCall(BVRestoreColumn(V,j,&v));
145: PetscCall(BVRestoreColumn(U,j+1,&x));
146: PetscCall(BVRestoreColumn(V,j+1,&y));
147: }
148: if (4*m > 100) PetscCall(PetscFree(hwork));
149: PetscCall(VecDestroy(&w));
150: PetscCall(VecDestroy(&f));
151: PetscCall(VecDestroy(&g));
152: PetscFunctionReturn(PETSC_SUCCESS);
153: }
155: static PetscErrorCode EPSComputeVectors_BSE_Shao(EPS eps)
156: {
157: Mat H;
158: Vec u1,v1;
159: BV U,V;
160: IS is[2];
161: PetscInt k;
162: PetscScalar lambda;
164: PetscFunctionBegin;
165: PetscCall(STGetMatrix(eps->st,0,&H));
166: PetscCall(MatNestGetISs(H,is,NULL));
167: PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
168: for (k=0; k<eps->nconv; k++) {
169: PetscCall(BVGetColumn(U,k,&u1));
170: PetscCall(BVGetColumn(V,k,&v1));
171: /* approx eigenvector is [ (eigr[k]*u1+v1)]
172: [conj(eigr[k]*u1-v1)] */
173: lambda = eps->eigr[k];
174: PetscCall(STBackTransform(eps->st,1,&lambda,&eps->eigi[k]));
175: PetscCall(VecAYPX(u1,lambda,v1));
176: PetscCall(VecAYPX(v1,-2.0,u1));
177: PetscCall(VecConjugate(v1));
178: PetscCall(BVRestoreColumn(U,k,&u1));
179: PetscCall(BVRestoreColumn(V,k,&v1));
180: }
181: PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
182: /* Normalize eigenvectors */
183: PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
184: PetscCall(BVNormalize(eps->V,NULL));
185: PetscFunctionReturn(PETSC_SUCCESS);
186: }
188: static PetscErrorCode Orthog_Gruning(Vec x,BV U,BV V,BV HU,BV HV,PetscInt j,PetscScalar *h,PetscScalar *c,PetscBool s,PetscBool *breakdown)
189: {
190: PetscInt i;
192: PetscFunctionBegin;
193: PetscCall(BVSetActiveColumns(U,0,j));
194: PetscCall(BVSetActiveColumns(HU,0,j));
195: if (s) {
196: PetscCall(BVSetActiveColumns(V,0,j));
197: PetscCall(BVSetActiveColumns(HV,0,j));
198: } else {
199: PetscCall(BVSetActiveColumns(V,0,j-1));
200: PetscCall(BVSetActiveColumns(HV,0,j-1));
201: }
202: #if defined(PETSC_USE_COMPLEX)
203: PetscCall(BVDotVecBegin(HU,x,c));
204: if (s || j>1) PetscCall(BVDotVecBegin(HV,x,c+j));
205: PetscCall(BVDotVecEnd(HU,x,c));
206: if (s || j>1) PetscCall(BVDotVecEnd(HV,x,c+j));
207: #else
208: if (s) PetscCall(BVDotVec(HU,x,c));
209: else PetscCall(BVDotVec(HV,x,c+j));
210: #endif
211: for (i=0; i<j; i++) {
212: if (s) { /* c1 = 2*real(HU^* x) ; c2 = 2*imag(HV^* x)*1i */
213: #if defined(PETSC_USE_COMPLEX)
214: c[i] = PetscRealPart(c[i]);
215: c[j+i] = PetscCMPLX(0.0,PetscImaginaryPart(c[j+i]));
216: #else
217: c[j+i] = 0.0;
218: #endif
219: } else { /* c1 = 2*imag(HU^* x)*1i ; c2 = 2*real(HV^* x) */
220: #if defined(PETSC_USE_COMPLEX)
221: c[i] = PetscCMPLX(0.0,PetscImaginaryPart(c[i]));
222: c[j+i] = PetscRealPart(c[j+i]);
223: #else
224: c[i] = 0.0;
225: #endif
226: }
227: }
228: /* x = x-U*c1-V*c2 */
229: #if defined(PETSC_USE_COMPLEX)
230: PetscCall(BVMultVec(U,-2.0,1.0,x,c));
231: PetscCall(BVMultVec(V,-2.0,1.0,x,c+j));
232: #else
233: if (s) PetscCall(BVMultVec(U,-2.0,1.0,x,c));
234: else PetscCall(BVMultVec(V,-2.0,1.0,x,c+j));
235: #endif
236: /* accumulate orthog coeffs into h */
237: for (i=0; i<2*j; i++) h[i] += 2*c[i];
238: PetscFunctionReturn(PETSC_SUCCESS);
239: }
241: /* Orthogonalize vector x against first j vectors in U and V */
242: static PetscErrorCode OrthogonalizeVector_Gruning(Vec x,BV U,BV V,BV HU,BV HV,PetscInt j,PetscReal *beta,PetscInt k,PetscScalar *h,PetscBool s,PetscBool *breakdown)
243: {
244: PetscInt l,i;
245: Vec u;
247: PetscFunctionBegin;
248: PetscCall(PetscArrayzero(h,4*j));
250: /* Local orthogonalization */
251: if (s) {
252: PetscCall(BVGetColumn(U,j-1,&u));
253: PetscCall(VecAXPY(x,-*beta,u));
254: PetscCall(BVRestoreColumn(U,j-1,&u));
255: h[j-1] = *beta;
256: } else {
257: l = j==k+1?0:j-2; /* 1st column to orthogonalize against */
258: for (i=l; i<j-1; i++) h[j+i] = beta[i];
259: /* x = x-V(:,l:j-2)*h(l:j-2) */
260: PetscCall(BVSetActiveColumns(V,l,j-1));
261: PetscCall(BVMultVec(V,-1.0,1.0,x,h+j+l));
262: }
264: /* Full orthogonalization */
265: PetscCall(Orthog_Gruning(x,U,V,HU,HV,j,h,h+2*j,s,breakdown));
266: PetscFunctionReturn(PETSC_SUCCESS);
267: }
269: static PetscErrorCode EPSBSELanczos_Gruning(EPS eps,BV U,BV V,BV HU,BV HV,PetscReal *beta1,PetscReal *beta2,PetscInt k,PetscInt *M,PetscBool *breakdown)
270: {
271: PetscInt j,m = *M;
272: Vec v,x,y,w,f,g,vecs[2];
273: Mat H;
274: IS is[2];
275: PetscReal nrm;
276: PetscScalar *hwork,lhwork[100],dot;
278: PetscFunctionBegin;
279: if (4*m > 100) PetscCall(PetscMalloc1(4*m,&hwork));
280: else hwork = lhwork;
281: PetscCall(STGetMatrix(eps->st,0,&H));
282: PetscCall(MatNestGetISs(H,is,NULL));
284: /* create work vectors */
285: PetscCall(BVGetColumn(V,0,&v));
286: PetscCall(VecDuplicate(v,&w));
287: vecs[0] = v;
288: vecs[1] = w;
289: PetscCall(VecCreateNest(PetscObjectComm((PetscObject)eps),2,is,vecs,&f));
290: PetscCall(VecCreateNest(PetscObjectComm((PetscObject)eps),2,is,vecs,&g));
291: PetscCall(BVRestoreColumn(V,0,&v));
293: /* Normalize initial vector */
294: if (k==0 && eps->nini==0) PetscCall(BVSetRandomColumn(eps->V,0));
296: /* y = Hmult(v1,1) */
297: PetscCall(BVGetColumn(U,k,&x));
298: PetscCall(BVGetColumn(HU,k,&y));
299: PetscCall(VecCopy(x,w));
300: PetscCall(VecConjugate(w));
301: PetscCall(VecNestSetSubVec(f,0,x));
302: PetscCall(VecNestSetSubVec(g,0,y));
303: PetscCall(STApply(eps->st,f,g));
304: /* nrm = sqrt(2*real(u1'*y)); */
305: PetscCall(VecDot(x,y,&dot));
306: nrm = PetscSqrtReal(PetscRealPart(2*dot));
307: /* U(:,j) = u1/nrm; */
308: /* HU(:,j) = y/nrm; */
309: PetscCall(VecScale(x,1.0/nrm));
310: PetscCall(VecScale(y,1.0/nrm));
311: PetscCall(BVRestoreColumn(U,k,&x));
312: PetscCall(BVRestoreColumn(HU,k,&y));
314: for (j=k;j<m;j++) {
315: /* j+1 columns (indexes 0 to j) have been computed */
316: PetscCall(BVGetColumn(HU,j,&x));
317: PetscCall(BVGetColumn(V,j,&v));
318: PetscCall(BVGetColumn(HV,j,&y));
319: PetscCall(VecCopy(x,v));
320: PetscCall(BVRestoreColumn(HU,j,&x));
321: /* v = Orthogonalize HU(:,j) */
322: PetscCall(OrthogonalizeVector_Gruning(v,U,V,HU,HV,j+1,beta2,k,hwork,PETSC_FALSE,breakdown));
323: /* y = Hmult(v,-1) */
324: PetscCall(VecCopy(v,w));
325: PetscCall(VecConjugate(w));
326: PetscCall(VecScale(w,-1.0));
327: PetscCall(VecNestSetSubVec(f,0,v));
328: PetscCall(VecNestSetSubVec(g,0,y));
329: PetscCall(STApply(eps->st,f,g));
330: /* beta = sqrt(2*real(v'*y)); */
331: PetscCall(VecDot(v,y,&dot));
332: beta1[j] = PetscSqrtReal(PetscRealPart(2*dot)); //FIXME Check beta != 0?
333: /* V(:,j) = v/beta1; */
334: /* HV(:,j) = y/beta1; */
335: PetscCall(VecScale(v,1.0/beta1[j]));
336: PetscCall(VecScale(y,1.0/beta1[j]));
337: PetscCall(BVRestoreColumn(V,j,&v));
338: PetscCall(BVRestoreColumn(HV,j,&y));
340: PetscCall(BVGetColumn(HV,j,&x));
341: PetscCall(BVGetColumn(U,j+1,&v));
342: PetscCall(BVGetColumn(HU,j+1,&y));
343: PetscCall(VecCopy(x,v));
344: PetscCall(BVRestoreColumn(HV,j,&x));
345: /* v = Orthogonalize HV(:,j) */
346: PetscCall(OrthogonalizeVector_Gruning(v,U,V,HU,HV,j+1,&beta1[j],k,hwork,PETSC_TRUE,breakdown));
347: /* y = Hmult(v,1) */
348: PetscCall(VecCopy(v,w));
349: PetscCall(VecConjugate(w));
350: PetscCall(VecNestSetSubVec(f,0,v));
351: PetscCall(VecNestSetSubVec(g,0,y));
352: PetscCall(STApply(eps->st,f,g));
353: /* beta = sqrt(2*real(v'*y)); */
354: PetscCall(VecDot(v,y,&dot));
355: beta2[j] = PetscSqrtReal(PetscRealPart(2*dot)); //FIXME Check beta != 0?
356: /* U(:,j) = v/beta2; */
357: /* HU(:,j) = y/beta2; */
358: PetscCall(VecScale(v,1.0/beta2[j]));
359: PetscCall(VecScale(y,1.0/beta2[j]));
360: PetscCall(BVRestoreColumn(U,j+1,&v));
361: PetscCall(BVRestoreColumn(HU,j+1,&y));
362: }
363: if (4*m > 100) PetscCall(PetscFree(hwork));
364: PetscCall(VecDestroy(&w));
365: PetscCall(VecDestroy(&f));
366: PetscCall(VecDestroy(&g));
367: PetscFunctionReturn(PETSC_SUCCESS);
368: }
370: static PetscErrorCode EPSComputeVectors_BSE_Gruning(EPS eps)
371: {
372: Mat H;
373: Vec u1,v1;
374: BV U,V;
375: IS is[2];
376: PetscInt k;
378: PetscFunctionBegin;
379: PetscCall(STGetMatrix(eps->st,0,&H));
380: PetscCall(MatNestGetISs(H,is,NULL));
381: PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
382: /* approx eigenvector [x1] is [ u1+v1 ]
383: [x2] [conj(u1)-conj(v1)] */
384: for (k=0; k<eps->nconv; k++) {
385: PetscCall(BVGetColumn(U,k,&u1));
386: PetscCall(BVGetColumn(V,k,&v1));
387: /* x1 = u1 + v1 */
388: PetscCall(VecAXPY(u1,1.0,v1));
389: /* x2 = conj(u1) - conj(v1) = conj(u1 - v1) = conj((u1 + v1) - 2*v1) */
390: PetscCall(VecAYPX(v1,-2.0,u1));
391: PetscCall(VecConjugate(v1));
392: PetscCall(BVRestoreColumn(U,k,&u1));
393: PetscCall(BVRestoreColumn(V,k,&v1));
394: }
395: PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
396: /* Normalize eigenvectors */
397: PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
398: PetscCall(BVNormalize(eps->V,NULL));
399: PetscFunctionReturn(PETSC_SUCCESS);
400: }
402: static PetscErrorCode Orthog_ProjectedBSE(Vec hx,Vec hy,BV X,BV Y,PetscInt j,PetscScalar *h,PetscScalar *c,PetscBool *breakdown)
403: {
404: PetscInt i;
405: Mat MX,MY,MXl,MYl;
406: Vec c1,c2,hxl,hyl,hz;
407: PetscScalar *cx1,*cx2;
408: PetscMPIInt len;
410: PetscFunctionBegin;
411: PetscCall(PetscMPIIntCast(j,&len));
412: PetscCall(BVSetActiveColumns(X,0,j));
413: PetscCall(BVSetActiveColumns(Y,0,j));
414: /* BVTDotVec does not exist yet, implemented via MatMult operations */
415: PetscCall(BVGetMat(X,&MX));
416: PetscCall(BVGetMat(Y,&MY));
417: PetscCall(MatDenseGetLocalMatrix(MX,&MXl));
418: PetscCall(MatDenseGetLocalMatrix(MY,&MYl));
419: PetscCall(MatCreateVecs(MXl,&c1,&hyl));
420: PetscCall(MatCreateVecs(MXl,&c2,&hxl));
422: /* c1 = X^* hx - Y^* hy
423: * c2 = -Y^T hx + X^T hy */
425: PetscCall(VecGetLocalVector(hx,hxl));
426: PetscCall(VecGetLocalVector(hy,hyl));
427: PetscCall(VecDuplicate(hx,&hz));
428: /* c1 = -(Y^* hy) */
429: PetscCall(MatMultHermitianTranspose(MYl,hyl,c1));
430: PetscCall(VecScale(c1,-1));
431: /* c1 = c1 + X^* hx */
432: PetscCall(MatMultHermitianTransposeAdd(MXl,hxl,c1,c1));
433: PetscCall(VecGetArray(c1,&cx1));
434: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE,cx1,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)hx)));
435: PetscCall(VecRestoreArray(c1,&cx1));
436: /* c2 = -(Y^T hx) */
437: PetscCall(MatMultTranspose(MYl,hxl,c2));
438: PetscCall(VecScale(c2,-1));
439: /* c2 = c2 + X^T hy */
440: PetscCall(MatMultTransposeAdd(MXl,hyl,c2,c2));
441: PetscCall(VecGetArray(c2,&cx2));
442: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE,cx2,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)hx)));
443: PetscCall(VecRestoreArray(c2,&cx2));
444: PetscCall(VecRestoreLocalVector(hx,hxl));
445: PetscCall(VecRestoreLocalVector(hy,hyl));
447: /* accumulate orthog coeffs into h */
448: PetscCall(VecGetArrayRead(c1,(const PetscScalar**)&cx1));
449: PetscCall(VecGetArrayRead(c2,(const PetscScalar**)&cx2));
450: for (i=0; i<j; i++) h[i] += cx1[i];
451: for (i=0; i<j; i++) h[i+j] += cx2[i];
452: PetscCall(VecRestoreArrayRead(c1,(const PetscScalar**)&cx1));
453: PetscCall(VecRestoreArrayRead(c2,(const PetscScalar**)&cx2));
455: /* u = hx - X c1 - conj(Y) c2 */
457: /* conj(Y) c2 */
458: PetscCall(VecConjugate(c2));
459: PetscCall(VecGetLocalVector(hz,hxl));
460: PetscCall(MatMult(MYl,c2,hxl));
461: PetscCall(VecConjugate(hxl));
462: /* X c1 */
463: PetscCall(MatMultAdd(MXl,c1,hxl,hxl));
464: PetscCall(VecRestoreLocalVector(hz,hxl));
465: PetscCall(VecAXPY(hx,-1,hz));
467: PetscCall(BVRestoreMat(X,&MX));
468: PetscCall(BVRestoreMat(Y,&MY));
469: PetscCall(VecDestroy(&c1));
470: PetscCall(VecDestroy(&c2));
471: PetscCall(VecDestroy(&hxl));
472: PetscCall(VecDestroy(&hyl));
473: PetscCall(VecDestroy(&hz));
474: PetscFunctionReturn(PETSC_SUCCESS);
475: }
477: /* Orthogonalize vector x against first j vectors in X and Y */
478: static PetscErrorCode OrthogonalizeVector_ProjectedBSE(Vec hx,Vec hy,BV X,BV Y,PetscInt j,PetscReal *beta,PetscInt k,PetscScalar *h,PetscBool *breakdown)
479: {
480: PetscInt l,i;
481: PetscScalar alpha,alpha1,alpha2;
482: Vec x,y;
484: PetscFunctionBegin;
485: PetscCall(PetscArrayzero(h,2*j));
487: /* Local orthogonalization */
488: l = j==k+1?0:j-2; /* 1st column to orthogonalize against */
489: /* alpha = X(:,j-1)'*hx-Y(:,j-1)'*hy */
490: PetscCall(BVGetColumn(X,j-1,&x));
491: PetscCall(BVGetColumn(Y,j-1,&y));
492: PetscCall(VecDotBegin(hx,x,&alpha1));
493: PetscCall(VecDotBegin(hy,y,&alpha2));
494: PetscCall(VecDotEnd(hx,x,&alpha1));
495: PetscCall(VecDotEnd(hy,y,&alpha2));
496: PetscCall(BVRestoreColumn(X,j-1,&x));
497: PetscCall(BVRestoreColumn(Y,j-1,&y));
498: alpha = alpha1-alpha2;
499: /* Store coeffs into h */
500: for (i=l; i<j-1; i++) h[i] = h[j+i] = beta[i];
501: h[j-1] = alpha;
502: h[2*j-1] = alpha-1.0;
503: /* Orthogonalize: hx = hx - X(:,l:j-1)*h1 - conj(Y(:,l:j-1))*h2 */
504: PetscCall(BVSetActiveColumns(X,l,j));
505: PetscCall(BVSetActiveColumns(Y,l,j));
506: PetscCall(BVMultVec(X,-1.0,1.0,hx,h+l));
507: PetscCall(VecConjugate(hx));
508: for (i=j+l; i<2*j; i++) h[j+i] = PetscConj(h[i]);
509: PetscCall(BVMultVec(Y,-1.0,1.0,hx,h+2*j+l));
510: PetscCall(VecConjugate(hx));
512: /* Full orthogonalization */
513: PetscCall(VecCopy(hx,hy));
514: PetscCall(VecConjugate(hy));
515: PetscCall(Orthog_ProjectedBSE(hx,hy,X,Y,j,h,h+2*j,breakdown));
516: PetscFunctionReturn(PETSC_SUCCESS);
517: }
519: static PetscErrorCode EPSBSELanczos_ProjectedBSE(EPS eps,BV X,BV Y,Vec v,PetscReal *alpha,PetscReal *beta,PetscInt k,PetscInt *M,PetscBool *breakdown)
520: {
521: PetscInt j,m = *M;
522: Vec u,x,y,w,f,g,vecs[2];
523: Mat H;
524: IS is[2];
525: PetscReal nrm;
526: PetscScalar *hwork,lhwork[100],gamma;
528: PetscFunctionBegin;
529: if (4*m > 100) PetscCall(PetscMalloc1(4*m,&hwork));
530: else hwork = lhwork;
531: PetscCall(STGetMatrix(eps->st,0,&H));
532: PetscCall(MatNestGetISs(H,is,NULL));
534: /* create work vectors */
535: PetscCall(BVGetColumn(Y,0,&u));
536: PetscCall(VecDuplicate(u,&w));
537: vecs[0] = u;
538: vecs[1] = w;
539: PetscCall(VecCreateNest(PetscObjectComm((PetscObject)eps),2,is,vecs,&f));
540: PetscCall(VecCreateNest(PetscObjectComm((PetscObject)eps),2,is,vecs,&g));
541: PetscCall(BVRestoreColumn(Y,0,&u));
543: /* Normalize initial vector */
544: if (k==0) {
545: if (eps->nini==0) PetscCall(BVSetRandomColumn(eps->V,0));
546: PetscCall(BVGetColumn(X,0,&x));
547: /* v = Hmult(u,1) */
548: PetscCall(BVGetColumn(Y,0,&y));
549: PetscCall(VecCopy(x,w));
550: PetscCall(VecConjugate(w));
551: PetscCall(VecNestSetSubVec(f,0,x));
552: PetscCall(VecNestSetSubVec(g,0,y));
553: PetscCall(STApply(eps->st,f,g));
554: /* nrm = sqrt(real(u'v)) */
555: PetscCall(VecDot(y,x,&gamma));
556: nrm = PetscSqrtReal(PetscRealPart(gamma));
557: /* u = u /(nrm*2) */
558: PetscCall(VecScale(x,1.0/(2.0*nrm)));
559: /* v = v /(nrm*2) */
560: PetscCall(VecScale(y,1.0/(2.0*nrm)));
561: PetscCall(VecCopy(y,v));
562: /* X(:,1) = (u+v) */
563: PetscCall(VecAXPY(x,1,y));
564: /* Y(:,1) = conj(u-v) */
565: PetscCall(VecAYPX(y,-2,x));
566: PetscCall(VecConjugate(y));
567: PetscCall(BVRestoreColumn(X,0,&x));
568: PetscCall(BVRestoreColumn(Y,0,&y));
569: }
571: for (j=k;j<m;j++) {
572: /* j+1 columns (indexes 0 to j) have been computed */
573: PetscCall(BVGetColumn(X,j+1,&x));
574: PetscCall(BVGetColumn(Y,j+1,&y));
575: /* u = Hmult(v,-1)*/
576: PetscCall(VecCopy(v,w));
577: PetscCall(VecConjugate(w));
578: PetscCall(VecScale(w,-1.0));
579: PetscCall(VecNestSetSubVec(f,0,v));
580: PetscCall(VecNestSetSubVec(g,0,x));
581: PetscCall(STApply(eps->st,f,g));
582: /* hx = (u+v) */
583: PetscCall(VecCopy(x,y));
584: PetscCall(VecAXPY(x,1,v));
585: /* hy = conj(u-v) */
586: PetscCall(VecAXPY(y,-1,v));
587: PetscCall(VecConjugate(y));
588: /* [u,cd] = orthog(hx,hy,X(:,1:j),Y(:,1:j),opt)*/
589: PetscCall(OrthogonalizeVector_ProjectedBSE(x,y,X,Y,j+1,beta,k,hwork,breakdown));
590: /* alpha(j) = real(cd(j))-1/2 */
591: alpha[j] = 2*(PetscRealPart(hwork[j]) - 0.5);
592: /* v = Hmult(u,1) */
593: PetscCall(VecCopy(x,w));
594: PetscCall(VecConjugate(w));
595: PetscCall(VecNestSetSubVec(f,0,x));
596: PetscCall(VecNestSetSubVec(g,0,y));
597: PetscCall(STApply(eps->st,f,g));
598: /* nrm = sqrt(real(u'*v)) */
599: /* beta(j) = nrm */
600: PetscCall(VecDot(x,y,&gamma));
601: beta[j] = 2.0*PetscSqrtReal(PetscRealPart(gamma));
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: /* Restore */
613: PetscCall(BVRestoreColumn(X,j+1,&x));
614: PetscCall(BVRestoreColumn(Y,j+1,&y));
615: }
616: if (4*m > 100) 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;
671: PetscInt nev=(eps->nev+1)/2;
673: PetscFunctionBegin;
674: PetscCheck((eps->problem_type==EPS_BSE),PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Problem type should be BSE");
675: EPSCheckUnsupportedCondition(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION | EPS_FEATURE_BALANCE,PETSC_TRUE," with BSE structure");
676: PetscCall(EPSSetDimensions_Default(eps,nev,&eps->ncv,&eps->mpd));
677: PetscCheck(eps->ncv<=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);
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 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: PetscCall(EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta,0.0,1.0,&k));
770: PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
771: nconv = k;
773: /* Update l */
774: if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
775: else l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
776: if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
777: if (l) PetscCall(PetscInfo(eps,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
779: if (eps->reason == EPS_CONVERGED_ITERATING) {
780: PetscCheck(!breakdown,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Breakdown in BSE Krylov-Schur (beta=%g)",(double)beta);
781: /* Prepare the Rayleigh quotient for restart */
782: PetscCall(DSTruncate(eps->ds,k+l,PETSC_FALSE));
783: }
784: /* Update the corresponding vectors
785: U(:,idx) = U*Q(:,idx), V(:,idx) = V*Q(:,idx) */
786: PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&Q));
787: PetscCall(BVMultInPlace(U,Q,eps->nconv,k+l));
788: PetscCall(BVMultInPlace(V,Q,eps->nconv,k+l));
789: PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&Q));
791: if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) PetscCall(BVCopyColumn(eps->V,nv,k+l));
792: eps->nconv = k;
793: for (k=0; k<eps->ncv; k++) eps->eigr[k] = PetscSqrtReal(PetscRealPart(eps->eigr[k]));
794: PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv));
795: }
797: eps->nev = nevsave;
799: PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
800: PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
801: PetscFunctionReturn(PETSC_SUCCESS);
802: }
804: /*
805: EPSConvergence_Gruning - convergence check based on SVDKrylovConvergence().
806: */
807: static PetscErrorCode EPSConvergence_Gruning(EPS eps,PetscBool getall,PetscInt kini,PetscInt nits,PetscInt *kout)
808: {
809: PetscInt k,marker,ld;
810: PetscReal *alpha,*beta,resnorm;
811: PetscBool extra;
813: PetscFunctionBegin;
814: *kout = 0;
815: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
816: PetscCall(DSGetExtraRow(eps->ds,&extra));
817: PetscCheck(extra,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Only implemented for DS with extra row");
818: marker = -1;
819: if (eps->trackall) getall = PETSC_TRUE;
820: PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&alpha));
821: beta = alpha + ld;
822: for (k=kini;k<kini+nits;k++) {
823: resnorm = PetscAbsReal(beta[k]);
824: PetscCall((*eps->converged)(eps,eps->eigr[k],eps->eigi[k],resnorm,&eps->errest[k],eps->convergedctx));
825: if (marker==-1 && eps->errest[k] >= eps->tol) marker = k;
826: if (marker!=-1 && !getall) break;
827: }
828: PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&alpha));
829: if (marker!=-1) k = marker;
830: *kout = k;
831: PetscFunctionReturn(PETSC_SUCCESS);
832: }
834: PetscErrorCode EPSSolve_KrylovSchur_BSE_Gruning(EPS eps)
835: {
836: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
837: PetscInt k,l,ld,nv,nconv=0,nevsave;
838: Mat H,Q,Z;
839: BV U,V,HU,HV;
840: IS is[2];
841: PetscReal *d1,*d2,beta;
842: PetscBool breakdown=PETSC_FALSE;
844: PetscFunctionBegin;
845: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
847: /* Extract matrix blocks */
848: PetscCall(STGetMatrix(eps->st,0,&H));
849: PetscCall(MatNestGetISs(H,is,NULL));
851: /* Get the split bases */
852: PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
854: /* Create HU and HV */
855: PetscCall(BVDuplicate(U,&HU));
856: PetscCall(BVDuplicate(V,&HV));
858: nevsave = eps->nev;
859: eps->nev = (eps->nev+1)/2;
860: l = 0;
862: /* Restart loop */
863: while (eps->reason == EPS_CONVERGED_ITERATING) {
864: eps->its++;
866: /* Compute an nv-step Lanczos factorization */
867: nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
868: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
869: PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&d1));
870: d2 = d1 + ld;
871: PetscCall(EPSBSELanczos_Gruning(eps,U,V,HU,HV,d1,d2,eps->nconv+l,&nv,&breakdown));
872: beta = d1[nv-1];
873: PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&d1));
875: /* Compute SVD */
876: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
877: PetscCall(DSSVDSetDimensions(eps->ds,nv));
878: PetscCall(DSSetState(eps->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
879: PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
881: PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
882: PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
883: PetscCall(DSUpdateExtraRow(eps->ds));
884: PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
886: /* Check convergence */
887: PetscCall(EPSConvergence_Gruning(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,&k));
888: PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
889: nconv = k;
891: /* Update l */
892: if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
893: else l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
894: if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
895: if (l) PetscCall(PetscInfo(eps,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
897: if (eps->reason == EPS_CONVERGED_ITERATING) {
898: PetscCheck(!breakdown,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Breakdown in BSE Krylov-Schur (beta=%g)",(double)beta);
899: /* Prepare the Rayleigh quotient for restart */
900: PetscCall(DSTruncate(eps->ds,k+l,PETSC_FALSE));
901: }
902: /* Update the corresponding vectors
903: U(:,idx) = U*Q(:,idx), V(:,idx) = V*Z(:,idx) */
904: PetscCall(DSGetMat(eps->ds,DS_MAT_U,&Q));
905: PetscCall(DSGetMat(eps->ds,DS_MAT_V,&Z));
906: PetscCall(BVMultInPlace(U,Z,eps->nconv,k+l));
907: PetscCall(BVMultInPlace(V,Q,eps->nconv,k+l));
908: PetscCall(BVMultInPlace(HU,Z,eps->nconv,k+l));
909: PetscCall(BVMultInPlace(HV,Q,eps->nconv,k+l));
910: PetscCall(DSRestoreMat(eps->ds,DS_MAT_U,&Q));
911: PetscCall(DSRestoreMat(eps->ds,DS_MAT_V,&Z));
913: if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) PetscCall(BVCopyColumn(U,nv,k+l));
914: eps->nconv = k;
915: PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv));
916: }
918: eps->nev = nevsave;
920: PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
921: PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
923: PetscCall(BVDestroy(&HU));
924: PetscCall(BVDestroy(&HV));
925: PetscFunctionReturn(PETSC_SUCCESS);
926: }
928: PetscErrorCode EPSSolve_KrylovSchur_BSE_ProjectedBSE(EPS eps)
929: {
930: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
931: PetscInt k,l,ld,nv,nconv=0,nevsave;
932: Mat H,Q;
933: Vec v;
934: BV U,V;
935: IS is[2];
936: PetscReal *a,*b,beta;
937: PetscBool breakdown=PETSC_FALSE;
939: PetscFunctionBegin;
940: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
942: /* Extract matrix blocks */
943: PetscCall(STGetMatrix(eps->st,0,&H));
944: PetscCall(MatNestGetISs(H,is,NULL));
946: /* Get the split bases */
947: PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
949: /* Vector v */
950: PetscCall(BVCreateVec(V,&v));
952: nevsave = eps->nev;
953: eps->nev = (eps->nev+1)/2;
954: l = 0;
956: /* Restart loop */
957: while (eps->reason == EPS_CONVERGED_ITERATING) {
958: eps->its++;
960: /* Compute an nv-step Lanczos factorization */
961: nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
962: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
963: PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&a));
964: b = a + ld;
965: PetscCall(EPSBSELanczos_ProjectedBSE(eps,U,V,v,a,b,eps->nconv+l,&nv,&breakdown));
966: beta = b[nv-1];
967: PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&a));
968: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
969: PetscCall(DSSetState(eps->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
970: PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
972: /* Solve projected problem */
973: PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
974: PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
975: PetscCall(DSUpdateExtraRow(eps->ds));
976: PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
978: /* Check convergence */
979: PetscCall(EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta,0.0,1.0,&k));
980: PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
981: nconv = k;
983: /* Update l */
984: if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
985: else l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
986: if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
987: if (l) PetscCall(PetscInfo(eps,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
989: if (eps->reason == EPS_CONVERGED_ITERATING) {
990: PetscCheck(!breakdown,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Breakdown in BSE Krylov-Schur (beta=%g)",(double)beta);
991: /* Prepare the Rayleigh quotient for restart */
992: PetscCall(DSTruncate(eps->ds,k+l,PETSC_FALSE));
993: }
994: /* Update the corresponding vectors
995: U(:,idx) = U*Q(:,idx), V(:,idx) = V*Q(:,idx) */
996: PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&Q));
997: PetscCall(BVMultInPlace(U,Q,eps->nconv,k+l));
998: PetscCall(BVMultInPlace(V,Q,eps->nconv,k+l));
999: PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&Q));
1001: if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) PetscCall(BVCopyColumn(eps->V,nv,k+l));
1002: eps->nconv = k;
1003: /* Obtain eigenvalues */
1004: for (k=0; k<eps->ncv; k++) eps->eigr[k] = PetscSqrtReal(PetscRealPart(eps->eigr[k]));
1005: PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv));
1006: }
1008: eps->nev = nevsave;
1010: PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
1011: PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
1013: PetscCall(VecDestroy(&v));
1014: PetscFunctionReturn(PETSC_SUCCESS);
1015: }