Actual source code: ks-lrep.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
10: /*
11: SLEPc eigensolver: "krylovschur"
13: Method: thick-restarted Lanczos for Linear Response eigenvalue problems
15: References:
17: [1] Z. Teng, R.-C. Li, "Convergence analysis of Lanczos-type methods for the
18: linear response eigenvalue problem", J. Comput. Appl. Math. 247, 2013.
20: [2] H.-X. Zhong, H. Xu, "Weighted Golub-Kahan-Lanczos bidiagonalization
21: algorithms", Elec. Trans. Numer. Anal. 47, 2017.
23: */
24: #include <slepc/private/epsimpl.h>
25: #include "krylovschur.h"
27: static PetscErrorCode Orthog_Teng(Vec x,BV U,BV V,PetscInt j,PetscScalar *h,PetscScalar *c)
28: {
29: PetscInt i;
31: PetscFunctionBegin;
32: PetscCall(BVSetActiveColumns(U,0,j));
33: PetscCall(BVSetActiveColumns(V,0,j));
34: /* c = U^* x */
35: PetscCall(BVDotVec(U,x,c));
36: /* x = x-V*c */
37: PetscCall(BVMultVec(V,-1.0,1.0,x,c));
38: /* accumulate orthog coeffs into h */
39: for (i=0;i<2*j;i++) h[i] += c[i];
40: PetscFunctionReturn(PETSC_SUCCESS);
41: }
43: /* Orthogonalize vector x against first j vectors in U and V
44: v is column j-1 of V */
45: static PetscErrorCode OrthogonalizeVector_Teng(Vec x,BV U,BV V,PetscInt j,Vec u,PetscReal *beta,PetscInt k,PetscScalar *h)
46: {
47: PetscReal alpha;
48: PetscInt i,l;
50: PetscFunctionBegin;
51: PetscCall(PetscArrayzero(h,2*j));
53: /* Local orthogonalization */
54: l = j==k+1?0:j-2; /* 1st column to orthogonalize against */
55: PetscCall(VecDotRealPart(x,u,&alpha));
56: for (i=l;i<j-1;i++) h[i] = beta[i];
57: h[j-1] = alpha;
58: /* x = x-V(:,l:j-1)*h(l:j-1) */
59: PetscCall(BVSetActiveColumns(V,l,j));
60: PetscCall(BVMultVec(V,-1.0,1.0,x,h+l));
62: /* Full orthogonalization */
63: PetscCall(Orthog_Teng(x,U,V,j,h,h+2*j));
64: PetscFunctionReturn(PETSC_SUCCESS);
65: }
67: static PetscErrorCode EPSLREPLanczos_Teng(EPS eps,Mat K,Mat M,BV U,BV V,PetscReal *alpha,PetscReal *beta,PetscInt k,PetscInt *min,PetscBool *breakdown)
68: {
69: PetscInt j,m = *min;
70: Vec u,v,uh,vh;
71: PetscReal beta0;
72: PetscScalar *hwork,lhwork[100],gamma;
73: PetscBool alloc=PETSC_FALSE;
75: PetscFunctionBegin;
76: if (4*m > 100) {
77: PetscCall(PetscMalloc1(4*m,&hwork));
78: alloc = PETSC_TRUE;
79: } else hwork = lhwork;
81: /* Normalize initial vector */
82: if (k==0) {
83: if (eps->nini==0) PetscCall(BVSetRandomColumn(V,0));
84: PetscCall(BVGetColumn(U,0,&u));
85: PetscCall(BVGetColumn(V,0,&v));
86: PetscCall(MatMult(M,v,u));
87: PetscCall(VecDot(u,v,&gamma));
88: beta0 = PetscSqrtReal(PetscRealPart(gamma));
89: if (beta0==0.0) {
90: if (breakdown) *breakdown = PETSC_TRUE;
91: *min = 1; m = 0;
92: } else {
93: PetscCall(VecScale(u,1.0/beta0));
94: PetscCall(VecScale(v,1.0/beta0));
95: }
96: PetscCall(BVRestoreColumn(U,0,&u));
97: PetscCall(BVRestoreColumn(V,0,&v));
98: }
100: for (j=k;j<m;j++) {
101: /* j+1 columns (indices 0 to j) have been computed */
102: PetscCall(BVGetColumn(U,j+1,&uh));
103: PetscCall(BVGetColumn(V,j+1,&vh));
104: PetscCall(BVGetColumn(U,j,&u));
105: PetscCall(MatMult(K,u,vh));
106: PetscCall(OrthogonalizeVector_Teng(vh,U,V,j+1,u,beta,k,hwork));
107: alpha[j] = PetscRealPart(hwork[j]);
108: PetscCall(MatMult(M,vh,uh));
109: PetscCall(VecDot(uh,vh,&gamma));
110: beta[j] = PetscSqrtReal(PetscRealPart(gamma));
111: if (beta[j]==0.0) {
112: if (breakdown) *breakdown = PETSC_TRUE;
113: *min = j+1; m = j;
114: } else {
115: PetscCall(VecScale(uh,1.0/beta[j]));
116: PetscCall(VecScale(vh,1.0/beta[j]));
117: }
118: PetscCall(BVRestoreColumn(U,j+1,&uh));
119: PetscCall(BVRestoreColumn(V,j+1,&vh));
120: PetscCall(BVRestoreColumn(U,j,&u));
121: }
122: if (alloc) PetscCall(PetscFree(hwork));
123: PetscFunctionReturn(PETSC_SUCCESS);
124: }
126: /* K-Orthogonalize vector vh against first j vectors in V using beta coeffs */
127: static PetscErrorCode OrthogonalizeVector_Zhong_v(Vec vh,BV V,PetscInt j,PetscReal *beta,PetscInt k,PetscScalar *h)
128: {
129: PetscInt i,l;
131: PetscFunctionBegin;
132: /* Local orthogonalization */
133: l = j==k?0:j-1; /* 1st column to orthogonalize against */
134: for (i=l;i<j;i++) h[i] = beta[i];
135: /* vh = vh-V[:,l:j-1]*h[l:j-1] */
136: PetscCall(BVSetActiveColumns(V,l,j));
137: PetscCall(BVMultVec(V,-1.0,1.0,vh,h+l));
138: PetscFunctionReturn(PETSC_SUCCESS);
139: }
141: /* M-Orthogonalize vector uh against first j vectors in U. Full orthog */
142: static PetscErrorCode Orthogonalize_Zhong_u(Vec uh,BV U,BV MU,PetscInt j,PetscScalar *h)
143: {
144: PetscFunctionBegin;
145: PetscCall(BVSetActiveColumns(U,0,j));
146: PetscCall(BVSetActiveColumns(MU,0,j));
147: /* h=MU'*uh */
148: PetscCall(BVDotVec(MU,uh,h));
149: /* uh=uh-U*h */
150: PetscCall(BVMultVec(U,-1.0,1.0,uh,h));
151: PetscFunctionReturn(PETSC_SUCCESS);
152: }
154: /* M-Orthogonalize vector uh against first j vectors in U. Local+full orthog */
155: static PetscErrorCode OrthogonalizeVector_Zhong_u(Vec uh,BV U,BV MU,PetscInt j,PetscReal *alpha,PetscInt k,PetscScalar *h)
156: {
157: Vec u;
159: PetscFunctionBegin;
160: /* Local orthogonalization: uh = uh-U[:,j-1]*alpha[j-1] */
161: PetscCall(BVGetColumn(U,j-1,&u));
162: PetscCall(VecAXPY(uh,-alpha[j-1],u));
163: PetscCall(BVRestoreColumn(U,j-1,&u));
164: /* Full orthogonalization */
165: PetscCall(Orthogonalize_Zhong_u(uh,U,MU,j,h));
166: PetscFunctionReturn(PETSC_SUCCESS);
167: }
169: static PetscErrorCode EPSLREPLanczos_Zhong(EPS eps,Mat K,Mat M,BV U,BV V,BV MU,PetscReal *alpha,PetscReal *beta,PetscInt k,PetscInt *min,PetscBool *breakdown)
170: {
171: PetscInt j,m = *min;
172: Vec uh,vh,x;
173: PetscReal beta0;
174: PetscScalar *hwork,lhwork[100],gamma;
175: PetscBool alloc=PETSC_FALSE;
177: PetscFunctionBegin;
178: if (m > 100) {
179: PetscCall(PetscMalloc1(m,&hwork));
180: alloc = PETSC_TRUE;
181: } else hwork = lhwork;
183: /* Normalize initial vector */
184: if (k==0) {
185: if (eps->nini==0) PetscCall(BVSetRandomColumn(U,0));
186: PetscCall(BVGetColumn(U,0,&uh));
187: PetscCall(BVGetColumn(MU,0,&vh));
188: PetscCall(MatMult(M,uh,vh));
189: PetscCall(VecDot(uh,vh,&gamma));
190: beta0 = PetscSqrtReal(PetscRealPart(gamma));
191: if (beta0==0.0) {
192: if (breakdown) *breakdown = PETSC_TRUE;
193: *min = 1; m = 0;
194: } else {
195: PetscCall(VecScale(uh,1.0/beta0));
196: PetscCall(VecScale(vh,1.0/beta0));
197: }
198: PetscCall(BVRestoreColumn(U,0,&uh));
199: PetscCall(BVRestoreColumn(MU,0,&vh));
200: }
202: for (j=k;j<m;j++) {
203: /* Compute column j of V, then column j+1 of U and MU */
204: PetscCall(BVGetColumn(U,j+1,&uh));
205: PetscCall(BVGetColumn(V,j,&vh));
206: PetscCall(BVGetColumn(MU,j,&x));
207: PetscCall(VecCopy(x,vh));
208: PetscCall(BVRestoreColumn(MU,j,&x));
209: PetscCall(OrthogonalizeVector_Zhong_v(vh,V,j,beta,k,hwork));
210: PetscCall(MatMult(K,vh,uh));
211: PetscCall(VecDot(uh,vh,&gamma));
212: alpha[j] = PetscSqrtReal(PetscRealPart(gamma));
213: if (alpha[j]==0.0) {
214: if (breakdown) *breakdown = PETSC_TRUE;
215: *min = j+1; m = j;
216: PetscCall(BVRestoreColumn(U,j+1,&uh));
217: } else {
218: PetscCall(VecScale(uh,1.0/alpha[j]));
219: PetscCall(VecScale(vh,1.0/alpha[j]));
220: }
221: PetscCall(BVRestoreColumn(V,j,&vh));
222: if (breakdown && *breakdown) continue;
224: PetscCall(OrthogonalizeVector_Zhong_u(uh,U,MU,j+1,alpha,k,hwork));
225: PetscCall(BVGetColumn(MU,j+1,&vh));
226: PetscCall(MatMult(M,uh,vh));
227: PetscCall(VecDot(uh,vh,&gamma));
228: beta[j] = PetscSqrtReal(PetscRealPart(gamma));
229: if (beta[j]==0.0) {
230: if (breakdown) *breakdown = PETSC_TRUE;
231: *min = j+1; m = j;
232: } else {
233: PetscCall(VecScale(uh,1.0/beta[j]));
234: PetscCall(VecScale(vh,1.0/beta[j]));
235: }
236: PetscCall(BVRestoreColumn(U,j+1,&uh));
237: PetscCall(BVRestoreColumn(MU,j+1,&vh));
238: }
239: if (alloc) PetscCall(PetscFree(hwork));
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: /*
244: EPSConvergence_Zhong - convergence check based on SVDKrylovConvergence().
245: FIXME: Code dulication. This is a copy of EPSConvergence_Gruning
246: */
247: static PetscErrorCode EPSConvergence_Zhong(EPS eps,PetscBool getall,PetscInt kini,PetscInt nits,PetscInt *kout)
248: {
249: PetscInt k,marker,ld;
250: PetscReal *alpha,*beta,resnorm;
251: PetscBool extra;
253: PetscFunctionBegin;
254: *kout = 0;
255: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
256: PetscCall(DSGetExtraRow(eps->ds,&extra));
257: PetscCheck(extra,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Only implemented for DS with extra row");
258: marker = -1;
259: if (eps->trackall) getall = PETSC_TRUE;
260: PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&alpha));
261: beta = alpha + ld;
262: for (k=kini;k<kini+nits;k++) {
263: resnorm = PetscAbsReal(beta[k]);
264: PetscCall((*eps->converged)(eps,eps->eigr[k],eps->eigi[k],resnorm,&eps->errest[k],eps->convergedctx));
265: if (marker==-1 && eps->errest[k] >= eps->tol) marker = k;
266: if (marker!=-1 && !getall) break;
267: }
268: PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&alpha));
269: if (marker!=-1) k = marker;
270: *kout = k;
271: PetscFunctionReturn(PETSC_SUCCESS);
272: }
274: static PetscErrorCode EPSUnreduceVectors(EPS eps,BV U,BV V)
275: {
276: PetscInt k;
277: Vec u,v,w;
279: PetscFunctionBegin;
280: /* The approximate eigenvector is [u+v; u-v], where [u; v] is the reduced eigenvector */
281: PetscCall(BVCreateVec(V,&w));
282: for (k=0;k<eps->nconv;k++) {
283: PetscCall(BVGetColumn(U,k,&u));
284: PetscCall(BVGetColumn(V,k,&v));
285: PetscCall(VecCopy(v,w));
286: PetscCall(VecCopy(u,v));
287: PetscCall(VecAXPY(u,1.0,w));
288: PetscCall(VecAXPY(v,-1.0,w));
289: PetscCall(BVRestoreColumn(U,k,&u));
290: PetscCall(BVRestoreColumn(V,k,&v));
291: }
292: PetscCall(VecDestroy(&w));
293: PetscFunctionReturn(PETSC_SUCCESS);
294: }
296: static PetscErrorCode EPSComputeVectors_LREP_Teng(EPS eps)
297: {
298: Mat H;
299: Vec v;
300: BV U,V;
301: IS is[2];
302: PetscInt k;
303: PetscScalar lambda;
304: PetscBool reduced;
306: PetscFunctionBegin;
307: PetscCall(STGetMatrix(eps->st,0,&H));
308: PetscCall(MatNestGetISs(H,is,NULL));
309: PetscCall(SlepcCheckMatLREPReduced(H,&reduced));
310: PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&V,&U));
311: for (k=0;k<eps->nconv;k++) {
312: PetscCall(BVGetColumn(V,k,&v));
313: /* approx eigenvector is [eigr[k]*v; u] */
314: lambda = eps->eigr[k];
315: PetscCall(STBackTransform(eps->st,1,&lambda,&eps->eigi[k]));
316: PetscCall(VecScale(v,lambda));
317: PetscCall(BVRestoreColumn(V,k,&v));
318: }
319: if (!reduced) PetscCall(EPSUnreduceVectors(eps,V,U));
320: PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&V,&U));
321: /* Normalize eigenvectors */
322: PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
323: PetscCall(BVNormalize(eps->V,NULL));
324: PetscFunctionReturn(PETSC_SUCCESS);
325: }
327: static PetscErrorCode EPSComputeVectors_LREP_Zhong(EPS eps)
328: {
329: Mat H;
330: BV U,V;
331: IS is[2];
332: PetscBool reduced;
334: PetscFunctionBegin;
335: PetscCall(STGetMatrix(eps->st,0,&H));
336: PetscCall(SlepcCheckMatLREPReduced(H,&reduced));
337: /* Approx eigenvector for the reduced form is [u; v] */
338: if (!reduced) {
339: PetscCall(MatNestGetISs(H,is,NULL));
340: PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
341: PetscCall(EPSUnreduceVectors(eps, U, V));
342: PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
343: }
344: /* Normalize eigenvectors */
345: PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
346: PetscCall(BVNormalize(eps->V,NULL));
347: PetscFunctionReturn(PETSC_SUCCESS);
348: }
350: PetscErrorCode EPSSetUp_KrylovSchur_LREP(EPS eps)
351: {
352: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
353: PetscBool flg;
355: PetscFunctionBegin;
356: PetscCheck((eps->problem_type==EPS_LREP),PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Problem type should be LREP");
357: EPSCheckUnsupportedCondition(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION | EPS_FEATURE_BALANCE,PETSC_TRUE," with LREP structure");
358: PetscCall(EPSSetDimensions_Default(eps,&eps->nev,&eps->ncv,&eps->mpd));
359: PetscCheck(eps->ncv<=eps->nev+eps->mpd,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must not be larger than nev+mpd");
360: if (eps->max_it==PETSC_DETERMINE) eps->max_it = PetscMax(100,2*eps->n/eps->ncv)*((eps->stop==EPS_STOP_THRESHOLD)?10:1);
362: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STSHIFT,&flg));
363: PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Krylov-Schur LREP only supports shift ST");
364: if (!eps->which) eps->which = EPS_SMALLEST_MAGNITUDE;
366: if (!ctx->keep) ctx->keep = 0.5;
367: PetscCall(STSetStructured(eps->st,PETSC_FALSE));
369: PetscCall(EPSAllocateSolution(eps,1));
370: switch (ctx->lrep) {
371: case EPS_KRYLOVSCHUR_LREP_TENG:
372: eps->ops->solve = EPSSolve_KrylovSchur_LREP_Teng;
373: eps->ops->computevectors = EPSComputeVectors_LREP_Teng;
374: PetscCall(DSSetType(eps->ds,DSHEP));
375: PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
376: PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
377: PetscCall(DSAllocate(eps->ds,eps->ncv+1));
378: break;
379: case EPS_KRYLOVSCHUR_LREP_ZHONG:
380: eps->ops->solve = EPSSolve_KrylovSchur_LREP_Zhong;
381: eps->ops->computevectors = EPSComputeVectors_LREP_Zhong;
382: PetscCall(DSSetType(eps->ds,DSSVD));
383: PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
384: PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
385: PetscCall(DSAllocate(eps->ds,eps->ncv+1));
386: break;
387: default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Unexpected error");
388: }
389: PetscFunctionReturn(PETSC_SUCCESS);
390: }
392: static PetscErrorCode EPSCreateReducedMats(Mat H,Mat *K,Mat *M)
393: {
394: PetscInt ma,na,Ma,Na;
395: Mat A,B;
396: const PetscScalar scal[] = { 1.0, -1.0 };
398: PetscFunctionBegin;
399: PetscCall(MatNestGetSubMat(H,0,0,&A));
400: PetscCall(MatNestGetSubMat(H,0,1,&B));
401: PetscCall(MatGetSize(A,&Ma,&Na));
402: PetscCall(MatGetLocalSize(A,&ma,&na));
403: /* K = A-B */
404: PetscCall(MatCreate(PetscObjectComm((PetscObject)A),K));
405: PetscCall(MatSetSizes(*K,ma,na,Ma,Na));
406: PetscCall(MatSetType(*K,MATCOMPOSITE));
407: PetscCall(MatCompositeAddMat(*K,A));
408: PetscCall(MatCompositeAddMat(*K,B));
409: PetscCall(MatAssemblyBegin(*K,MAT_FINAL_ASSEMBLY));
410: PetscCall(MatAssemblyEnd(*K,MAT_FINAL_ASSEMBLY));
411: PetscCall(MatCompositeSetScalings(*K,scal));
412: /* M = A+B */
413: PetscCall(MatCreate(PetscObjectComm((PetscObject)A),M));
414: PetscCall(MatSetSizes(*M,ma,na,Ma,Na));
415: PetscCall(MatSetType(*M,MATCOMPOSITE));
416: PetscCall(MatCompositeAddMat(*M,A));
417: PetscCall(MatCompositeAddMat(*M,B));
418: PetscCall(MatAssemblyBegin(*M,MAT_FINAL_ASSEMBLY));
419: PetscCall(MatAssemblyEnd(*M,MAT_FINAL_ASSEMBLY));
420: PetscFunctionReturn(PETSC_SUCCESS);
421: }
423: PetscErrorCode EPSSolve_KrylovSchur_LREP_Teng(EPS eps)
424: {
425: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
426: PetscInt i,k,l,ld,nv,nconv=0,nevsave;
427: Mat H,Q,K,M;
428: BV U,V;
429: IS is[2];
430: PetscReal *a,*b,beta;
431: PetscBool reduced,breakdown=PETSC_FALSE;
433: PetscFunctionBegin;
434: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
436: /* Extract matrix blocks */
437: PetscCall(STGetMatrix(eps->st,0,&H));
438: PetscCall(MatNestGetISs(H,is,NULL));
439: PetscCall(SlepcCheckMatLREPReduced(H,&reduced));
440: if (reduced) {
441: PetscCall(MatNestGetSubMat(H,0,1,&K));
442: PetscCall(MatNestGetSubMat(H,1,0,&M));
443: } else PetscCall(EPSCreateReducedMats(H,&K,&M));
445: /* Get the split bases */
446: PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&V,&U));
448: nevsave = eps->nev;
449: eps->nev = (eps->nev+1)/2;
450: l = 0;
452: /* Restart loop */
453: while (eps->reason == EPS_CONVERGED_ITERATING) {
454: eps->its++;
456: /* Compute an nv-step Lanczos factorization */
457: nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
458: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
459: PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&a));
460: b = a + ld;
461: PetscCall(EPSLREPLanczos_Teng(eps,K,M,U,V,a,b,eps->nconv+l,&nv,&breakdown));
462: beta = b[nv-1];
463: PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&a));
464: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
465: PetscCall(DSSetState(eps->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
466: PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
468: /* Solve projected problem */
469: PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
470: PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
471: PetscCall(DSUpdateExtraRow(eps->ds));
472: PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
474: /* Check convergence */
475: for (i=0;i<nv;i++) eps->eigr[i] = PetscSqrtReal(PetscRealPart(eps->eigr[i]));
476: PetscCall(EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta,0.0,1.0,&k));
477: EPSSetCtxThreshold(eps,eps->eigr,eps->eigi,eps->errest,k,nv);
478: PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
479: nconv = k;
481: /* Update l */
482: if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
483: else l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
484: if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
485: if (l) PetscCall(PetscInfo(eps,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
487: if (eps->reason == EPS_CONVERGED_ITERATING) {
488: PetscCheck(!breakdown,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Breakdown in LREP Krylov-Schur (beta=%g)",(double)beta);
489: /* Prepare the Rayleigh quotient for restart */
490: PetscCall(DSTruncate(eps->ds,k+l,PETSC_FALSE));
491: }
492: /* Update the corresponding vectors
493: U(:,idx) = U*Q(:,idx), V(:,idx) = V*Q(:,idx) */
494: PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&Q));
495: PetscCall(BVMultInPlace(U,Q,eps->nconv,k+l));
496: PetscCall(BVMultInPlace(V,Q,eps->nconv,k+l));
497: PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&Q));
499: if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
500: PetscCall(BVCopyColumn(eps->V,nv,k+l));
501: if (eps->stop==EPS_STOP_THRESHOLD && nv-k<5) { /* reallocate */
502: eps->ncv = eps->mpd+k;
503: PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&V,&U));
504: PetscCall(EPSReallocateSolution(eps,eps->ncv+1));
505: PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&V,&U));
506: for (i=nv;i<eps->ncv;i++) eps->perm[i] = i;
507: PetscCall(DSReallocate(eps->ds,eps->ncv+1));
508: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
509: }
510: }
511: eps->nconv = k;
512: PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv));
513: }
515: eps->nev = nevsave;
517: PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
518: PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&V,&U));
519: if (!reduced) {
520: PetscCall(MatDestroy(&K));
521: PetscCall(MatDestroy(&M));
522: }
523: PetscFunctionReturn(PETSC_SUCCESS);
524: }
526: PetscErrorCode EPSSolve_KrylovSchur_LREP_Zhong(EPS eps)
527: {
528: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
529: PetscInt i,k,l,ld,nv,nconv=0,nevsave;
530: Mat H,Q,Z,K,M;
531: BV U,V,MU;
532: IS is[2];
533: PetscReal *a,*b,beta;
534: PetscBool reduced,breakdown=PETSC_FALSE;
536: PetscFunctionBegin;
537: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
539: /* Extract matrix blocks */
540: PetscCall(STGetMatrix(eps->st,0,&H));
541: PetscCall(MatNestGetISs(H,is,NULL));
542: PetscCall(SlepcCheckMatLREPReduced(H,&reduced));
543: if (reduced) {
544: PetscCall(MatNestGetSubMat(H,0,1,&K));
545: PetscCall(MatNestGetSubMat(H,1,0,&M));
546: } else PetscCall(EPSCreateReducedMats(H,&K,&M));
548: /* Get the split bases */
549: PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
551: /* Create MU */
552: PetscCall(BVDuplicate(U,&MU));
554: nevsave = eps->nev;
555: eps->nev = (eps->nev+1)/2;
556: l = 0;
558: /* Restart loop */
559: while (eps->reason == EPS_CONVERGED_ITERATING) {
560: eps->its++;
562: /* Compute an nv-step Lanczos factorization */
563: nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
564: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
565: PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&a));
566: b = a + ld;
567: PetscCall(EPSLREPLanczos_Zhong(eps,K,M,U,V,MU,a,b,eps->nconv+l,&nv,&breakdown));
568: beta = b[nv-1];
569: PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&a));
570: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
571: PetscCall(DSSVDSetDimensions(eps->ds,nv));
572: PetscCall(DSSetState(eps->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
573: PetscCall(BVSetActiveColumns(U,eps->nconv,nv));
574: PetscCall(BVSetActiveColumns(V,eps->nconv,nv));
576: /* Solve projected problem */
577: PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
578: PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
579: PetscCall(DSUpdateExtraRow(eps->ds));
580: PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
582: /* Check convergence */
583: PetscCall(EPSConvergence_Zhong(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,&k));
584: EPSSetCtxThreshold(eps,eps->eigr,eps->eigi,eps->errest,k,nv);
585: PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
586: nconv = k;
588: /* Update l */
589: if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
590: else l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
591: if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
592: if (l) PetscCall(PetscInfo(eps,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
594: if (eps->reason == EPS_CONVERGED_ITERATING) {
595: PetscCheck(!breakdown,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Breakdown in LREP Krylov-Schur (beta=%g)",(double)beta);
596: /* Prepare the Rayleigh quotient for restart */
597: PetscCall(DSTruncate(eps->ds,k+l,PETSC_FALSE));
598: }
599: /* Update the corresponding vectors
600: U(:,idx) = U*Q(:,idx), MU(:,idx) = MU*Q(:,idx), V(:,idx) = V*Z(:,idx), */
601: PetscCall(DSGetMat(eps->ds,DS_MAT_U,&Z));
602: PetscCall(DSGetMat(eps->ds,DS_MAT_V,&Q));
603: PetscCall(BVMultInPlace(U,Q,eps->nconv,k+l));
604: PetscCall(BVMultInPlace(MU,Q,eps->nconv,k+l));
605: PetscCall(BVMultInPlace(V,Z,eps->nconv,k+l));
606: PetscCall(DSRestoreMat(eps->ds,DS_MAT_U,&Z));
607: PetscCall(DSRestoreMat(eps->ds,DS_MAT_V,&Q));
609: if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
610: PetscCall(BVCopyColumn(U,nv,k+l));
611: PetscCall(BVCopyColumn(MU,nv,k+l));
612: if (eps->stop==EPS_STOP_THRESHOLD && nv-k<5) { /* reallocate */
613: eps->ncv = eps->mpd+k;
614: PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
615: PetscCall(EPSReallocateSolution(eps,eps->ncv+1));
616: PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
617: PetscCall(BVResize(MU,eps->ncv+1,PETSC_TRUE));
618: for (i=nv;i<eps->ncv;i++) eps->perm[i] = i;
619: PetscCall(DSReallocate(eps->ds,eps->ncv+1));
620: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
621: }
622: }
624: eps->nconv = k;
625: PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv));
626: }
628: eps->nev = nevsave;
630: PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
631: PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
632: PetscCall(BVDestroy(&MU));
633: if (!reduced) {
634: PetscCall(MatDestroy(&K));
635: PetscCall(MatDestroy(&M));
636: }
637: PetscFunctionReturn(PETSC_SUCCESS);
638: }