Actual source code: qarnoldi.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, 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 quadratic eigensolver: "qarnoldi"

 13:    Method: Q-Arnoldi

 15:    Algorithm:

 17:        Quadratic Arnoldi with Krylov-Schur type restart.

 19:    References:

 21:        [1] K. Meerbergen, "The Quadratic Arnoldi method for the solution
 22:            of the quadratic eigenvalue problem", SIAM J. Matrix Anal.
 23:            Appl. 30(4):1462-1482, 2008.
 24: */

 26: #include <slepc/private/pepimpl.h>    /*I "slepcpep.h" I*/
 27: #include <petscblaslapack.h>

 29: typedef struct {
 30:   PetscReal keep;         /* restart parameter */
 31:   PetscBool lock;         /* locking/non-locking variant */
 32: } PEP_QARNOLDI;

 34: PetscErrorCode PEPSetUp_QArnoldi(PEP pep)
 35: {
 37:   PEP_QARNOLDI   *ctx = (PEP_QARNOLDI*)pep->data;
 38:   PetscBool      shift,sinv,flg;

 41:   pep->lineariz = PETSC_TRUE;
 42:   PEPSetDimensions_Default(pep,pep->nev,&pep->ncv,&pep->mpd);
 43:   if (!ctx->lock && pep->mpd<pep->ncv) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Should not use mpd parameter in non-locking variant");
 44:   if (!pep->max_it) pep->max_it = PetscMax(100,4*pep->n/pep->ncv);

 46:   PetscObjectTypeCompare((PetscObject)pep->st,STSHIFT,&shift);
 47:   PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv);
 48:   if (!shift && !sinv) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Only STSHIFT and STSINVERT spectral transformations can be used");
 49:   if (!pep->which) {
 50:     if (sinv) pep->which = PEP_TARGET_MAGNITUDE;
 51:     else pep->which = PEP_LARGEST_MAGNITUDE;
 52:   }
 53:   if (pep->which==PEP_ALL) SETERRQ(PetscObjectComm((PetscObject)pep),1,"Wrong value of pep->which");

 55:   if (pep->nmat!=3) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Solver only available for quadratic problems");
 56:   if (pep->basis!=PEP_BASIS_MONOMIAL) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Solver not implemented for non-monomial bases");
 57:   STGetTransform(pep->st,&flg);
 58:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Solver requires the ST transformation flag set, see STSetTransform()");

 60:   /* set default extraction */
 61:   if (!pep->extract) {
 62:     pep->extract = PEP_EXTRACT_NONE;
 63:   }
 64:   if (pep->extract!=PEP_EXTRACT_NONE) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Solver does not support requested extraction");

 66:   if (!ctx->keep) ctx->keep = 0.5;

 68:   PEPAllocateSolution(pep,0);
 69:   PEPSetWorkVecs(pep,4);

 71:   DSSetType(pep->ds,DSNHEP);
 72:   DSSetExtraRow(pep->ds,PETSC_TRUE);
 73:   DSAllocate(pep->ds,pep->ncv+1);

 75:   return(0);
 76: }

 78: PetscErrorCode PEPExtractVectors_QArnoldi(PEP pep)
 79: {
 81:   PetscInt       i,k=pep->nconv,ldds;
 82:   PetscScalar    *X,*pX0;
 83:   Mat            X0;

 86:   if (pep->nconv==0) return(0);
 87:   DSGetLeadingDimension(pep->ds,&ldds);
 88:   DSVectors(pep->ds,DS_MAT_X,NULL,NULL);
 89:   DSGetArray(pep->ds,DS_MAT_X,&X);

 91:   /* update vectors V = V*X */
 92:   MatCreateSeqDense(PETSC_COMM_SELF,k,k,NULL,&X0);
 93:   MatDenseGetArray(X0,&pX0);
 94:   for (i=0;i<k;i++) {
 95:     PetscMemcpy(pX0+i*k,X+i*ldds,k*sizeof(PetscScalar));
 96:   }
 97:   MatDenseRestoreArray(X0,&pX0);
 98:   BVMultInPlace(pep->V,X0,0,k);
 99:   MatDestroy(&X0);
100:   DSRestoreArray(pep->ds,DS_MAT_X,&X);
101:   return(0);
102: }

