Actual source code: stoar.c

slepc-3.9.2 2018-07-02
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2018, 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: "stoar"

 13:    Method: S-TOAR

 15:    Algorithm:

 17:        Symmetric Two-Level Orthogonal Arnoldi.

 19:    References:

 21:        [1] C. Campos and J.E. Roman, "Restarted Q-Arnoldi-type methods
 22:            exploiting symmetry in quadratic eigenvalue problems", BIT
 23:            Numer. Math. 56(4):1213-1236, 2016.
 24: */

 26: #include <slepc/private/pepimpl.h>         /*I "slepcpep.h" I*/
 27: #include "../src/pep/impls/krylov/pepkrylov.h"
 28: #include <slepcblaslapack.h>

 30: static PetscBool  cited = PETSC_FALSE;
 31: static const char citation[] =
 32:   "@Article{slepc-stoar,\n"
 33:   "   author = \"C. Campos and J. E. Roman\",\n"
 34:   "   title = \"Restarted {Q-Arnoldi-type} methods exploiting symmetry in quadratic eigenvalue problems\",\n"
 35:   "   journal = \"{BIT} Numer. Math.\",\n"
 36:   "   volume = \"56\",\n"
 37:   "   number = \"4\",\n"
 38:   "   pages = \"1213--1236\",\n"
 39:   "   year = \"2016,\"\n"
 40:   "   doi = \"https://doi.org/10.1007/s10543-016-0601-5\"\n"
 41:   "}\n";


 44: PetscErrorCode MatMult_Func(Mat A,Vec x,Vec y)
 45: {
 47:   ShellMatCtx    *ctx;

 50:   MatShellGetContext(A,&ctx);
 51:   MatMult(ctx->A,x,y);
 52:   VecScale(y,ctx->scal);
 53:   return(0);
 54: }

 56: PetscErrorCode MatDestroy_Func(Mat A)
 57: {
 58:   ShellMatCtx    *ctx;

 62:   MatShellGetContext(A,(void**)&ctx);
 63:   PetscFree(ctx);
 64:   return(0);
 65: }

 67: PetscErrorCode PEPSetUp_STOAR(PEP pep)
 68: {
 69:   PetscErrorCode    ierr;
 70:   PetscBool         shift,sinv,flg;
 71:   PEP_TOAR          *ctx = (PEP_TOAR*)pep->data;
 72:   PetscInt          ld;
 73:   PetscReal         eta;
 74:   BVOrthogType      otype;
 75:   BVOrthogBlockType obtype;

 78:   if (pep->problem_type!=PEP_HERMITIAN && pep->problem_type!=PEP_HYPERBOLIC) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Requested method is only available for Hermitian (or hyperbolic) problems");
 79:   if (pep->nmat!=3) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Solver only available for quadratic problems");
 80:   if (pep->basis!=PEP_BASIS_MONOMIAL) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Solver not implemented for non-monomial bases");
 81:   /* spectrum slicing requires special treatment of default values */
 82:   if (pep->which==PEP_ALL) {
 83:     if (pep->stopping!=PEPStoppingBasic) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Spectrum slicing does not support user-defined stopping test");
 84:     pep->ops->solve = PEPSolve_STOAR_QSlice;
 85:     pep->ops->extractvectors = NULL;
 86:     pep->ops->setdefaultst   = NULL;
 87:     PEPSetUp_STOAR_QSlice(pep);
 88:   } else {
 89:     PEPSetDimensions_Default(pep,pep->nev,&pep->ncv,&pep->mpd);
 90:     if (!ctx->lock && pep->mpd<pep->ncv) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Should not use mpd parameter in non-locking variant");
 91:     if (!pep->max_it) pep->max_it = PetscMax(100,2*(pep->nmat-1)*pep->n/pep->ncv);
 92:     pep->ops->solve = PEPSolve_STOAR;
 93:     ld   = pep->ncv+2;
 94:     DSSetType(pep->ds,DSGHIEP);
 95:     DSSetCompact(pep->ds,PETSC_TRUE);
 96:     DSAllocate(pep->ds,ld);
 97:     STGetTransform(pep->st,&flg);
 98:     if (!flg) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Solver requires the ST transformation flag set, see STSetTransform()");
 99:   }

101:   pep->lineariz = PETSC_TRUE;
102:   PetscObjectTypeCompare((PetscObject)pep->st,STSHIFT,&shift);
103:   PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv);
104:   if (!shift && !sinv) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Only STSHIFT and STSINVERT spectral transformations can be used");
105:   if (!pep->which) {
106:     if (sinv) pep->which = PEP_TARGET_MAGNITUDE;
107:     else pep->which = PEP_LARGEST_MAGNITUDE;
108:   }


