Actual source code: ks-bse.c

slepc-3.23.1 2025-05-01
Report Typos and Errors
  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,PetscBool *breakdown)
 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,PetscBool *breakdown)
 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,breakdown));
 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:   PetscContainer container;
105:   SlepcMatStruct mctx;

107:   PetscFunctionBegin;
108:   if (4*m > 100) PetscCall(PetscMalloc1(4*m,&hwork));
109:   else hwork = lhwork;
110:   PetscCall(STGetMatrix(eps->st,0,&H));
111:   PetscCall(MatNestGetISs(H,is,NULL));
112:   PetscCall(PetscObjectQuery((PetscObject)H,"SlepcMatStruct",(PetscObject*)&container));
113:   PetscCall(PetscContainerGetPointer(container,(void**)&mctx));

115:   /* create work vectors */
116:   PetscCall(BVGetColumn(V,0,&v));
117:   PetscCall(VecDuplicate(v,&w));
118:   vecs[0] = v;
119:   vecs[1] = w;
120:   PetscCall(VecCreateNest(PetscObjectComm((PetscObject)eps),2,is,vecs,&f));
121:   PetscCall(VecCreateNest(PetscObjectComm((PetscObject)eps),2,is,vecs,&g));
122:   PetscCall(BVRestoreColumn(V,0,&v));

124:   /* Normalize initial vector */
125:   if (k==0) {
126:     if (eps->nini==0) PetscCall(BVSetRandomColumn(eps->V,0));
127:     PetscCall(BVGetColumn(U,0,&x));
128:     PetscCall(BVGetColumn(V,0,&y));
129:     PetscCall(VecNestSetSubVec(f,0,x));
130:     PetscCall(VecNestSetSubVec(g,0,y));
131:     mctx->s = 1.0;
132:     PetscCall(STApply(eps->st,f,g));
133:     PetscCall(VecDot(y,x,&gamma));
134:     nrm = PetscSqrtReal(PetscRealPart(gamma));
135:     PetscCall(VecScale(x,1.0/nrm));
136:     PetscCall(VecScale(y,1.0/nrm));
137:     PetscCall(BVRestoreColumn(U,0,&x));
138:     PetscCall(BVRestoreColumn(V,0,&y));
139:   }

141:   for (j=k;j<m;j++) {
142:     /* j+1 columns (indexes 0 to j) have been computed */
143:     PetscCall(BVGetColumn(V,j,&v));
144:     PetscCall(BVGetColumn(U,j+1,&x));
145:     PetscCall(BVGetColumn(V,j+1,&y));
146:     PetscCall(VecNestSetSubVec(f,0,v));
147:     PetscCall(VecNestSetSubVec(g,0,x));
148:     mctx->s = -1.0;
149:     PetscCall(STApply(eps->st,f,g));
150:     PetscCall(OrthogonalizeVector_Shao(x,U,V,j+1,v,beta,k,hwork,breakdown));
151:     alpha[j] = PetscRealPart(hwork[j]);
152:     PetscCall(VecNestSetSubVec(f,0,x));
153:     PetscCall(VecNestSetSubVec(g,0,y));
154:     mctx->s = 1.0;
155:     PetscCall(STApply(eps->st,f,g));
156:     PetscCall(VecDot(x,y,&gamma));
157:     beta[j] = PetscSqrtReal(PetscRealPart(gamma));
158:     PetscCall(VecScale(x,1.0/beta[j]));
159:     PetscCall(VecScale(y,1.0/beta[j]));
160:     PetscCall(BVRestoreColumn(V,j,&v));
161:     PetscCall(BVRestoreColumn(U,j+1,&x));
162:     PetscCall(BVRestoreColumn(V,j+1,&y));
163:   }
164:   if (4*m > 100) PetscCall(PetscFree(hwork));
165:   PetscCall(VecDestroy(&w));
166:   PetscCall(VecDestroy(&f));
167:   PetscCall(VecDestroy(&g));
168:   PetscFunctionReturn(PETSC_SUCCESS);
169: }

171: static PetscErrorCode EPSComputeVectors_BSE_Shao(EPS eps)
172: {
173:   Mat         H;
174:   Vec         u1,v1;
175:   BV          U,V;
176:   IS          is[2];
177:   PetscInt    k;
178:   PetscScalar lambda;

180:   PetscFunctionBegin;
181:   PetscCall(STGetMatrix(eps->st,0,&H));
182:   PetscCall(MatNestGetISs(H,is,NULL));
183:   PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
184:   for (k=0; k<eps->nconv; k++) {
185:     PetscCall(BVGetColumn(U,k,&u1));
186:     PetscCall(BVGetColumn(V,k,&v1));
187:     /* approx eigenvector is [    (eigr[k]*u1+v1)]
188:                              [conj(eigr[k]*u1-v1)]  */
189:     lambda = eps->eigr[k];
190:     PetscCall(STBackTransform(eps->st,1,&lambda,&eps->eigi[k]));
191:     PetscCall(VecAYPX(u1,lambda,v1));
192:     PetscCall(VecAYPX(v1,-2.0,u1));
193:     PetscCall(VecConjugate(v1));
194:     PetscCall(BVRestoreColumn(U,k,&u1));
195:     PetscCall(BVRestoreColumn(V,k,&v1));
196:   }
197:   PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
198:   /* Normalize eigenvectors */
199:   PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
200:   PetscCall(BVNormalize(eps->V,NULL));
201:   PetscFunctionReturn(PETSC_SUCCESS);
202: }