104: /*
105:   Compute a step of Classical Gram-Schmidt orthogonalization
106: */
107: static PetscErrorCode PEPQArnoldiCGS(PEP pep,PetscScalar *H,PetscBLASInt ldh,PetscScalar *h,PetscBLASInt j,BV V,Vec t,Vec v,Vec w,PetscReal *onorm,PetscReal *norm,PetscScalar *work)
108: {
110:   PetscBLASInt   ione = 1,j_1 = j+1;
111:   PetscReal      x,y;
112:   PetscScalar    dot,one = 1.0,zero = 0.0;

115:   /* compute norm of v and w */
116:   if (onorm) {
117:     VecNorm(v,NORM_2,&x);
118:     VecNorm(w,NORM_2,&y);
119:     *onorm = PetscSqrtReal(x*x+y*y);
120:   }

122:   /* orthogonalize: compute h */
123:   BVDotVec(V,v,h);
124:   BVDotVec(V,w,work);
125:   if (j>0)
126:     PetscStackCallBLAS("BLASgemv",BLASgemv_("C",&j_1,&j,&one,H,&ldh,work,&ione,&one,h,&ione));
127:   VecDot(w,t,&dot);
128:   h[j] += dot;

130:   /* orthogonalize: update v and w */
131:   BVMultVec(V,-1.0,1.0,v,h);
132:   if (j>0) {
133:     PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&j_1,&j,&one,H,&ldh,h,&ione,&zero,work,&ione));
134:     BVMultVec(V,-1.0,1.0,w,work);
135:   }
136:   VecAXPY(w,-h[j],t);

138:   /* compute norm of v and w */
139:   if (norm) {
140:     VecNorm(v,NORM_2,&x);
141:     VecNorm(w,NORM_2,&y);
142:     *norm = PetscSqrtReal(x*x+y*y);
143:   }
144:   return(0);
145: }

147: /*
148:   Compute a run of Q-Arnoldi iterations
149: */
150: static PetscErrorCode PEPQArnoldi(PEP pep,PetscScalar *H,PetscInt ldh,PetscInt k,PetscInt *M,Vec v,Vec w,PetscReal *beta,PetscBool *breakdown,PetscScalar *work)
151: {
152:   PetscErrorCode     ierr;
153:   PetscInt           i,j,l,m = *M;
154:   Vec                t = pep->work[2],u = pep->work[3];
155:   BVOrthogRefineType refinement;
156:   PetscReal          norm=0.0,onorm,eta;
157:   PetscScalar        *c = work + m;

160:   BVGetOrthogonalization(pep->V,NULL,&refinement,&eta,NULL);
161:   BVInsertVec(pep->V,k,v);
162:   for (j=k;j<m;j++) {
163:     /* apply operator */
164:     VecCopy(w,t);
165:     if (pep->Dr) {
166:       VecPointwiseMult(v,v,pep->Dr);
167:     }
168:     STMatMult(pep->st,0,v,u);
169:     VecCopy(t,v);
170:     if (pep->Dr) {
171:       VecPointwiseMult(t,t,pep->Dr);
172:     }
173:     STMatMult(pep->st,1,t,w);
174:     VecAXPY(u,pep->sfactor,w);
175:     STMatSolve(pep->st,u,w);
176:     VecScale(w,-1.0/(pep->sfactor*pep->sfactor));
177:     if (pep->Dr) {
178:       VecPointwiseDivide(w,w,pep->Dr);
179:     }
180:     VecCopy(v,t);
181:     BVSetActiveColumns(pep->V,0,j+1);

183:     /* orthogonalize */
184:     switch (refinement) {
185:       case BV_ORTHOG_REFINE_NEVER:
186:         PEPQArnoldiCGS(pep,H,ldh,H+ldh*j,j,pep->V,t,v,w,NULL,&norm,work);
187:         *breakdown = PETSC_FALSE;
188:         break;
189:       case BV_ORTHOG_REFINE_ALWAYS:
190:         PEPQArnoldiCGS(pep,H,ldh,H+ldh*j,j,pep->V,t,v,w,NULL,NULL,work);
191:         PEPQArnoldiCGS(pep,H,ldh,c,j,pep->V,t,v,w,&onorm,&norm,work);
192:         for (i=0;i<=j;i++) H[ldh*j+i] += c[i];
193:         if (norm < eta * onorm) *breakdown = PETSC_TRUE;
194:         else *breakdown = PETSC_FALSE;
195:         break;
196:       case BV_ORTHOG_REFINE_IFNEEDED:
197:         PEPQArnoldiCGS(pep,H,ldh,H+ldh*j,j,pep->V,t,v,w,&onorm,&norm,work);
198:         /* ||q|| < eta ||h|| */
199:         l = 1;
200:         while (l<3 && norm < eta * onorm) {
201:           l++;
202:           onorm = norm;
203:           PEPQArnoldiCGS(pep,H,ldh,c,j,pep->V,t,v,w,NULL,&norm,work);
204:           for (i=0;i<=j;i++) H[ldh*j+i] += c[i];
205:         }
206:         if (norm < eta * onorm) *breakdown = PETSC_TRUE;
207:         else *breakdown = PETSC_FALSE;
208:         break;
209:     }
210:     VecScale(v,1.0/norm);
211:     VecScale(w,1.0/norm);

213:     H[j+1+ldh*j] = norm;
214:     if (j<m-1) {
215:       BVInsertVec(pep->V,j+1,v);
216:     }
217:   }
218:   *beta = norm;
219:   return(0);
220: }

