Actual source code: ptoar.c

slepc-3.18.2 2023-01-26
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: */
 10: /*
 11:    SLEPc 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:     PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv);
 78:     if (sinv) PEPEvaluateBasis(pep,pep->target,0,pep->solvematcoeffs,NULL);
 79:     else {
 80:       for (i=0;i<pep->nmat-1;i++) pep->solvematcoeffs[i] = 0.0;
 81:       pep->solvematcoeffs[pep->nmat-1] = 1.0;
 82:     }
 83:   }
 84:   BVDestroy(&ctx->V);
 85:   BVCreateTensor(pep->V,pep->nmat-1,&ctx->V);
 86:   return 0;
 87: }

 89: /*
 90:   Extend the TOAR basis by applying the matrix operator
 91:   over a vector which is decomposed in the TOAR way
 92:   Input:
 93:     - pbc: array containing the polynomial basis coefficients
 94:     - S,V: define the latest Arnoldi vector (nv vectors in V)
 95:   Output:
 96:     - t: new vector extending the TOAR basis
 97:     - r: temporary coefficients to compute the TOAR coefficients
 98:          for the new Arnoldi vector
 99:   Workspace: t_ (two vectors)
100: */
101: 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_)
102: {
103:   PetscInt       nmat=pep->nmat,deg=nmat-1,k,j,off=0,lss;
104:   Vec            v=t_[0],ve=t_[1],q=t_[2];
105:   PetscScalar    alpha=1.0,*ss,a;
106:   PetscReal      *ca=pep->pbc,*cb=pep->pbc+nmat,*cg=pep->pbc+2*nmat;
107:   PetscBool      flg;

109:   BVSetActiveColumns(pep->V,0,nv);
110:   STGetTransform(pep->st,&flg);
111:   if (sinvert) {
112:     for (j=0;j<nv;j++) {
113:       if (deg>1) r[lr+j] = S[j]/ca[0];
114:       if (deg>2) r[2*lr+j] = (S[ls+j]+(sigma-cb[1])*r[lr+j])/ca[1];
115:     }
116:     for (k=2;k<deg-1;k++) {
117:       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];
118:     }
119:     k = deg-1;
120:     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];
121:     ss = r; lss = lr; off = 1; alpha = -1.0; a = pep->sfactor;
122:   } else {
123:     ss = S; lss = ls; off = 0; alpha = -ca[deg-1]; a = 1.0;
124:   }
125:   BVMultVec(V,1.0,0.0,v,ss+off*lss);
126:   if (PetscUnlikely(pep->Dr)) { /* balancing */
127:     VecPointwiseMult(v,v,pep->Dr);
128:   }
129:   STMatMult(pep->st,off,v,q);
130:   VecScale(q,a);
131:   for (j=1+off;j<deg+off-1;j++) {
132:     BVMultVec(V,1.0,0.0,v,ss+j*lss);
133:     if (PetscUnlikely(pep->Dr)) VecPointwiseMult(v,v,pep->Dr);
134:     STMatMult(pep->st,j,v,t);
135:     a *= pep->sfactor;
136:     VecAXPY(q,a,t);
137:   }
138:   if (sinvert) {
139:     BVMultVec(V,1.0,0.0,v,ss);
140:     if (PetscUnlikely(pep->Dr)) VecPointwiseMult(v,v,pep->Dr);
141:     STMatMult(pep->st,deg,v,t);
142:     a *= pep->sfactor;
143:     VecAXPY(q,a,t);
144:   } else {
145:     BVMultVec(V,1.0,0.0,ve,ss+(deg-1)*lss);
146:     if (PetscUnlikely(pep->Dr)) VecPointwiseMult(ve,ve,pep->Dr);
147:     a *= pep->sfactor;
148:     STMatMult(pep->st,deg-1,ve,t);
149:     VecAXPY(q,a,t);
150:     a *= pep->sfactor;
151:   }
152:   if (flg || !sinvert) alpha /= a;
153:   STMatSolve(pep->st,q,t);
154:   VecScale(t,alpha);
155:   if (!sinvert) {
156:     VecAXPY(t,cg[deg-1],v);
157:     VecAXPY(t,cb[deg-1],ve);
158:   }
159:   if (PetscUnlikely(pep->Dr)) VecPointwiseDivide(t,t,pep->Dr);
160:   return 0;
161: }