204: static PetscErrorCode Orthog_Gruning(Vec x,BV U,BV V,BV HU,BV HV,PetscInt j,PetscScalar *h,PetscScalar *c,PetscBool *breakdown)
205: {
206:   PetscInt i;

208:   PetscFunctionBegin;
209:   PetscCall(BVSetActiveColumns(U,0,j));
210:   PetscCall(BVSetActiveColumns(HU,0,j));
211:   PetscCall(BVSetActiveColumns(V,0,j));
212:   PetscCall(BVSetActiveColumns(HV,0,j));
213: #if defined(PETSC_USE_COMPLEX)
214:   PetscCall(BVDotVecBegin(HU,x,c));
215:   PetscCall(BVDotVecBegin(HV,x,c+j));
216:   PetscCall(BVDotVecEnd(HU,x,c));
217:   PetscCall(BVDotVecEnd(HV,x,c+j));
218: #else
219:   PetscCall(BVDotVec(HU,x,c));
220: #endif
221:   for (i=0; i<j; i++) {
222:     /* c1 = 2*real(HU^* x) ; c2 = 2*imag(HV^* x)*1i */
223: #if defined(PETSC_USE_COMPLEX)
224:     c[i] = PetscRealPart(c[i]);
225:     c[j+i] = PetscCMPLX(0.0,PetscImaginaryPart(c[j+i]));
226: #else
227:     c[j+i] = 0.0;
228: #endif
229:   }
230:   /* x = x-U*c1-V*c2 */
231: #if defined(PETSC_USE_COMPLEX)
232:   PetscCall(BVMultVec(U,-2.0,1.0,x,c));
233:   PetscCall(BVMultVec(V,-2.0,1.0,x,c+j));
234: #else
235:   PetscCall(BVMultVec(U,-2.0,1.0,x,c));
236: #endif
237:   /* accumulate orthog coeffs into h */
238:   for (i=0; i<2*j; i++) h[i] += 2*c[i];
239:   PetscFunctionReturn(PETSC_SUCCESS);
240: }

242: /* Orthogonalize vector x against first j vectors in U and V */
243: 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)
244: {
245:   PetscInt l,i;
246:   Vec      u;

248:   PetscFunctionBegin;
249:   PetscCall(PetscArrayzero(h,4*j));

251:   if (s) {
252:     /* Local orthogonalization */
253:     PetscCall(BVGetColumn(U,j-1,&u));
254:     PetscCall(VecAXPY(x,-*beta,u));
255:     PetscCall(BVRestoreColumn(U,j-1,&u));
256:     h[j-1] = *beta;
257:     /* Full orthogonalization */
258:     PetscCall(Orthog_Gruning(x,U,V,HU,HV,j,h,h+2*j,breakdown));
259:   } else {
260:     /* Local orthogonalization */
261:     l = j==k+1?0:j-2;  /* 1st column to orthogonalize against */
262:     for (i=l; i<j-1; i++) h[j+i] = beta[i];
263:     /* x = x-V(:,l:j-2)*h(l:j-2) */
264:     PetscCall(BVSetActiveColumns(V,l,j-1));
265:     PetscCall(BVMultVec(V,-1.0,1.0,x,h+j+l));
266:   }
267:   PetscFunctionReturn(PETSC_SUCCESS);
268: }

