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