163: /*
164:   Compute TOAR coefficients of the blocks of the new Arnoldi vector computed
165: */
166: static PetscErrorCode PEPTOARCoefficients(PEP pep,PetscBool sinvert,PetscScalar sigma,PetscInt nv,PetscScalar *S,PetscInt ls,PetscScalar *r,PetscInt lr,PetscScalar *x)
167: {
168:   PetscInt    k,j,nmat=pep->nmat,d=nmat-1;
169:   PetscReal   *ca=pep->pbc,*cb=pep->pbc+nmat,*cg=pep->pbc+2*nmat;
170:   PetscScalar t=1.0,tp=0.0,tt;

172:   if (sinvert) {
173:     for (k=1;k<d;k++) {
174:       tt = t;
175:       t = ((sigma-cb[k-1])*t-cg[k-1]*tp)/ca[k-1]; /* k-th basis polynomial */
176:       tp = tt;
177:       for (j=0;j<=nv;j++) r[k*lr+j] += t*x[j];
178:     }
179:   } else {
180:     for (j=0;j<=nv;j++) r[j] = (cb[0]-sigma)*S[j]+ca[0]*S[ls+j];
181:     for (k=1;k<d-1;k++) {
182:       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];
183:     }
184:     if (sigma!=0.0) for (j=0;j<=nv;j++) r[(d-1)*lr+j] -= sigma*S[(d-1)*ls+j];
185:   }
186:   return 0;
187: }

189: /*
190:   Compute a run of Arnoldi iterations dim(work)=ld
191: */
192: static PetscErrorCode PEPTOARrun(PEP pep,PetscScalar sigma,Mat A,PetscInt k,PetscInt *M,PetscReal *beta,PetscBool *breakdown,Vec *t_)
193: {
194:   PEP_TOAR       *ctx = (PEP_TOAR*)pep->data;
195:   PetscInt       j,m=*M,deg=pep->nmat-1,ld;
196:   PetscInt       ldh,lds,nqt,l;
197:   Vec            t;
198:   PetscReal      norm=0.0;
199:   PetscBool      flg,sinvert=PETSC_FALSE,lindep;
200:   PetscScalar    *H,*x,*S;
201:   Mat            MS;

203:   *beta = 0.0;
204:   MatDenseGetArray(A,&H);
205:   MatDenseGetLDA(A,&ldh);
206:   BVTensorGetFactors(ctx->V,NULL,&MS);
207:   MatDenseGetArray(MS,&S);
208:   BVGetSizes(pep->V,NULL,NULL,&ld);
209:   lds = ld*deg;
210:   BVGetActiveColumns(pep->V,&l,&nqt);
211:   STGetTransform(pep->st,&flg);
212:   if (!flg) {
213:     /* spectral transformation handled by the solver */
214:     PetscObjectTypeCompareAny((PetscObject)pep->st,&flg,STSINVERT,STSHIFT,"");
216:     PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinvert);
217:   }
218:   BVSetActiveColumns(ctx->V,0,m);
219:   for (j=k;j<m;j++) {
220:     /* apply operator */
221:     BVGetColumn(pep->V,nqt,&t);
222:     PEPTOARExtendBasis(pep,sinvert,sigma,S+j*lds,ld,nqt,pep->V,t,S+(j+1)*lds,ld,t_);
223:     BVRestoreColumn(pep->V,nqt,&t);

225:     /* orthogonalize */
226:     if (sinvert) x = S+(j+1)*lds;
227:     else x = S+(deg-1)*ld+(j+1)*lds;
228:     BVOrthogonalizeColumn(pep->V,nqt,x,&norm,&lindep);
229:     if (!lindep) {
230:       x[nqt] = norm;
231:       BVScaleColumn(pep->V,nqt,1.0/norm);
232:       nqt++;
233:     }

235:     PEPTOARCoefficients(pep,sinvert,sigma,nqt-1,S+j*lds,ld,S+(j+1)*lds,ld,x);

237:     /* level-2 orthogonalization */
238:     BVOrthogonalizeColumn(ctx->V,j+1,H+j*ldh,&norm,breakdown);
239:     H[j+1+ldh*j] = norm;
240:     if (PetscUnlikely(*breakdown)) {
241:       *M = j+1;
242:       break;
243:     }
244:     BVScaleColumn(ctx->V,j+1,1.0/norm);
245:     BVSetActiveColumns(pep->V,l,nqt);
246:   }
247:   *beta = norm;
248:   BVSetActiveColumns(ctx->V,0,*M);
249:   MatDenseRestoreArray(MS,&S);
250:   BVTensorRestoreFactors(ctx->V,NULL,&MS);
251:   MatDenseRestoreArray(A,&H);
252:   return 0;
253: }