270: static PetscErrorCode EPSBSELanczos_Gruning(EPS eps,BV U,BV V,BV HU,BV HV,PetscReal *beta1,PetscReal *beta2,PetscInt k,PetscInt *M,PetscBool *breakdown)
271: {
272:   PetscInt       j,m = *M;
273:   Vec            v,x,y,w,f,g,vecs[2];
274:   Mat            H;
275:   IS             is[2];
276:   PetscReal      nrm;
277:   PetscScalar    *hwork,lhwork[100],dot;
278:   PetscContainer container;
279:   SlepcMatStruct mctx;

281:   PetscFunctionBegin;
282:   if (4*m > 100) PetscCall(PetscMalloc1(4*m,&hwork));
283:   else hwork = lhwork;
284:   PetscCall(STGetMatrix(eps->st,0,&H));
285:   PetscCall(MatNestGetISs(H,is,NULL));
286:   PetscCall(PetscObjectQuery((PetscObject)H,"SlepcMatStruct",(PetscObject*)&container));
287:   PetscCall(PetscContainerGetPointer(container,(void**)&mctx));

289:   /* create work vectors */
290:   PetscCall(BVGetColumn(V,0,&v));
291:   PetscCall(VecDuplicate(v,&w));
292:   vecs[0] = v;
293:   vecs[1] = w;
294:   PetscCall(VecCreateNest(PetscObjectComm((PetscObject)eps),2,is,vecs,&f));
295:   PetscCall(VecCreateNest(PetscObjectComm((PetscObject)eps),2,is,vecs,&g));
296:   PetscCall(BVRestoreColumn(V,0,&v));

298:   /* Normalize initial vector */
299:   if (k==0) {
300:     if (eps->nini==0) PetscCall(BVSetRandomColumn(eps->V,0));
301:     /* y = Hmult(v1,1) */
302:     PetscCall(BVGetColumn(U,k,&x));
303:     PetscCall(BVGetColumn(HU,k,&y));
304:     PetscCall(VecNestSetSubVec(f,0,x));
305:     PetscCall(VecNestSetSubVec(g,0,y));
306:     mctx->s = 1.0;
307:     PetscCall(STApply(eps->st,f,g));
308:     /* nrm = sqrt(2*real(u1'*y)); */
309:     PetscCall(VecDot(x,y,&dot));
310:     nrm = PetscSqrtReal(PetscRealPart(2*dot));
311:     /* U(:,j) = u1/nrm; */
312:     /* HU(:,j) = y/nrm; */
313:     PetscCall(VecScale(x,1.0/nrm));
314:     PetscCall(VecScale(y,1.0/nrm));
315:     PetscCall(BVRestoreColumn(U,k,&x));
316:     PetscCall(BVRestoreColumn(HU,k,&y));
317:   }

319:   for (j=k;j<m;j++) {
320:     /* j+1 columns (indexes 0 to j) have been computed */
321:     PetscCall(BVGetColumn(HU,j,&x));
322:     PetscCall(BVGetColumn(V,j,&v));
323:     PetscCall(BVGetColumn(HV,j,&y));
324:     PetscCall(VecCopy(x,v));
325:     PetscCall(BVRestoreColumn(HU,j,&x));
326:     /* v = Orthogonalize HU(:,j) */
327:     PetscCall(OrthogonalizeVector_Gruning(v,U,V,HU,HV,j+1,beta2,k,hwork,PETSC_FALSE,breakdown));
328:     /* y = Hmult(v,-1) */
329:     PetscCall(VecNestSetSubVec(f,0,v));
330:     PetscCall(VecNestSetSubVec(g,0,y));
331:     mctx->s = -1.0;
332:     PetscCall(STApply(eps->st,f,g));
333:     /* beta = sqrt(2*real(v'*y)); */
334:     PetscCall(VecDot(v,y,&dot));
335:     beta1[j] = PetscSqrtReal(PetscRealPart(2*dot)); //FIXME Check beta != 0?
336:     /* V(:,j) = v/beta1; */
337:     /* HV(:,j) = y/beta1; */
338:     PetscCall(VecScale(v,1.0/beta1[j]));
339:     PetscCall(VecScale(y,1.0/beta1[j]));
340:     PetscCall(BVRestoreColumn(V,j,&v));
341:     PetscCall(BVRestoreColumn(HV,j,&y));

343:     PetscCall(BVGetColumn(HV,j,&x));
344:     PetscCall(BVGetColumn(U,j+1,&v));
345:     PetscCall(BVGetColumn(HU,j+1,&y));
346:     PetscCall(VecCopy(x,v));
347:     PetscCall(BVRestoreColumn(HV,j,&x));
348:     /* v = Orthogonalize HV(:,j) */
349:     PetscCall(OrthogonalizeVector_Gruning(v,U,V,HU,HV,j+1,&beta1[j],k,hwork,PETSC_TRUE,breakdown));
350:     /* y = Hmult(v,1) */
351:     PetscCall(VecNestSetSubVec(f,0,v));
352:     PetscCall(VecNestSetSubVec(g,0,y));
353:     mctx->s = 1.0;
354:     PetscCall(STApply(eps->st,f,g));
355:     /* beta = sqrt(2*real(v'*y)); */
356:     PetscCall(VecDot(v,y,&dot));
357:     beta2[j] = PetscSqrtReal(PetscRealPart(2*dot)); //FIXME Check beta != 0?
358:     /* U(:,j) = v/beta2; */
359:     /* HU(:,j) = y/beta2; */
360:     PetscCall(VecScale(v,1.0/beta2[j]));
361:     PetscCall(VecScale(y,1.0/beta2[j]));
362:     PetscCall(BVRestoreColumn(U,j+1,&v));
363:     PetscCall(BVRestoreColumn(HU,j+1,&y));
364:   }
365:   if (4*m > 100) PetscCall(PetscFree(hwork));
366:   PetscCall(VecDestroy(&w));
367:   PetscCall(VecDestroy(&f));
368:   PetscCall(VecDestroy(&g));
369:   PetscFunctionReturn(PETSC_SUCCESS);
370: }

372: static PetscErrorCode EPSComputeVectors_BSE_Gruning(EPS eps)
373: {
374:   Mat         H;
375:   Vec         u1,v1;
376:   BV          U,V;
377:   IS          is[2];
378:   PetscInt    k;

380:   PetscFunctionBegin;
381:   PetscCall(STGetMatrix(eps->st,0,&H));
382:   PetscCall(MatNestGetISs(H,is,NULL));
383:   PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
384:   /* approx eigenvector [x1] is [     u1+v1       ]
385:                         [x2]    [conj(u1)-conj(v1)]  */
386:   for (k=0; k<eps->nconv; k++) {
387:     PetscCall(BVGetColumn(U,k,&u1));
388:     PetscCall(BVGetColumn(V,k,&v1));
389:     /* x1 = u1 + v1 */
390:     PetscCall(VecAXPY(u1,1.0,v1));
391:     /* x2 = conj(u1) - conj(v1) = conj(u1 - v1) = conj((u1 + v1) - 2*v1) */
392:     PetscCall(VecAYPX(v1,-2.0,u1));
393:     PetscCall(VecConjugate(v1));
394:     PetscCall(BVRestoreColumn(U,k,&u1));
395:     PetscCall(BVRestoreColumn(V,k,&v1));
396:   }
397:   PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
398:   /* Normalize eigenvectors */
399:   PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
400:   PetscCall(BVNormalize(eps->V,NULL));
401:   PetscFunctionReturn(PETSC_SUCCESS);
402: }

404: /* Full orthogonalization of vector [hx, conj(hx)] against first j vectors in X and Y */
405: static PetscErrorCode Orthog_ProjectedBSE(Vec hx,BV X,BV Y,PetscInt j,PetscScalar *h,PetscScalar *c,PetscBool *breakdown)
406: {
407:   PetscInt i;

409:   PetscFunctionBegin;
410:   PetscCall(BVSetActiveColumns(X,0,j));
411:   PetscCall(BVSetActiveColumns(Y,0,j));
412:   /* c1 = X^* hx */
413:   PetscCall(BVDotVec(X,hx,c));
414:   /* c2 = Y^* conj(hx) */
415:   PetscCall(VecConjugate(hx));
416:   PetscCall(BVDotVec(Y,hx,c+j));
417:   /* c = c1 - c2 */
418:   for (i=0;i<j;i++) c[i] -= c[i+j];
419:   /* hx = hx - conj(Y*c) */
420:   PetscCall(BVMultVec(Y,-1.0,1.0,hx,c));
421:   PetscCall(VecConjugate(hx));
422:   /* hx = hx - X*c */
423:   PetscCall(BVMultVec(X,-1.0,1.0,hx,c));
424:   /* accumulate orthog coeffs into h */
425:   for (i=0;i<j;i++) h[i] += c[i];
426:   PetscFunctionReturn(PETSC_SUCCESS);
427: }