222: PetscErrorCode PEPSolve_QArnoldi(PEP pep)
223: {
225:   PEP_QARNOLDI   *ctx = (PEP_QARNOLDI*)pep->data;
226:   PetscInt       j,k,l,lwork,nv,ld,newn,nconv;
227:   Vec            v=pep->work[0],w=pep->work[1];
228:   Mat            Q;
229:   PetscScalar    *S,*work;
230:   PetscReal      beta=0.0,norm,x,y;
231:   PetscBool      breakdown=PETSC_FALSE,sinv;

234:   DSGetLeadingDimension(pep->ds,&ld);
235:   lwork = 7*pep->ncv;
236:   PetscMalloc1(lwork,&work);
237:   PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv);
238:   RGPushScale(pep->rg,sinv?pep->sfactor:1.0/pep->sfactor);
239:   STScaleShift(pep->st,sinv?pep->sfactor:1.0/pep->sfactor);

241:   /* Get the starting Arnoldi vector */
242:   for (j=0;j<2;j++) {
243:     if (j>=pep->nini) {
244:       BVSetRandomColumn(pep->V,j);
245:     }
246:   }
247:   BVCopyVec(pep->V,0,v);
248:   BVCopyVec(pep->V,1,w);
249:   VecNorm(v,NORM_2,&x);
250:   VecNorm(w,NORM_2,&y);
251:   norm = PetscSqrtReal(x*x+y*y);
252:   VecScale(v,1.0/norm);
253:   VecScale(w,1.0/norm);

255:   /* clean projected matrix (including the extra-arrow) */
256:   DSGetArray(pep->ds,DS_MAT_A,&S);
257:   PetscMemzero(S,ld*ld*sizeof(PetscScalar));
258:   DSRestoreArray(pep->ds,DS_MAT_A,&S);
259:   
260:    /* Restart loop */
261:   l = 0;
262:   while (pep->reason == PEP_CONVERGED_ITERATING) {
263:     pep->its++;

265:     /* Compute an nv-step Arnoldi factorization */
266:     nv = PetscMin(pep->nconv+pep->mpd,pep->ncv);
267:     DSGetArray(pep->ds,DS_MAT_A,&S);
268:     PEPQArnoldi(pep,S,ld,pep->nconv+l,&nv,v,w,&beta,&breakdown,work);
269:     DSRestoreArray(pep->ds,DS_MAT_A,&S);
270:     DSSetDimensions(pep->ds,nv,0,pep->nconv,pep->nconv+l);
271:     if (l==0) {
272:       DSSetState(pep->ds,DS_STATE_INTERMEDIATE);
273:     } else {
274:       DSSetState(pep->ds,DS_STATE_RAW);
275:     }
276:     BVSetActiveColumns(pep->V,pep->nconv,nv);

278:     /* Solve projected problem */
279:     DSSolve(pep->ds,pep->eigr,pep->eigi);
280:     DSSort(pep->ds,pep->eigr,pep->eigi,NULL,NULL,NULL);
281:     DSUpdateExtraRow(pep->ds);
282:     DSSynchronize(pep->ds,pep->eigr,pep->eigi);

284:     /* Check convergence */
285:     PEPKrylovConvergence(pep,PETSC_FALSE,pep->nconv,nv-pep->nconv,beta,&k);
286:     (*pep->stopping)(pep,pep->its,pep->max_it,k,pep->nev,&pep->reason,pep->stoppingctx);
287:     nconv = k;

289:     /* Update l */
290:     if (pep->reason != PEP_CONVERGED_ITERATING || breakdown) l = 0;
291:     else l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
292:     if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */

294:     if (pep->reason == PEP_CONVERGED_ITERATING) {
295:       if (breakdown) {
296:         /* Stop if breakdown */
297:         PetscInfo2(pep,"Breakdown Quadratic Arnoldi method (it=%D norm=%g)\n",pep->its,(double)beta);
298:         pep->reason = PEP_DIVERGED_BREAKDOWN;
299:       } else {
300:         /* Prepare the Rayleigh quotient for restart */
301:         DSTruncate(pep->ds,k+l);
302:         DSGetDimensions(pep->ds,&newn,NULL,NULL,NULL,NULL);
303:         l = newn-k;
304:       }
305:     }
306:     /* Update the corresponding vectors V(:,idx) = V*Q(:,idx) */
307:     DSGetMat(pep->ds,DS_MAT_Q,&Q);
308:     BVMultInPlace(pep->V,Q,pep->nconv,k+l);
309:     MatDestroy(&Q);

311:     pep->nconv = k;
312:     PEPMonitor(pep,pep->its,nconv,pep->eigr,pep->eigi,pep->errest,nv);
313:   }
314:   BVSetActiveColumns(pep->V,0,pep->nconv);
315:   for (j=0;j<pep->nconv;j++) {
316:     pep->eigr[j] *= pep->sfactor;
317:     pep->eigi[j] *= pep->sfactor;
318:   }