255: /*
256:   Computes T_j = phi_idx(T). In T_j and T_p are phi_{idx-1}(T)
257:    and phi_{idx-2}(T) respectively or null if idx=0,1.
258:    Tp and Tj are input/output arguments
259: */
260: static PetscErrorCode PEPEvaluateBasisM(PEP pep,PetscInt k,PetscScalar *T,PetscInt ldt,PetscInt idx,PetscScalar **Tp,PetscScalar **Tj)
261: {
262:   PetscInt       i;
263:   PetscReal      *ca,*cb,*cg;
264:   PetscScalar    *pt,g,a;
265:   PetscBLASInt   k_,ldt_;

267:   if (idx==0) {
268:     PetscArrayzero(*Tj,k*k);
269:     PetscArrayzero(*Tp,k*k);
270:     for (i=0;i<k;i++) (*Tj)[i+i*k] = 1.0;
271:   } else {
272:     PetscBLASIntCast(ldt,&ldt_);
273:     PetscBLASIntCast(k,&k_);
274:     ca = pep->pbc; cb = pep->pbc+pep->nmat; cg = pep->pbc+2*pep->nmat;
275:     for (i=0;i<k;i++) T[i*ldt+i] -= cb[idx-1];
276:     a = 1/ca[idx-1];
277:     g = (idx==1)?0.0:-cg[idx-1]/ca[idx-1];
278:     PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&a,T,&ldt_,*Tj,&k_,&g,*Tp,&k_));
279:     pt = *Tj; *Tj = *Tp; *Tp = pt;
280:     for (i=0;i<k;i++) T[i*ldt+i] += cb[idx-1];
281:   }
282:   return 0;
283: }

285: static PetscErrorCode PEPExtractInvariantPair(PEP pep,PetscScalar sigma,PetscInt sr,PetscInt k,PetscScalar *S,PetscInt ld,PetscInt deg,Mat HH)
286: {
287:   PetscInt       i,j,jj,ldh,lds,ldt,d=pep->nmat-1,idxcpy=0;
288:   PetscScalar    *H,*At,*Bt,*Hj,*Hp,*T,sone=1.0,g,a,*pM,*work;
289:   PetscBLASInt   k_,sr_,lds_,ldh_,info,*p,lwork,ldt_;
290:   PetscBool      transf=PETSC_FALSE,flg;
291:   PetscReal      norm,maxnrm,*rwork;
292:   BV             *R,Y;
293:   Mat            M,*A;

295:   if (k==0) return 0;
296:   MatDenseGetArray(HH,&H);
297:   MatDenseGetLDA(HH,&ldh);
298:   lds = deg*ld;
299:   PetscCalloc6(k,&p,sr*k,&At,k*k,&Bt,k*k,&Hj,k*k,&Hp,sr*k,&work);
300:   PetscBLASIntCast(sr,&sr_);
301:   PetscBLASIntCast(k,&k_);
302:   PetscBLASIntCast(lds,&lds_);
303:   PetscBLASIntCast(ldh,&ldh_);
304:   STGetTransform(pep->st,&flg);
305:   if (!flg) {
306:     PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&flg);
307:     if (flg || sigma!=0.0) transf=PETSC_TRUE;
308:   }
309:   if (transf) {
310:     PetscMalloc1(k*k,&T);
311:     ldt = k;
312:     for (i=0;i<k;i++) PetscArraycpy(T+k*i,H+i*ldh,k);
313:     if (flg) {
314:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
315:       PetscCallBLAS("LAPACKgetrf",LAPACKgetrf_(&k_,&k_,T,&k_,p,&info));
316:       SlepcCheckLapackInfo("getrf",info);
317:       PetscBLASIntCast(sr*k,&lwork);
318:       PetscCallBLAS("LAPACKgetri",LAPACKgetri_(&k_,T,&k_,p,work,&lwork,&info));
319:       SlepcCheckLapackInfo("getri",info);
320:       PetscFPTrapPop();
321:     }
322:     if (sigma!=0.0) for (i=0;i<k;i++) T[i+k*i] += sigma;
323:   } else {
324:     T = H; ldt = ldh;
325:   }
326:   PetscBLASIntCast(ldt,&ldt_);
327:   switch (pep->extract) {
328:   case PEP_EXTRACT_NONE:
329:     break;
330:   case PEP_EXTRACT_NORM:
331:     if (pep->basis == PEP_BASIS_MONOMIAL) {
332:       PetscBLASIntCast(ldt,&ldt_);
333:       PetscMalloc1(k,&rwork);
334:       norm = LAPACKlange_("F",&k_,&k_,T,&ldt_,rwork);
335:       PetscFree(rwork);
336:       if (norm>1.0) idxcpy = d-1;
337:     } else {
338:       PetscBLASIntCast(ldt,&ldt_);
339:       PetscMalloc1(k,&rwork);
340:       maxnrm = 0.0;
341:       for (i=0;i<pep->nmat-1;i++) {
342:         PEPEvaluateBasisM(pep,k,T,ldt,i,&Hp,&Hj);
343:         norm = LAPACKlange_("F",&k_,&k_,Hj,&k_,rwork);
344:         if (norm > maxnrm) {
345:           idxcpy = i;
346:           maxnrm = norm;
347:         }
348:       }
349:       PetscFree(rwork);
350:     }
351:     if (idxcpy>0) {
352:       /* copy block idxcpy of S to the first one */
353:       for (j=0;j<k;j++) PetscArraycpy(S+j*lds,S+idxcpy*ld+j*lds,sr);
354:     }
355:     break;
356:   case PEP_EXTRACT_RESIDUAL:
357:     STGetTransform(pep->st,&flg);
358:     if (flg) {
359:       PetscMalloc1(pep->nmat,&A);
360:       for (i=0;i<pep->nmat;i++) STGetMatrixTransformed(pep->st,i,A+i);
361:     } else A = pep->A;
362:     PetscMalloc1(pep->nmat-1,&R);
363:     for (i=0;i<pep->nmat-1;i++) BVDuplicateResize(pep->V,k,R+i);
364:     BVDuplicateResize(pep->V,sr,&Y);
365:     MatCreateSeqDense(PETSC_COMM_SELF,sr,k,NULL,&M);
366:     g = 0.0; a = 1.0;
367:     BVSetActiveColumns(pep->V,0,sr);
368:     for (j=0;j<pep->nmat;j++) {
369:       BVMatMult(pep->V,A[j],Y);
370:       PEPEvaluateBasisM(pep,k,T,ldt,i,&Hp,&Hj);
371:       for (i=0;i<pep->nmat-1;i++) {
372:         PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&sr_,&k_,&k_,&a,S+i*ld,&lds_,Hj,&k_,&g,At,&sr_));
373:         MatDenseGetArray(M,&pM);
374:         for (jj=0;jj<k;jj++) PetscArraycpy(pM+jj*sr,At+jj*sr,sr);
375:         MatDenseRestoreArray(M,&pM);
376:         BVMult(R[i],1.0,(i==0)?0.0:1.0,Y,M);
377:       }
378:     }