429: /* Orthogonalize vector [hx; hy] against first j vectors in X and Y
430:    The result is a vector [u; conj(u)]. Vector hx is overwritten with u. */
431: static PetscErrorCode OrthogonalizeVector_ProjectedBSE(Vec hx,Vec hy,BV X,BV Y,PetscInt j,PetscReal *beta,PetscInt k,PetscScalar *h,PetscBool *breakdown)
432: {
433:   PetscInt    l,i;
434:   PetscScalar alpha,alpha1,alpha2;
435:   Vec         x,y;

437:   PetscFunctionBegin;
438:   PetscCall(PetscArrayzero(h,2*j));

440:   /* Local orthogonalization */
441:   l = j==k+1?0:j-2;  /* 1st column to orthogonalize against */
442:   /* alpha = X(:,j-1)'*hx-Y(:,j-1)'*hy */
443:   PetscCall(BVGetColumn(X,j-1,&x));
444:   PetscCall(BVGetColumn(Y,j-1,&y));
445:   PetscCall(VecDotBegin(hx,x,&alpha1));
446:   PetscCall(VecDotBegin(hy,y,&alpha2));
447:   PetscCall(VecDotEnd(hx,x,&alpha1));
448:   PetscCall(VecDotEnd(hy,y,&alpha2));
449:   PetscCall(BVRestoreColumn(X,j-1,&x));
450:   PetscCall(BVRestoreColumn(Y,j-1,&y));
451:   alpha = alpha1-alpha2;
452:   /* Store coeffs into h */
453:   for (i=l; i<j-1; i++) h[i] = h[j+i] = beta[i]/2.0;
454:   h[j-1] = alpha;

456:   /* Orthogonalize: hx = hx - X(:,l:j-1)*h1 - conj(Y(:,l:j-1))*h2 */
457:   /* hx = hx - X(:,l:j-1)*h1 */
458:   PetscCall(BVSetActiveColumns(X,l,j));
459:   PetscCall(BVSetActiveColumns(Y,l,j));
460:   PetscCall(BVMultVec(X,-1.0,1.0,hx,h+l));
461:   /* hx = conj(hx) */
462:   PetscCall(VecConjugate(hx));
463:   /* hx = hx - Y(:,l:j-1)*conj(h2) */
464:   h[2*j-1] = PetscConj(alpha-1.0);
465:   PetscCall(BVMultVec(Y,-1.0,1.0,hx,h+j+l));
466:   h[2*j-1] = alpha-1.0;
467:   /* hx = conj(hx) */
468:   PetscCall(VecConjugate(hx));

470:   /* Full orthogonalization */
471:   PetscCall(Orthog_ProjectedBSE(hx,X,Y,j,h,h+2*j,breakdown));
472:   PetscFunctionReturn(PETSC_SUCCESS);
473: }