111:   PEPAllocateSolution(pep,2);
112:   PEPSetWorkVecs(pep,4);
113:   BVDestroy(&ctx->V);
114:   BVCreateTensor(pep->V,pep->nmat-1,&ctx->V);
115:   BVGetOrthogonalization(pep->V,&otype,NULL,&eta,&obtype);
116:   BVSetOrthogonalization(ctx->V,otype,BV_ORTHOG_REFINE_ALWAYS,eta,obtype);
117:   return(0);
118: }

120: /*
121:   Compute a run of Lanczos iterations. dim(work)=(ctx->ld)*4
122: */
123: static PetscErrorCode PEPSTOARrun(PEP pep,PetscReal *a,PetscReal *b,PetscReal *omega,PetscInt k,PetscInt *M,PetscBool *breakdown,PetscBool *symmlost,Vec *t_)
124: {
126:   PEP_TOAR       *ctx = (PEP_TOAR*)pep->data;
127:   PetscInt       i,j,m=*M,l,lock;
128:   PetscInt       lds,d,ld,offq,nqt;
129:   Vec            v=t_[0],t=t_[1],q=t_[2];
130:   PetscReal      norm,sym=0.0,fro=0.0,*f;
131:   PetscScalar    *y,*S;
132:   PetscBLASInt   j_,one=1;
133:   PetscBool      lindep;
134:   Mat            MS;

137:   PetscMalloc1(*M,&y);
138:   BVGetSizes(pep->V,NULL,NULL,&ld);
139:   BVTensorGetDegree(ctx->V,&d);
140:   BVGetActiveColumns(pep->V,&lock,&nqt);
141:   lds = d*ld;
142:   offq = ld;
143:   *breakdown = PETSC_FALSE; /* ----- */
144:   DSGetDimensions(pep->ds,NULL,NULL,&l,NULL,NULL);
145:   BVSetActiveColumns(ctx->V,0,m);
146:   BVSetActiveColumns(pep->V,0,nqt);
147:   for (j=k;j<m;j++) {
148:     /* apply operator */
149:     BVTensorGetFactors(ctx->V,NULL,&MS);
150:     MatDenseGetArray(MS,&S);
151:     BVGetColumn(pep->V,nqt,&t);
152:     BVMultVec(pep->V,1.0,0.0,v,S+j*lds);
153:     STMatMult(pep->st,0,v,t);
154:     BVMultVec(pep->V,1.0,0.0,v,S+offq+j*lds);
155:     STMatMult(pep->st,1,v,q);
156:     VecAXPY(q,pep->sfactor,t);
157:     STMatSolve(pep->st,q,t);
158:     VecScale(t,-1.0/(pep->sfactor*pep->sfactor));
159:     BVRestoreColumn(pep->V,nqt,&t);

161:     /* orthogonalize */
162:     BVOrthogonalizeColumn(pep->V,nqt,S+offq+(j+1)*lds,&norm,&lindep);
163:     if (!lindep) {
164:       *(S+offq+(j+1)*lds+nqt) = norm;
165:       BVScaleColumn(pep->V,nqt,1.0/norm);
166:       nqt++;
167:     }
168:     for (i=0;i<=nqt-1;i++) *(S+(j+1)*lds+i) = *(S+offq+j*lds+i);
169:     BVSetActiveColumns(pep->V,0,nqt);
170:     MatDenseRestoreArray(MS,&S);
171:     BVTensorRestoreFactors(ctx->V,NULL,&MS);

173:     /* level-2 orthogonalization */
174:     BVOrthogonalizeColumn(ctx->V,j+1,y,&norm,&lindep);
175:     a[j] = PetscRealPart(y[j]);
176:     omega[j+1] = (norm > 0)?1.0:-1.0;
177:     BVScaleColumn(ctx->V,j+1,1.0/norm);
178:     b[j] = PetscAbsReal(norm);

180:     /* check symmetry */
181:     DSGetArrayReal(pep->ds,DS_MAT_T,&f);
182:     if (j==k) {
183:       for (i=l;i<j-1;i++) y[i] = PetscAbsScalar(y[i])-PetscAbsReal(f[2*ld+i]);
184:       for (i=0;i<l;i++) y[i] = 0.0;
185:     }
186:     DSRestoreArrayReal(pep->ds,DS_MAT_T,&f);
187:     if (j>0) y[j-1] = PetscAbsScalar(y[j-1])-PetscAbsReal(b[j-1]);
188:     PetscBLASIntCast(j,&j_);
189:     sym = SlepcAbs(BLASnrm2_(&j_,y,&one),sym);
190:     fro = SlepcAbs(fro,SlepcAbs(a[j],b[j]));
191:     if (j>0) fro = SlepcAbs(fro,b[j-1]);
192:     if (sym/fro>PetscMax(PETSC_SQRT_MACHINE_EPSILON,10*pep->tol)) {
193:       *symmlost = PETSC_TRUE;
194:       *M=j;
195:       break;
196:     }
197:   }
198:   BVSetActiveColumns(pep->V,lock,nqt);
199:   BVSetActiveColumns(ctx->V,0,*M);
200:   PetscFree(y);
201:   return(0);
202: }