320:   STScaleShift(pep->st,sinv?1.0/pep->sfactor:pep->sfactor);
321:   RGPopScale(pep->rg);

323:   /* truncate Schur decomposition and change the state to raw so that
324:      DSVectors() computes eigenvectors from scratch */
325:   DSSetDimensions(pep->ds,pep->nconv,0,0,0);
326:   DSSetState(pep->ds,DS_STATE_RAW);
327:   PetscFree(work);
328:   return(0);
329: }

331: static PetscErrorCode PEPQArnoldiSetRestart_QArnoldi(PEP pep,PetscReal keep)
332: {
333:   PEP_QARNOLDI *ctx = (PEP_QARNOLDI*)pep->data;

336:   if (keep==PETSC_DEFAULT) ctx->keep = 0.5;
337:   else {
338:     if (keep<0.1 || keep>0.9) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"The keep argument must be in the range [0.1,0.9]");
339:     ctx->keep = keep;
340:   }
341:   return(0);
342: }

344: /*@
345:    PEPQArnoldiSetRestart - Sets the restart parameter for the Q-Arnoldi
346:    method, in particular the proportion of basis vectors that must be kept
347:    after restart.

349:    Logically Collective on PEP

351:    Input Parameters:
352: +  pep  - the eigenproblem solver context
353: -  keep - the number of vectors to be kept at restart

355:    Options Database Key:
356: .  -pep_qarnoldi_restart - Sets the restart parameter

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

361:    Level: advanced

363: .seealso: PEPQArnoldiGetRestart()
364: @*/
365: PetscErrorCode PEPQArnoldiSetRestart(PEP pep,PetscReal keep)
366: {

372:   PetscTryMethod(pep,"PEPQArnoldiSetRestart_C",(PEP,PetscReal),(pep,keep));
373:   return(0);
374: }

376: static PetscErrorCode PEPQArnoldiGetRestart_QArnoldi(PEP pep,PetscReal *keep)
377: {
378:   PEP_QARNOLDI *ctx = (PEP_QARNOLDI*)pep->data;

381:   *keep = ctx->keep;
382:   return(0);
383: }

385: /*@
386:    PEPQArnoldiGetRestart - Gets the restart parameter used in the Q-Arnoldi method.

388:    Not Collective

390:    Input Parameter:
391: .  pep - the eigenproblem solver context

393:    Output Parameter:
394: .  keep - the restart parameter

396:    Level: advanced

398: .seealso: PEPQArnoldiSetRestart()
399: @*/
400: PetscErrorCode PEPQArnoldiGetRestart(PEP pep,PetscReal *keep)
401: {

407:   PetscUseMethod(pep,"PEPQArnoldiGetRestart_C",(PEP,PetscReal*),(pep,keep));
408:   return(0);
409: }

411: static PetscErrorCode PEPQArnoldiSetLocking_QArnoldi(PEP pep,PetscBool lock)
412: {
413:   PEP_QARNOLDI *ctx = (PEP_QARNOLDI*)pep->data;

416:   ctx->lock = lock;
417:   return(0);
418: }

420: /*@
421:    PEPQArnoldiSetLocking - Choose between locking and non-locking variants of
422:    the Q-Arnoldi method.

424:    Logically Collective on PEP

426:    Input Parameters:
427: +  pep  - the eigenproblem solver context
428: -  lock - true if the locking variant must be selected

430:    Options Database Key:
431: .  -pep_qarnoldi_locking - Sets the locking flag

433:    Notes:
434:    The default is to keep all directions in the working subspace even if
435:    already converged to working accuracy (the non-locking variant).
436:    This behaviour can be changed so that converged eigenpairs are locked
437:    when the method restarts.

439:    Note that the default behaviour is the opposite to Krylov solvers in EPS.

441:    Level: advanced

443: .seealso: PEPQArnoldiGetLocking()
444: @*/
445: PetscErrorCode PEPQArnoldiSetLocking(PEP pep,PetscBool lock)
446: {

452:   PetscTryMethod(pep,"PEPQArnoldiSetLocking_C",(PEP,PetscBool),(pep,lock));
453:   return(0);
454: }

