Actual source code: ptoar.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 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: static PetscErrorCode PEPSetUp_TOAR(PEP pep)
51: {
52: PEP_TOAR *ctx = (PEP_TOAR*)pep->data;
53: PetscBool sinv,flg;
54: PetscInt i;
56: PetscFunctionBegin;
57: PEPCheckShiftSinvert(pep);
58: PetscCall(PEPSetDimensions_Default(pep,pep->nev,&pep->ncv,&pep->mpd));
59: PetscCheck(ctx->lock || pep->mpd>=pep->ncv,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Should not use mpd parameter in non-locking variant");
60: if (pep->max_it==PETSC_DETERMINE) pep->max_it = PetscMax(100,2*(pep->nmat-1)*pep->n/pep->ncv);
61: if (!pep->which) PetscCall(PEPSetWhichEigenpairs_Default(pep));
62: PetscCheck(pep->which!=PEP_ALL,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"This solver does not support computing all eigenvalues");
63: if (pep->problem_type!=PEP_GENERAL) PetscCall(PetscInfo(pep,"Problem type ignored, performing a non-symmetric linearization\n"));
65: if (!ctx->keep) ctx->keep = 0.5;
67: PetscCall(PEPAllocateSolution(pep,pep->nmat-1));
68: PetscCall(PEPSetWorkVecs(pep,3));
69: PetscCall(DSSetType(pep->ds,DSNHEP));
70: PetscCall(DSSetExtraRow(pep->ds,PETSC_TRUE));
71: PetscCall(DSAllocate(pep->ds,pep->ncv+1));
73: PetscCall(PEPBasisCoefficients(pep,pep->pbc));
74: PetscCall(STGetTransform(pep->st,&flg));
75: if (!flg) {
76: PetscCall(PetscFree(pep->solvematcoeffs));
77: PetscCall(PetscMalloc1(pep->nmat,&pep->solvematcoeffs));
78: PetscCall(PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv));
79: if (sinv) PetscCall(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: PetscCall(BVDestroy(&ctx->V));
86: PetscCall(BVCreateTensor(pep->V,pep->nmat-1,&ctx->V));
87: PetscFunctionReturn(PETSC_SUCCESS);
88: }
90: /*
91: Extend the TOAR basis by applying 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: PetscFunctionBegin;
111: PetscCall(BVSetActiveColumns(pep->V,0,nv));
112: PetscCall(STGetTransform(pep->st,&flg));
113: if (sinvert) {
114: for (j=0;j<nv;j++) {
115: if (deg>1) r[lr+j] = S[j]/ca[0];
116: if (deg>2) r[2*lr+j] = (S[ls+j]+(sigma-cb[1])*r[lr+j])/ca[1];
117: }
118: for (k=2;k<deg-1;k++) {
119: 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];
120: }
121: k = deg-1;
122: 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];
123: ss = r; lss = lr; off = 1; alpha = -1.0; a = pep->sfactor;
124: } else {
125: ss = S; lss = ls; off = 0; alpha = -ca[deg-1]; a = 1.0;
126: }
127: PetscCall(BVMultVec(V,1.0,0.0,v,ss+off*lss));
128: if (PetscUnlikely(pep->Dr)) { /* balancing */
129: PetscCall(VecPointwiseMult(v,v,pep->Dr));
130: }
131: PetscCall(STMatMult(pep->st,off,v,q));
132: PetscCall(VecScale(q,a));
133: for (j=1+off;j<deg+off-1;j++) {
134: PetscCall(BVMultVec(V,1.0,0.0,v,ss+j*lss));
135: if (PetscUnlikely(pep->Dr)) PetscCall(VecPointwiseMult(v,v,pep->Dr));
136: PetscCall(STMatMult(pep->st,j,v,t));
137: a *= pep->sfactor;
138: PetscCall(VecAXPY(q,a,t));
139: }
140: if (sinvert) {
141: PetscCall(BVMultVec(V,1.0,0.0,v,ss));
142: if (PetscUnlikely(pep->Dr)) PetscCall(VecPointwiseMult(v,v,pep->Dr));
143: PetscCall(STMatMult(pep->st,deg,v,t));
144: a *= pep->sfactor;
145: PetscCall(VecAXPY(q,a,t));
146: } else {
147: PetscCall(BVMultVec(V,1.0,0.0,ve,ss+(deg-1)*lss));
148: if (PetscUnlikely(pep->Dr)) PetscCall(VecPointwiseMult(ve,ve,pep->Dr));
149: a *= pep->sfactor;
150: PetscCall(STMatMult(pep->st,deg-1,ve,t));
151: PetscCall(VecAXPY(q,a,t));
152: a *= pep->sfactor;
153: }
154: if (flg || !sinvert) alpha /= a;
155: PetscCall(STMatSolve(pep->st,q,t));
156: PetscCall(VecScale(t,alpha));
157: if (!sinvert) {
158: PetscCall(VecAXPY(t,cg[deg-1],v));
159: PetscCall(VecAXPY(t,cb[deg-1],ve));
160: }
161: if (PetscUnlikely(pep->Dr)) PetscCall(VecPointwiseDivide(t,t,pep->Dr));
162: PetscFunctionReturn(PETSC_SUCCESS);
163: }
165: /*
166: Compute TOAR coefficients of the blocks of the new Arnoldi vector computed
167: */
168: static PetscErrorCode PEPTOARCoefficients(PEP pep,PetscBool sinvert,PetscScalar sigma,PetscInt nv,PetscScalar *S,PetscInt ls,PetscScalar *r,PetscInt lr,PetscScalar *x)
169: {
170: PetscInt k,j,nmat=pep->nmat,d=nmat-1;
171: PetscReal *ca=pep->pbc,*cb=pep->pbc+nmat,*cg=pep->pbc+2*nmat;
172: PetscScalar t=1.0,tp=0.0,tt;
174: PetscFunctionBegin;
175: if (sinvert) {
176: for (k=1;k<d;k++) {
177: tt = t;
178: t = ((sigma-cb[k-1])*t-cg[k-1]*tp)/ca[k-1]; /* k-th basis polynomial */
179: tp = tt;
180: for (j=0;j<=nv;j++) r[k*lr+j] += t*x[j];
181: }
182: } else {
183: for (j=0;j<=nv;j++) r[j] = (cb[0]-sigma)*S[j]+ca[0]*S[ls+j];
184: for (k=1;k<d-1;k++) {
185: 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];
186: }
187: if (sigma!=0.0) for (j=0;j<=nv;j++) r[(d-1)*lr+j] -= sigma*S[(d-1)*ls+j];
188: }
189: PetscFunctionReturn(PETSC_SUCCESS);
190: }
192: /*
193: Compute a run of Arnoldi iterations dim(work)=ld
194: */
195: static PetscErrorCode PEPTOARrun(PEP pep,PetscScalar sigma,Mat A,PetscInt k,PetscInt *M,PetscReal *beta,PetscBool *breakdown,Vec *t_)
196: {
197: PEP_TOAR *ctx = (PEP_TOAR*)pep->data;
198: PetscInt j,m=*M,deg=pep->nmat-1,ld;
199: PetscInt ldh,lds,nqt,l;
200: Vec t;
201: PetscReal norm=0.0;
202: PetscBool flg,sinvert=PETSC_FALSE,lindep;
203: PetscScalar *H,*x,*S;
204: Mat MS;
206: PetscFunctionBegin;
207: *beta = 0.0;
208: PetscCall(MatDenseGetArray(A,&H));
209: PetscCall(MatDenseGetLDA(A,&ldh));
210: PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
211: PetscCall(MatDenseGetArray(MS,&S));
212: PetscCall(BVGetSizes(pep->V,NULL,NULL,&ld));
213: lds = ld*deg;
214: PetscCall(BVGetActiveColumns(pep->V,&l,&nqt));
215: PetscCall(STGetTransform(pep->st,&flg));
216: if (!flg) {
217: /* spectral transformation handled by the solver */
218: PetscCall(PetscObjectTypeCompareAny((PetscObject)pep->st,&flg,STSINVERT,STSHIFT,""));
219: PetscCheck(flg,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"ST type not supported for TOAR without transforming matrices");
220: PetscCall(PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinvert));
221: }
222: PetscCall(BVSetActiveColumns(ctx->V,0,m));
223: for (j=k;j<m;j++) {
224: /* apply operator */
225: PetscCall(BVGetColumn(pep->V,nqt,&t));
226: PetscCall(PEPTOARExtendBasis(pep,sinvert,sigma,S+j*lds,ld,nqt,pep->V,t,S+(j+1)*lds,ld,t_));
227: PetscCall(BVRestoreColumn(pep->V,nqt,&t));
229: /* orthogonalize */
230: if (sinvert) x = S+(j+1)*lds;
231: else x = S+(deg-1)*ld+(j+1)*lds;
232: PetscCall(BVOrthogonalizeColumn(pep->V,nqt,x,&norm,&lindep));
233: if (!lindep) {
234: x[nqt] = norm;
235: PetscCall(BVScaleColumn(pep->V,nqt,1.0/norm));
236: nqt++;
237: }
239: PetscCall(PEPTOARCoefficients(pep,sinvert,sigma,nqt-1,S+j*lds,ld,S+(j+1)*lds,ld,x));
241: /* level-2 orthogonalization */
242: PetscCall(BVOrthogonalizeColumn(ctx->V,j+1,H+j*ldh,&norm,breakdown));
243: H[j+1+ldh*j] = norm;
244: if (PetscUnlikely(*breakdown)) {
245: *M = j+1;
246: break;
247: }
248: PetscCall(BVScaleColumn(ctx->V,j+1,1.0/norm));
249: PetscCall(BVSetActiveColumns(pep->V,l,nqt));
250: }
251: *beta = norm;
252: PetscCall(BVSetActiveColumns(ctx->V,0,*M));
253: PetscCall(MatDenseRestoreArray(MS,&S));
254: PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
255: PetscCall(MatDenseRestoreArray(A,&H));
256: PetscFunctionReturn(PETSC_SUCCESS);
257: }
259: /*
260: Computes T_j = phi_idx(T). In T_j and T_p are phi_{idx-1}(T)
261: and phi_{idx-2}(T) respectively or null if idx=0,1.
262: Tp and Tj are input/output arguments
263: */
264: static PetscErrorCode PEPEvaluateBasisM(PEP pep,PetscInt k,PetscScalar *T,PetscInt ldt,PetscInt idx,PetscScalar **Tp,PetscScalar **Tj)
265: {
266: PetscInt i;
267: PetscReal *ca,*cb,*cg;
268: PetscScalar *pt,g,a;
269: PetscBLASInt k_,ldt_;
271: PetscFunctionBegin;
272: if (idx==0) {
273: PetscCall(PetscArrayzero(*Tj,k*k));
274: PetscCall(PetscArrayzero(*Tp,k*k));
275: for (i=0;i<k;i++) (*Tj)[i+i*k] = 1.0;
276: } else {
277: PetscCall(PetscBLASIntCast(ldt,&ldt_));
278: PetscCall(PetscBLASIntCast(k,&k_));
279: ca = pep->pbc; cb = pep->pbc+pep->nmat; cg = pep->pbc+2*pep->nmat;
280: for (i=0;i<k;i++) T[i*ldt+i] -= cb[idx-1];
281: a = 1/ca[idx-1];
282: g = (idx==1)?0.0:-cg[idx-1]/ca[idx-1];
283: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&a,T,&ldt_,*Tj,&k_,&g,*Tp,&k_));
284: pt = *Tj; *Tj = *Tp; *Tp = pt;
285: for (i=0;i<k;i++) T[i*ldt+i] += cb[idx-1];
286: }
287: PetscFunctionReturn(PETSC_SUCCESS);
288: }
290: static PetscErrorCode PEPExtractInvariantPair(PEP pep,PetscScalar sigma,PetscInt sr,PetscInt k,PetscScalar *S,PetscInt ld,PetscInt deg,Mat HH)
291: {
292: PetscInt i,j,jj,ldh,lds,ldt,d=pep->nmat-1,idxcpy=0;
293: PetscScalar *H,*At,*Bt,*Hj,*Hp,*T,sone=1.0,g,a,*pM,*work;
294: PetscBLASInt k_,sr_,lds_,ldh_,*p,lwork,ldt_;
295: PetscBool transf=PETSC_FALSE,flg;
296: PetscReal norm,maxnrm,*rwork;
297: BV *R,Y;
298: Mat M,*A;
300: PetscFunctionBegin;
301: if (k==0) PetscFunctionReturn(PETSC_SUCCESS);
302: PetscCall(MatDenseGetArray(HH,&H));
303: PetscCall(MatDenseGetLDA(HH,&ldh));
304: lds = deg*ld;
305: PetscCall(PetscCalloc6(k,&p,sr*k,&At,k*k,&Bt,k*k,&Hj,k*k,&Hp,sr*k,&work));
306: PetscCall(PetscBLASIntCast(sr,&sr_));
307: PetscCall(PetscBLASIntCast(k,&k_));
308: PetscCall(PetscBLASIntCast(lds,&lds_));
309: PetscCall(PetscBLASIntCast(ldh,&ldh_));
310: PetscCall(STGetTransform(pep->st,&flg));
311: if (!flg) {
312: PetscCall(PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&flg));
313: if (flg || sigma!=0.0) transf=PETSC_TRUE;
314: }
315: if (transf) {
316: PetscCall(PetscMalloc1(k*k,&T));
317: ldt = k;
318: for (i=0;i<k;i++) PetscCall(PetscArraycpy(T+k*i,H+i*ldh,k));
319: if (flg) {
320: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
321: PetscCallLAPACKInfo("LAPACKgetrf",LAPACKgetrf_(&k_,&k_,T,&k_,p,&info));
322: PetscCall(PetscBLASIntCast(sr*k,&lwork));
323: PetscCallLAPACKInfo("LAPACKgetri",LAPACKgetri_(&k_,T,&k_,p,work,&lwork,&info));
324: PetscCall(PetscFPTrapPop());
325: }
326: if (sigma!=0.0) for (i=0;i<k;i++) T[i+k*i] += sigma;
327: } else {
328: T = H; ldt = ldh;
329: }
330: PetscCall(PetscBLASIntCast(ldt,&ldt_));
331: switch (pep->extract) {
332: case PEP_EXTRACT_NONE:
333: break;
334: case PEP_EXTRACT_NORM:
335: if (pep->basis == PEP_BASIS_MONOMIAL) {
336: PetscCall(PetscBLASIntCast(ldt,&ldt_));
337: PetscCall(PetscMalloc1(k,&rwork));
338: norm = LAPACKlange_("F",&k_,&k_,T,&ldt_,rwork);
339: PetscCall(PetscFree(rwork));
340: if (norm>1.0) idxcpy = d-1;
341: } else {
342: PetscCall(PetscBLASIntCast(ldt,&ldt_));
343: PetscCall(PetscMalloc1(k,&rwork));
344: maxnrm = 0.0;
345: for (i=0;i<pep->nmat-1;i++) {
346: PetscCall(PEPEvaluateBasisM(pep,k,T,ldt,i,&Hp,&Hj));
347: norm = LAPACKlange_("F",&k_,&k_,Hj,&k_,rwork);
348: if (norm > maxnrm) {
349: idxcpy = i;
350: maxnrm = norm;
351: }
352: }
353: PetscCall(PetscFree(rwork));
354: }
355: if (idxcpy>0) {
356: /* copy block idxcpy of S to the first one */
357: for (j=0;j<k;j++) PetscCall(PetscArraycpy(S+j*lds,S+idxcpy*ld+j*lds,sr));
358: }
359: break;
360: case PEP_EXTRACT_RESIDUAL:
361: PetscCall(STGetTransform(pep->st,&flg));
362: if (flg) {
363: PetscCall(PetscMalloc1(pep->nmat,&A));
364: for (i=0;i<pep->nmat;i++) PetscCall(STGetMatrixTransformed(pep->st,i,A+i));
365: } else A = pep->A;
366: PetscCall(PetscMalloc1(pep->nmat-1,&R));
367: for (i=0;i<pep->nmat-1;i++) PetscCall(BVDuplicateResize(pep->V,k,R+i));
368: PetscCall(BVDuplicateResize(pep->V,sr,&Y));
369: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,sr,k,NULL,&M));
370: g = 0.0; a = 1.0;
371: PetscCall(BVSetActiveColumns(pep->V,0,sr));
372: for (j=0;j<pep->nmat;j++) {
373: PetscCall(BVMatMult(pep->V,A[j],Y));
374: PetscCall(PEPEvaluateBasisM(pep,k,T,ldt,i,&Hp,&Hj));
375: for (i=0;i<pep->nmat-1;i++) {
376: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&sr_,&k_,&k_,&a,S+i*ld,&lds_,Hj,&k_,&g,At,&sr_));
377: PetscCall(MatDenseGetArray(M,&pM));
378: for (jj=0;jj<k;jj++) PetscCall(PetscArraycpy(pM+jj*sr,At+jj*sr,sr));
379: PetscCall(MatDenseRestoreArray(M,&pM));
380: PetscCall(BVMult(R[i],1.0,(i==0)?0.0:1.0,Y,M));
381: }
382: }
384: /* frobenius norm */
385: maxnrm = 0.0;
386: for (i=0;i<pep->nmat-1;i++) {
387: PetscCall(BVNorm(R[i],NORM_FROBENIUS,&norm));
388: if (maxnrm > norm) {
389: maxnrm = norm;
390: idxcpy = i;
391: }
392: }
393: if (idxcpy>0) {
394: /* copy block idxcpy of S to the first one */
395: for (j=0;j<k;j++) PetscCall(PetscArraycpy(S+j*lds,S+idxcpy*ld+j*lds,sr));
396: }
397: if (flg) PetscCall(PetscFree(A));
398: for (i=0;i<pep->nmat-1;i++) PetscCall(BVDestroy(&R[i]));
399: PetscCall(PetscFree(R));
400: PetscCall(BVDestroy(&Y));
401: PetscCall(MatDestroy(&M));
402: break;
403: case PEP_EXTRACT_STRUCTURED:
404: for (j=0;j<k;j++) Bt[j+j*k] = 1.0;
405: for (j=0;j<sr;j++) {
406: for (i=0;i<k;i++) At[j*k+i] = PetscConj(S[i*lds+j]);
407: }
408: PetscCall(PEPEvaluateBasisM(pep,k,T,ldt,0,&Hp,&Hj));
409: for (i=1;i<deg;i++) {
410: PetscCall(PEPEvaluateBasisM(pep,k,T,ldt,i,&Hp,&Hj));
411: PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&k_,&sr_,&k_,&sone,Hj,&k_,S+i*ld,&lds_,&sone,At,&k_));
412: PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&k_,&k_,&k_,&sone,Hj,&k_,Hj,&k_,&sone,Bt,&k_));
413: }
414: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
415: PetscCallLAPACKInfo("LAPACKgesv",LAPACKgesv_(&k_,&sr_,Bt,&k_,p,At,&k_,&info));
416: PetscCall(PetscFPTrapPop());
417: for (j=0;j<sr;j++) {
418: for (i=0;i<k;i++) S[i*lds+j] = PetscConj(At[j*k+i]);
419: }
420: break;
421: }
422: if (transf) PetscCall(PetscFree(T));
423: PetscCall(PetscFree6(p,At,Bt,Hj,Hp,work));
424: PetscCall(MatDenseRestoreArray(HH,&H));
425: PetscFunctionReturn(PETSC_SUCCESS);
426: }
428: static PetscErrorCode PEPSolve_TOAR(PEP pep)
429: {
430: PEP_TOAR *ctx = (PEP_TOAR*)pep->data;
431: PetscInt i,j,k,l,nv=0,ld,lds,nq=0,nconv=0;
432: PetscInt nmat=pep->nmat,deg=nmat-1;
433: PetscScalar *S,sigma;
434: PetscReal beta;
435: PetscBool breakdown=PETSC_FALSE,flg,falselock=PETSC_FALSE,sinv=PETSC_FALSE;
436: Mat H,MS,MQ;
438: PetscFunctionBegin;
439: PetscCall(PetscCitationsRegister(citation,&cited));
440: if (ctx->lock) {
441: /* undocumented option to use a cheaper locking instead of the true locking */
442: PetscCall(PetscOptionsGetBool(NULL,NULL,"-pep_toar_falselocking",&falselock,NULL));
443: }
444: PetscCall(STGetShift(pep->st,&sigma));
446: /* update polynomial basis coefficients */
447: PetscCall(STGetTransform(pep->st,&flg));
448: if (pep->sfactor!=1.0) {
449: for (i=0;i<nmat;i++) {
450: pep->pbc[nmat+i] /= pep->sfactor;
451: pep->pbc[2*nmat+i] /= pep->sfactor*pep->sfactor;
452: }
453: if (!flg) {
454: pep->target /= pep->sfactor;
455: PetscCall(RGPushScale(pep->rg,1.0/pep->sfactor));
456: PetscCall(STScaleShift(pep->st,1.0/pep->sfactor));
457: sigma /= pep->sfactor;
458: } else {
459: PetscCall(PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv));
460: pep->target = sinv?pep->target*pep->sfactor:pep->target/pep->sfactor;
461: PetscCall(RGPushScale(pep->rg,sinv?pep->sfactor:1.0/pep->sfactor));
462: PetscCall(STScaleShift(pep->st,sinv?pep->sfactor:1.0/pep->sfactor));
463: }
464: }
466: if (flg) sigma = 0.0;
468: /* clean projected matrix (including the extra-arrow) */
469: PetscCall(DSSetDimensions(pep->ds,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_DETERMINE));
470: PetscCall(DSGetMat(pep->ds,DS_MAT_A,&H));
471: PetscCall(MatZeroEntries(H));
472: PetscCall(DSRestoreMat(pep->ds,DS_MAT_A,&H));
474: /* Get the starting Arnoldi vector */
475: PetscCall(BVTensorBuildFirstColumn(ctx->V,pep->nini));
477: /* restart loop */
478: l = 0;
479: while (pep->reason == PEP_CONVERGED_ITERATING) {
480: pep->its++;
482: /* compute an nv-step Lanczos factorization */
483: nv = PetscMax(PetscMin(nconv+pep->mpd,pep->ncv),nv);
484: PetscCall(DSGetMat(pep->ds,DS_MAT_A,&H));
485: PetscCall(PEPTOARrun(pep,sigma,H,pep->nconv+l,&nv,&beta,&breakdown,pep->work));
486: PetscCall(DSRestoreMat(pep->ds,DS_MAT_A,&H));
487: PetscCall(DSSetDimensions(pep->ds,nv,pep->nconv,pep->nconv+l));
488: PetscCall(DSSetState(pep->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
489: PetscCall(BVSetActiveColumns(ctx->V,pep->nconv,nv));
491: /* solve projected problem */
492: PetscCall(DSSolve(pep->ds,pep->eigr,pep->eigi));
493: PetscCall(DSSort(pep->ds,pep->eigr,pep->eigi,NULL,NULL,NULL));
494: PetscCall(DSUpdateExtraRow(pep->ds));
495: PetscCall(DSSynchronize(pep->ds,pep->eigr,pep->eigi));
497: /* check convergence */
498: PetscCall(PEPKrylovConvergence(pep,PETSC_FALSE,pep->nconv,nv-pep->nconv,beta,&k));
499: PetscCall((*pep->stopping)(pep,pep->its,pep->max_it,k,pep->nev,&pep->reason,pep->stoppingctx));
501: /* update l */
502: if (pep->reason != PEP_CONVERGED_ITERATING || breakdown) l = 0;
503: else {
504: l = (nv==k)?0:PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
505: PetscCall(DSGetTruncateSize(pep->ds,k,nv,&l));
506: if (!breakdown) {
507: /* prepare the Rayleigh quotient for restart */
508: PetscCall(DSTruncate(pep->ds,k+l,PETSC_FALSE));
509: }
510: }
511: nconv = k;
512: if (!ctx->lock && pep->reason == PEP_CONVERGED_ITERATING && !breakdown) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
513: if (l) PetscCall(PetscInfo(pep,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
515: /* update S */
516: PetscCall(DSGetMat(pep->ds,DS_MAT_Q,&MQ));
517: PetscCall(BVMultInPlace(ctx->V,MQ,pep->nconv,k+l));
518: PetscCall(DSRestoreMat(pep->ds,DS_MAT_Q,&MQ));
520: /* copy last column of S */
521: PetscCall(BVCopyColumn(ctx->V,nv,k+l));
523: if (PetscUnlikely(breakdown && pep->reason == PEP_CONVERGED_ITERATING)) {
524: /* stop if breakdown */
525: PetscCall(PetscInfo(pep,"Breakdown TOAR method (it=%" PetscInt_FMT " norm=%g)\n",pep->its,(double)beta));
526: pep->reason = PEP_DIVERGED_BREAKDOWN;
527: }
528: if (pep->reason != PEP_CONVERGED_ITERATING) l--;
529: /* truncate S */
530: PetscCall(BVGetActiveColumns(pep->V,NULL,&nq));
531: if (k+l+deg<=nq) {
532: PetscCall(BVSetActiveColumns(ctx->V,pep->nconv,k+l+1));
533: if (!falselock && ctx->lock) PetscCall(BVTensorCompress(ctx->V,k-pep->nconv));
534: else PetscCall(BVTensorCompress(ctx->V,0));
535: }
536: pep->nconv = k;
537: PetscCall(PEPMonitor(pep,pep->its,nconv,pep->eigr,pep->eigi,pep->errest,nv));
538: }
539: if (pep->nconv>0) {
540: /* {V*S_nconv^i}_{i=0}^{d-1} has rank nconv instead of nconv+d-1. Force zeros in each S_nconv^i block */
541: PetscCall(BVSetActiveColumns(ctx->V,0,pep->nconv));
542: PetscCall(BVGetActiveColumns(pep->V,NULL,&nq));
543: PetscCall(BVSetActiveColumns(pep->V,0,nq));
544: if (nq>pep->nconv) {
545: PetscCall(BVTensorCompress(ctx->V,pep->nconv));
546: PetscCall(BVSetActiveColumns(pep->V,0,pep->nconv));
547: nq = pep->nconv;
548: }
550: /* perform Newton refinement if required */
551: if (pep->refine==PEP_REFINE_MULTIPLE && pep->rits>0) {
552: /* extract invariant pair */
553: PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
554: PetscCall(MatDenseGetArray(MS,&S));
555: PetscCall(DSGetMat(pep->ds,DS_MAT_A,&H));
556: PetscCall(BVGetSizes(pep->V,NULL,NULL,&ld));
557: lds = deg*ld;
558: PetscCall(PEPExtractInvariantPair(pep,sigma,nq,pep->nconv,S,ld,deg,H));
559: PetscCall(DSRestoreMat(pep->ds,DS_MAT_A,&H));
560: PetscCall(DSSetDimensions(pep->ds,pep->nconv,0,0));
561: PetscCall(DSSetState(pep->ds,DS_STATE_RAW));
562: PetscCall(PEPNewtonRefinement_TOAR(pep,sigma,&pep->rits,NULL,pep->nconv,S,lds));
563: PetscCall(DSSolve(pep->ds,pep->eigr,pep->eigi));
564: PetscCall(DSSort(pep->ds,pep->eigr,pep->eigi,NULL,NULL,NULL));
565: PetscCall(DSSynchronize(pep->ds,pep->eigr,pep->eigi));
566: PetscCall(DSGetMat(pep->ds,DS_MAT_Q,&MQ));
567: PetscCall(BVMultInPlace(ctx->V,MQ,0,pep->nconv));
568: PetscCall(DSRestoreMat(pep->ds,DS_MAT_Q,&MQ));
569: PetscCall(MatDenseRestoreArray(MS,&S));
570: PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
571: }
572: }
573: PetscCall(STGetTransform(pep->st,&flg));
574: if (pep->refine!=PEP_REFINE_MULTIPLE || pep->rits==0) {
575: if (!flg) PetscTryTypeMethod(pep,backtransform);
576: if (pep->sfactor!=1.0) {
577: for (j=0;j<pep->nconv;j++) {
578: pep->eigr[j] *= pep->sfactor;
579: pep->eigi[j] *= pep->sfactor;
580: }
581: /* restore original values */
582: for (i=0;i<pep->nmat;i++) {
583: pep->pbc[pep->nmat+i] *= pep->sfactor;
584: pep->pbc[2*pep->nmat+i] *= pep->sfactor*pep->sfactor;
585: }
586: }
587: }
588: /* restore original values */
589: if (!flg) {
590: pep->target *= pep->sfactor;
591: PetscCall(STScaleShift(pep->st,pep->sfactor));
592: } else {
593: PetscCall(STScaleShift(pep->st,sinv?1.0/pep->sfactor:pep->sfactor));
594: pep->target = sinv?pep->target/pep->sfactor:pep->target*pep->sfactor;
595: }
596: if (pep->sfactor!=1.0) PetscCall(RGPopScale(pep->rg));
598: /* change the state to raw so that DSVectors() computes eigenvectors from scratch */
599: PetscCall(DSSetDimensions(pep->ds,pep->nconv,0,0));
600: PetscCall(DSSetState(pep->ds,DS_STATE_RAW));
601: PetscFunctionReturn(PETSC_SUCCESS);
602: }
604: static PetscErrorCode PEPTOARSetRestart_TOAR(PEP pep,PetscReal keep)
605: {
606: PEP_TOAR *ctx = (PEP_TOAR*)pep->data;
608: PetscFunctionBegin;
609: if (keep==(PetscReal)PETSC_DEFAULT || keep==(PetscReal)PETSC_DECIDE) ctx->keep = 0.5;
610: else {
611: PetscCheck(keep>=0.1 && keep<=0.9,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"The keep argument must be in the range [0.1,0.9]");
612: ctx->keep = keep;
613: }
614: PetscFunctionReturn(PETSC_SUCCESS);
615: }
617: /*@
618: PEPTOARSetRestart - Sets the restart parameter for the TOAR
619: method, in particular the proportion of basis vectors that must be kept
620: after restart.
622: Logically Collective
624: Input Parameters:
625: + pep - the polynomial eigensolver context
626: - keep - the number of vectors to be kept at restart
628: Options Database Key:
629: . -pep_toar_restart keep - sets the restart parameter
631: Note:
632: Allowed values are in the range [0.1,0.9]. The default is 0.5.
634: Level: advanced
636: .seealso: [](ch:pep), `PEPTOAR`, `PEPTOARGetRestart()`
637: @*/
638: PetscErrorCode PEPTOARSetRestart(PEP pep,PetscReal keep)
639: {
640: PetscFunctionBegin;
643: PetscTryMethod(pep,"PEPTOARSetRestart_C",(PEP,PetscReal),(pep,keep));
644: PetscFunctionReturn(PETSC_SUCCESS);
645: }
647: static PetscErrorCode PEPTOARGetRestart_TOAR(PEP pep,PetscReal *keep)
648: {
649: PEP_TOAR *ctx = (PEP_TOAR*)pep->data;
651: PetscFunctionBegin;
652: *keep = ctx->keep;
653: PetscFunctionReturn(PETSC_SUCCESS);
654: }
656: /*@
657: PEPTOARGetRestart - Gets the restart parameter used in the TOAR method.
659: Not Collective
661: Input Parameter:
662: . pep - the polynomial eigensolver context
664: Output Parameter:
665: . keep - the restart parameter
667: Level: advanced
669: .seealso: [](ch:pep), `PEPTOAR`, `PEPTOARSetRestart()`
670: @*/
671: PetscErrorCode PEPTOARGetRestart(PEP pep,PetscReal *keep)
672: {
673: PetscFunctionBegin;
675: PetscAssertPointer(keep,2);
676: PetscUseMethod(pep,"PEPTOARGetRestart_C",(PEP,PetscReal*),(pep,keep));
677: PetscFunctionReturn(PETSC_SUCCESS);
678: }
680: static PetscErrorCode PEPTOARSetLocking_TOAR(PEP pep,PetscBool lock)
681: {
682: PEP_TOAR *ctx = (PEP_TOAR*)pep->data;
684: PetscFunctionBegin;
685: ctx->lock = lock;
686: PetscFunctionReturn(PETSC_SUCCESS);
687: }
689: /*@
690: PEPTOARSetLocking - Choose between locking and non-locking variants of
691: the TOAR method.
693: Logically Collective
695: Input Parameters:
696: + pep - the polynomial eigensolver context
697: - lock - `PETSC_TRUE` if the locking variant must be selected
699: Options Database Key:
700: . -pep_toar_locking (true|false) - sets the locking flag
702: Note:
703: The default is to lock converged eigenpairs when the method restarts.
704: This behavior can be changed so that all directions are kept in the
705: working subspace even if already converged to working accuracy (the
706: non-locking variant).
708: Level: advanced
710: .seealso: [](ch:pep), `PEPTOAR`, `PEPTOARGetLocking()`
711: @*/
712: PetscErrorCode PEPTOARSetLocking(PEP pep,PetscBool lock)
713: {
714: PetscFunctionBegin;
717: PetscTryMethod(pep,"PEPTOARSetLocking_C",(PEP,PetscBool),(pep,lock));
718: PetscFunctionReturn(PETSC_SUCCESS);
719: }
721: static PetscErrorCode PEPTOARGetLocking_TOAR(PEP pep,PetscBool *lock)
722: {
723: PEP_TOAR *ctx = (PEP_TOAR*)pep->data;
725: PetscFunctionBegin;
726: *lock = ctx->lock;
727: PetscFunctionReturn(PETSC_SUCCESS);
728: }
730: /*@
731: PEPTOARGetLocking - Gets the locking flag used in the TOAR method.
733: Not Collective
735: Input Parameter:
736: . pep - the polynomial eigensolver context
738: Output Parameter:
739: . lock - the locking flag
741: Level: advanced
743: .seealso: [](ch:pep), `PEPTOAR`, `PEPTOARSetLocking()`
744: @*/
745: PetscErrorCode PEPTOARGetLocking(PEP pep,PetscBool *lock)
746: {
747: PetscFunctionBegin;
749: PetscAssertPointer(lock,2);
750: PetscUseMethod(pep,"PEPTOARGetLocking_C",(PEP,PetscBool*),(pep,lock));
751: PetscFunctionReturn(PETSC_SUCCESS);
752: }
754: static PetscErrorCode PEPSetFromOptions_TOAR(PEP pep,PetscOptionItems PetscOptionsObject)
755: {
756: PetscBool flg,lock;
757: PetscReal keep;
759: PetscFunctionBegin;
760: PetscOptionsHeadBegin(PetscOptionsObject,"PEP TOAR Options");
762: PetscCall(PetscOptionsReal("-pep_toar_restart","Proportion of vectors kept after restart","PEPTOARSetRestart",0.5,&keep,&flg));
763: if (flg) PetscCall(PEPTOARSetRestart(pep,keep));
765: PetscCall(PetscOptionsBool("-pep_toar_locking","Choose between locking and non-locking variants","PEPTOARSetLocking",PETSC_FALSE,&lock,&flg));
766: if (flg) PetscCall(PEPTOARSetLocking(pep,lock));
768: PetscOptionsHeadEnd();
769: PetscFunctionReturn(PETSC_SUCCESS);
770: }
772: static PetscErrorCode PEPView_TOAR(PEP pep,PetscViewer viewer)
773: {
774: PEP_TOAR *ctx = (PEP_TOAR*)pep->data;
775: PetscBool isascii;
777: PetscFunctionBegin;
778: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
779: if (isascii) {
780: PetscCall(PetscViewerASCIIPrintf(viewer," %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep)));
781: PetscCall(PetscViewerASCIIPrintf(viewer," using the %slocking variant\n",ctx->lock?"":"non-"));
782: }
783: PetscFunctionReturn(PETSC_SUCCESS);
784: }
786: static PetscErrorCode PEPDestroy_TOAR(PEP pep)
787: {
788: PEP_TOAR *ctx = (PEP_TOAR*)pep->data;
790: PetscFunctionBegin;
791: PetscCall(BVDestroy(&ctx->V));
792: PetscCall(PetscFree(pep->data));
793: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPTOARSetRestart_C",NULL));
794: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPTOARGetRestart_C",NULL));
795: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPTOARSetLocking_C",NULL));
796: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPTOARGetLocking_C",NULL));
797: PetscFunctionReturn(PETSC_SUCCESS);
798: }
800: /*MC
801: PEPTOAR - PEPTOAR = "toar" - The Two-level Orthogonal Arnoldi method (TOAR)
802: for polynomial eigenvalue problems.
804: Notes:
805: This is the default solver, and is recommended in most situations.
807: It implements the Two-level Orthogonal Arnoldi procedure {cite:p}`Lu16`,
808: which carries out an implicit linearization and operates with an
809: orthogonal Krylov basis stored in a compact form $V = (I \otimes U) S$.
810: The details of the SLEPc implementation can be found in {cite:p}`Cam16a`.
812: Level: beginner
814: .seealso: [](ch:pep), `PEP`, `PEPType`, `PEPSetType()`
815: M*/
816: SLEPC_EXTERN PetscErrorCode PEPCreate_TOAR(PEP pep)
817: {
818: PEP_TOAR *ctx;
820: PetscFunctionBegin;
821: PetscCall(PetscNew(&ctx));
822: pep->data = (void*)ctx;
824: pep->lineariz = PETSC_TRUE;
825: ctx->lock = PETSC_TRUE;
827: pep->ops->solve = PEPSolve_TOAR;
828: pep->ops->setup = PEPSetUp_TOAR;
829: pep->ops->setfromoptions = PEPSetFromOptions_TOAR;
830: pep->ops->destroy = PEPDestroy_TOAR;
831: pep->ops->view = PEPView_TOAR;
832: pep->ops->backtransform = PEPBackTransform_Default;
833: pep->ops->computevectors = PEPComputeVectors_Default;
834: pep->ops->extractvectors = PEPExtractVectors_TOAR;
836: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPTOARSetRestart_C",PEPTOARSetRestart_TOAR));
837: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPTOARGetRestart_C",PEPTOARGetRestart_TOAR));
838: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPTOARSetLocking_C",PEPTOARSetLocking_TOAR));
839: PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPTOARGetLocking_C",PEPTOARGetLocking_TOAR));
840: PetscFunctionReturn(PETSC_SUCCESS);
841: }