204: #if 0
205: static PetscErrorCode PEPSTOARpreKConvergence(PEP pep,PetscInt nv,PetscReal *norm,Vec *w)
206: {
208:   PEP_TOAR      *ctx = (PEP_TOAR*)pep->data;
209:   PetscBLASInt   n_,one=1;
210:   PetscInt       lds=2*ctx->ld;
211:   PetscReal      t1,t2;
212:   PetscScalar    *S=ctx->S;

215:   PetscBLASIntCast(nv+2,&n_);
216:   t1 = BLASnrm2_(&n_,S+nv*2*ctx->ld,&one);
217:   t2 = BLASnrm2_(&n_,S+(nv*2+1)*ctx->ld,&one);
218:   *norm = SlepcAbs(t1,t2);
219:   BVSetActiveColumns(pep->V,0,nv+2);
220:   BVMultVec(pep->V,1.0,0.0,w[1],S+nv*lds);
221:   STMatMult(pep->st,0,w[1],w[2]);
222:   VecNorm(w[2],NORM_2,&t1);
223:   BVMultVec(pep->V,1.0,0.0,w[1],S+ctx->ld+nv*lds);
224:   STMatMult(pep->st,2,w[1],w[2]);
225:   VecNorm(w[2],NORM_2,&t2);
226:   t2 *= pep->sfactor*pep->sfactor;
227:   *norm = PetscMax(*norm,SlepcAbs(t1,t2));
228:   return(0);
229: }
230: #endif

232: PetscErrorCode PEPSolve_STOAR(PEP pep)
233: {
235:   PEP_TOAR       *ctx = (PEP_TOAR*)pep->data;
236:   PetscInt       j,k,l,nv=0,ld,ldds,t,nq=0,m,n;
237:   PetscInt       nconv=0,deg=pep->nmat-1;
238:   PetscScalar    *Q,*om,scal[2];
239:   PetscReal      beta,norm=1.0,*omega,*a,*b,*r;
240:   PetscBool      breakdown,symmlost=PETSC_FALSE,sinv,falselock=PETSC_TRUE;
241:   Mat            MQ,A,pA[4],As[2],D[2];
242:   Vec            vomega;
243:   ShellMatCtx    *ctxMat[2];

246:   PetscCitationsRegister(citation,&cited);
247:   STGetMatrixTransformed(pep->st,2,&D[1]); /* M */
248:   MatGetLocalSize(D[1],&m,&n);
249:   STGetMatrixTransformed(pep->st,0,&D[0]); /* K */
250:   scal[0] = -1.0; scal[1] = pep->sfactor*pep->sfactor;
251:   for (j=0;j<2;j++) {
252:     PetscNew(ctxMat+j);
253:     (ctxMat[j])->A = D[j]; (ctxMat[j])->scal = scal[j];
254:     MatCreateShell(PetscObjectComm((PetscObject)pep),m,n,PETSC_DETERMINE,PETSC_DETERMINE,ctxMat[j],&As[j]);
255:     MatShellSetOperation(As[j],MATOP_MULT,(void(*)())MatMult_Func);
256:     MatShellSetOperation(As[j],MATOP_DESTROY,(void(*)())MatDestroy_Func);
257:   }
258:   pA[0] = As[0]; pA[1] = pA[2] = NULL; pA[3] = As[1];
259:   MatCreateNest(PetscObjectComm((PetscObject)pep),2,NULL,2,NULL,pA,&A);
260:   for (j=0;j<2;j++) { MatDestroy(&As[j]); }
261:   BVSetMatrix(ctx->V,A,PETSC_TRUE);
262:   MatDestroy(&A);
263:   if (ctx->lock) {
264:     PetscOptionsGetBool(NULL,NULL,"-pep_stoar_falselocking",&falselock,NULL);
265:   }
266:   BVGetSizes(pep->V,NULL,NULL,&ld);
267:   PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv);
268:   RGPushScale(pep->rg,sinv?pep->sfactor:1.0/pep->sfactor);
269:   STScaleShift(pep->st,sinv?pep->sfactor:1.0/pep->sfactor);