380:     /* frobenius norm */
381:     maxnrm = 0.0;
382:     for (i=0;i<pep->nmat-1;i++) {
383:       BVNorm(R[i],NORM_FROBENIUS,&norm);
384:       if (maxnrm > norm) {
385:         maxnrm = norm;
386:         idxcpy = i;
387:       }
388:     }
389:     if (idxcpy>0) {
390:       /* copy block idxcpy of S to the first one */
391:       for (j=0;j<k;j++) PetscArraycpy(S+j*lds,S+idxcpy*ld+j*lds,sr);
392:     }
393:     if (flg) PetscFree(A);
394:     for (i=0;i<pep->nmat-1;i++) BVDestroy(&R[i]);
395:     PetscFree(R);
396:     BVDestroy(&Y);
397:     MatDestroy(&M);
398:     break;
399:   case PEP_EXTRACT_STRUCTURED:
400:     for (j=0;j<k;j++) Bt[j+j*k] = 1.0;
401:     for (j=0;j<sr;j++) {
402:       for (i=0;i<k;i++) At[j*k+i] = PetscConj(S[i*lds+j]);
403:     }
404:     PEPEvaluateBasisM(pep,k,T,ldt,0,&Hp,&Hj);
405:     for (i=1;i<deg;i++) {
406:       PEPEvaluateBasisM(pep,k,T,ldt,i,&Hp,&Hj);
407:       PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&k_,&sr_,&k_,&sone,Hj,&k_,S+i*ld,&lds_,&sone,At,&k_));
408:       PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&k_,&k_,&k_,&sone,Hj,&k_,Hj,&k_,&sone,Bt,&k_));
409:     }
410:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
411:     PetscCallBLAS("LAPACKgesv",LAPACKgesv_(&k_,&sr_,Bt,&k_,p,At,&k_,&info));
412:     PetscFPTrapPop();
413:     SlepcCheckLapackInfo("gesv",info);
414:     for (j=0;j<sr;j++) {
415:       for (i=0;i<k;i++) S[i*lds+j] = PetscConj(At[j*k+i]);
416:     }
417:     break;
418:   }
419:   if (transf) PetscFree(T);
420:   PetscFree6(p,At,Bt,Hj,Hp,work);
421:   MatDenseRestoreArray(HH,&H);
422:   return 0;
423: }