456: static PetscErrorCode PEPQArnoldiGetLocking_QArnoldi(PEP pep,PetscBool *lock)
457: {
458:   PEP_QARNOLDI *ctx = (PEP_QARNOLDI*)pep->data;

461:   *lock = ctx->lock;
462:   return(0);
463: }

465: /*@
466:    PEPQArnoldiGetLocking - Gets the locking flag used in the Q-Arnoldi method.

468:    Not Collective

470:    Input Parameter:
471: .  pep - the eigenproblem solver context

473:    Output Parameter:
474: .  lock - the locking flag

476:    Level: advanced

478: .seealso: PEPQArnoldiSetLocking()
479: @*/
480: PetscErrorCode PEPQArnoldiGetLocking(PEP pep,PetscBool *lock)
481: {

487:   PetscUseMethod(pep,"PEPQArnoldiGetLocking_C",(PEP,PetscBool*),(pep,lock));
488:   return(0);
489: }

491: PetscErrorCode PEPSetFromOptions_QArnoldi(PetscOptionItems *PetscOptionsObject,PEP pep)
492: {
494:   PetscBool      flg,lock;
495:   PetscReal      keep;

498:   PetscOptionsHead(PetscOptionsObject,"PEP Q-Arnoldi Options");

500:     PetscOptionsReal("-pep_qarnoldi_restart","Proportion of vectors kept after restart","PEPQArnoldiSetRestart",0.5,&keep,&flg);
501:     if (flg) { PEPQArnoldiSetRestart(pep,keep); }

503:     PetscOptionsBool("-pep_qarnoldi_locking","Choose between locking and non-locking variants","PEPQArnoldiSetLocking",PETSC_FALSE,&lock,&flg);
504:     if (flg) { PEPQArnoldiSetLocking(pep,lock); }

506:   PetscOptionsTail();
507:   return(0);
508: }

510: PetscErrorCode PEPView_QArnoldi(PEP pep,PetscViewer viewer)
511: {
513:   PEP_QARNOLDI   *ctx = (PEP_QARNOLDI*)pep->data;
514:   PetscBool      isascii;

517:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
518:   if (isascii) {
519:     PetscViewerASCIIPrintf(viewer,"  %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep));
520:     PetscViewerASCIIPrintf(viewer,"  using the %slocking variant\n",ctx->lock?"":"non-");
521:   }
522:   return(0);
523: }

525: PetscErrorCode PEPDestroy_QArnoldi(PEP pep)
526: {

530:   PetscFree(pep->data);
531:   PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiSetRestart_C",NULL);
532:   PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiGetRestart_C",NULL);
533:   PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiSetLocking_C",NULL);
534:   PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiGetLocking_C",NULL);
535:   return(0);
536: }

538: SLEPC_EXTERN PetscErrorCode PEPCreate_QArnoldi(PEP pep)
539: {
540:   PEP_QARNOLDI   *ctx;

544:   PetscNewLog(pep,&ctx);
545:   pep->data = (void*)ctx;
546:   ctx->lock = PETSC_TRUE;

548:   pep->ops->solve          = PEPSolve_QArnoldi;
549:   pep->ops->setup          = PEPSetUp_QArnoldi;
550:   pep->ops->setfromoptions = PEPSetFromOptions_QArnoldi;
551:   pep->ops->destroy        = PEPDestroy_QArnoldi;
552:   pep->ops->view           = PEPView_QArnoldi;
553:   pep->ops->backtransform  = PEPBackTransform_Default;
554:   pep->ops->computevectors = PEPComputeVectors_Default;
555:   pep->ops->extractvectors = PEPExtractVectors_QArnoldi;
556:   pep->ops->setdefaultst   = PEPSetDefaultST_Transform;

558:   PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiSetRestart_C",PEPQArnoldiSetRestart_QArnoldi);
559:   PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiGetRestart_C",PEPQArnoldiGetRestart_QArnoldi);
560:   PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiSetLocking_C",PEPQArnoldiSetLocking_QArnoldi);
561:   PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiGetLocking_C",PEPQArnoldiGetLocking_QArnoldi);
562:   return(0);
563: }