271:   /* Get the starting Arnoldi vector */
272:   BVTensorBuildFirstColumn(ctx->V,pep->nini);
273:   DSGetArrayReal(pep->ds,DS_MAT_D,&omega);
274:   VecCreateSeq(PETSC_COMM_SELF,1,&vomega);
275:   BVSetActiveColumns(ctx->V,0,1);
276:   BVGetSignature(ctx->V,vomega);
277:   VecGetArray(vomega,&om);
278:   omega[0] = PetscRealPart(om[0]);
279:   VecRestoreArray(vomega,&om);
280:   DSRestoreArrayReal(pep->ds,DS_MAT_D,&omega);
281:   VecDestroy(&vomega);

283:   /* Restart loop */
284:   l = 0;
285:   DSGetLeadingDimension(pep->ds,&ldds);
286:   while (pep->reason == PEP_CONVERGED_ITERATING) {
287:     pep->its++;
288:     DSGetArrayReal(pep->ds,DS_MAT_T,&a);
289:     b = a+ldds;
290:     DSGetArrayReal(pep->ds,DS_MAT_D,&omega);

292:     /* Compute an nv-step Lanczos factorization */
293:     nv = PetscMin(pep->nconv+pep->mpd,pep->ncv);
294:     PEPSTOARrun(pep,a,b,omega,pep->nconv+l,&nv,&breakdown,&symmlost,pep->work);
295:     beta = b[nv-1];
296:     if (symmlost && nv==pep->nconv+l) {
297:       pep->reason = PEP_DIVERGED_SYMMETRY_LOST;
298:       pep->nconv = nconv;
299:       if (falselock || !ctx->lock) {
300:        BVSetActiveColumns(ctx->V,0,pep->nconv);
301:        BVTensorCompress(ctx->V,0);
302:       }
303:       break;
304:     }
305:     DSRestoreArrayReal(pep->ds,DS_MAT_T,&a);
306:     DSRestoreArrayReal(pep->ds,DS_MAT_D,&omega);
307:     DSSetDimensions(pep->ds,nv,0,pep->nconv,pep->nconv+l);
308:     if (l==0) {
309:       DSSetState(pep->ds,DS_STATE_INTERMEDIATE);
310:     } else {
311:       DSSetState(pep->ds,DS_STATE_RAW);
312:     }

314:     /* Solve projected problem */
315:     DSSolve(pep->ds,pep->eigr,pep->eigi);
316:     DSSort(pep->ds,pep->eigr,pep->eigi,NULL,NULL,NULL);
317:     DSSynchronize(pep->ds,pep->eigr,pep->eigi);

319:     /* Check convergence */
320:     /* PEPSTOARpreKConvergence(pep,nv,&norm,pep->work);*/
321:     norm = 1.0;
322:     DSGetDimensions(pep->ds,NULL,NULL,NULL,NULL,&t);
323:     PEPKrylovConvergence(pep,PETSC_FALSE,pep->nconv,t-pep->nconv,PetscAbsReal(beta)*norm,&k);
324:     (*pep->stopping)(pep,pep->its,pep->max_it,k,pep->nev,&pep->reason,pep->stoppingctx);

326:     /* Update l */
327:     if (pep->reason != PEP_CONVERGED_ITERATING || breakdown) l = 0;
328:     else {
329:       l = PetscMax(1,(PetscInt)((nv-k)/2));
330:       l = PetscMin(l,t);
331:       if (!breakdown) {
332:         DSGetArrayReal(pep->ds,DS_MAT_T,&a);
333:         if (*(a+ldds+k+l-1)!=0) {
334:           if (k+l<nv-1) l = l+1;
335:           else l = l-1;
336:         }
337:         /* Prepare the Rayleigh quotient for restart */
338:         DSGetArray(pep->ds,DS_MAT_Q,&Q);
339:         DSGetArrayReal(pep->ds,DS_MAT_D,&omega);
340:         r = a + 2*ldds;
341:         for (j=k;j<k+l;j++) {
342:           r[j] = PetscRealPart(Q[nv-1+j*ldds]*beta);
343:         }
344:         b = a+ldds;
345:         b[k+l-1] = r[k+l-1];
346:         omega[k+l] = omega[nv];
347:         DSRestoreArray(pep->ds,DS_MAT_Q,&Q);
348:         DSRestoreArrayReal(pep->ds,DS_MAT_T,&a);
349:         DSRestoreArrayReal(pep->ds,DS_MAT_D,&omega);
350:       }
351:     }
352:     nconv = k;
353:     if (!ctx->lock && pep->reason == PEP_CONVERGED_ITERATING && !breakdown) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */

355:     /* Update S */
356:     DSGetMat(pep->ds,DS_MAT_Q,&MQ);
357:     BVMultInPlace(ctx->V,MQ,pep->nconv,k+l);
358:     MatDestroy(&MQ);

360:     /* Copy last column of S */
361:     BVCopyColumn(ctx->V,nv,k+l);
362:     DSGetArrayReal(pep->ds,DS_MAT_D,&omega);
363:     VecCreateSeq(PETSC_COMM_SELF,k+l,&vomega);
364:     VecGetArray(vomega,&om);
365:     for (j=0;j<k+l;j++) om[j] = omega[j];
366:     VecRestoreArray(vomega,&om);
367:     BVSetActiveColumns(ctx->V,0,k+l);
368:     BVSetSignature(ctx->V,vomega);
369:     VecDestroy(&vomega);
370:     DSRestoreArrayReal(pep->ds,DS_MAT_D,&omega);

372:     if (breakdown && pep->reason == PEP_CONVERGED_ITERATING) {
373:       /* stop if breakdown */
374:       PetscInfo2(pep,"Breakdown TOAR method (it=%D norm=%g)\n",pep->its,(double)beta);
375:       pep->reason = PEP_DIVERGED_BREAKDOWN;
376:     }
377:     if (pep->reason != PEP_CONVERGED_ITERATING) l--; 
378:     BVGetActiveColumns(pep->V,NULL,&nq);
379:     if (k+l+deg<=nq) {
380:       BVSetActiveColumns(ctx->V,pep->nconv,k+l+1);
381:       if (!falselock && ctx->lock) {
382:         BVTensorCompress(ctx->V,k-pep->nconv);
383:       } else {
384:         BVTensorCompress(ctx->V,0);
385:       }
386:     }
387:     pep->nconv = k;
388:     PEPMonitor(pep,pep->its,nconv,pep->eigr,pep->eigi,pep->errest,nv);
389:   }