425: PetscErrorCode PEPSolve_TOAR(PEP pep)
426: {
427:   PEP_TOAR       *ctx = (PEP_TOAR*)pep->data;
428:   PetscInt       i,j,k,l,nv=0,ld,lds,nq=0,nconv=0;
429:   PetscInt       nmat=pep->nmat,deg=nmat-1;
430:   PetscScalar    *S,sigma;
431:   PetscReal      beta;
432:   PetscBool      breakdown=PETSC_FALSE,flg,falselock=PETSC_FALSE,sinv=PETSC_FALSE;
433:   Mat            H,MS,MQ;

435:   PetscCitationsRegister(citation,&cited);
436:   if (ctx->lock) {
437:     /* undocumented option to use a cheaper locking instead of the true locking */
438:     PetscOptionsGetBool(NULL,NULL,"-pep_toar_falselocking",&falselock,NULL);
439:   }
440:   STGetShift(pep->st,&sigma);

442:   /* update polynomial basis coefficients */
443:   STGetTransform(pep->st,&flg);
444:   if (pep->sfactor!=1.0) {
445:     for (i=0;i<nmat;i++) {
446:       pep->pbc[nmat+i] /= pep->sfactor;
447:       pep->pbc[2*nmat+i] /= pep->sfactor*pep->sfactor;
448:     }
449:     if (!flg) {
450:       pep->target /= pep->sfactor;
451:       RGPushScale(pep->rg,1.0/pep->sfactor);
452:       STScaleShift(pep->st,1.0/pep->sfactor);
453:       sigma /= pep->sfactor;
454:     } else {
455:       PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv);
456:       pep->target = sinv?pep->target*pep->sfactor:pep->target/pep->sfactor;
457:       RGPushScale(pep->rg,sinv?pep->sfactor:1.0/pep->sfactor);
458:       STScaleShift(pep->st,sinv?pep->sfactor:1.0/pep->sfactor);
459:     }
460:   }

462:   if (flg) sigma = 0.0;