475: static PetscErrorCode EPSBSELanczos_ProjectedBSE(EPS eps,BV X,BV Y,Vec v,PetscReal *alpha,PetscReal *beta,PetscInt k,PetscInt *M,PetscBool *breakdown)
476: {
477:   PetscInt       j,m = *M;
478:   Vec            u,x,y,w,f,g,vecs[2];
479:   Mat            H;
480:   IS             is[2];
481:   PetscReal      nrm;
482:   PetscScalar    *hwork,lhwork[100],gamma;
483:   PetscContainer container;
484:   SlepcMatStruct mctx;

486:   PetscFunctionBegin;
487:   if (4*m > 100) PetscCall(PetscMalloc1(4*m,&hwork));
488:   else hwork = lhwork;
489:   PetscCall(STGetMatrix(eps->st,0,&H));
490:   PetscCall(MatNestGetISs(H,is,NULL));
491:   PetscCall(PetscObjectQuery((PetscObject)H,"SlepcMatStruct",(PetscObject*)&container));
492:   PetscCall(PetscContainerGetPointer(container,(void**)&mctx));

494:   /* create work vectors */
495:   PetscCall(BVGetColumn(Y,0,&u));
496:   PetscCall(VecDuplicate(u,&w));
497:   vecs[0] = u;
498:   vecs[1] = w;
499:   PetscCall(VecCreateNest(PetscObjectComm((PetscObject)eps),2,is,vecs,&f));
500:   PetscCall(VecCreateNest(PetscObjectComm((PetscObject)eps),2,is,vecs,&g));
501:   PetscCall(BVRestoreColumn(Y,0,&u));

503:   /* Normalize initial vector */
504:   if (k==0) {
505:     if (eps->nini==0) PetscCall(BVSetRandomColumn(eps->V,0));
506:     PetscCall(BVGetColumn(X,0,&x));
507:     /* v = Hmult(u,1) */
508:     PetscCall(BVGetColumn(Y,0,&y));
509:     PetscCall(VecNestSetSubVec(f,0,x));
510:     PetscCall(VecNestSetSubVec(g,0,y));
511:     mctx->s = 1.0;
512:     PetscCall(STApply(eps->st,f,g));
513:     /* nrm = sqrt(real(u'*v)) */
514:     PetscCall(VecDot(y,x,&gamma));
515:     nrm = PetscSqrtReal(PetscRealPart(gamma));
516:     /* u = u /(nrm*2) */
517:     PetscCall(VecScale(x,1.0/(2.0*nrm)));
518:     /* v = v /(nrm*2) */
519:     PetscCall(VecScale(y,1.0/(2.0*nrm)));
520:     PetscCall(VecCopy(y,v));
521:     /* X(:,1) = (u+v) */
522:     PetscCall(VecAXPY(x,1,y));
523:     /* Y(:,1) = conj(u-v) */
524:     PetscCall(VecAYPX(y,-2,x));
525:     PetscCall(VecConjugate(y));
526:     PetscCall(BVRestoreColumn(X,0,&x));
527:     PetscCall(BVRestoreColumn(Y,0,&y));
528:   }

530:   for (j=k;j<m;j++) {
531:     /* j+1 columns (indexes 0 to j) have been computed */
532:     PetscCall(BVGetColumn(X,j+1,&x));
533:     PetscCall(BVGetColumn(Y,j+1,&y));
534:     /* u = Hmult(v,-1)*/
535:     PetscCall(VecNestSetSubVec(f,0,v));
536:     PetscCall(VecNestSetSubVec(g,0,x));
537:     mctx->s = -1.0;
538:     PetscCall(STApply(eps->st,f,g));
539:     /* hx = (u+v) */
540:     PetscCall(VecCopy(x,y));
541:     PetscCall(VecAXPY(x,1,v));
542:     /* hy = conj(u-v) */
543:     PetscCall(VecAXPY(y,-1,v));
544:     PetscCall(VecConjugate(y));
545:     /* [u,cd] = orthog(hx,hy,X(:,1:j),Y(:,1:j),opt)*/
546:     PetscCall(OrthogonalizeVector_ProjectedBSE(x,y,X,Y,j+1,beta,k,hwork,breakdown));
547:     /* alpha(j) = 2*(real(cd(j))-1/2) */
548:     alpha[j] = 2*(PetscRealPart(hwork[j]) - 0.5);
549:     /* v = Hmult(u,1) */
550:     PetscCall(VecNestSetSubVec(f,0,x));
551:     PetscCall(VecNestSetSubVec(g,0,y));
552:     mctx->s = 1.0;
553:     PetscCall(STApply(eps->st,f,g));
554:     /* nrm = sqrt(real(u'*v)) */
555:     /* beta(j) = 2*nrm */
556:     PetscCall(VecDot(x,y,&gamma));
557:     beta[j] = 2.0*PetscSqrtReal(PetscRealPart(gamma));
558:     /* u = u/(nrm*2) */
559:     PetscCall(VecScale(x,1.0/beta[j]));
560:     /* v = v/(nrm*2) */
561:     PetscCall(VecScale(y,1.0/beta[j]));
562:     PetscCall(VecCopy(y,v));
563:     /* X(:,j+1) = (u+v) */
564:     PetscCall(VecAXPY(x,1,y));
565:     /* Y(:,j+1) = conj(u-v) */
566:     PetscCall(VecAYPX(y,-2,x));
567:     PetscCall(VecConjugate(y));
568:     /* Restore */
569:     PetscCall(BVRestoreColumn(X,j+1,&x));
570:     PetscCall(BVRestoreColumn(Y,j+1,&y));
571:   }
572:   if (4*m > 100) PetscCall(PetscFree(hwork));
573:   PetscCall(VecDestroy(&w));
574:   PetscCall(VecDestroy(&f));
575:   PetscCall(VecDestroy(&g));
576:   PetscFunctionReturn(PETSC_SUCCESS);
577: }

579: static PetscErrorCode EPSComputeVectors_BSE_ProjectedBSE(EPS eps)
580: {
581:   Mat         H;
582:   Vec         x1,y1,cx1,cy1;
583:   BV          X,Y;
584:   IS          is[2];
585:   PetscInt    k;
586:   PetscScalar delta1,delta2,lambda;

588:   PetscFunctionBegin;
589:   PetscCall(STGetMatrix(eps->st,0,&H));
590:   PetscCall(MatNestGetISs(H,is,NULL));
591:   PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&X,&Y));
592:   PetscCall(BVCreateVec(X,&cx1));
593:   PetscCall(BVCreateVec(Y,&cy1));
594:   for (k=0; k<eps->nconv; k++) {
595:     /* approx eigenvector is [ delta1*x1 + delta2*conj(y1) ]
596:                              [ delta1*y1 + delta2*conj(x1) ] */
597:     lambda = eps->eigr[k];
598:     PetscCall(STBackTransform(eps->st,1,&lambda,&eps->eigi[k]));
599:     delta1 = lambda+1.0;
600:     delta2 = lambda-1.0;
601:     PetscCall(BVGetColumn(X,k,&x1));
602:     PetscCall(BVGetColumn(Y,k,&y1));
603:     PetscCall(VecCopy(x1,cx1));
604:     PetscCall(VecCopy(y1,cy1));
605:     PetscCall(VecConjugate(cx1));
606:     PetscCall(VecConjugate(cy1));
607:     PetscCall(VecScale(x1,delta1));
608:     PetscCall(VecScale(y1,delta1));
609:     PetscCall(VecAXPY(x1,delta2,cy1));
610:     PetscCall(VecAXPY(y1,delta2,cx1));
611:     PetscCall(BVRestoreColumn(X,k,&x1));
612:     PetscCall(BVRestoreColumn(Y,k,&y1));
613:   }
614:   PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&X,&Y));
615:   PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
616:   /* Normalize eigenvector */
617:   PetscCall(BVNormalize(eps->V,NULL));
618:   PetscCall(VecDestroy(&cx1));
619:   PetscCall(VecDestroy(&cy1));
620:   PetscFunctionReturn(PETSC_SUCCESS);
621: }

