Actual source code: qarnoldi.c
slepc-3.14.1 2020-12-08
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2020, 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>
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 flg;
41: PEPCheckQuadratic(pep);
42: PEPCheckShiftSinvert(pep);
43: PEPSetDimensions_Default(pep,pep->nev,&pep->ncv,&pep->mpd);
44: if (!ctx->lock && pep->mpd<pep->ncv) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Should not use mpd parameter in non-locking variant");
45: if (pep->max_it==PETSC_DEFAULT) pep->max_it = PetscMax(100,4*pep->n/pep->ncv);
46: if (!pep->which) { PEPSetWhichEigenpairs_Default(pep); }
47: if (pep->which==PEP_ALL) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"This solver does not support computing all eigenvalues");
49: STGetTransform(pep->st,&flg);
50: if (!flg) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Solver requires the ST transformation flag set, see STSetTransform()");
52: /* set default extraction */
53: if (!pep->extract) pep->extract = PEP_EXTRACT_NONE;
54: PEPCheckUnsupported(pep,PEP_FEATURE_NONMONOMIAL | PEP_FEATURE_EXTRACT);
56: if (!ctx->keep) ctx->keep = 0.5;
58: PEPAllocateSolution(pep,0);
59: PEPSetWorkVecs(pep,4);
61: DSSetType(pep->ds,DSNHEP);
62: DSSetExtraRow(pep->ds,PETSC_TRUE);
63: DSAllocate(pep->ds,pep->ncv+1);
65: return(0);
66: }
68: PetscErrorCode PEPExtractVectors_QArnoldi(PEP pep)
69: {
71: PetscInt i,k=pep->nconv,ldds;
72: PetscScalar *X,*pX0;
73: Mat X0;
76: if (pep->nconv==0) return(0);
77: DSGetLeadingDimension(pep->ds,&ldds);
78: DSVectors(pep->ds,DS_MAT_X,NULL,NULL);
79: DSGetArray(pep->ds,DS_MAT_X,&X);
81: /* update vectors V = V*X */
82: MatCreateSeqDense(PETSC_COMM_SELF,k,k,NULL,&X0);
83: MatDenseGetArray(X0,&pX0);
84: for (i=0;i<k;i++) {
85: PetscArraycpy(pX0+i*k,X+i*ldds,k);
86: }
87: MatDenseRestoreArray(X0,&pX0);
88: BVMultInPlace(pep->V,X0,0,k);
89: MatDestroy(&X0);
90: DSRestoreArray(pep->ds,DS_MAT_X,&X);
91: return(0);
92: }
94: /*
95: Compute a step of Classical Gram-Schmidt orthogonalization
96: */
97: 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)
98: {
100: PetscBLASInt ione = 1,j_1 = j+1;
101: PetscReal x,y;
102: PetscScalar dot,one = 1.0,zero = 0.0;
105: /* compute norm of v and w */
106: if (onorm) {
107: VecNorm(v,NORM_2,&x);
108: VecNorm(w,NORM_2,&y);
109: *onorm = PetscSqrtReal(x*x+y*y);
110: }
112: /* orthogonalize: compute h */
113: BVDotVec(V,v,h);
114: BVDotVec(V,w,work);
115: if (j>0)
116: PetscStackCallBLAS("BLASgemv",BLASgemv_("C",&j_1,&j,&one,H,&ldh,work,&ione,&one,h,&ione));
117: VecDot(w,t,&dot);
118: h[j] += dot;
120: /* orthogonalize: update v and w */
121: BVMultVec(V,-1.0,1.0,v,h);
122: if (j>0) {
123: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&j_1,&j,&one,H,&ldh,h,&ione,&zero,work,&ione));
124: BVMultVec(V,-1.0,1.0,w,work);
125: }
126: VecAXPY(w,-h[j],t);
128: /* compute norm of v and w */
129: if (norm) {
130: VecNorm(v,NORM_2,&x);
131: VecNorm(w,NORM_2,&y);
132: *norm = PetscSqrtReal(x*x+y*y);
133: }
134: return(0);
135: }
137: /*
138: Compute a run of Q-Arnoldi iterations
139: */
140: static PetscErrorCode PEPQArnoldi(PEP pep,PetscScalar *H,PetscInt ldh,PetscInt k,PetscInt *M,Vec v,Vec w,PetscReal *beta,PetscBool *breakdown,PetscScalar *work)
141: {
142: PetscErrorCode ierr;
143: PetscInt i,j,l,m = *M;
144: Vec t = pep->work[2],u = pep->work[3];
145: BVOrthogRefineType refinement;
146: PetscReal norm=0.0,onorm,eta;
147: PetscScalar *c = work + m;
150: BVGetOrthogonalization(pep->V,NULL,&refinement,&eta,NULL);
151: BVInsertVec(pep->V,k,v);
152: for (j=k;j<m;j++) {
153: /* apply operator */
154: VecCopy(w,t);
155: if (pep->Dr) {
156: VecPointwiseMult(v,v,pep->Dr);
157: }
158: STMatMult(pep->st,0,v,u);
159: VecCopy(t,v);
160: if (pep->Dr) {
161: VecPointwiseMult(t,t,pep->Dr);
162: }
163: STMatMult(pep->st,1,t,w);
164: VecAXPY(u,pep->sfactor,w);
165: STMatSolve(pep->st,u,w);
166: VecScale(w,-1.0/(pep->sfactor*pep->sfactor));
167: if (pep->Dr) {
168: VecPointwiseDivide(w,w,pep->Dr);
169: }
170: VecCopy(v,t);
171: BVSetActiveColumns(pep->V,0,j+1);
173: /* orthogonalize */
174: switch (refinement) {
175: case BV_ORTHOG_REFINE_NEVER:
176: PEPQArnoldiCGS(pep,H,ldh,H+ldh*j,j,pep->V,t,v,w,NULL,&norm,work);
177: *breakdown = PETSC_FALSE;
178: break;
179: case BV_ORTHOG_REFINE_ALWAYS:
180: PEPQArnoldiCGS(pep,H,ldh,H+ldh*j,j,pep->V,t,v,w,NULL,NULL,work);
181: PEPQArnoldiCGS(pep,H,ldh,c,j,pep->V,t,v,w,&onorm,&norm,work);
182: for (i=0;i<=j;i++) H[ldh*j+i] += c[i];
183: if (norm < eta * onorm) *breakdown = PETSC_TRUE;
184: else *breakdown = PETSC_FALSE;
185: break;
186: case BV_ORTHOG_REFINE_IFNEEDED:
187: PEPQArnoldiCGS(pep,H,ldh,H+ldh*j,j,pep->V,t,v,w,&onorm,&norm,work);
188: /* ||q|| < eta ||h|| */
189: l = 1;
190: while (l<3 && norm < eta * onorm) {
191: l++;
192: onorm = norm;
193: PEPQArnoldiCGS(pep,H,ldh,c,j,pep->V,t,v,w,NULL,&norm,work);
194: for (i=0;i<=j;i++) H[ldh*j+i] += c[i];
195: }
196: if (norm < eta * onorm) *breakdown = PETSC_TRUE;
197: else *breakdown = PETSC_FALSE;
198: break;
199: }
200: VecScale(v,1.0/norm);
201: VecScale(w,1.0/norm);
203: H[j+1+ldh*j] = norm;
204: if (j<m-1) {
205: BVInsertVec(pep->V,j+1,v);
206: }
207: }
208: *beta = norm;
209: return(0);
210: }
212: PetscErrorCode PEPSolve_QArnoldi(PEP pep)
213: {
215: PEP_QARNOLDI *ctx = (PEP_QARNOLDI*)pep->data;
216: PetscInt j,k,l,lwork,nv,ld,newn,nconv;
217: Vec v=pep->work[0],w=pep->work[1];
218: Mat Q;
219: PetscScalar *S,*work;
220: PetscReal beta=0.0,norm,x,y;
221: PetscBool breakdown=PETSC_FALSE,sinv;
224: DSGetLeadingDimension(pep->ds,&ld);
225: lwork = 7*pep->ncv;
226: PetscMalloc1(lwork,&work);
227: PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv);
228: RGPushScale(pep->rg,sinv?pep->sfactor:1.0/pep->sfactor);
229: STScaleShift(pep->st,sinv?pep->sfactor:1.0/pep->sfactor);
231: /* Get the starting Arnoldi vector */
232: for (j=0;j<2;j++) {
233: if (j>=pep->nini) {
234: BVSetRandomColumn(pep->V,j);
235: }
236: }
237: BVCopyVec(pep->V,0,v);
238: BVCopyVec(pep->V,1,w);
239: VecNorm(v,NORM_2,&x);
240: VecNorm(w,NORM_2,&y);
241: norm = PetscSqrtReal(x*x+y*y);
242: VecScale(v,1.0/norm);
243: VecScale(w,1.0/norm);
245: /* clean projected matrix (including the extra-arrow) */
246: DSGetArray(pep->ds,DS_MAT_A,&S);
247: PetscArrayzero(S,ld*ld);
248: DSRestoreArray(pep->ds,DS_MAT_A,&S);
250: /* Restart loop */
251: l = 0;
252: while (pep->reason == PEP_CONVERGED_ITERATING) {
253: pep->its++;
255: /* Compute an nv-step Arnoldi factorization */
256: nv = PetscMin(pep->nconv+pep->mpd,pep->ncv);
257: DSGetArray(pep->ds,DS_MAT_A,&S);
258: PEPQArnoldi(pep,S,ld,pep->nconv+l,&nv,v,w,&beta,&breakdown,work);
259: DSRestoreArray(pep->ds,DS_MAT_A,&S);
260: DSSetDimensions(pep->ds,nv,0,pep->nconv,pep->nconv+l);
261: if (l==0) {
262: DSSetState(pep->ds,DS_STATE_INTERMEDIATE);
263: } else {
264: DSSetState(pep->ds,DS_STATE_RAW);
265: }
266: BVSetActiveColumns(pep->V,pep->nconv,nv);
268: /* Solve projected problem */
269: DSSolve(pep->ds,pep->eigr,pep->eigi);
270: DSSort(pep->ds,pep->eigr,pep->eigi,NULL,NULL,NULL);
271: DSUpdateExtraRow(pep->ds);
272: DSSynchronize(pep->ds,pep->eigr,pep->eigi);
274: /* Check convergence */
275: PEPKrylovConvergence(pep,PETSC_FALSE,pep->nconv,nv-pep->nconv,beta,&k);
276: (*pep->stopping)(pep,pep->its,pep->max_it,k,pep->nev,&pep->reason,pep->stoppingctx);
277: nconv = k;
279: /* Update l */
280: if (pep->reason != PEP_CONVERGED_ITERATING || breakdown) l = 0;
281: else l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
282: if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
284: if (pep->reason == PEP_CONVERGED_ITERATING) {
285: if (breakdown) {
286: /* Stop if breakdown */
287: PetscInfo2(pep,"Breakdown Quadratic Arnoldi method (it=%D norm=%g)\n",pep->its,(double)beta);
288: pep->reason = PEP_DIVERGED_BREAKDOWN;
289: } else {
290: /* Prepare the Rayleigh quotient for restart */
291: DSTruncate(pep->ds,k+l);
292: DSGetDimensions(pep->ds,&newn,NULL,NULL,NULL,NULL);
293: l = newn-k;
294: }
295: }
296: /* Update the corresponding vectors V(:,idx) = V*Q(:,idx) */
297: DSGetMat(pep->ds,DS_MAT_Q,&Q);
298: BVMultInPlace(pep->V,Q,pep->nconv,k+l);
299: MatDestroy(&Q);
301: pep->nconv = k;
302: PEPMonitor(pep,pep->its,nconv,pep->eigr,pep->eigi,pep->errest,nv);
303: }
304: BVSetActiveColumns(pep->V,0,pep->nconv);
305: for (j=0;j<pep->nconv;j++) {
306: pep->eigr[j] *= pep->sfactor;
307: pep->eigi[j] *= pep->sfactor;
308: }
310: STScaleShift(pep->st,sinv?1.0/pep->sfactor:pep->sfactor);
311: RGPopScale(pep->rg);
313: /* truncate Schur decomposition and change the state to raw so that
314: DSVectors() computes eigenvectors from scratch */
315: DSSetDimensions(pep->ds,pep->nconv,0,0,0);
316: DSSetState(pep->ds,DS_STATE_RAW);
317: PetscFree(work);
318: return(0);
319: }
321: static PetscErrorCode PEPQArnoldiSetRestart_QArnoldi(PEP pep,PetscReal keep)
322: {
323: PEP_QARNOLDI *ctx = (PEP_QARNOLDI*)pep->data;
326: if (keep==PETSC_DEFAULT) ctx->keep = 0.5;
327: else {
328: 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]");
329: ctx->keep = keep;
330: }
331: return(0);
332: }
334: /*@
335: PEPQArnoldiSetRestart - Sets the restart parameter for the Q-Arnoldi
336: method, in particular the proportion of basis vectors that must be kept
337: after restart.
339: Logically Collective on pep
341: Input Parameters:
342: + pep - the eigenproblem solver context
343: - keep - the number of vectors to be kept at restart
345: Options Database Key:
346: . -pep_qarnoldi_restart - Sets the restart parameter
348: Notes:
349: Allowed values are in the range [0.1,0.9]. The default is 0.5.
351: Level: advanced
353: .seealso: PEPQArnoldiGetRestart()
354: @*/
355: PetscErrorCode PEPQArnoldiSetRestart(PEP pep,PetscReal keep)
356: {
362: PetscTryMethod(pep,"PEPQArnoldiSetRestart_C",(PEP,PetscReal),(pep,keep));
363: return(0);
364: }
366: static PetscErrorCode PEPQArnoldiGetRestart_QArnoldi(PEP pep,PetscReal *keep)
367: {
368: PEP_QARNOLDI *ctx = (PEP_QARNOLDI*)pep->data;
371: *keep = ctx->keep;
372: return(0);
373: }
375: /*@
376: PEPQArnoldiGetRestart - Gets the restart parameter used in the Q-Arnoldi method.
378: Not Collective
380: Input Parameter:
381: . pep - the eigenproblem solver context
383: Output Parameter:
384: . keep - the restart parameter
386: Level: advanced
388: .seealso: PEPQArnoldiSetRestart()
389: @*/
390: PetscErrorCode PEPQArnoldiGetRestart(PEP pep,PetscReal *keep)
391: {
397: PetscUseMethod(pep,"PEPQArnoldiGetRestart_C",(PEP,PetscReal*),(pep,keep));
398: return(0);
399: }
401: static PetscErrorCode PEPQArnoldiSetLocking_QArnoldi(PEP pep,PetscBool lock)
402: {
403: PEP_QARNOLDI *ctx = (PEP_QARNOLDI*)pep->data;
406: ctx->lock = lock;
407: return(0);
408: }
410: /*@
411: PEPQArnoldiSetLocking - Choose between locking and non-locking variants of
412: the Q-Arnoldi method.
414: Logically Collective on pep
416: Input Parameters:
417: + pep - the eigenproblem solver context
418: - lock - true if the locking variant must be selected
420: Options Database Key:
421: . -pep_qarnoldi_locking - Sets the locking flag
423: Notes:
424: The default is to keep all directions in the working subspace even if
425: already converged to working accuracy (the non-locking variant).
426: This behaviour can be changed so that converged eigenpairs are locked
427: when the method restarts.
429: Note that the default behaviour is the opposite to Krylov solvers in EPS.
431: Level: advanced
433: .seealso: PEPQArnoldiGetLocking()
434: @*/
435: PetscErrorCode PEPQArnoldiSetLocking(PEP pep,PetscBool lock)
436: {
442: PetscTryMethod(pep,"PEPQArnoldiSetLocking_C",(PEP,PetscBool),(pep,lock));
443: return(0);
444: }
446: static PetscErrorCode PEPQArnoldiGetLocking_QArnoldi(PEP pep,PetscBool *lock)
447: {
448: PEP_QARNOLDI *ctx = (PEP_QARNOLDI*)pep->data;
451: *lock = ctx->lock;
452: return(0);
453: }
455: /*@
456: PEPQArnoldiGetLocking - Gets the locking flag used in the Q-Arnoldi method.
458: Not Collective
460: Input Parameter:
461: . pep - the eigenproblem solver context
463: Output Parameter:
464: . lock - the locking flag
466: Level: advanced
468: .seealso: PEPQArnoldiSetLocking()
469: @*/
470: PetscErrorCode PEPQArnoldiGetLocking(PEP pep,PetscBool *lock)
471: {
477: PetscUseMethod(pep,"PEPQArnoldiGetLocking_C",(PEP,PetscBool*),(pep,lock));
478: return(0);
479: }
481: PetscErrorCode PEPSetFromOptions_QArnoldi(PetscOptionItems *PetscOptionsObject,PEP pep)
482: {
484: PetscBool flg,lock;
485: PetscReal keep;
488: PetscOptionsHead(PetscOptionsObject,"PEP Q-Arnoldi Options");
490: PetscOptionsReal("-pep_qarnoldi_restart","Proportion of vectors kept after restart","PEPQArnoldiSetRestart",0.5,&keep,&flg);
491: if (flg) { PEPQArnoldiSetRestart(pep,keep); }
493: PetscOptionsBool("-pep_qarnoldi_locking","Choose between locking and non-locking variants","PEPQArnoldiSetLocking",PETSC_FALSE,&lock,&flg);
494: if (flg) { PEPQArnoldiSetLocking(pep,lock); }
496: PetscOptionsTail();
497: return(0);
498: }
500: PetscErrorCode PEPView_QArnoldi(PEP pep,PetscViewer viewer)
501: {
503: PEP_QARNOLDI *ctx = (PEP_QARNOLDI*)pep->data;
504: PetscBool isascii;
507: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
508: if (isascii) {
509: PetscViewerASCIIPrintf(viewer," %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep));
510: PetscViewerASCIIPrintf(viewer," using the %slocking variant\n",ctx->lock?"":"non-");
511: }
512: return(0);
513: }
515: PetscErrorCode PEPDestroy_QArnoldi(PEP pep)
516: {
520: PetscFree(pep->data);
521: PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiSetRestart_C",NULL);
522: PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiGetRestart_C",NULL);
523: PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiSetLocking_C",NULL);
524: PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiGetLocking_C",NULL);
525: return(0);
526: }
528: SLEPC_EXTERN PetscErrorCode PEPCreate_QArnoldi(PEP pep)
529: {
530: PEP_QARNOLDI *ctx;
534: PetscNewLog(pep,&ctx);
535: pep->data = (void*)ctx;
537: pep->lineariz = PETSC_TRUE;
538: ctx->lock = PETSC_TRUE;
540: pep->ops->solve = PEPSolve_QArnoldi;
541: pep->ops->setup = PEPSetUp_QArnoldi;
542: pep->ops->setfromoptions = PEPSetFromOptions_QArnoldi;
543: pep->ops->destroy = PEPDestroy_QArnoldi;
544: pep->ops->view = PEPView_QArnoldi;
545: pep->ops->backtransform = PEPBackTransform_Default;
546: pep->ops->computevectors = PEPComputeVectors_Default;
547: pep->ops->extractvectors = PEPExtractVectors_QArnoldi;
548: pep->ops->setdefaultst = PEPSetDefaultST_Transform;
550: PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiSetRestart_C",PEPQArnoldiSetRestart_QArnoldi);
551: PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiGetRestart_C",PEPQArnoldiGetRestart_QArnoldi);
552: PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiSetLocking_C",PEPQArnoldiSetLocking_QArnoldi);
553: PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiGetLocking_C",PEPQArnoldiGetLocking_QArnoldi);
554: return(0);
555: }