464:   /* clean projected matrix (including the extra-arrow) */
465:   DSSetDimensions(pep->ds,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
466:   DSGetMat(pep->ds,DS_MAT_A,&H);
467:   MatZeroEntries(H);
468:   DSRestoreMat(pep->ds,DS_MAT_A,&H);

470:   /* Get the starting Arnoldi vector */
471:   BVTensorBuildFirstColumn(ctx->V,pep->nini);

473:   /* restart loop */
474:   l = 0;
475:   while (pep->reason == PEP_CONVERGED_ITERATING) {
476:     pep->its++;

478:     /* compute an nv-step Lanczos factorization */
479:     nv = PetscMax(PetscMin(nconv+pep->mpd,pep->ncv),nv);
480:     DSGetMat(pep->ds,DS_MAT_A,&H);
481:     PEPTOARrun(pep,sigma,H,pep->nconv+l,&nv,&beta,&breakdown,pep->work);
482:     DSRestoreMat(pep->ds,DS_MAT_A,&H);
483:     DSSetDimensions(pep->ds,nv,pep->nconv,pep->nconv+l);
484:     DSSetState(pep->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE);
485:     BVSetActiveColumns(ctx->V,pep->nconv,nv);

487:     /* solve projected problem */
488:     DSSolve(pep->ds,pep->eigr,pep->eigi);
489:     DSSort(pep->ds,pep->eigr,pep->eigi,NULL,NULL,NULL);
490:     DSUpdateExtraRow(pep->ds);
491:     DSSynchronize(pep->ds,pep->eigr,pep->eigi);

493:     /* check convergence */
494:     PEPKrylovConvergence(pep,PETSC_FALSE,pep->nconv,nv-pep->nconv,beta,&k);
495:     (*pep->stopping)(pep,pep->its,pep->max_it,k,pep->nev,&pep->reason,pep->stoppingctx);

497:     /* update l */
498:     if (pep->reason != PEP_CONVERGED_ITERATING || breakdown) l = 0;
499:     else {
500:       l = (nv==k)?0:PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
501:       DSGetTruncateSize(pep->ds,k,nv,&l);
502:       if (!breakdown) {
503:         /* prepare the Rayleigh quotient for restart */
504:         DSTruncate(pep->ds,k+l,PETSC_FALSE);
505:       }
506:     }
507:     nconv = k;
508:     if (!ctx->lock && pep->reason == PEP_CONVERGED_ITERATING && !breakdown) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
509:     if (l) PetscInfo(pep,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l);

511:     /* update S */
512:     DSGetMat(pep->ds,DS_MAT_Q,&MQ);
513:     BVMultInPlace(ctx->V,MQ,pep->nconv,k+l);
514:     DSRestoreMat(pep->ds,DS_MAT_Q,&MQ);

516:     /* copy last column of S */
517:     BVCopyColumn(ctx->V,nv,k+l);

519:     if (PetscUnlikely(breakdown && pep->reason == PEP_CONVERGED_ITERATING)) {
520:       /* stop if breakdown */
521:       PetscInfo(pep,"Breakdown TOAR method (it=%" PetscInt_FMT " norm=%g)\n",pep->its,(double)beta);
522:       pep->reason = PEP_DIVERGED_BREAKDOWN;
523:     }
524:     if (pep->reason != PEP_CONVERGED_ITERATING) l--;
525:     /* truncate S */
526:     BVGetActiveColumns(pep->V,NULL,&nq);
527:     if (k+l+deg<=nq) {
528:       BVSetActiveColumns(ctx->V,pep->nconv,k+l+1);
529:       if (!falselock && ctx->lock) BVTensorCompress(ctx->V,k-pep->nconv);
530:       else BVTensorCompress(ctx->V,0);
531:     }
532:     pep->nconv = k;
533:     PEPMonitor(pep,pep->its,nconv,pep->eigr,pep->eigi,pep->errest,nv);
534:   }
535:   if (pep->nconv>0) {
536:     /* {V*S_nconv^i}_{i=0}^{d-1} has rank nconv instead of nconv+d-1. Force zeros in each S_nconv^i block */
537:     BVSetActiveColumns(ctx->V,0,pep->nconv);
538:     BVGetActiveColumns(pep->V,NULL,&nq);
539:     BVSetActiveColumns(pep->V,0,nq);
540:     if (nq>pep->nconv) {
541:       BVTensorCompress(ctx->V,pep->nconv);
542:       BVSetActiveColumns(pep->V,0,pep->nconv);
543:       nq = pep->nconv;
544:     }

546:     /* perform Newton refinement if required */
547:     if (pep->refine==PEP_REFINE_MULTIPLE && pep->rits>0) {
548:       /* extract invariant pair */
549:       BVTensorGetFactors(ctx->V,NULL,&MS);
550:       MatDenseGetArray(MS,&S);
551:       DSGetMat(pep->ds,DS_MAT_A,&H);
552:       BVGetSizes(pep->V,NULL,NULL,&ld);
553:       lds = deg*ld;
554:       PEPExtractInvariantPair(pep,sigma,nq,pep->nconv,S,ld,deg,H);
555:       DSRestoreMat(pep->ds,DS_MAT_A,&H);
556:       DSSetDimensions(pep->ds,pep->nconv,0,0);
557:       DSSetState(pep->ds,DS_STATE_RAW);
558:       PEPNewtonRefinement_TOAR(pep,sigma,&pep->rits,NULL,pep->nconv,S,lds);
559:       DSSolve(pep->ds,pep->eigr,pep->eigi);
560:       DSSort(pep->ds,pep->eigr,pep->eigi,NULL,NULL,NULL);
561:       DSSynchronize(pep->ds,pep->eigr,pep->eigi);
562:       DSGetMat(pep->ds,DS_MAT_Q,&MQ);
563:       BVMultInPlace(ctx->V,MQ,0,pep->nconv);
564:       DSRestoreMat(pep->ds,DS_MAT_Q,&MQ);
565:       MatDenseRestoreArray(MS,&S);
566:       BVTensorRestoreFactors(ctx->V,NULL,&MS);
567:     }
568:   }
569:   STGetTransform(pep->st,&flg);
570:   if (pep->refine!=PEP_REFINE_MULTIPLE || pep->rits==0) {
571:     if (!flg) PetscTryTypeMethod(pep,backtransform);
572:     if (pep->sfactor!=1.0) {
573:       for (j=0;j<pep->nconv;j++) {
574:         pep->eigr[j] *= pep->sfactor;
575:         pep->eigi[j] *= pep->sfactor;
576:       }
577:       /* restore original values */
578:       for (i=0;i<pep->nmat;i++) {
579:         pep->pbc[pep->nmat+i] *= pep->sfactor;
580:         pep->pbc[2*pep->nmat+i] *= pep->sfactor*pep->sfactor;
581:       }
582:     }
583:   }
584:   /* restore original values */
585:   if (!flg) {
586:     pep->target *= pep->sfactor;
587:     STScaleShift(pep->st,pep->sfactor);
588:   } else {
589:     STScaleShift(pep->st,sinv?1.0/pep->sfactor:pep->sfactor);
590:     pep->target = (sinv)?pep->target/pep->sfactor:pep->target*pep->sfactor;
591:   }
592:   if (pep->sfactor!=1.0) RGPopScale(pep->rg);

594:   /* change the state to raw so that DSVectors() computes eigenvectors from scratch */
595:   DSSetDimensions(pep->ds,pep->nconv,0,0);
596:   DSSetState(pep->ds,DS_STATE_RAW);
597:   return 0;
598: }

