Actual source code: ptoar.c
slepc-3.17.2 2022-08-09
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 polynomial eigensolver: "toar"
13: Method: TOAR
15: Algorithm:
17: Two-Level Orthogonal Arnoldi.
19: References:
21: [1] Y. Su, J. Zhang and Z. Bai, "A compact Arnoldi algorithm for
22: polynomial eigenvalue problems", talk presented at RANMEP 2008.
24: [2] C. Campos and J.E. Roman, "Parallel Krylov solvers for the
25: polynomial eigenvalue problem in SLEPc", SIAM J. Sci. Comput.
26: 38(5):S385-S411, 2016.
28: [3] D. Lu, Y. Su and Z. Bai, "Stability analysis of the two-level
29: orthogonal Arnoldi procedure", SIAM J. Matrix Anal. App.
30: 37(1):195-214, 2016.
31: */
33: #include <slepc/private/pepimpl.h>
34: #include "../src/pep/impls/krylov/pepkrylov.h"
35: #include <slepcblaslapack.h>
37: static PetscBool cited = PETSC_FALSE;
38: static const char citation[] =
39: "@Article{slepc-pep,\n"
40: " author = \"C. Campos and J. E. Roman\",\n"
41: " title = \"Parallel {Krylov} solvers for the polynomial eigenvalue problem in {SLEPc}\",\n"
42: " journal = \"{SIAM} J. Sci. Comput.\",\n"
43: " volume = \"38\",\n"
44: " number = \"5\",\n"
45: " pages = \"S385--S411\",\n"
46: " year = \"2016,\"\n"
47: " doi = \"https://doi.org/10.1137/15M1022458\"\n"
48: "}\n";
50: PetscErrorCode PEPSetUp_TOAR(PEP pep)
51: {
52: PEP_TOAR *ctx = (PEP_TOAR*)pep->data;
53: PetscBool sinv,flg;
54: PetscInt i;
56: PEPCheckShiftSinvert(pep);
57: PEPSetDimensions_Default(pep,pep->nev,&pep->ncv,&pep->mpd);
59: if (pep->max_it==PETSC_DEFAULT) pep->max_it = PetscMax(100,2*(pep->nmat-1)*pep->n/pep->ncv);
60: if (!pep->which) PEPSetWhichEigenpairs_Default(pep);
62: if (pep->problem_type!=PEP_GENERAL) PetscInfo(pep,"Problem type ignored, performing a non-symmetric linearization\n");
64: if (!ctx->keep) ctx->keep = 0.5;
66: PEPAllocateSolution(pep,pep->nmat-1);
67: PEPSetWorkVecs(pep,3);
68: DSSetType(pep->ds,DSNHEP);
69: DSSetExtraRow(pep->ds,PETSC_TRUE);
70: DSAllocate(pep->ds,pep->ncv+1);
72: PEPBasisCoefficients(pep,pep->pbc);
73: STGetTransform(pep->st,&flg);
74: if (!flg) {
75: PetscFree(pep->solvematcoeffs);
76: PetscMalloc1(pep->nmat,&pep->solvematcoeffs);
77: PetscLogObjectMemory((PetscObject)pep,pep->nmat*sizeof(PetscScalar));
78: PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv);
79: if (sinv) PEPEvaluateBasis(pep,pep->target,0,pep->solvematcoeffs,NULL);
80: else {
81: for (i=0;i<pep->nmat-1;i++) pep->solvematcoeffs[i] = 0.0;
82: pep->solvematcoeffs[pep->nmat-1] = 1.0;
83: }
84: }
85: BVDestroy(&ctx->V);
86: BVCreateTensor(pep->V,pep->nmat-1,&ctx->V);
87: PetscFunctionReturn(0);
88: }
90: /*
91: Extend the TOAR basis by applying the the matrix operator
92: over a vector which is decomposed in the TOAR way
93: Input:
94: - pbc: array containing the polynomial basis coefficients
95: - S,V: define the latest Arnoldi vector (nv vectors in V)
96: Output:
97: - t: new vector extending the TOAR basis
98: - r: temporary coefficients to compute the TOAR coefficients
99: for the new Arnoldi vector
100: Workspace: t_ (two vectors)
101: */
102: static PetscErrorCode PEPTOARExtendBasis(PEP pep,PetscBool sinvert,PetscScalar sigma,PetscScalar *S,PetscInt ls,PetscInt nv,BV V,Vec t,PetscScalar *r,PetscInt lr,Vec *t_)
103: {
104: PetscInt nmat=pep->nmat,deg=nmat-1,k,j,off=0,lss;
105: Vec v=t_[0],ve=t_[1],q=t_[2];
106: PetscScalar alpha=1.0,*ss,a;
107: PetscReal *ca=pep->pbc,*cb=pep->pbc+nmat,*cg=pep->pbc+2*nmat;
108: PetscBool flg;
110: BVSetActiveColumns(pep->V,0,nv);
111: STGetTransform(pep->st,&flg);
112: if (sinvert) {
113: for (j=0;j<nv;j++) {
114: if (deg>1) r[lr+j] = S[j]/ca[0];
115: if (deg>2) r[2*lr+j] = (S[ls+j]+(sigma-cb[1])*r[lr+j])/ca[1];
116: }
117: for (k=2;k<deg-1;k++) {
118: for (j=0;j<nv;j++) r[(k+1)*lr+j] = (S[k*ls+j]+(sigma-cb[k])*r[k*lr+j]-cg[k]*r[(k-1)*lr+j])/ca[k];
119: }
120: k = deg-1;
121: for (j=0;j<nv;j++) r[j] = (S[k*ls+j]+(sigma-cb[k])*r[k*lr+j]-cg[k]*r[(k-1)*lr+j])/ca[k];
122: ss = r; lss = lr; off = 1; alpha = -1.0; a = pep->sfactor;
123: } else {
124: ss = S; lss = ls; off = 0; alpha = -ca[deg-1]; a = 1.0;
125: }
126: BVMultVec(V,1.0,0.0,v,ss+off*lss);
127: if (PetscUnlikely(pep->Dr)) { /* balancing */
128: VecPointwiseMult(v,v,pep->Dr);
129: }
130: STMatMult(pep->st,off,v,q);
131: VecScale(q,a);
132: for (j=1+off;j<deg+off-1;j++) {
133: BVMultVec(V,1.0,0.0,v,ss+j*lss);
134: if (PetscUnlikely(pep->Dr)) VecPointwiseMult(v,v,pep->Dr);
135: STMatMult(pep->st,j,v,t);
136: a *= pep->sfactor;
137: VecAXPY(q,a,t);
138: }
139: if (sinvert) {
140: BVMultVec(V,1.0,0.0,v,ss);
141: if (PetscUnlikely(pep->Dr)) VecPointwiseMult(v,v,pep->Dr);
142: STMatMult(pep->st,deg,v,t);
143: a *= pep->sfactor;
144: VecAXPY(q,a,t);
145: } else {
146: BVMultVec(V,1.0,0.0,ve,ss+(deg-1)*lss);
147: if (PetscUnlikely(pep->Dr)) VecPointwiseMult(ve,ve,pep->Dr);
148: a *= pep->sfactor;
149: STMatMult(pep->st,deg-1,ve,t);
150: VecAXPY(q,a,t);
151: a *= pep->sfactor;
152: }
153: if (flg || !sinvert) alpha /= a;
154: STMatSolve(pep->st,q,t);
155: VecScale(t,alpha);
156: if (!sinvert) {
157: VecAXPY(t,cg[deg-1],v);
158: VecAXPY(t,cb[deg-1],ve);
159: }
160: if (PetscUnlikely(pep->Dr)) VecPointwiseDivide(t,t,pep->Dr);
161: PetscFunctionReturn(0);
162: }
164: /*
165: Compute TOAR coefficients of the blocks of the new Arnoldi vector computed
166: */
167: static PetscErrorCode PEPTOARCoefficients(PEP pep,PetscBool sinvert,PetscScalar sigma,PetscInt nv,PetscScalar *S,PetscInt ls,PetscScalar *r,PetscInt lr,PetscScalar *x)
168: {
169: PetscInt k,j,nmat=pep->nmat,d=nmat-1;
170: PetscReal *ca=pep->pbc,*cb=pep->pbc+nmat,*cg=pep->pbc+2*nmat;
171: PetscScalar t=1.0,tp=0.0,tt;
173: if (sinvert) {
174: for (k=1;k<d;k++) {
175: tt = t;
176: t = ((sigma-cb[k-1])*t-cg[k-1]*tp)/ca[k-1]; /* k-th basis polynomial */
177: tp = tt;
178: for (j=0;j<=nv;j++) r[k*lr+j] += t*x[j];
179: }
180: } else {
181: for (j=0;j<=nv;j++) r[j] = (cb[0]-sigma)*S[j]+ca[0]*S[ls+j];
182: for (k=1;k<d-1;k++) {
183: for (j=0;j<=nv;j++) r[k*lr+j] = (cb[k]-sigma)*S[k*ls+j]+ca[k]*S[(k+1)*ls+j]+cg[k]*S[(k-1)*ls+j];
184: }
185: if (sigma!=0.0) for (j=0;j<=nv;j++) r[(d-1)*lr+j] -= sigma*S[(d-1)*ls+j];
186: }
187: PetscFunctionReturn(0);
188: }
190: /*
191: Compute a run of Arnoldi iterations dim(work)=ld
192: */
193: static PetscErrorCode PEPTOARrun(PEP pep,PetscScalar sigma,PetscScalar *H,PetscInt ldh,PetscInt k,PetscInt *M,PetscBool *breakdown,Vec *t_)
194: {
195: PEP_TOAR *ctx = (PEP_TOAR*)pep->data;
196: PetscInt j,m=*M,deg=pep->nmat-1,ld;
197: PetscInt lds,nqt,l;
198: Vec t;
199: PetscReal norm;
200: PetscBool flg,sinvert=PETSC_FALSE,lindep;
201: PetscScalar *x,*S;
202: Mat MS;
204: BVTensorGetFactors(ctx->V,NULL,&MS);
205: MatDenseGetArray(MS,&S);
206: BVGetSizes(pep->V,NULL,NULL,&ld);
207: lds = ld*deg;
208: BVGetActiveColumns(pep->V,&l,&nqt);
209: STGetTransform(pep->st,&flg);
210: if (!flg) {
211: /* spectral transformation handled by the solver */
212: PetscObjectTypeCompareAny((PetscObject)pep->st,&flg,STSINVERT,STSHIFT,"");
214: PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinvert);
215: }
216: BVSetActiveColumns(ctx->V,0,m);
217: for (j=k;j<m;j++) {
218: /* apply operator */
219: BVGetColumn(pep->V,nqt,&t);
220: PEPTOARExtendBasis(pep,sinvert,sigma,S+j*lds,ld,nqt,pep->V,t,S+(j+1)*lds,ld,t_);
221: BVRestoreColumn(pep->V,nqt,&t);
223: /* orthogonalize */
224: if (sinvert) x = S+(j+1)*lds;
225: else x = S+(deg-1)*ld+(j+1)*lds;
226: BVOrthogonalizeColumn(pep->V,nqt,x,&norm,&lindep);
227: if (!lindep) {
228: x[nqt] = norm;
229: BVScaleColumn(pep->V,nqt,1.0/norm);
230: nqt++;
231: }
233: PEPTOARCoefficients(pep,sinvert,sigma,nqt-1,S+j*lds,ld,S+(j+1)*lds,ld,x);
235: /* level-2 orthogonalization */
236: BVOrthogonalizeColumn(ctx->V,j+1,H+j*ldh,&norm,breakdown);
237: H[j+1+ldh*j] = norm;
238: if (PetscUnlikely(*breakdown)) {
239: *M = j+1;
240: break;
241: }
242: BVScaleColumn(ctx->V,j+1,1.0/norm);
243: BVSetActiveColumns(pep->V,l,nqt);
244: }
245: BVSetActiveColumns(ctx->V,0,*M);
246: MatDenseRestoreArray(MS,&S);
247: BVTensorRestoreFactors(ctx->V,NULL,&MS);
248: PetscFunctionReturn(0);
249: }
251: /*
252: Computes T_j = phi_idx(T). In T_j and T_p are phi_{idx-1}(T)
253: and phi_{idx-2}(T) respectively or null if idx=0,1.
254: Tp and Tj are input/output arguments
255: */
256: static PetscErrorCode PEPEvaluateBasisM(PEP pep,PetscInt k,PetscScalar *T,PetscInt ldt,PetscInt idx,PetscScalar **Tp,PetscScalar **Tj)
257: {
258: PetscInt i;
259: PetscReal *ca,*cb,*cg;
260: PetscScalar *pt,g,a;
261: PetscBLASInt k_,ldt_;
263: if (idx==0) {
264: PetscArrayzero(*Tj,k*k);
265: PetscArrayzero(*Tp,k*k);
266: for (i=0;i<k;i++) (*Tj)[i+i*k] = 1.0;
267: } else {
268: PetscBLASIntCast(ldt,&ldt_);
269: PetscBLASIntCast(k,&k_);
270: ca = pep->pbc; cb = pep->pbc+pep->nmat; cg = pep->pbc+2*pep->nmat;
271: for (i=0;i<k;i++) T[i*ldt+i] -= cb[idx-1];
272: a = 1/ca[idx-1];
273: g = (idx==1)?0.0:-cg[idx-1]/ca[idx-1];
274: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&a,T,&ldt_,*Tj,&k_,&g,*Tp,&k_));
275: pt = *Tj; *Tj = *Tp; *Tp = pt;
276: for (i=0;i<k;i++) T[i*ldt+i] += cb[idx-1];
277: }
278: PetscFunctionReturn(0);
279: }
281: static PetscErrorCode PEPExtractInvariantPair(PEP pep,PetscScalar sigma,PetscInt sr,PetscInt k,PetscScalar *S,PetscInt ld,PetscInt deg,PetscScalar *H,PetscInt ldh)
282: {
283: PetscInt i,j,jj,lds,ldt,d=pep->nmat-1,idxcpy=0;
284: PetscScalar *At,*Bt,*Hj,*Hp,*T,sone=1.0,g,a,*pM,*work;
285: PetscBLASInt k_,sr_,lds_,ldh_,info,*p,lwork,ldt_;
286: PetscBool transf=PETSC_FALSE,flg;
287: PetscReal norm,maxnrm,*rwork;
288: BV *R,Y;
289: Mat M,*A;
291: if (k==0) PetscFunctionReturn(0);
292: lds = deg*ld;
293: PetscCalloc6(k,&p,sr*k,&At,k*k,&Bt,k*k,&Hj,k*k,&Hp,sr*k,&work);
294: PetscBLASIntCast(sr,&sr_);
295: PetscBLASIntCast(k,&k_);
296: PetscBLASIntCast(lds,&lds_);
297: PetscBLASIntCast(ldh,&ldh_);
298: STGetTransform(pep->st,&flg);
299: if (!flg) {
300: PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&flg);
301: if (flg || sigma!=0.0) transf=PETSC_TRUE;
302: }
303: if (transf) {
304: PetscMalloc1(k*k,&T);
305: ldt = k;
306: for (i=0;i<k;i++) PetscArraycpy(T+k*i,H+i*ldh,k);
307: if (flg) {
308: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
309: PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&k_,&k_,T,&k_,p,&info));
310: SlepcCheckLapackInfo("getrf",info);
311: PetscBLASIntCast(sr*k,&lwork);
312: PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&k_,T,&k_,p,work,&lwork,&info));
313: SlepcCheckLapackInfo("getri",info);
314: PetscFPTrapPop();
315: }
316: if (sigma!=0.0) for (i=0;i<k;i++) T[i+k*i] += sigma;
317: } else {
318: T = H; ldt = ldh;
319: }
320: PetscBLASIntCast(ldt,&ldt_);
321: switch (pep->extract) {
322: case PEP_EXTRACT_NONE:
323: break;
324: case PEP_EXTRACT_NORM:
325: if (pep->basis == PEP_BASIS_MONOMIAL) {
326: PetscBLASIntCast(ldt,&ldt_);
327: PetscMalloc1(k,&rwork);
328: norm = LAPACKlange_("F",&k_,&k_,T,&ldt_,rwork);
329: PetscFree(rwork);
330: if (norm>1.0) idxcpy = d-1;
331: } else {
332: PetscBLASIntCast(ldt,&ldt_);
333: PetscMalloc1(k,&rwork);
334: maxnrm = 0.0;
335: for (i=0;i<pep->nmat-1;i++) {
336: PEPEvaluateBasisM(pep,k,T,ldt,i,&Hp,&Hj);
337: norm = LAPACKlange_("F",&k_,&k_,Hj,&k_,rwork);
338: if (norm > maxnrm) {
339: idxcpy = i;
340: maxnrm = norm;
341: }
342: }
343: PetscFree(rwork);
344: }
345: if (idxcpy>0) {
346: /* copy block idxcpy of S to the first one */
347: for (j=0;j<k;j++) PetscArraycpy(S+j*lds,S+idxcpy*ld+j*lds,sr);
348: }
349: break;
350: case PEP_EXTRACT_RESIDUAL:
351: STGetTransform(pep->st,&flg);
352: if (flg) {
353: PetscMalloc1(pep->nmat,&A);
354: for (i=0;i<pep->nmat;i++) STGetMatrixTransformed(pep->st,i,A+i);
355: } else A = pep->A;
356: PetscMalloc1(pep->nmat-1,&R);
357: for (i=0;i<pep->nmat-1;i++) BVDuplicateResize(pep->V,k,R+i);
358: BVDuplicateResize(pep->V,sr,&Y);
359: MatCreateSeqDense(PETSC_COMM_SELF,sr,k,NULL,&M);
360: g = 0.0; a = 1.0;
361: BVSetActiveColumns(pep->V,0,sr);
362: for (j=0;j<pep->nmat;j++) {
363: BVMatMult(pep->V,A[j],Y);
364: PEPEvaluateBasisM(pep,k,T,ldt,i,&Hp,&Hj);
365: for (i=0;i<pep->nmat-1;i++) {
366: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&sr_,&k_,&k_,&a,S+i*ld,&lds_,Hj,&k_,&g,At,&sr_));
367: MatDenseGetArray(M,&pM);
368: for (jj=0;jj<k;jj++) PetscArraycpy(pM+jj*sr,At+jj*sr,sr);
369: MatDenseRestoreArray(M,&pM);
370: BVMult(R[i],1.0,(i==0)?0.0:1.0,Y,M);
371: }
372: }
374: /* frobenius norm */
375: maxnrm = 0.0;
376: for (i=0;i<pep->nmat-1;i++) {
377: BVNorm(R[i],NORM_FROBENIUS,&norm);
378: if (maxnrm > norm) {
379: maxnrm = norm;
380: idxcpy = i;
381: }
382: }
383: if (idxcpy>0) {
384: /* copy block idxcpy of S to the first one */
385: for (j=0;j<k;j++) PetscArraycpy(S+j*lds,S+idxcpy*ld+j*lds,sr);
386: }
387: if (flg) PetscFree(A);
388: for (i=0;i<pep->nmat-1;i++) BVDestroy(&R[i]);
389: PetscFree(R);
390: BVDestroy(&Y);
391: MatDestroy(&M);
392: break;
393: case PEP_EXTRACT_STRUCTURED:
394: for (j=0;j<k;j++) Bt[j+j*k] = 1.0;
395: for (j=0;j<sr;j++) {
396: for (i=0;i<k;i++) At[j*k+i] = PetscConj(S[i*lds+j]);
397: }
398: PEPEvaluateBasisM(pep,k,T,ldt,0,&Hp,&Hj);
399: for (i=1;i<deg;i++) {
400: PEPEvaluateBasisM(pep,k,T,ldt,i,&Hp,&Hj);
401: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","C",&k_,&sr_,&k_,&sone,Hj,&k_,S+i*ld,&lds_,&sone,At,&k_));
402: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","C",&k_,&k_,&k_,&sone,Hj,&k_,Hj,&k_,&sone,Bt,&k_));
403: }
404: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
405: PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&k_,&sr_,Bt,&k_,p,At,&k_,&info));
406: PetscFPTrapPop();
407: SlepcCheckLapackInfo("gesv",info);
408: for (j=0;j<sr;j++) {
409: for (i=0;i<k;i++) S[i*lds+j] = PetscConj(At[j*k+i]);
410: }
411: break;
412: }
413: if (transf) PetscFree(T);
414: PetscFree6(p,At,Bt,Hj,Hp,work);
415: PetscFunctionReturn(0);
416: }
418: PetscErrorCode PEPSolve_TOAR(PEP pep)
419: {
420: PEP_TOAR *ctx = (PEP_TOAR*)pep->data;
421: PetscInt i,j,k,l,nv=0,ld,lds,ldds,nq=0,nconv=0;
422: PetscInt nmat=pep->nmat,deg=nmat-1;
423: PetscScalar *S,*H,sigma;
424: PetscReal beta;
425: PetscBool breakdown=PETSC_FALSE,flg,falselock=PETSC_FALSE,sinv=PETSC_FALSE;
426: Mat MS,MQ;
428: PetscCitationsRegister(citation,&cited);
429: if (ctx->lock) {
430: /* undocumented option to use a cheaper locking instead of the true locking */
431: PetscOptionsGetBool(NULL,NULL,"-pep_toar_falselocking",&falselock,NULL);
432: }
433: DSGetLeadingDimension(pep->ds,&ldds);
434: STGetShift(pep->st,&sigma);
436: /* update polynomial basis coefficients */
437: STGetTransform(pep->st,&flg);
438: if (pep->sfactor!=1.0) {
439: for (i=0;i<nmat;i++) {
440: pep->pbc[nmat+i] /= pep->sfactor;
441: pep->pbc[2*nmat+i] /= pep->sfactor*pep->sfactor;
442: }
443: if (!flg) {
444: pep->target /= pep->sfactor;
445: RGPushScale(pep->rg,1.0/pep->sfactor);
446: STScaleShift(pep->st,1.0/pep->sfactor);
447: sigma /= pep->sfactor;
448: } else {
449: PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv);
450: pep->target = sinv?pep->target*pep->sfactor:pep->target/pep->sfactor;
451: RGPushScale(pep->rg,sinv?pep->sfactor:1.0/pep->sfactor);
452: STScaleShift(pep->st,sinv?pep->sfactor:1.0/pep->sfactor);
453: }
454: }
456: if (flg) sigma = 0.0;
458: /* clean projected matrix (including the extra-arrow) */
459: DSGetArray(pep->ds,DS_MAT_A,&H);
460: PetscArrayzero(H,ldds*ldds);
461: DSRestoreArray(pep->ds,DS_MAT_A,&H);
463: /* Get the starting Arnoldi vector */
464: BVTensorBuildFirstColumn(ctx->V,pep->nini);
466: /* restart loop */
467: l = 0;
468: while (pep->reason == PEP_CONVERGED_ITERATING) {
469: pep->its++;
471: /* compute an nv-step Lanczos factorization */
472: nv = PetscMax(PetscMin(nconv+pep->mpd,pep->ncv),nv);
473: DSGetArray(pep->ds,DS_MAT_A,&H);
474: PEPTOARrun(pep,sigma,H,ldds,pep->nconv+l,&nv,&breakdown,pep->work);
475: beta = PetscAbsScalar(H[(nv-1)*ldds+nv]);
476: DSRestoreArray(pep->ds,DS_MAT_A,&H);
477: DSSetDimensions(pep->ds,nv,pep->nconv,pep->nconv+l);
478: DSSetState(pep->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE);
479: BVSetActiveColumns(ctx->V,pep->nconv,nv);
481: /* solve projected problem */
482: DSSolve(pep->ds,pep->eigr,pep->eigi);
483: DSSort(pep->ds,pep->eigr,pep->eigi,NULL,NULL,NULL);
484: DSUpdateExtraRow(pep->ds);
485: DSSynchronize(pep->ds,pep->eigr,pep->eigi);
487: /* check convergence */
488: PEPKrylovConvergence(pep,PETSC_FALSE,pep->nconv,nv-pep->nconv,beta,&k);
489: (*pep->stopping)(pep,pep->its,pep->max_it,k,pep->nev,&pep->reason,pep->stoppingctx);
491: /* update l */
492: if (pep->reason != PEP_CONVERGED_ITERATING || breakdown) l = 0;
493: else {
494: l = (nv==k)?0:PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
495: DSGetTruncateSize(pep->ds,k,nv,&l);
496: if (!breakdown) {
497: /* prepare the Rayleigh quotient for restart */
498: DSTruncate(pep->ds,k+l,PETSC_FALSE);
499: }
500: }
501: nconv = k;
502: if (!ctx->lock && pep->reason == PEP_CONVERGED_ITERATING && !breakdown) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
503: if (l) PetscInfo(pep,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l);
505: /* update S */
506: DSGetMat(pep->ds,DS_MAT_Q,&MQ);
507: BVMultInPlace(ctx->V,MQ,pep->nconv,k+l);
508: MatDestroy(&MQ);
510: /* copy last column of S */
511: BVCopyColumn(ctx->V,nv,k+l);
513: if (PetscUnlikely(breakdown && pep->reason == PEP_CONVERGED_ITERATING)) {
514: /* stop if breakdown */
515: PetscInfo(pep,"Breakdown TOAR method (it=%" PetscInt_FMT " norm=%g)\n",pep->its,(double)beta);
516: pep->reason = PEP_DIVERGED_BREAKDOWN;
517: }
518: if (pep->reason != PEP_CONVERGED_ITERATING) l--;
519: /* truncate S */
520: BVGetActiveColumns(pep->V,NULL,&nq);
521: if (k+l+deg<=nq) {
522: BVSetActiveColumns(ctx->V,pep->nconv,k+l+1);
523: if (!falselock && ctx->lock) BVTensorCompress(ctx->V,k-pep->nconv);
524: else BVTensorCompress(ctx->V,0);
525: }
526: pep->nconv = k;
527: PEPMonitor(pep,pep->its,nconv,pep->eigr,pep->eigi,pep->errest,nv);
528: }
529: if (pep->nconv>0) {
530: /* {V*S_nconv^i}_{i=0}^{d-1} has rank nconv instead of nconv+d-1. Force zeros in each S_nconv^i block */
531: BVSetActiveColumns(ctx->V,0,pep->nconv);
532: BVGetActiveColumns(pep->V,NULL,&nq);
533: BVSetActiveColumns(pep->V,0,nq);
534: if (nq>pep->nconv) {
535: BVTensorCompress(ctx->V,pep->nconv);
536: BVSetActiveColumns(pep->V,0,pep->nconv);
537: nq = pep->nconv;
538: }
540: /* perform Newton refinement if required */
541: if (pep->refine==PEP_REFINE_MULTIPLE && pep->rits>0) {
542: /* extract invariant pair */
543: BVTensorGetFactors(ctx->V,NULL,&MS);
544: MatDenseGetArray(MS,&S);
545: DSGetArray(pep->ds,DS_MAT_A,&H);
546: BVGetSizes(pep->V,NULL,NULL,&ld);
547: lds = deg*ld;
548: PEPExtractInvariantPair(pep,sigma,nq,pep->nconv,S,ld,deg,H,ldds);
549: DSRestoreArray(pep->ds,DS_MAT_A,&H);
550: DSSetDimensions(pep->ds,pep->nconv,0,0);
551: DSSetState(pep->ds,DS_STATE_RAW);
552: PEPNewtonRefinement_TOAR(pep,sigma,&pep->rits,NULL,pep->nconv,S,lds);
553: DSSolve(pep->ds,pep->eigr,pep->eigi);
554: DSSort(pep->ds,pep->eigr,pep->eigi,NULL,NULL,NULL);
555: DSSynchronize(pep->ds,pep->eigr,pep->eigi);
556: DSGetMat(pep->ds,DS_MAT_Q,&MQ);
557: BVMultInPlace(ctx->V,MQ,0,pep->nconv);
558: MatDestroy(&MQ);
559: MatDenseRestoreArray(MS,&S);
560: BVTensorRestoreFactors(ctx->V,NULL,&MS);
561: }
562: }
563: STGetTransform(pep->st,&flg);
564: if (pep->refine!=PEP_REFINE_MULTIPLE || pep->rits==0) {
565: if (!flg && pep->ops->backtransform) (*pep->ops->backtransform)(pep);
566: if (pep->sfactor!=1.0) {
567: for (j=0;j<pep->nconv;j++) {
568: pep->eigr[j] *= pep->sfactor;
569: pep->eigi[j] *= pep->sfactor;
570: }
571: /* restore original values */
572: for (i=0;i<pep->nmat;i++) {
573: pep->pbc[pep->nmat+i] *= pep->sfactor;
574: pep->pbc[2*pep->nmat+i] *= pep->sfactor*pep->sfactor;
575: }
576: }
577: }
578: /* restore original values */
579: if (!flg) {
580: pep->target *= pep->sfactor;
581: STScaleShift(pep->st,pep->sfactor);
582: } else {
583: STScaleShift(pep->st,sinv?1.0/pep->sfactor:pep->sfactor);
584: pep->target = (sinv)?pep->target/pep->sfactor:pep->target*pep->sfactor;
585: }
586: if (pep->sfactor!=1.0) RGPopScale(pep->rg);
588: /* change the state to raw so that DSVectors() computes eigenvectors from scratch */
589: DSSetDimensions(pep->ds,pep->nconv,0,0);
590: DSSetState(pep->ds,DS_STATE_RAW);
591: PetscFunctionReturn(0);
592: }
594: static PetscErrorCode PEPTOARSetRestart_TOAR(PEP pep,PetscReal keep)
595: {
596: PEP_TOAR *ctx = (PEP_TOAR*)pep->data;
598: if (keep==PETSC_DEFAULT) ctx->keep = 0.5;
599: else {
601: ctx->keep = keep;
602: }
603: PetscFunctionReturn(0);
604: }
606: /*@
607: PEPTOARSetRestart - Sets the restart parameter for the TOAR
608: method, in particular the proportion of basis vectors that must be kept
609: after restart.
611: Logically Collective on pep
613: Input Parameters:
614: + pep - the eigenproblem solver context
615: - keep - the number of vectors to be kept at restart
617: Options Database Key:
618: . -pep_toar_restart - Sets the restart parameter
620: Notes:
621: Allowed values are in the range [0.1,0.9]. The default is 0.5.
623: Level: advanced
625: .seealso: PEPTOARGetRestart()
626: @*/
627: PetscErrorCode PEPTOARSetRestart(PEP pep,PetscReal keep)
628: {
631: PetscTryMethod(pep,"PEPTOARSetRestart_C",(PEP,PetscReal),(pep,keep));
632: PetscFunctionReturn(0);
633: }
635: static PetscErrorCode PEPTOARGetRestart_TOAR(PEP pep,PetscReal *keep)
636: {
637: PEP_TOAR *ctx = (PEP_TOAR*)pep->data;
639: *keep = ctx->keep;
640: PetscFunctionReturn(0);
641: }
643: /*@
644: PEPTOARGetRestart - Gets the restart parameter used in the TOAR method.
646: Not Collective
648: Input Parameter:
649: . pep - the eigenproblem solver context
651: Output Parameter:
652: . keep - the restart parameter
654: Level: advanced
656: .seealso: PEPTOARSetRestart()
657: @*/
658: PetscErrorCode PEPTOARGetRestart(PEP pep,PetscReal *keep)
659: {
662: PetscUseMethod(pep,"PEPTOARGetRestart_C",(PEP,PetscReal*),(pep,keep));
663: PetscFunctionReturn(0);
664: }
666: static PetscErrorCode PEPTOARSetLocking_TOAR(PEP pep,PetscBool lock)
667: {
668: PEP_TOAR *ctx = (PEP_TOAR*)pep->data;
670: ctx->lock = lock;
671: PetscFunctionReturn(0);
672: }
674: /*@
675: PEPTOARSetLocking - Choose between locking and non-locking variants of
676: the TOAR method.
678: Logically Collective on pep
680: Input Parameters:
681: + pep - the eigenproblem solver context
682: - lock - true if the locking variant must be selected
684: Options Database Key:
685: . -pep_toar_locking - Sets the locking flag
687: Notes:
688: The default is to lock converged eigenpairs when the method restarts.
689: This behaviour can be changed so that all directions are kept in the
690: working subspace even if already converged to working accuracy (the
691: non-locking variant).
693: Level: advanced
695: .seealso: PEPTOARGetLocking()
696: @*/
697: PetscErrorCode PEPTOARSetLocking(PEP pep,PetscBool lock)
698: {
701: PetscTryMethod(pep,"PEPTOARSetLocking_C",(PEP,PetscBool),(pep,lock));
702: PetscFunctionReturn(0);
703: }
705: static PetscErrorCode PEPTOARGetLocking_TOAR(PEP pep,PetscBool *lock)
706: {
707: PEP_TOAR *ctx = (PEP_TOAR*)pep->data;
709: *lock = ctx->lock;
710: PetscFunctionReturn(0);
711: }
713: /*@
714: PEPTOARGetLocking - Gets the locking flag used in the TOAR method.
716: Not Collective
718: Input Parameter:
719: . pep - the eigenproblem solver context
721: Output Parameter:
722: . lock - the locking flag
724: Level: advanced
726: .seealso: PEPTOARSetLocking()
727: @*/
728: PetscErrorCode PEPTOARGetLocking(PEP pep,PetscBool *lock)
729: {
732: PetscUseMethod(pep,"PEPTOARGetLocking_C",(PEP,PetscBool*),(pep,lock));
733: PetscFunctionReturn(0);
734: }
736: PetscErrorCode PEPSetFromOptions_TOAR(PetscOptionItems *PetscOptionsObject,PEP pep)
737: {
738: PetscBool flg,lock;
739: PetscReal keep;
741: PetscOptionsHead(PetscOptionsObject,"PEP TOAR Options");
743: PetscOptionsReal("-pep_toar_restart","Proportion of vectors kept after restart","PEPTOARSetRestart",0.5,&keep,&flg);
744: if (flg) PEPTOARSetRestart(pep,keep);
746: PetscOptionsBool("-pep_toar_locking","Choose between locking and non-locking variants","PEPTOARSetLocking",PETSC_FALSE,&lock,&flg);
747: if (flg) PEPTOARSetLocking(pep,lock);
749: PetscOptionsTail();
750: PetscFunctionReturn(0);
751: }
753: PetscErrorCode PEPView_TOAR(PEP pep,PetscViewer viewer)
754: {
755: PEP_TOAR *ctx = (PEP_TOAR*)pep->data;
756: PetscBool isascii;
758: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
759: if (isascii) {
760: PetscViewerASCIIPrintf(viewer," %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep));
761: PetscViewerASCIIPrintf(viewer," using the %slocking variant\n",ctx->lock?"":"non-");
762: }
763: PetscFunctionReturn(0);
764: }
766: PetscErrorCode PEPDestroy_TOAR(PEP pep)
767: {
768: PEP_TOAR *ctx = (PEP_TOAR*)pep->data;
770: BVDestroy(&ctx->V);
771: PetscFree(pep->data);
772: PetscObjectComposeFunction((PetscObject)pep,"PEPTOARSetRestart_C",NULL);
773: PetscObjectComposeFunction((PetscObject)pep,"PEPTOARGetRestart_C",NULL);
774: PetscObjectComposeFunction((PetscObject)pep,"PEPTOARSetLocking_C",NULL);
775: PetscObjectComposeFunction((PetscObject)pep,"PEPTOARGetLocking_C",NULL);
776: PetscFunctionReturn(0);
777: }
779: SLEPC_EXTERN PetscErrorCode PEPCreate_TOAR(PEP pep)
780: {
781: PEP_TOAR *ctx;
783: PetscNewLog(pep,&ctx);
784: pep->data = (void*)ctx;
786: pep->lineariz = PETSC_TRUE;
787: ctx->lock = PETSC_TRUE;
789: pep->ops->solve = PEPSolve_TOAR;
790: pep->ops->setup = PEPSetUp_TOAR;
791: pep->ops->setfromoptions = PEPSetFromOptions_TOAR;
792: pep->ops->destroy = PEPDestroy_TOAR;
793: pep->ops->view = PEPView_TOAR;
794: pep->ops->backtransform = PEPBackTransform_Default;
795: pep->ops->computevectors = PEPComputeVectors_Default;
796: pep->ops->extractvectors = PEPExtractVectors_TOAR;
798: PetscObjectComposeFunction((PetscObject)pep,"PEPTOARSetRestart_C",PEPTOARSetRestart_TOAR);
799: PetscObjectComposeFunction((PetscObject)pep,"PEPTOARGetRestart_C",PEPTOARGetRestart_TOAR);
800: PetscObjectComposeFunction((PetscObject)pep,"PEPTOARSetLocking_C",PEPTOARSetLocking_TOAR);
801: PetscObjectComposeFunction((PetscObject)pep,"PEPTOARGetLocking_C",PEPTOARGetLocking_TOAR);
802: PetscFunctionReturn(0);
803: }