623: PetscErrorCode EPSSetUp_KrylovSchur_BSE(EPS eps)
624: {
625:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
626:   PetscBool       flg,sinvert;
627:   PetscInt        nev;

629:   PetscFunctionBegin;
630:   PetscCheck((eps->problem_type==EPS_BSE),PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Problem type should be BSE");
631:   EPSCheckUnsupportedCondition(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION | EPS_FEATURE_BALANCE,PETSC_TRUE," with BSE structure");
632:   PetscCall(PetscCitationsRegister(citation,&cited));
633:   if (eps->nev==0 && eps->stop!=EPS_STOP_THRESHOLD) eps->nev = 1;
634:   nev = (eps->nev+1)/2;
635:   PetscCall(EPSSetDimensions_Default(eps,&nev,&eps->ncv,&eps->mpd));
636:   PetscCheck(eps->ncv<=nev+eps->mpd,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must not be larger than nev+mpd");
637:   if (eps->max_it==PETSC_DETERMINE) eps->max_it = PetscMax(100,2*eps->n/eps->ncv)*((eps->stop==EPS_STOP_THRESHOLD)?10:1);

639:   PetscCall(PetscObjectTypeCompareAny((PetscObject)eps->st,&flg,STSINVERT,STSHIFT,""));
640:   PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Krylov-Schur BSE only supports shift and shift-and-invert ST");
641:   PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STSINVERT,&sinvert));
642:   if (!eps->which) {
643:     if (sinvert) eps->which = EPS_TARGET_MAGNITUDE;
644:     else eps->which = EPS_SMALLEST_MAGNITUDE;
645:   }

647:   if (!ctx->keep) ctx->keep = 0.5;
648:   PetscCall(STSetStructured(eps->st,PETSC_TRUE));

650:   PetscCall(EPSAllocateSolution(eps,1));
651:   switch (ctx->bse) {
652:     case EPS_KRYLOVSCHUR_BSE_SHAO:
653:       eps->ops->solve = EPSSolve_KrylovSchur_BSE_Shao;
654:       eps->ops->computevectors = EPSComputeVectors_BSE_Shao;
655:       PetscCall(DSSetType(eps->ds,DSHEP));
656:       PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
657:       PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
658:       PetscCall(DSAllocate(eps->ds,eps->ncv+1));
659:       break;
660:     case EPS_KRYLOVSCHUR_BSE_GRUNING:
661:       eps->ops->solve = EPSSolve_KrylovSchur_BSE_Gruning;
662:       eps->ops->computevectors = EPSComputeVectors_BSE_Gruning;
663:       PetscCall(DSSetType(eps->ds,DSSVD));
664:       PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
665:       PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
666:       PetscCall(DSAllocate(eps->ds,eps->ncv+1));
667:       break;
668:     case EPS_KRYLOVSCHUR_BSE_PROJECTEDBSE:
669:       eps->ops->solve = EPSSolve_KrylovSchur_BSE_ProjectedBSE;
670:       eps->ops->computevectors = EPSComputeVectors_BSE_ProjectedBSE;
671:       PetscCall(DSSetType(eps->ds,DSHEP));
672:       PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
673:       PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
674:       PetscCall(DSAllocate(eps->ds,eps->ncv+1));
675:       break;
676:     default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Unexpected error");
677:   }
678:   PetscFunctionReturn(PETSC_SUCCESS);
679: }

681: PetscErrorCode EPSSolve_KrylovSchur_BSE_Shao(EPS eps)
682: {
683:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
684:   PetscInt        i,k,l,ld,nv,nconv=0,nevsave;
685:   Mat             H,Q;
686:   BV              U,V;
687:   IS              is[2];
688:   PetscReal       *a,*b,beta;
689:   PetscBool       breakdown=PETSC_FALSE;

691:   PetscFunctionBegin;
692:   PetscCall(DSGetLeadingDimension(eps->ds,&ld));

694:   /* Extract matrix blocks */
695:   PetscCall(STGetMatrix(eps->st,0,&H));
696:   PetscCall(MatNestGetISs(H,is,NULL));

698:   /* Get the split bases */
699:   PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));

701:   nevsave  = eps->nev;
702:   eps->nev = (eps->nev+1)/2;
703:   l = 0;

705:   /* Restart loop */
706:   while (eps->reason == EPS_CONVERGED_ITERATING) {
707:     eps->its++;

709:     /* Compute an nv-step Lanczos factorization */
710:     nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
711:     PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
712:     PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&a));
713:     b = a + ld;
714:     PetscCall(EPSBSELanczos_Shao(eps,U,V,a,b,eps->nconv+l,&nv,&breakdown));
715:     beta = b[nv-1];
716:     PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&a));
717:     PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
718:     PetscCall(DSSetState(eps->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
719:     PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));

721:     /* Solve projected problem */
722:     PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
723:     PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
724:     PetscCall(DSUpdateExtraRow(eps->ds));
725:     PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));

727:     /* Check convergence */
728:     for (i=0;i<nv;i++) eps->eigr[i] = PetscSqrtReal(PetscRealPart(eps->eigr[i]));
729:     PetscCall(EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta,0.0,1.0,&k));
730:     EPSSetCtxThreshold(eps,eps->eigr,eps->eigi,k);
731:     PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
732:     nconv = k;

734:     /* Update l */
735:     if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
736:     else l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
737:     if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
738:     if (l) PetscCall(PetscInfo(eps,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));