391:   if (pep->nconv>0) {
392:     BVSetActiveColumns(ctx->V,0,pep->nconv);
393:     BVGetActiveColumns(pep->V,NULL,&nq);
394:     BVSetActiveColumns(pep->V,0,nq);
395:     if (nq>pep->nconv) {
396:       BVTensorCompress(ctx->V,pep->nconv);
397:       BVSetActiveColumns(pep->V,0,pep->nconv);
398:     }
399:     for (j=0;j<pep->nconv;j++) {
400:       pep->eigr[j] *= pep->sfactor;
401:       pep->eigi[j] *= pep->sfactor;
402:     }
403:   }
404:   STScaleShift(pep->st,sinv?1.0/pep->sfactor:pep->sfactor);
405:   RGPopScale(pep->rg);

407:   /* truncate Schur decomposition and change the state to raw so that
408:      DSVectors() computes eigenvectors from scratch */
409:   DSSetDimensions(pep->ds,pep->nconv,0,0,0);
410:   DSSetState(pep->ds,DS_STATE_RAW);
411:   return(0);
412: }

414: PetscErrorCode PEPSetFromOptions_STOAR(PetscOptionItems *PetscOptionsObject,PEP pep)
415: {
417:   PetscBool      flg,lock,b,f1,f2,f3;
418:   PetscInt       i,j,k;
419:   PEP_TOAR *ctx = (PEP_TOAR*)pep->data;

422:   PetscOptionsHead(PetscOptionsObject,"PEP STOAR Options");

424:   PetscOptionsBool("-pep_stoar_locking","Choose between locking and non-locking variants","PEPSTOARSetLocking",PETSC_FALSE,&lock,&flg);
425:   if (flg) { PEPSTOARSetLocking(pep,lock); }

427:   b = ctx->detect;
428:   PetscOptionsBool("-pep_stoar_detect_zeros","Check zeros during factorizations at interval boundaries","PEPSTOARSetDetectZeros",ctx->detect,&b,&flg);
429:   if (flg) { PEPSTOARSetDetectZeros(pep,b); }

431:   i = 1;
432:   j = k = PETSC_DECIDE;
433:   PetscOptionsInt("-pep_stoar_nev","Number of eigenvalues to compute in each subsolve (only for spectrum slicing)","PEPSTOARSetDimensions",20,&i,&f1);
434:   PetscOptionsInt("-pep_stoar_ncv","Number of basis vectors in each subsolve (only for spectrum slicing)","PEPSTOARSetDimensions",40,&j,&f2);
435:   PetscOptionsInt("-pep_stoar_mpd","Maximum dimension of projected problem in each subsolve (only for spectrum slicing)","PEPSTOARSetDimensions",40,&k,&f3);
436:   if (f1 || f2 || f3) { PEPSTOARSetDimensions(pep,i,j,k); }

438:   PetscOptionsTail();
439:   return(0);
440: }

442: static PetscErrorCode PEPSTOARSetLocking_STOAR(PEP pep,PetscBool lock)
443: {
444:   PEP_TOAR *ctx = (PEP_TOAR*)pep->data;

447:   ctx->lock = lock;
448:   return(0);
449: }

