Actual source code: qarnoldi.c
slepc-3.18.0 2022-10-01
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 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: {
36: PEP_QARNOLDI *ctx = (PEP_QARNOLDI*)pep->data;
37: PetscBool flg;
39: PEPCheckQuadratic(pep);
40: PEPCheckShiftSinvert(pep);
41: PEPSetDimensions_Default(pep,pep->nev,&pep->ncv,&pep->mpd);
43: if (pep->max_it==PETSC_DEFAULT) pep->max_it = PetscMax(100,4*pep->n/pep->ncv);
44: if (!pep->which) PEPSetWhichEigenpairs_Default(pep);
47: STGetTransform(pep->st,&flg);
50: /* set default extraction */
51: if (!pep->extract) pep->extract = PEP_EXTRACT_NONE;
52: PEPCheckUnsupported(pep,PEP_FEATURE_NONMONOMIAL | PEP_FEATURE_EXTRACT);
54: if (!ctx->keep) ctx->keep = 0.5;
56: PEPAllocateSolution(pep,0);
57: PEPSetWorkVecs(pep,4);
59: DSSetType(pep->ds,DSNHEP);
60: DSSetExtraRow(pep->ds,PETSC_TRUE);
61: DSAllocate(pep->ds,pep->ncv+1);
63: return 0;
64: }
66: PetscErrorCode PEPExtractVectors_QArnoldi(PEP pep)
67: {
68: PetscInt k=pep->nconv;
69: Mat X,X0;
71: if (pep->nconv==0) return 0;
72: DSVectors(pep->ds,DS_MAT_X,NULL,NULL);
74: /* update vectors V = V*X */
75: DSGetMat(pep->ds,DS_MAT_X,&X);
76: MatDenseGetSubMatrix(X,0,k,0,k,&X0);
77: BVMultInPlace(pep->V,X0,0,k);
78: MatDenseRestoreSubMatrix(X,&X0);
79: DSRestoreMat(pep->ds,DS_MAT_X,&X);
80: return 0;
81: }
83: /*
84: Compute a step of Classical Gram-Schmidt orthogonalization
85: */
86: 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)
87: {
88: PetscBLASInt ione = 1,j_1 = j+1;
89: PetscReal x,y;
90: PetscScalar dot,one = 1.0,zero = 0.0;
92: /* compute norm of v and w */
93: if (onorm) {
94: VecNorm(v,NORM_2,&x);
95: VecNorm(w,NORM_2,&y);
96: *onorm = SlepcAbs(x,y);
97: }
99: /* orthogonalize: compute h */
100: BVDotVec(V,v,h);
101: BVDotVec(V,w,work);
102: if (j>0) PetscCallBLAS("BLASgemv",BLASgemv_("C",&j_1,&j,&one,H,&ldh,work,&ione,&one,h,&ione));
103: VecDot(w,t,&dot);
104: h[j] += dot;
106: /* orthogonalize: update v and w */
107: BVMultVec(V,-1.0,1.0,v,h);
108: if (j>0) {
109: PetscCallBLAS("BLASgemv",BLASgemv_("N",&j_1,&j,&one,H,&ldh,h,&ione,&zero,work,&ione));
110: BVMultVec(V,-1.0,1.0,w,work);
111: }
112: VecAXPY(w,-h[j],t);
114: /* compute norm of v and w */
115: if (norm) {
116: VecNorm(v,NORM_2,&x);
117: VecNorm(w,NORM_2,&y);
118: *norm = SlepcAbs(x,y);
119: }
120: return 0;
121: }
123: /*
124: Compute a run of Q-Arnoldi iterations
125: */
126: static PetscErrorCode PEPQArnoldi(PEP pep,Mat A,PetscInt k,PetscInt *M,Vec v,Vec w,PetscReal *beta,PetscBool *breakdown,PetscScalar *work)
127: {
128: PetscInt i,j,l,m = *M,ldh;
129: Vec t = pep->work[2],u = pep->work[3];
130: BVOrthogRefineType refinement;
131: PetscReal norm=0.0,onorm,eta;
132: PetscScalar *H,*c = work + m;
134: *beta = 0.0;
135: MatDenseGetArray(A,&H);
136: MatDenseGetLDA(A,&ldh);
137: BVGetOrthogonalization(pep->V,NULL,&refinement,&eta,NULL);
138: BVInsertVec(pep->V,k,v);
139: for (j=k;j<m;j++) {
140: /* apply operator */
141: VecCopy(w,t);
142: if (pep->Dr) VecPointwiseMult(v,v,pep->Dr);
143: STMatMult(pep->st,0,v,u);
144: VecCopy(t,v);
145: if (pep->Dr) VecPointwiseMult(t,t,pep->Dr);
146: STMatMult(pep->st,1,t,w);
147: VecAXPY(u,pep->sfactor,w);
148: STMatSolve(pep->st,u,w);
149: VecScale(w,-1.0/(pep->sfactor*pep->sfactor));
150: if (pep->Dr) VecPointwiseDivide(w,w,pep->Dr);
151: VecCopy(v,t);
152: BVSetActiveColumns(pep->V,0,j+1);
154: /* orthogonalize */
155: switch (refinement) {
156: case BV_ORTHOG_REFINE_NEVER:
157: PEPQArnoldiCGS(pep,H,ldh,H+ldh*j,j,pep->V,t,v,w,NULL,&norm,work);
158: *breakdown = PETSC_FALSE;
159: break;
160: case BV_ORTHOG_REFINE_ALWAYS:
161: PEPQArnoldiCGS(pep,H,ldh,H+ldh*j,j,pep->V,t,v,w,NULL,NULL,work);
162: PEPQArnoldiCGS(pep,H,ldh,c,j,pep->V,t,v,w,&onorm,&norm,work);
163: for (i=0;i<=j;i++) H[ldh*j+i] += c[i];
164: if (norm < eta * onorm) *breakdown = PETSC_TRUE;
165: else *breakdown = PETSC_FALSE;
166: break;
167: case BV_ORTHOG_REFINE_IFNEEDED:
168: PEPQArnoldiCGS(pep,H,ldh,H+ldh*j,j,pep->V,t,v,w,&onorm,&norm,work);
169: /* ||q|| < eta ||h|| */
170: l = 1;
171: while (l<3 && norm < eta * onorm) {
172: l++;
173: onorm = norm;
174: PEPQArnoldiCGS(pep,H,ldh,c,j,pep->V,t,v,w,NULL,&norm,work);
175: for (i=0;i<=j;i++) H[ldh*j+i] += c[i];
176: }
177: if (norm < eta * onorm) *breakdown = PETSC_TRUE;
178: else *breakdown = PETSC_FALSE;
179: break;
180: }
181: VecScale(v,1.0/norm);
182: VecScale(w,1.0/norm);
184: H[j+1+ldh*j] = norm;
185: if (j<m-1) BVInsertVec(pep->V,j+1,v);
186: }
187: *beta = norm;
188: MatDenseRestoreArray(A,&H);
189: return 0;
190: }
192: PetscErrorCode PEPSolve_QArnoldi(PEP pep)
193: {
194: PEP_QARNOLDI *ctx = (PEP_QARNOLDI*)pep->data;
195: PetscInt j,k,l,lwork,nv,nconv;
196: Vec v=pep->work[0],w=pep->work[1];
197: Mat Q,S;
198: PetscScalar *work;
199: PetscReal beta,norm,x,y;
200: PetscBool breakdown=PETSC_FALSE,sinv;
202: lwork = 7*pep->ncv;
203: PetscMalloc1(lwork,&work);
204: PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv);
205: RGPushScale(pep->rg,sinv?pep->sfactor:1.0/pep->sfactor);
206: STScaleShift(pep->st,sinv?pep->sfactor:1.0/pep->sfactor);
208: /* Get the starting Arnoldi vector */
209: for (j=0;j<2;j++) {
210: if (j>=pep->nini) BVSetRandomColumn(pep->V,j);
211: }
212: BVCopyVec(pep->V,0,v);
213: BVCopyVec(pep->V,1,w);
214: VecNorm(v,NORM_2,&x);
215: VecNorm(w,NORM_2,&y);
216: norm = SlepcAbs(x,y);
217: VecScale(v,1.0/norm);
218: VecScale(w,1.0/norm);
220: /* clean projected matrix (including the extra-arrow) */
221: DSSetDimensions(pep->ds,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
222: DSGetMat(pep->ds,DS_MAT_A,&S);
223: MatZeroEntries(S);
224: DSRestoreMat(pep->ds,DS_MAT_A,&S);
226: /* Restart loop */
227: l = 0;
228: while (pep->reason == PEP_CONVERGED_ITERATING) {
229: pep->its++;
231: /* Compute an nv-step Arnoldi factorization */
232: nv = PetscMin(pep->nconv+pep->mpd,pep->ncv);
233: DSGetMat(pep->ds,DS_MAT_A,&S);
234: PEPQArnoldi(pep,S,pep->nconv+l,&nv,v,w,&beta,&breakdown,work);
235: DSRestoreMat(pep->ds,DS_MAT_A,&S);
236: DSSetDimensions(pep->ds,nv,pep->nconv,pep->nconv+l);
237: DSSetState(pep->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE);
238: BVSetActiveColumns(pep->V,pep->nconv,nv);
240: /* Solve projected problem */
241: DSSolve(pep->ds,pep->eigr,pep->eigi);
242: DSSort(pep->ds,pep->eigr,pep->eigi,NULL,NULL,NULL);
243: DSUpdateExtraRow(pep->ds);
244: DSSynchronize(pep->ds,pep->eigr,pep->eigi);
246: /* Check convergence */
247: PEPKrylovConvergence(pep,PETSC_FALSE,pep->nconv,nv-pep->nconv,beta,&k);
248: (*pep->stopping)(pep,pep->its,pep->max_it,k,pep->nev,&pep->reason,pep->stoppingctx);
249: nconv = k;
251: /* Update l */
252: if (pep->reason != PEP_CONVERGED_ITERATING || breakdown) l = 0;
253: else {
254: l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
255: DSGetTruncateSize(pep->ds,k,nv,&l);
256: }
257: if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
258: if (l) PetscInfo(pep,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l);
260: if (pep->reason == PEP_CONVERGED_ITERATING) {
261: if (PetscUnlikely(breakdown)) {
262: /* Stop if breakdown */
263: PetscInfo(pep,"Breakdown Quadratic Arnoldi method (it=%" PetscInt_FMT " norm=%g)\n",pep->its,(double)beta);
264: pep->reason = PEP_DIVERGED_BREAKDOWN;
265: } else {
266: /* Prepare the Rayleigh quotient for restart */
267: DSTruncate(pep->ds,k+l,PETSC_FALSE);
268: }
269: }
270: /* Update the corresponding vectors V(:,idx) = V*Q(:,idx) */
271: DSGetMat(pep->ds,DS_MAT_Q,&Q);
272: BVMultInPlace(pep->V,Q,pep->nconv,k+l);
273: DSRestoreMat(pep->ds,DS_MAT_Q,&Q);
275: pep->nconv = k;
276: PEPMonitor(pep,pep->its,nconv,pep->eigr,pep->eigi,pep->errest,nv);
277: }
278: BVSetActiveColumns(pep->V,0,pep->nconv);
279: for (j=0;j<pep->nconv;j++) {
280: pep->eigr[j] *= pep->sfactor;
281: pep->eigi[j] *= pep->sfactor;
282: }
284: STScaleShift(pep->st,sinv?1.0/pep->sfactor:pep->sfactor);
285: RGPopScale(pep->rg);
287: DSTruncate(pep->ds,pep->nconv,PETSC_TRUE);
288: PetscFree(work);
289: return 0;
290: }
292: static PetscErrorCode PEPQArnoldiSetRestart_QArnoldi(PEP pep,PetscReal keep)
293: {
294: PEP_QARNOLDI *ctx = (PEP_QARNOLDI*)pep->data;
296: if (keep==PETSC_DEFAULT) ctx->keep = 0.5;
297: else {
299: ctx->keep = keep;
300: }
301: return 0;
302: }
304: /*@
305: PEPQArnoldiSetRestart - Sets the restart parameter for the Q-Arnoldi
306: method, in particular the proportion of basis vectors that must be kept
307: after restart.
309: Logically Collective on pep
311: Input Parameters:
312: + pep - the eigenproblem solver context
313: - keep - the number of vectors to be kept at restart
315: Options Database Key:
316: . -pep_qarnoldi_restart - Sets the restart parameter
318: Notes:
319: Allowed values are in the range [0.1,0.9]. The default is 0.5.
321: Level: advanced
323: .seealso: PEPQArnoldiGetRestart()
324: @*/
325: PetscErrorCode PEPQArnoldiSetRestart(PEP pep,PetscReal keep)
326: {
329: PetscTryMethod(pep,"PEPQArnoldiSetRestart_C",(PEP,PetscReal),(pep,keep));
330: return 0;
331: }
333: static PetscErrorCode PEPQArnoldiGetRestart_QArnoldi(PEP pep,PetscReal *keep)
334: {
335: PEP_QARNOLDI *ctx = (PEP_QARNOLDI*)pep->data;
337: *keep = ctx->keep;
338: return 0;
339: }
341: /*@
342: PEPQArnoldiGetRestart - Gets the restart parameter used in the Q-Arnoldi method.
344: Not Collective
346: Input Parameter:
347: . pep - the eigenproblem solver context
349: Output Parameter:
350: . keep - the restart parameter
352: Level: advanced
354: .seealso: PEPQArnoldiSetRestart()
355: @*/
356: PetscErrorCode PEPQArnoldiGetRestart(PEP pep,PetscReal *keep)
357: {
360: PetscUseMethod(pep,"PEPQArnoldiGetRestart_C",(PEP,PetscReal*),(pep,keep));
361: return 0;
362: }
364: static PetscErrorCode PEPQArnoldiSetLocking_QArnoldi(PEP pep,PetscBool lock)
365: {
366: PEP_QARNOLDI *ctx = (PEP_QARNOLDI*)pep->data;
368: ctx->lock = lock;
369: return 0;
370: }
372: /*@
373: PEPQArnoldiSetLocking - Choose between locking and non-locking variants of
374: the Q-Arnoldi method.
376: Logically Collective on pep
378: Input Parameters:
379: + pep - the eigenproblem solver context
380: - lock - true if the locking variant must be selected
382: Options Database Key:
383: . -pep_qarnoldi_locking - Sets the locking flag
385: Notes:
386: The default is to lock converged eigenpairs when the method restarts.
387: This behaviour can be changed so that all directions are kept in the
388: working subspace even if already converged to working accuracy (the
389: non-locking variant).
391: Level: advanced
393: .seealso: PEPQArnoldiGetLocking()
394: @*/
395: PetscErrorCode PEPQArnoldiSetLocking(PEP pep,PetscBool lock)
396: {
399: PetscTryMethod(pep,"PEPQArnoldiSetLocking_C",(PEP,PetscBool),(pep,lock));
400: return 0;
401: }
403: static PetscErrorCode PEPQArnoldiGetLocking_QArnoldi(PEP pep,PetscBool *lock)
404: {
405: PEP_QARNOLDI *ctx = (PEP_QARNOLDI*)pep->data;
407: *lock = ctx->lock;
408: return 0;
409: }
411: /*@
412: PEPQArnoldiGetLocking - Gets the locking flag used in the Q-Arnoldi method.
414: Not Collective
416: Input Parameter:
417: . pep - the eigenproblem solver context
419: Output Parameter:
420: . lock - the locking flag
422: Level: advanced
424: .seealso: PEPQArnoldiSetLocking()
425: @*/
426: PetscErrorCode PEPQArnoldiGetLocking(PEP pep,PetscBool *lock)
427: {
430: PetscUseMethod(pep,"PEPQArnoldiGetLocking_C",(PEP,PetscBool*),(pep,lock));
431: return 0;
432: }
434: PetscErrorCode PEPSetFromOptions_QArnoldi(PEP pep,PetscOptionItems *PetscOptionsObject)
435: {
436: PetscBool flg,lock;
437: PetscReal keep;
439: PetscOptionsHeadBegin(PetscOptionsObject,"PEP Q-Arnoldi Options");
441: PetscOptionsReal("-pep_qarnoldi_restart","Proportion of vectors kept after restart","PEPQArnoldiSetRestart",0.5,&keep,&flg);
442: if (flg) PEPQArnoldiSetRestart(pep,keep);
444: PetscOptionsBool("-pep_qarnoldi_locking","Choose between locking and non-locking variants","PEPQArnoldiSetLocking",PETSC_FALSE,&lock,&flg);
445: if (flg) PEPQArnoldiSetLocking(pep,lock);
447: PetscOptionsHeadEnd();
448: return 0;
449: }
451: PetscErrorCode PEPView_QArnoldi(PEP pep,PetscViewer viewer)
452: {
453: PEP_QARNOLDI *ctx = (PEP_QARNOLDI*)pep->data;
454: PetscBool isascii;
456: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
457: if (isascii) {
458: PetscViewerASCIIPrintf(viewer," %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep));
459: PetscViewerASCIIPrintf(viewer," using the %slocking variant\n",ctx->lock?"":"non-");
460: }
461: return 0;
462: }
464: PetscErrorCode PEPDestroy_QArnoldi(PEP pep)
465: {
466: PetscFree(pep->data);
467: PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiSetRestart_C",NULL);
468: PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiGetRestart_C",NULL);
469: PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiSetLocking_C",NULL);
470: PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiGetLocking_C",NULL);
471: return 0;
472: }
474: SLEPC_EXTERN PetscErrorCode PEPCreate_QArnoldi(PEP pep)
475: {
476: PEP_QARNOLDI *ctx;
478: PetscNew(&ctx);
479: pep->data = (void*)ctx;
481: pep->lineariz = PETSC_TRUE;
482: ctx->lock = PETSC_TRUE;
484: pep->ops->solve = PEPSolve_QArnoldi;
485: pep->ops->setup = PEPSetUp_QArnoldi;
486: pep->ops->setfromoptions = PEPSetFromOptions_QArnoldi;
487: pep->ops->destroy = PEPDestroy_QArnoldi;
488: pep->ops->view = PEPView_QArnoldi;
489: pep->ops->backtransform = PEPBackTransform_Default;
490: pep->ops->computevectors = PEPComputeVectors_Default;
491: pep->ops->extractvectors = PEPExtractVectors_QArnoldi;
492: pep->ops->setdefaultst = PEPSetDefaultST_Transform;
494: PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiSetRestart_C",PEPQArnoldiSetRestart_QArnoldi);
495: PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiGetRestart_C",PEPQArnoldiGetRestart_QArnoldi);
496: PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiSetLocking_C",PEPQArnoldiSetLocking_QArnoldi);
497: PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiGetLocking_C",PEPQArnoldiGetLocking_QArnoldi);
498: return 0;
499: }