Actual source code: qarnoldi.c
slepc-3.11.2 2019-07-30
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: }