451: /*@
452:    PEPSTOARSetLocking - Choose between locking and non-locking variants of
453:    the STOAR method.

455:    Logically Collective on PEP

457:    Input Parameters:
458: +  pep  - the eigenproblem solver context
459: -  lock - true if the locking variant must be selected

461:    Options Database Key:
462: .  -pep_stoar_locking - Sets the locking flag

464:    Notes:
465:    The default is to lock converged eigenpairs when the method restarts.
466:    This behaviour can be changed so that all directions are kept in the
467:    working subspace even if already converged to working accuracy (the
468:    non-locking variant).

470:    Level: advanced

472: .seealso: PEPSTOARGetLocking()
473: @*/
474: PetscErrorCode PEPSTOARSetLocking(PEP pep,PetscBool lock)
475: {

481:   PetscTryMethod(pep,"PEPSTOARSetLocking_C",(PEP,PetscBool),(pep,lock));
482:   return(0);
483: }

485: static PetscErrorCode PEPSTOARGetLocking_STOAR(PEP pep,PetscBool *lock)
486: {
487:   PEP_TOAR *ctx = (PEP_TOAR*)pep->data;

490:   *lock = ctx->lock;
491:   return(0);
492: }

494: /*@
495:    PEPSTOARGetLocking - Gets the locking flag used in the STOAR method.

497:    Not Collective

499:    Input Parameter:
500: .  pep - the eigenproblem solver context

502:    Output Parameter:
503: .  lock - the locking flag

505:    Level: advanced

507: .seealso: PEPSTOARSetLocking()
508: @*/
509: PetscErrorCode PEPSTOARGetLocking(PEP pep,PetscBool *lock)
510: {

516:   PetscUseMethod(pep,"PEPSTOARGetLocking_C",(PEP,PetscBool*),(pep,lock));
517:   return(0);
518: }

520: static PetscErrorCode PEPSTOARGetInertias_STOAR(PEP pep,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
521: {
523:   PetscInt       i,numsh;
524:   PEP_TOAR       *ctx = (PEP_TOAR*)pep->data;
525:   PEP_SR         sr = ctx->sr;

528:   if (!pep->state) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONGSTATE,"Must call PEPSetUp() first");
529:   if (!ctx->sr) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see PEPSetInterval()");
530:   switch (pep->state) {
531:   case PEP_STATE_INITIAL:
532:     break;
533:   case PEP_STATE_SETUP:
534:     if (n) *n = 2;
535:     if (shifts) {
536:       PetscMalloc1(2,shifts);
537:       (*shifts)[0] = pep->inta;
538:       (*shifts)[1] = pep->intb;
539:     }
540:     if (inertias) {
541:       PetscMalloc1(2,inertias);
542:       (*inertias)[0] = (sr->dir==1)?sr->inertia0:sr->inertia1;
543:       (*inertias)[1] = (sr->dir==1)?sr->inertia1:sr->inertia0;
544:     }
545:     break;
546:   case PEP_STATE_SOLVED:
547:   case PEP_STATE_EIGENVECTORS:
548:     numsh = ctx->nshifts;
549:     if (n) *n = numsh;
550:     if (shifts) {
551:       PetscMalloc1(numsh,shifts);
552:       for (i=0;i<numsh;i++) (*shifts)[i] = ctx->shifts[i];
553:     }
554:     if (inertias) {
555:       PetscMalloc1(numsh,inertias);
556:       for (i=0;i<numsh;i++) (*inertias)[i] = ctx->inertias[i];
557:     }
558:     break;
559:   }
560:   return(0);
561: }