600: static PetscErrorCode PEPTOARSetRestart_TOAR(PEP pep,PetscReal keep)
601: {
602:   PEP_TOAR *ctx = (PEP_TOAR*)pep->data;

604:   if (keep==PETSC_DEFAULT) ctx->keep = 0.5;
605:   else {
607:     ctx->keep = keep;
608:   }
609:   return 0;
610: }

612: /*@
613:    PEPTOARSetRestart - Sets the restart parameter for the TOAR
614:    method, in particular the proportion of basis vectors that must be kept
615:    after restart.

617:    Logically Collective on pep

619:    Input Parameters:
620: +  pep  - the eigenproblem solver context
621: -  keep - the number of vectors to be kept at restart

623:    Options Database Key:
624: .  -pep_toar_restart - Sets the restart parameter

626:    Notes:
627:    Allowed values are in the range [0.1,0.9]. The default is 0.5.

629:    Level: advanced

631: .seealso: PEPTOARGetRestart()
632: @*/
633: PetscErrorCode PEPTOARSetRestart(PEP pep,PetscReal keep)
634: {
637:   PetscTryMethod(pep,"PEPTOARSetRestart_C",(PEP,PetscReal),(pep,keep));
638:   return 0;
639: }

641: static PetscErrorCode PEPTOARGetRestart_TOAR(PEP pep,PetscReal *keep)
642: {
643:   PEP_TOAR *ctx = (PEP_TOAR*)pep->data;

645:   *keep = ctx->keep;
646:   return 0;
647: }

649: /*@
650:    PEPTOARGetRestart - Gets the restart parameter used in the TOAR method.

652:    Not Collective

654:    Input Parameter:
655: .  pep - the eigenproblem solver context

657:    Output Parameter:
658: .  keep - the restart parameter

660:    Level: advanced

662: .seealso: PEPTOARSetRestart()
663: @*/
664: PetscErrorCode PEPTOARGetRestart(PEP pep,PetscReal *keep)
665: {
668:   PetscUseMethod(pep,"PEPTOARGetRestart_C",(PEP,PetscReal*),(pep,keep));
669:   return 0;
670: }

672: static PetscErrorCode PEPTOARSetLocking_TOAR(PEP pep,PetscBool lock)
673: {
674:   PEP_TOAR *ctx = (PEP_TOAR*)pep->data;

676:   ctx->lock = lock;
677:   return 0;
678: }

680: /*@
681:    PEPTOARSetLocking - Choose between locking and non-locking variants of
682:    the TOAR method.

684:    Logically Collective on pep

686:    Input Parameters:
687: +  pep  - the eigenproblem solver context
688: -  lock - true if the locking variant must be selected

690:    Options Database Key:
691: .  -pep_toar_locking - Sets the locking flag

693:    Notes:
694:    The default is to lock converged eigenpairs when the method restarts.
695:    This behaviour can be changed so that all directions are kept in the
696:    working subspace even if already converged to working accuracy (the
697:    non-locking variant).

699:    Level: advanced

701: .seealso: PEPTOARGetLocking()
702: @*/
703: PetscErrorCode PEPTOARSetLocking(PEP pep,PetscBool lock)
704: {
707:   PetscTryMethod(pep,"PEPTOARSetLocking_C",(PEP,PetscBool),(pep,lock));
708:   return 0;
709: }

