Actual source code: ks-bse.c
slepc-3.22.2 2024-12-02
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) {
295: if (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));
313: }
315: for (j=k;j<m;j++) {
316: /* j+1 columns (indexes 0 to j) have been computed */
317: PetscCall(BVGetColumn(HU,j,&x));
318: PetscCall(BVGetColumn(V,j,&v));
319: PetscCall(BVGetColumn(HV,j,&y));
320: PetscCall(VecCopy(x,v));
321: PetscCall(BVRestoreColumn(HU,j,&x));
322: /* v = Orthogonalize HU(:,j) */
323: PetscCall(OrthogonalizeVector_Gruning(v,U,V,HU,HV,j+1,beta2,k,hwork,PETSC_FALSE,breakdown));
324: /* y = Hmult(v,-1) */
325: PetscCall(VecCopy(v,w));
326: PetscCall(VecConjugate(w));
327: PetscCall(VecScale(w,-1.0));
328: PetscCall(VecNestSetSubVec(f,0,v));
329: PetscCall(VecNestSetSubVec(g,0,y));
330: PetscCall(STApply(eps->st,f,g));
331: /* beta = sqrt(2*real(v'*y)); */
332: PetscCall(VecDot(v,y,&dot));
333: beta1[j] = PetscSqrtReal(PetscRealPart(2*dot)); //FIXME Check beta != 0?
334: /* V(:,j) = v/beta1; */
335: /* HV(:,j) = y/beta1; */
336: PetscCall(VecScale(v,1.0/beta1[j]));
337: PetscCall(VecScale(y,1.0/beta1[j]));
338: PetscCall(BVRestoreColumn(V,j,&v));
339: PetscCall(BVRestoreColumn(HV,j,&y));
341: PetscCall(BVGetColumn(HV,j,&x));
342: PetscCall(BVGetColumn(U,j+1,&v));
343: PetscCall(BVGetColumn(HU,j+1,&y));
344: PetscCall(VecCopy(x,v));
345: PetscCall(BVRestoreColumn(HV,j,&x));
346: /* v = Orthogonalize HV(:,j) */
347: PetscCall(OrthogonalizeVector_Gruning(v,U,V,HU,HV,j+1,&beta1[j],k,hwork,PETSC_TRUE,breakdown));
348: /* y = Hmult(v,1) */
349: PetscCall(VecCopy(v,w));
350: PetscCall(VecConjugate(w));
351: PetscCall(VecNestSetSubVec(f,0,v));
352: PetscCall(VecNestSetSubVec(g,0,y));
353: PetscCall(STApply(eps->st,f,g));
354: /* beta = sqrt(2*real(v'*y)); */
355: PetscCall(VecDot(v,y,&dot));
356: beta2[j] = PetscSqrtReal(PetscRealPart(2*dot)); //FIXME Check beta != 0?
357: /* U(:,j) = v/beta2; */
358: /* HU(:,j) = y/beta2; */
359: PetscCall(VecScale(v,1.0/beta2[j]));
360: PetscCall(VecScale(y,1.0/beta2[j]));
361: PetscCall(BVRestoreColumn(U,j+1,&v));
362: PetscCall(BVRestoreColumn(HU,j+1,&y));
363: }
364: if (4*m > 100) PetscCall(PetscFree(hwork));
365: PetscCall(VecDestroy(&w));
366: PetscCall(VecDestroy(&f));
367: PetscCall(VecDestroy(&g));
368: PetscFunctionReturn(PETSC_SUCCESS);
369: }
371: static PetscErrorCode EPSComputeVectors_BSE_Gruning(EPS eps)
372: {
373: Mat H;
374: Vec u1,v1;
375: BV U,V;
376: IS is[2];
377: PetscInt k;
379: PetscFunctionBegin;
380: PetscCall(STGetMatrix(eps->st,0,&H));
381: PetscCall(MatNestGetISs(H,is,NULL));
382: PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
383: /* approx eigenvector [x1] is [ u1+v1 ]
384: [x2] [conj(u1)-conj(v1)] */
385: for (k=0; k<eps->nconv; k++) {
386: PetscCall(BVGetColumn(U,k,&u1));
387: PetscCall(BVGetColumn(V,k,&v1));
388: /* x1 = u1 + v1 */
389: PetscCall(VecAXPY(u1,1.0,v1));
390: /* x2 = conj(u1) - conj(v1) = conj(u1 - v1) = conj((u1 + v1) - 2*v1) */
391: PetscCall(VecAYPX(v1,-2.0,u1));
392: PetscCall(VecConjugate(v1));
393: PetscCall(BVRestoreColumn(U,k,&u1));
394: PetscCall(BVRestoreColumn(V,k,&v1));
395: }
396: PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
397: /* Normalize eigenvectors */
398: PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
399: PetscCall(BVNormalize(eps->V,NULL));
400: PetscFunctionReturn(PETSC_SUCCESS);
401: }
403: static PetscErrorCode Orthog_ProjectedBSE(Vec hx,Vec hy,BV X,BV Y,PetscInt j,PetscScalar *h,PetscScalar *c,PetscBool *breakdown)
404: {
405: PetscInt i;
406: Mat MX,MY,MXl,MYl;
407: Vec c1,c2,hxl,hyl,hz;
408: PetscScalar *cx1,*cx2;
409: PetscMPIInt len;
411: PetscFunctionBegin;
412: PetscCall(PetscMPIIntCast(j,&len));
413: PetscCall(BVSetActiveColumns(X,0,j));
414: PetscCall(BVSetActiveColumns(Y,0,j));
415: /* BVTDotVec does not exist yet, implemented via MatMult operations */
416: PetscCall(BVGetMat(X,&MX));
417: PetscCall(BVGetMat(Y,&MY));
418: PetscCall(MatDenseGetLocalMatrix(MX,&MXl));
419: PetscCall(MatDenseGetLocalMatrix(MY,&MYl));
420: PetscCall(MatCreateVecs(MXl,&c1,&hyl));
421: PetscCall(MatCreateVecs(MXl,&c2,&hxl));
423: /* c1 = X^* hx - Y^* hy
424: * c2 = -Y^T hx + X^T hy */
426: PetscCall(VecGetLocalVector(hx,hxl));
427: PetscCall(VecGetLocalVector(hy,hyl));
428: PetscCall(VecDuplicate(hx,&hz));
429: /* c1 = -(Y^* hy) */
430: PetscCall(MatMultHermitianTranspose(MYl,hyl,c1));
431: PetscCall(VecScale(c1,-1));
432: /* c1 = c1 + X^* hx */
433: PetscCall(MatMultHermitianTransposeAdd(MXl,hxl,c1,c1));
434: PetscCall(VecGetArray(c1,&cx1));
435: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE,cx1,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)hx)));
436: PetscCall(VecRestoreArray(c1,&cx1));
437: /* c2 = -(Y^T hx) */
438: PetscCall(MatMultTranspose(MYl,hxl,c2));
439: PetscCall(VecScale(c2,-1));
440: /* c2 = c2 + X^T hy */
441: PetscCall(MatMultTransposeAdd(MXl,hyl,c2,c2));
442: PetscCall(VecGetArray(c2,&cx2));
443: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE,cx2,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)hx)));
444: PetscCall(VecRestoreArray(c2,&cx2));
445: PetscCall(VecRestoreLocalVector(hx,hxl));
446: PetscCall(VecRestoreLocalVector(hy,hyl));
448: /* accumulate orthog coeffs into h */
449: PetscCall(VecGetArrayRead(c1,(const PetscScalar**)&cx1));
450: PetscCall(VecGetArrayRead(c2,(const PetscScalar**)&cx2));
451: for (i=0; i<j; i++) h[i] += cx1[i];
452: for (i=0; i<j; i++) h[i+j] += cx2[i];
453: PetscCall(VecRestoreArrayRead(c1,(const PetscScalar**)&cx1));
454: PetscCall(VecRestoreArrayRead(c2,(const PetscScalar**)&cx2));
456: /* u = hx - X c1 - conj(Y) c2 */
458: /* conj(Y) c2 */
459: PetscCall(VecConjugate(c2));
460: PetscCall(VecGetLocalVector(hz,hxl));
461: PetscCall(MatMult(MYl,c2,hxl));
462: PetscCall(VecConjugate(hxl));
463: /* X c1 */
464: PetscCall(MatMultAdd(MXl,c1,hxl,hxl));
465: PetscCall(VecRestoreLocalVector(hz,hxl));
466: PetscCall(VecAXPY(hx,-1,hz));
468: PetscCall(BVRestoreMat(X,&MX));
469: PetscCall(BVRestoreMat(Y,&MY));
470: PetscCall(VecDestroy(&c1));
471: PetscCall(VecDestroy(&c2));
472: PetscCall(VecDestroy(&hxl));
473: PetscCall(VecDestroy(&hyl));
474: PetscCall(VecDestroy(&hz));
475: PetscFunctionReturn(PETSC_SUCCESS);
476: }
478: /* Orthogonalize vector x against first j vectors in X and Y */
479: static PetscErrorCode OrthogonalizeVector_ProjectedBSE(Vec hx,Vec hy,BV X,BV Y,PetscInt j,PetscReal *beta,PetscInt k,PetscScalar *h,PetscBool *breakdown)
480: {
481: PetscInt l,i;
482: PetscScalar alpha,alpha1,alpha2;
483: Vec x,y;
485: PetscFunctionBegin;
486: PetscCall(PetscArrayzero(h,2*j));
488: /* Local orthogonalization */
489: l = j==k+1?0:j-2; /* 1st column to orthogonalize against */
490: /* alpha = X(:,j-1)'*hx-Y(:,j-1)'*hy */
491: PetscCall(BVGetColumn(X,j-1,&x));
492: PetscCall(BVGetColumn(Y,j-1,&y));
493: PetscCall(VecDotBegin(hx,x,&alpha1));
494: PetscCall(VecDotBegin(hy,y,&alpha2));
495: PetscCall(VecDotEnd(hx,x,&alpha1));
496: PetscCall(VecDotEnd(hy,y,&alpha2));
497: PetscCall(BVRestoreColumn(X,j-1,&x));
498: PetscCall(BVRestoreColumn(Y,j-1,&y));
499: alpha = alpha1-alpha2;
500: /* Store coeffs into h */
501: for (i=l; i<j-1; i++) h[i] = h[j+i] = beta[i];
502: h[j-1] = alpha;
503: h[2*j-1] = alpha-1.0;
504: /* Orthogonalize: hx = hx - X(:,l:j-1)*h1 - conj(Y(:,l:j-1))*h2 */
505: PetscCall(BVSetActiveColumns(X,l,j));
506: PetscCall(BVSetActiveColumns(Y,l,j));
507: PetscCall(BVMultVec(X,-1.0,1.0,hx,h+l));
508: PetscCall(VecConjugate(hx));
509: for (i=j+l; i<2*j; i++) h[j+i] = PetscConj(h[i]);
510: PetscCall(BVMultVec(Y,-1.0,1.0,hx,h+2*j+l));
511: PetscCall(VecConjugate(hx));
513: /* Full orthogonalization */
514: PetscCall(VecCopy(hx,hy));
515: PetscCall(VecConjugate(hy));
516: PetscCall(Orthog_ProjectedBSE(hx,hy,X,Y,j,h,h+2*j,breakdown));
517: PetscFunctionReturn(PETSC_SUCCESS);
518: }
520: static PetscErrorCode EPSBSELanczos_ProjectedBSE(EPS eps,BV X,BV Y,Vec v,PetscReal *alpha,PetscReal *beta,PetscInt k,PetscInt *M,PetscBool *breakdown)
521: {
522: PetscInt j,m = *M;
523: Vec u,x,y,w,f,g,vecs[2];
524: Mat H;
525: IS is[2];
526: PetscReal nrm;
527: PetscScalar *hwork,lhwork[100],gamma;
529: PetscFunctionBegin;
530: if (4*m > 100) PetscCall(PetscMalloc1(4*m,&hwork));
531: else hwork = lhwork;
532: PetscCall(STGetMatrix(eps->st,0,&H));
533: PetscCall(MatNestGetISs(H,is,NULL));
535: /* create work vectors */
536: PetscCall(BVGetColumn(Y,0,&u));
537: PetscCall(VecDuplicate(u,&w));
538: vecs[0] = u;
539: vecs[1] = w;
540: PetscCall(VecCreateNest(PetscObjectComm((PetscObject)eps),2,is,vecs,&f));
541: PetscCall(VecCreateNest(PetscObjectComm((PetscObject)eps),2,is,vecs,&g));
542: PetscCall(BVRestoreColumn(Y,0,&u));
544: /* Normalize initial vector */
545: if (k==0) {
546: if (eps->nini==0) PetscCall(BVSetRandomColumn(eps->V,0));
547: PetscCall(BVGetColumn(X,0,&x));
548: /* v = Hmult(u,1) */
549: PetscCall(BVGetColumn(Y,0,&y));
550: PetscCall(VecCopy(x,w));
551: PetscCall(VecConjugate(w));
552: PetscCall(VecNestSetSubVec(f,0,x));
553: PetscCall(VecNestSetSubVec(g,0,y));
554: PetscCall(STApply(eps->st,f,g));
555: /* nrm = sqrt(real(u'v)) */
556: PetscCall(VecDot(y,x,&gamma));
557: nrm = PetscSqrtReal(PetscRealPart(gamma));
558: /* u = u /(nrm*2) */
559: PetscCall(VecScale(x,1.0/(2.0*nrm)));
560: /* v = v /(nrm*2) */
561: PetscCall(VecScale(y,1.0/(2.0*nrm)));
562: PetscCall(VecCopy(y,v));
563: /* X(:,1) = (u+v) */
564: PetscCall(VecAXPY(x,1,y));
565: /* Y(:,1) = conj(u-v) */
566: PetscCall(VecAYPX(y,-2,x));
567: PetscCall(VecConjugate(y));
568: PetscCall(BVRestoreColumn(X,0,&x));
569: PetscCall(BVRestoreColumn(Y,0,&y));
570: }
572: for (j=k;j<m;j++) {
573: /* j+1 columns (indexes 0 to j) have been computed */
574: PetscCall(BVGetColumn(X,j+1,&x));
575: PetscCall(BVGetColumn(Y,j+1,&y));
576: /* u = Hmult(v,-1)*/
577: PetscCall(VecCopy(v,w));
578: PetscCall(VecConjugate(w));
579: PetscCall(VecScale(w,-1.0));
580: PetscCall(VecNestSetSubVec(f,0,v));
581: PetscCall(VecNestSetSubVec(g,0,x));
582: PetscCall(STApply(eps->st,f,g));
583: /* hx = (u+v) */
584: PetscCall(VecCopy(x,y));
585: PetscCall(VecAXPY(x,1,v));
586: /* hy = conj(u-v) */
587: PetscCall(VecAXPY(y,-1,v));
588: PetscCall(VecConjugate(y));
589: /* [u,cd] = orthog(hx,hy,X(:,1:j),Y(:,1:j),opt)*/
590: PetscCall(OrthogonalizeVector_ProjectedBSE(x,y,X,Y,j+1,beta,k,hwork,breakdown));
591: /* alpha(j) = real(cd(j))-1/2 */
592: alpha[j] = 2*(PetscRealPart(hwork[j]) - 0.5);
593: /* v = Hmult(u,1) */
594: PetscCall(VecCopy(x,w));
595: PetscCall(VecConjugate(w));
596: PetscCall(VecNestSetSubVec(f,0,x));
597: PetscCall(VecNestSetSubVec(g,0,y));
598: PetscCall(STApply(eps->st,f,g));
599: /* nrm = sqrt(real(u'*v)) */
600: /* beta(j) = nrm */
601: PetscCall(VecDot(x,y,&gamma));
602: beta[j] = 2.0*PetscSqrtReal(PetscRealPart(gamma));
603: /* u = u/(nrm*2) */
604: PetscCall(VecScale(x,1.0/beta[j]));
605: /* v = v/(nrm*2) */
606: PetscCall(VecScale(y,1.0/beta[j]));
607: PetscCall(VecCopy(y,v));
608: /* X(:,j+1) = (u+v) */
609: PetscCall(VecAXPY(x,1,y));
610: /* Y(:,j+1) = conj(u-v) */
611: PetscCall(VecAYPX(y,-2,x));
612: PetscCall(VecConjugate(y));
613: /* Restore */
614: PetscCall(BVRestoreColumn(X,j+1,&x));
615: PetscCall(BVRestoreColumn(Y,j+1,&y));
616: }
617: if (4*m > 100) PetscCall(PetscFree(hwork));
618: PetscCall(VecDestroy(&w));
619: PetscCall(VecDestroy(&f));
620: PetscCall(VecDestroy(&g));
621: PetscFunctionReturn(PETSC_SUCCESS);
622: }
624: static PetscErrorCode EPSComputeVectors_BSE_ProjectedBSE(EPS eps)
625: {
626: Mat H;
627: Vec x1,y1,cx1,cy1;
628: BV X,Y;
629: IS is[2];
630: PetscInt k;
631: PetscScalar delta1,delta2,lambda;
633: PetscFunctionBegin;
634: PetscCall(STGetMatrix(eps->st,0,&H));
635: PetscCall(MatNestGetISs(H,is,NULL));
636: PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&X,&Y));
637: PetscCall(BVCreateVec(X,&cx1));
638: PetscCall(BVCreateVec(Y,&cy1));
639: for (k=0; k<eps->nconv; k++) {
640: /* approx eigenvector is [ delta1*x1 + delta2*conj(y1) ]
641: [ delta1*y1 + delta2*conj(x1) ] */
642: lambda = eps->eigr[k];
643: PetscCall(STBackTransform(eps->st,1,&lambda,&eps->eigi[k]));
644: delta1 = lambda+1.0;
645: delta2 = lambda-1.0;
646: PetscCall(BVGetColumn(X,k,&x1));
647: PetscCall(BVGetColumn(Y,k,&y1));
648: PetscCall(VecCopy(x1,cx1));
649: PetscCall(VecCopy(y1,cy1));
650: PetscCall(VecConjugate(cx1));
651: PetscCall(VecConjugate(cy1));
652: PetscCall(VecScale(x1,delta1));
653: PetscCall(VecScale(y1,delta1));
654: PetscCall(VecAXPY(x1,delta2,cy1));
655: PetscCall(VecAXPY(y1,delta2,cx1));
656: PetscCall(BVRestoreColumn(X,k,&x1));
657: PetscCall(BVRestoreColumn(Y,k,&y1));
658: }
659: PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&X,&Y));
660: PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
661: /* Normalize eigenvector */
662: PetscCall(BVNormalize(eps->V,NULL));
663: PetscCall(VecDestroy(&cx1));
664: PetscCall(VecDestroy(&cy1));
665: PetscFunctionReturn(PETSC_SUCCESS);
666: }
668: PetscErrorCode EPSSetUp_KrylovSchur_BSE(EPS eps)
669: {
670: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
671: PetscBool flg,sinvert;
672: PetscInt nev=(eps->nev+1)/2;
674: PetscFunctionBegin;
675: PetscCheck((eps->problem_type==EPS_BSE),PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Problem type should be BSE");
676: EPSCheckUnsupportedCondition(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION | EPS_FEATURE_BALANCE,PETSC_TRUE," with BSE structure");
677: PetscCall(EPSSetDimensions_Default(eps,nev,&eps->ncv,&eps->mpd));
678: PetscCheck(eps->ncv<=nev+eps->mpd,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must not be larger than nev+mpd");
679: if (eps->max_it==PETSC_DETERMINE) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
681: PetscCall(PetscObjectTypeCompareAny((PetscObject)eps->st,&flg,STSINVERT,STSHIFT,""));
682: PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Krylov-Schur BSE only supports shift and shift-and-invert ST");
683: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STSINVERT,&sinvert));
684: if (!eps->which) {
685: if (sinvert) eps->which = EPS_TARGET_MAGNITUDE;
686: else eps->which = EPS_SMALLEST_MAGNITUDE;
687: }
689: if (!ctx->keep) ctx->keep = 0.5;
690: PetscCall(STSetStructured(eps->st,PETSC_TRUE));
692: PetscCall(EPSAllocateSolution(eps,1));
693: switch (ctx->bse) {
694: case EPS_KRYLOVSCHUR_BSE_SHAO:
695: eps->ops->solve = EPSSolve_KrylovSchur_BSE_Shao;
696: eps->ops->computevectors = EPSComputeVectors_BSE_Shao;
697: PetscCall(DSSetType(eps->ds,DSHEP));
698: PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
699: PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
700: PetscCall(DSAllocate(eps->ds,eps->ncv+1));
701: break;
702: case EPS_KRYLOVSCHUR_BSE_GRUNING:
703: eps->ops->solve = EPSSolve_KrylovSchur_BSE_Gruning;
704: eps->ops->computevectors = EPSComputeVectors_BSE_Gruning;
705: PetscCall(DSSetType(eps->ds,DSSVD));
706: PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
707: PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
708: PetscCall(DSAllocate(eps->ds,eps->ncv+1));
709: break;
710: case EPS_KRYLOVSCHUR_BSE_PROJECTEDBSE:
711: eps->ops->solve = EPSSolve_KrylovSchur_BSE_ProjectedBSE;
712: eps->ops->computevectors = EPSComputeVectors_BSE_ProjectedBSE;
713: PetscCall(DSSetType(eps->ds,DSHEP));
714: PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
715: PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
716: PetscCall(DSAllocate(eps->ds,eps->ncv+1));
717: break;
718: default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Unexpected error");
719: }
720: PetscFunctionReturn(PETSC_SUCCESS);
721: }
723: PetscErrorCode EPSSolve_KrylovSchur_BSE_Shao(EPS eps)
724: {
725: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
726: PetscInt i,k,l,ld,nv,nconv=0,nevsave;
727: Mat H,Q;
728: BV U,V;
729: IS is[2];
730: PetscReal *a,*b,beta;
731: PetscBool breakdown=PETSC_FALSE;
733: PetscFunctionBegin;
734: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
736: /* Extract matrix blocks */
737: PetscCall(STGetMatrix(eps->st,0,&H));
738: PetscCall(MatNestGetISs(H,is,NULL));
740: /* Get the split bases */
741: PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
743: nevsave = eps->nev;
744: eps->nev = (eps->nev+1)/2;
745: l = 0;
747: /* Restart loop */
748: while (eps->reason == EPS_CONVERGED_ITERATING) {
749: eps->its++;
751: /* Compute an nv-step Lanczos factorization */
752: nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
753: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
754: PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&a));
755: b = a + ld;
756: PetscCall(EPSBSELanczos_Shao(eps,U,V,a,b,eps->nconv+l,&nv,&breakdown));
757: beta = b[nv-1];
758: PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&a));
759: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
760: PetscCall(DSSetState(eps->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
761: PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
763: /* Solve projected problem */
764: PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
765: PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
766: PetscCall(DSUpdateExtraRow(eps->ds));
767: PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
769: /* Check convergence */
770: for (i=0;i<eps->ncv;i++) eps->eigr[i] = PetscSqrtReal(PetscRealPart(eps->eigr[i]));
771: PetscCall(EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta,0.0,1.0,&k));
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) PetscCall(BVCopyColumn(eps->V,nv,k+l));
794: eps->nconv = k;
795: PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv));
796: }
798: eps->nev = nevsave;
800: PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
801: PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
802: PetscFunctionReturn(PETSC_SUCCESS);
803: }
805: /*
806: EPSConvergence_Gruning - convergence check based on SVDKrylovConvergence().
807: */
808: static PetscErrorCode EPSConvergence_Gruning(EPS eps,PetscBool getall,PetscInt kini,PetscInt nits,PetscInt *kout)
809: {
810: PetscInt k,marker,ld;
811: PetscReal *alpha,*beta,resnorm;
812: PetscBool extra;
814: PetscFunctionBegin;
815: *kout = 0;
816: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
817: PetscCall(DSGetExtraRow(eps->ds,&extra));
818: PetscCheck(extra,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Only implemented for DS with extra row");
819: marker = -1;
820: if (eps->trackall) getall = PETSC_TRUE;
821: PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&alpha));
822: beta = alpha + ld;
823: for (k=kini;k<kini+nits;k++) {
824: resnorm = PetscAbsReal(beta[k]);
825: PetscCall((*eps->converged)(eps,eps->eigr[k],eps->eigi[k],resnorm,&eps->errest[k],eps->convergedctx));
826: if (marker==-1 && eps->errest[k] >= eps->tol) marker = k;
827: if (marker!=-1 && !getall) break;
828: }
829: PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&alpha));
830: if (marker!=-1) k = marker;
831: *kout = k;
832: PetscFunctionReturn(PETSC_SUCCESS);
833: }
835: PetscErrorCode EPSSolve_KrylovSchur_BSE_Gruning(EPS eps)
836: {
837: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
838: PetscInt k,l,ld,nv,nconv=0,nevsave;
839: Mat H,Q,Z;
840: BV U,V,HU,HV;
841: IS is[2];
842: PetscReal *d1,*d2,beta;
843: PetscBool breakdown=PETSC_FALSE;
845: PetscFunctionBegin;
846: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
848: /* Extract matrix blocks */
849: PetscCall(STGetMatrix(eps->st,0,&H));
850: PetscCall(MatNestGetISs(H,is,NULL));
852: /* Get the split bases */
853: PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
855: /* Create HU and HV */
856: PetscCall(BVDuplicate(U,&HU));
857: PetscCall(BVDuplicate(V,&HV));
859: nevsave = eps->nev;
860: eps->nev = (eps->nev+1)/2;
861: l = 0;
863: /* Restart loop */
864: while (eps->reason == EPS_CONVERGED_ITERATING) {
865: eps->its++;
867: /* Compute an nv-step Lanczos factorization */
868: nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
869: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
870: PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&d1));
871: d2 = d1 + ld;
872: PetscCall(EPSBSELanczos_Gruning(eps,U,V,HU,HV,d1,d2,eps->nconv+l,&nv,&breakdown));
873: beta = d1[nv-1];
874: PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&d1));
876: /* Compute SVD */
877: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
878: PetscCall(DSSVDSetDimensions(eps->ds,nv));
879: PetscCall(DSSetState(eps->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
880: PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
882: PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
883: PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
884: PetscCall(DSUpdateExtraRow(eps->ds));
885: PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
887: /* Check convergence */
888: PetscCall(EPSConvergence_Gruning(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,&k));
889: PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
890: nconv = k;
892: /* Update l */
893: if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
894: else l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
895: if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
896: if (l) PetscCall(PetscInfo(eps,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
898: if (eps->reason == EPS_CONVERGED_ITERATING) {
899: PetscCheck(!breakdown,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Breakdown in BSE Krylov-Schur (beta=%g)",(double)beta);
900: /* Prepare the Rayleigh quotient for restart */
901: PetscCall(DSTruncate(eps->ds,k+l,PETSC_FALSE));
902: }
903: /* Update the corresponding vectors
904: U(:,idx) = U*Q(:,idx), V(:,idx) = V*Z(:,idx) */
905: PetscCall(DSGetMat(eps->ds,DS_MAT_U,&Q));
906: PetscCall(DSGetMat(eps->ds,DS_MAT_V,&Z));
907: PetscCall(BVMultInPlace(U,Z,eps->nconv,k+l));
908: PetscCall(BVMultInPlace(V,Q,eps->nconv,k+l));
909: PetscCall(BVMultInPlace(HU,Z,eps->nconv,k+l));
910: PetscCall(BVMultInPlace(HV,Q,eps->nconv,k+l));
911: PetscCall(DSRestoreMat(eps->ds,DS_MAT_U,&Q));
912: PetscCall(DSRestoreMat(eps->ds,DS_MAT_V,&Z));
914: if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
915: PetscCall(BVCopyColumn(U,nv,k+l));
916: PetscCall(BVCopyColumn(HU,nv,k+l));
917: }
918: eps->nconv = k;
919: PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv));
920: }
922: eps->nev = nevsave;
924: PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
925: PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
927: PetscCall(BVDestroy(&HU));
928: PetscCall(BVDestroy(&HV));
929: PetscFunctionReturn(PETSC_SUCCESS);
930: }
932: PetscErrorCode EPSSolve_KrylovSchur_BSE_ProjectedBSE(EPS eps)
933: {
934: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
935: PetscInt i,k,l,ld,nv,nconv=0,nevsave;
936: Mat H,Q;
937: Vec v;
938: BV U,V;
939: IS is[2];
940: PetscReal *a,*b,beta;
941: PetscBool breakdown=PETSC_FALSE;
943: PetscFunctionBegin;
944: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
946: /* Extract matrix blocks */
947: PetscCall(STGetMatrix(eps->st,0,&H));
948: PetscCall(MatNestGetISs(H,is,NULL));
950: /* Get the split bases */
951: PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
953: /* Vector v */
954: PetscCall(BVCreateVec(V,&v));
956: nevsave = eps->nev;
957: eps->nev = (eps->nev+1)/2;
958: l = 0;
960: /* Restart loop */
961: while (eps->reason == EPS_CONVERGED_ITERATING) {
962: eps->its++;
964: /* Compute an nv-step Lanczos factorization */
965: nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
966: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
967: PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&a));
968: b = a + ld;
969: PetscCall(EPSBSELanczos_ProjectedBSE(eps,U,V,v,a,b,eps->nconv+l,&nv,&breakdown));
970: beta = b[nv-1];
971: PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&a));
972: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
973: PetscCall(DSSetState(eps->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
974: PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
976: /* Solve projected problem */
977: PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
978: PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
979: PetscCall(DSUpdateExtraRow(eps->ds));
980: PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
982: /* Check convergence */
983: for (i=0;i<eps->ncv;i++) eps->eigr[i] = PetscSqrtReal(PetscRealPart(eps->eigr[i]));
984: PetscCall(EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta,0.0,1.0,&k));
985: PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
986: nconv = k;
988: /* Update l */
989: if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
990: else l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
991: if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
992: if (l) PetscCall(PetscInfo(eps,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
994: if (eps->reason == EPS_CONVERGED_ITERATING) {
995: PetscCheck(!breakdown,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Breakdown in BSE Krylov-Schur (beta=%g)",(double)beta);
996: /* Prepare the Rayleigh quotient for restart */
997: PetscCall(DSTruncate(eps->ds,k+l,PETSC_FALSE));
998: }
999: /* Update the corresponding vectors
1000: U(:,idx) = U*Q(:,idx), V(:,idx) = V*Q(:,idx) */
1001: PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&Q));
1002: PetscCall(BVMultInPlace(U,Q,eps->nconv,k+l));
1003: PetscCall(BVMultInPlace(V,Q,eps->nconv,k+l));
1004: PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&Q));
1006: if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) PetscCall(BVCopyColumn(eps->V,nv,k+l));
1007: eps->nconv = k;
1008: PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv));
1009: }
1011: eps->nev = nevsave;
1013: PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
1014: PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
1016: PetscCall(VecDestroy(&v));
1017: PetscFunctionReturn(PETSC_SUCCESS);
1018: }