563: /*@C
564:    PEPSTOARGetInertias - Gets the values of the shifts and their
565:    corresponding inertias in case of doing spectrum slicing for a
566:    computational interval.

568:    Not Collective

570:    Input Parameter:
571: .  pep - the eigenproblem solver context

573:    Output Parameters:
574: +  n        - number of shifts, including the endpoints of the interval
575: .  shifts   - the values of the shifts used internally in the solver
576: -  inertias - the values of the inertia in each shift

578:    Notes:
579:    If called after PEPSolve(), all shifts used internally by the solver are
580:    returned (including both endpoints and any intermediate ones). If called
581:    before PEPSolve() and after PEPSetUp() then only the information of the
582:    endpoints of subintervals is available.

584:    This function is only available for spectrum slicing runs.

586:    The returned arrays should be freed by the user. Can pass NULL in any of
587:    the two arrays if not required.

589:    Fortran Notes:
590:    The calling sequence from Fortran is
591: .vb
592:    PEPSTOARGetInertias(pep,n,shifts,inertias,ierr)
593:    integer n
594:    double precision shifts(*)
595:    integer inertias(*)
596: .ve
597:    The arrays should be at least of length n. The value of n can be determined
598:    by an initial call
599: .vb
600:    PEPSTOARGetInertias(pep,n,PETSC_NULL_REAL,PETSC_NULL_INTEGER,ierr)
601: .ve

603:    Level: advanced

605: .seealso: PEPSetInterval()
606: @*/
607: PetscErrorCode PEPSTOARGetInertias(PEP pep,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
608: {

614:   PetscUseMethod(pep,"PEPSTOARGetInertias_C",(PEP,PetscInt*,PetscReal**,PetscInt**),(pep,n,shifts,inertias));
615:   return(0);
616: }

618: static PetscErrorCode PEPSTOARSetDetectZeros_STOAR(PEP pep,PetscBool detect)
619: {
620:   PEP_TOAR *ctx = (PEP_TOAR*)pep->data;

623:   ctx->detect = detect;
624:   pep->state  = PEP_STATE_INITIAL;
625:   return(0);
626: }

628: /*@
629:    PEPSTOARSetDetectZeros - Sets a flag to enforce detection of
630:    zeros during the factorizations throughout the spectrum slicing computation.

632:    Logically Collective on PEP

634:    Input Parameters:
635: +  pep    - the eigenproblem solver context
636: -  detect - check for zeros

638:    Options Database Key:
639: .  -pep_stoar_detect_zeros - Check for zeros; this takes an optional
640:    bool value (0/1/no/yes/true/false)

642:    Notes:
643:    A zero in the factorization indicates that a shift coincides with an eigenvalue.

645:    This flag is turned off by default, and may be necessary in some cases.
646:    This feature currently requires an external package for factorizations
647:    with support for zero detection, e.g. MUMPS.

649:    Level: advanced

651: .seealso: PEPSTOARSetPartitions(), PEPSetInterval()
652: @*/
653: PetscErrorCode PEPSTOARSetDetectZeros(PEP pep,PetscBool detect)
654: {

660:   PetscTryMethod(pep,"PEPSTOARSetDetectZeros_C",(PEP,PetscBool),(pep,detect));
661:   return(0);
662: }

664: static PetscErrorCode PEPSTOARGetDetectZeros_STOAR(PEP pep,PetscBool *detect)
665: {
666:   PEP_TOAR *ctx = (PEP_TOAR*)pep->data;

669:   *detect = ctx->detect;
670:   return(0);
671: }

673: /*@
674:    PEPSTOARGetDetectZeros - Gets the flag that enforces zero detection
675:    in spectrum slicing.

677:    Not Collective

679:    Input Parameter:
680: .  pep - the eigenproblem solver context

682:    Output Parameter:
683: .  detect - whether zeros detection is enforced during factorizations

685:    Level: advanced

687: .seealso: PEPSTOARSetDetectZeros()
688: @*/
689: PetscErrorCode PEPSTOARGetDetectZeros(PEP pep,PetscBool *detect)
690: {

696:   PetscUseMethod(pep,"PEPSTOARGetDetectZeros_C",(PEP,PetscBool*),(pep,detect));
697:   return(0);
698: }

700: static PetscErrorCode PEPSTOARSetDimensions_STOAR(PEP pep,PetscInt nev,PetscInt ncv,PetscInt mpd)
701: {
702:   PEP_TOAR *ctx = (PEP_TOAR*)pep->data;

705:   if (nev<1) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of nev. Must be > 0");
706:   ctx->nev = nev;
707:   if (ncv == PETSC_DECIDE || ncv == PETSC_DEFAULT) {
708:     ctx->ncv = 0;
709:   } else {
710:     if (ncv<1) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of ncv. Must be > 0");
711:     ctx->ncv = ncv;
712:   }
713:   if (mpd == PETSC_DECIDE || mpd == PETSC_DEFAULT) {
714:     ctx->mpd = 0;
715:   } else {
716:     if (mpd<1) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of mpd. Must be > 0");
717:     ctx->mpd = mpd;
718:   }
719:   pep->state = PEP_STATE_INITIAL;
720:   return(0);
721: }

723: /*@
724:    PEPSTOARSetDimensions - Sets the dimensions used for each subsolve
725:    step in case of doing spectrum slicing for a computational interval.
726:    The meaning of the parameters is the same as in PEPSetDimensions().

728:    Logically Collective on PEP

730:    Input Parameters:
731: +  pep - the eigenproblem solver context
732: .  nev - number of eigenvalues to compute
733: .  ncv - the maximum dimension of the subspace to be used by the subsolve
734: -  mpd - the maximum dimension allowed for the projected problem

736:    Options Database Key:
737: +  -eps_stoar_nev <nev> - Sets the number of eigenvalues
738: .  -eps_stoar_ncv <ncv> - Sets the dimension of the subspace
739: -  -eps_stoar_mpd <mpd> - Sets the maximum projected dimension

741:    Level: advanced

743: .seealso: PEPSTOARGetDimensions(), PEPSetDimensions(), PEPSetInterval()
744: @*/
745: PetscErrorCode PEPSTOARSetDimensions(PEP pep,PetscInt nev,PetscInt ncv,PetscInt mpd)
746: {

754:   PetscTryMethod(pep,"PEPSTOARSetDimensions_C",(PEP,PetscInt,PetscInt,PetscInt),(pep,nev,ncv,mpd));
755:   return(0);
756: }

758: static PetscErrorCode PEPSTOARGetDimensions_STOAR(PEP pep,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
759: {
760:   PEP_TOAR *ctx = (PEP_TOAR*)pep->data;

763:   if (nev) *nev = ctx->nev;
764:   if (ncv) *ncv = ctx->ncv;
765:   if (mpd) *mpd = ctx->mpd;
766:   return(0);
767: }

769: /*@
770:    PEPSTOARGetDimensions - Gets the dimensions used for each subsolve
771:    step in case of doing spectrum slicing for a computational interval.

773:    Not Collective

775:    Input Parameter:
776: .  pep - the eigenproblem solver context

778:    Output Parameters:
779: +  nev - number of eigenvalues to compute
780: .  ncv - the maximum dimension of the subspace to be used by the subsolve
781: -  mpd - the maximum dimension allowed for the projected problem

783:    Level: advanced

785: .seealso: PEPSTOARSetDimensions()
786: @*/
787: PetscErrorCode PEPSTOARGetDimensions(PEP pep,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
788: {

793:   PetscUseMethod(pep,"PEPSTOARGetDimensions_C",(PEP,PetscInt*,PetscInt*,PetscInt*),(pep,nev,ncv,mpd));
794:   return(0);
795: }

797: PetscErrorCode PEPView_STOAR(PEP pep,PetscViewer viewer)
798: {
800:   PEP_TOAR      *ctx = (PEP_TOAR*)pep->data;
801:   PetscBool      isascii;

804:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
805:   if (isascii) {
806:     PetscViewerASCIIPrintf(viewer,"  using the %slocking variant\n",ctx->lock?"":"non-");
807:   }
808:   return(0);
809: }

811: PetscErrorCode PEPReset_STOAR(PEP pep)
812: {

816:   if (pep->which==PEP_ALL) {
817:     PEPReset_STOAR_QSlice(pep);
818:   }
819:   return(0);
820: }

822: PetscErrorCode PEPDestroy_STOAR(PEP pep)
823: {
825:   PEP_TOAR       *ctx = (PEP_TOAR*)pep->data;

828:   BVDestroy(&ctx->V);
829:   PetscFree(pep->data);
830:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLocking_C",NULL);
831:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLocking_C",NULL);
832:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDetectZeros_C",NULL);
833:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDetectZeros_C",NULL);
834:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetInertias_C",NULL);
835:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDimensions_C",NULL);
836:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDimensions_C",NULL);
837:   return(0);
838: }

840: PETSC_EXTERN PetscErrorCode PEPCreate_STOAR(PEP pep)
841: {
843:   PEP_TOAR      *ctx;

846:   PetscNewLog(pep,&ctx);
847:   pep->data = (void*)ctx;
848:   ctx->lock = PETSC_TRUE;

850:   pep->ops->setup          = PEPSetUp_STOAR;
851:   pep->ops->setfromoptions = PEPSetFromOptions_STOAR;
852:   pep->ops->destroy        = PEPDestroy_STOAR;
853:   pep->ops->view           = PEPView_STOAR;
854:   pep->ops->backtransform  = PEPBackTransform_Default;
855:   pep->ops->computevectors = PEPComputeVectors_Default;
856:   pep->ops->extractvectors = PEPExtractVectors_TOAR;
857:   pep->ops->setdefaultst   = PEPSetDefaultST_Transform;
858:   pep->ops->reset          = PEPReset_STOAR;

860:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLocking_C",PEPSTOARSetLocking_STOAR);
861:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLocking_C",PEPSTOARGetLocking_STOAR);
862:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDetectZeros_C",PEPSTOARSetDetectZeros_STOAR);
863:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDetectZeros_C",PEPSTOARGetDetectZeros_STOAR);
864:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetInertias_C",PEPSTOARGetInertias_STOAR);
865:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDimensions_C",PEPSTOARGetDimensions_STOAR);
866:   PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDimensions_C",PEPSTOARSetDimensions_STOAR);
867:   return(0);
868: }