740:     if (eps->reason == EPS_CONVERGED_ITERATING) {
741:       PetscCheck(!breakdown,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Breakdown in BSE Krylov-Schur (beta=%g)",(double)beta);
742:       /* Prepare the Rayleigh quotient for restart */
743:       PetscCall(DSTruncate(eps->ds,k+l,PETSC_FALSE));
744:     }
745:     /* Update the corresponding vectors
746:        U(:,idx) = U*Q(:,idx),  V(:,idx) = V*Q(:,idx) */
747:     PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&Q));
748:     PetscCall(BVMultInPlace(U,Q,eps->nconv,k+l));
749:     PetscCall(BVMultInPlace(V,Q,eps->nconv,k+l));
750:     PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&Q));

752:     if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
753:       PetscCall(BVCopyColumn(eps->V,nv,k+l));
754:       if (eps->stop==EPS_STOP_THRESHOLD && nv-k<5) {  /* reallocate */
755:         eps->ncv = eps->mpd+k;
756:         PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
757:         PetscCall(EPSReallocateSolution(eps,eps->ncv+1));
758:         PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
759:         for (i=nv;i<eps->ncv;i++) eps->perm[i] = i;
760:         PetscCall(DSReallocate(eps->ds,eps->ncv+1));
761:         PetscCall(DSGetLeadingDimension(eps->ds,&ld));
762:       }
763:     }
764:     eps->nconv = k;
765:     PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv));
766:   }

768:   eps->nev = nevsave;

770:   PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
771:   PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
772:   PetscFunctionReturn(PETSC_SUCCESS);
773: }

775: /*
776:    EPSConvergence_Gruning - convergence check based on SVDKrylovConvergence().
777: */
778: static PetscErrorCode EPSConvergence_Gruning(EPS eps,PetscBool getall,PetscInt kini,PetscInt nits,PetscInt *kout)
779: {
780:   PetscInt       k,marker,ld;
781:   PetscReal      *alpha,*beta,resnorm;
782:   PetscBool      extra;

784:   PetscFunctionBegin;
785:   *kout = 0;
786:   PetscCall(DSGetLeadingDimension(eps->ds,&ld));
787:   PetscCall(DSGetExtraRow(eps->ds,&extra));
788:   PetscCheck(extra,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Only implemented for DS with extra row");
789:   marker = -1;
790:   if (eps->trackall) getall = PETSC_TRUE;
791:   PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&alpha));
792:   beta = alpha + ld;
793:   for (k=kini;k<kini+nits;k++) {
794:     resnorm = PetscAbsReal(beta[k]);
795:     PetscCall((*eps->converged)(eps,eps->eigr[k],eps->eigi[k],resnorm,&eps->errest[k],eps->convergedctx));
796:     if (marker==-1 && eps->errest[k] >= eps->tol) marker = k;
797:     if (marker!=-1 && !getall) break;
798:   }
799:   PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&alpha));
800:   if (marker!=-1) k = marker;
801:   *kout = k;
802:   PetscFunctionReturn(PETSC_SUCCESS);
803: }

805: PetscErrorCode EPSSolve_KrylovSchur_BSE_Gruning(EPS eps)
806: {
807:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
808:   PetscInt        i,k,l,ld,nv,nconv=0,nevsave;
809:   Mat             H,Q,Z;
810:   BV              U,V,HU,HV;
811:   IS              is[2];
812:   PetscReal       *d1,*d2,beta;
813:   PetscBool       breakdown=PETSC_FALSE;

815:   PetscFunctionBegin;
816:   PetscCall(DSGetLeadingDimension(eps->ds,&ld));

818:   /* Extract matrix blocks */
819:   PetscCall(STGetMatrix(eps->st,0,&H));
820:   PetscCall(MatNestGetISs(H,is,NULL));

822:   /* Get the split bases */
823:   PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));

825:   /* Create HU and HV */
826:   PetscCall(BVDuplicate(U,&HU));
827:   PetscCall(BVDuplicate(V,&HV));

829:   nevsave  = eps->nev;
830:   eps->nev = (eps->nev+1)/2;
831:   l = 0;

833:   /* Restart loop */
834:   while (eps->reason == EPS_CONVERGED_ITERATING) {
835:     eps->its++;

837:     /* Compute an nv-step Lanczos factorization */
838:     nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
839:     PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
840:     PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&d1));
841:     d2 = d1 + ld;
842:     PetscCall(EPSBSELanczos_Gruning(eps,U,V,HU,HV,d1,d2,eps->nconv+l,&nv,&breakdown));
843:     beta = d1[nv-1];
844:     PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&d1));

846:     /* Compute SVD */
847:     PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
848:     PetscCall(DSSVDSetDimensions(eps->ds,nv));
849:     PetscCall(DSSetState(eps->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
850:     PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));

852:     PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
853:     PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
854:     PetscCall(DSUpdateExtraRow(eps->ds));
855:     PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));

857:     /* Check convergence */
858:     PetscCall(EPSConvergence_Gruning(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,&k));
859:     EPSSetCtxThreshold(eps,eps->eigr,eps->eigi,k);
860:     PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
861:     nconv = k;

863:     /* Update l */
864:     if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
865:     else l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
866:     if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
867:     if (l) PetscCall(PetscInfo(eps,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));