711: static PetscErrorCode PEPTOARGetLocking_TOAR(PEP pep,PetscBool *lock)
712: {
713:   PEP_TOAR *ctx = (PEP_TOAR*)pep->data;

715:   *lock = ctx->lock;
716:   return 0;
717: }

719: /*@
720:    PEPTOARGetLocking - Gets the locking flag used in the TOAR method.

722:    Not Collective

724:    Input Parameter:
725: .  pep - the eigenproblem solver context

727:    Output Parameter:
728: .  lock - the locking flag

730:    Level: advanced

732: .seealso: PEPTOARSetLocking()
733: @*/
734: PetscErrorCode PEPTOARGetLocking(PEP pep,PetscBool *lock)
735: {
738:   PetscUseMethod(pep,"PEPTOARGetLocking_C",(PEP,PetscBool*),(pep,lock));
739:   return 0;
740: }

742: PetscErrorCode PEPSetFromOptions_TOAR(PEP pep,PetscOptionItems *PetscOptionsObject)
743: {
744:   PetscBool      flg,lock;
745:   PetscReal      keep;

747:   PetscOptionsHeadBegin(PetscOptionsObject,"PEP TOAR Options");

749:     PetscOptionsReal("-pep_toar_restart","Proportion of vectors kept after restart","PEPTOARSetRestart",0.5,&keep,&flg);
750:     if (flg) PEPTOARSetRestart(pep,keep);

752:     PetscOptionsBool("-pep_toar_locking","Choose between locking and non-locking variants","PEPTOARSetLocking",PETSC_FALSE,&lock,&flg);
753:     if (flg) PEPTOARSetLocking(pep,lock);

755:   PetscOptionsHeadEnd();
756:   return 0;
757: }

759: PetscErrorCode PEPView_TOAR(PEP pep,PetscViewer viewer)
760: {
761:   PEP_TOAR       *ctx = (PEP_TOAR*)pep->data;
762:   PetscBool      isascii;

764:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
765:   if (isascii) {
766:     PetscViewerASCIIPrintf(viewer,"  %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep));
767:     PetscViewerASCIIPrintf(viewer,"  using the %slocking variant\n",ctx->lock?"":"non-");
768:   }
769:   return 0;
770: }

772: PetscErrorCode PEPDestroy_TOAR(PEP pep)
773: {
774:   PEP_TOAR       *ctx = (PEP_TOAR*)pep->data;

776:   BVDestroy(&ctx->V);
777:   PetscFree(pep->data);
778:   PetscObjectComposeFunction((PetscObject)pep,"PEPTOARSetRestart_C",NULL);
779:   PetscObjectComposeFunction((PetscObject)pep,"PEPTOARGetRestart_C",NULL);
780:   PetscObjectComposeFunction((PetscObject)pep,"PEPTOARSetLocking_C",NULL);
781:   PetscObjectComposeFunction((PetscObject)pep,"PEPTOARGetLocking_C",NULL);
782:   return 0;
783: }

785: SLEPC_EXTERN PetscErrorCode PEPCreate_TOAR(PEP pep)
786: {
787:   PEP_TOAR       *ctx;

789:   PetscNew(&ctx);
790:   pep->data = (void*)ctx;

792:   pep->lineariz = PETSC_TRUE;
793:   ctx->lock     = PETSC_TRUE;

795:   pep->ops->solve          = PEPSolve_TOAR;
796:   pep->ops->setup          = PEPSetUp_TOAR;
797:   pep->ops->setfromoptions = PEPSetFromOptions_TOAR;
798:   pep->ops->destroy        = PEPDestroy_TOAR;
799:   pep->ops->view           = PEPView_TOAR;
800:   pep->ops->backtransform  = PEPBackTransform_Default;
801:   pep->ops->computevectors = PEPComputeVectors_Default;
802:   pep->ops->extractvectors = PEPExtractVectors_TOAR;

804:   PetscObjectComposeFunction((PetscObject)pep,"PEPTOARSetRestart_C",PEPTOARSetRestart_TOAR);
805:   PetscObjectComposeFunction((PetscObject)pep,"PEPTOARGetRestart_C",PEPTOARGetRestart_TOAR);
806:   PetscObjectComposeFunction((PetscObject)pep,"PEPTOARSetLocking_C",PEPTOARSetLocking_TOAR);
807:   PetscObjectComposeFunction((PetscObject)pep,"PEPTOARGetLocking_C",PEPTOARGetLocking_TOAR);
808:   return 0;
809: }