Actual source code: ks-bse.c

slepc-3.22.1 2024-10-28
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: */
 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: }