869:     if (eps->reason == EPS_CONVERGED_ITERATING) {
870:       PetscCheck(!breakdown,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Breakdown in BSE Krylov-Schur (beta=%g)",(double)beta);
871:       /* Prepare the Rayleigh quotient for restart */
872:       PetscCall(DSTruncate(eps->ds,k+l,PETSC_FALSE));
873:     }
874:     /* Update the corresponding vectors
875:        U(:,idx) = U*Q(:,idx),  V(:,idx) = V*Z(:,idx) */
876:     PetscCall(DSGetMat(eps->ds,DS_MAT_U,&Q));
877:     PetscCall(DSGetMat(eps->ds,DS_MAT_V,&Z));
878:     PetscCall(BVMultInPlace(U,Z,eps->nconv,k+l));
879:     PetscCall(BVMultInPlace(V,Q,eps->nconv,k+l));
880:     PetscCall(BVMultInPlace(HU,Z,eps->nconv,k+l));
881:     PetscCall(BVMultInPlace(HV,Q,eps->nconv,k+l));
882:     PetscCall(DSRestoreMat(eps->ds,DS_MAT_U,&Q));
883:     PetscCall(DSRestoreMat(eps->ds,DS_MAT_V,&Z));

885:     if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
886:       PetscCall(BVCopyColumn(U,nv,k+l));
887:       PetscCall(BVCopyColumn(HU,nv,k+l));
888:       if (eps->stop==EPS_STOP_THRESHOLD && nv-k<5) {  /* reallocate */
889:         eps->ncv = eps->mpd+k;
890:         PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
891:         PetscCall(EPSReallocateSolution(eps,eps->ncv+1));
892:         PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
893:         PetscCall(BVResize(HU,eps->ncv+1,PETSC_TRUE));
894:         PetscCall(BVResize(HV,eps->ncv+1,PETSC_TRUE));
895:         for (i=nv;i<eps->ncv;i++) { eps->perm[i] = i; eps->eigi[i] = 0.0; }
896:         PetscCall(DSReallocate(eps->ds,eps->ncv+1));
897:         PetscCall(DSGetLeadingDimension(eps->ds,&ld));
898:       }
899:     }
900:     eps->nconv = k;
901:     PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv));
902:   }

904:   eps->nev = nevsave;

906:   PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
907:   PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));

909:   PetscCall(BVDestroy(&HU));
910:   PetscCall(BVDestroy(&HV));
911:   PetscFunctionReturn(PETSC_SUCCESS);
912: }

914: PetscErrorCode EPSSolve_KrylovSchur_BSE_ProjectedBSE(EPS eps)
915: {
916:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
917:   PetscInt        i,k,l,ld,nv,nconv=0,nevsave;
918:   Mat             H,Q;
919:   Vec             v;
920:   BV              U,V;
921:   IS              is[2];
922:   PetscReal       *a,*b,beta;
923:   PetscBool       breakdown=PETSC_FALSE;

925:   PetscFunctionBegin;
926:   PetscCall(DSGetLeadingDimension(eps->ds,&ld));

928:   /* Extract matrix blocks */
929:   PetscCall(STGetMatrix(eps->st,0,&H));
930:   PetscCall(MatNestGetISs(H,is,NULL));

932:   /* Get the split bases */
933:   PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));

935:   /* Vector v */
936:   PetscCall(BVCreateVec(V,&v));

938:   nevsave  = eps->nev;
939:   eps->nev = (eps->nev+1)/2;
940:   l = 0;

942:   /* Restart loop */
943:   while (eps->reason == EPS_CONVERGED_ITERATING) {
944:     eps->its++;

946:     /* Compute an nv-step Lanczos factorization */
947:     nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
948:     PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
949:     PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&a));
950:     b = a + ld;
951:     PetscCall(EPSBSELanczos_ProjectedBSE(eps,U,V,v,a,b,eps->nconv+l,&nv,&breakdown));
952:     beta = b[nv-1];
953:     PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&a));
954:     PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
955:     PetscCall(DSSetState(eps->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
956:     PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));

958:     /* Solve projected problem */
959:     PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
960:     PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
961:     PetscCall(DSUpdateExtraRow(eps->ds));
962:     PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));

964:     /* Check convergence */
965:     for (i=0;i<nv;i++) eps->eigr[i] = PetscSqrtReal(PetscRealPart(eps->eigr[i]));
966:     PetscCall(EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta,0.0,1.0,&k));
967:     EPSSetCtxThreshold(eps,eps->eigr,eps->eigi,k);
968:     PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
969:     nconv = k;

971:     /* Update l */
972:     if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
973:     else l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
974:     if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
975:     if (l) PetscCall(PetscInfo(eps,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));

977:     if (eps->reason == EPS_CONVERGED_ITERATING) {
978:       PetscCheck(!breakdown,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Breakdown in BSE Krylov-Schur (beta=%g)",(double)beta);
979:       /* Prepare the Rayleigh quotient for restart */
980:       PetscCall(DSTruncate(eps->ds,k+l,PETSC_FALSE));
981:     }
982:     /* Update the corresponding vectors
983:        U(:,idx) = U*Q(:,idx),  V(:,idx) = V*Q(:,idx) */
984:     PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&Q));
985:     PetscCall(BVMultInPlace(U,Q,eps->nconv,k+l));
986:     PetscCall(BVMultInPlace(V,Q,eps->nconv,k+l));
987:     PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&Q));

989:     if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
990:       PetscCall(BVCopyColumn(eps->V,nv,k+l));
991:       if (eps->stop==EPS_STOP_THRESHOLD && nv-k<5) {  /* reallocate */
992:         eps->ncv = eps->mpd+k;
993:         PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
994:         PetscCall(EPSReallocateSolution(eps,eps->ncv+1));
995:         PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
996:         for (i=nv;i<eps->ncv;i++) eps->perm[i] = i;
997:         PetscCall(DSReallocate(eps->ds,eps->ncv+1));
998:         PetscCall(DSGetLeadingDimension(eps->ds,&ld));
999:       }
1000:     }
1001:     eps->nconv = k;
1002:     PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv));
1003:   }

1005:   eps->nev = nevsave;

1007:   PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
1008:   PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));

1010:   PetscCall(VecDestroy(&v));
1011:   PetscFunctionReturn(PETSC_SUCCESS);
1012: }