Actual source code: lyapii.c
slepc-3.16.2 2022-02-01
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2021, 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 eigensolver: "lyapii"
13: Method: Lyapunov inverse iteration
15: Algorithm:
17: Lyapunov inverse iteration using LME solvers
19: References:
21: [1] H.C. Elman and M. Wu, "Lyapunov inverse iteration for computing a
22: few rightmost eigenvalues of large generalized eigenvalue problems",
23: SIAM J. Matrix Anal. Appl. 34(4):1685-1707, 2013.
25: [2] K. Meerbergen and A. Spence, "Inverse iteration for purely imaginary
26: eigenvalues with application to the detection of Hopf bifurcations in
27: large-scale problems", SIAM J. Matrix Anal. Appl. 31:1982-1999, 2010.
28: */
30: #include <slepc/private/epsimpl.h>
31: #include <slepcblaslapack.h>
33: typedef struct {
34: LME lme; /* Lyapunov solver */
35: DS ds; /* used to compute the SVD for compression */
36: PetscInt rkl; /* prescribed rank for the Lyapunov solver */
37: PetscInt rkc; /* the compressed rank, cannot be larger than rkl */
38: } EPS_LYAPII;
40: typedef struct {
41: Mat S; /* the operator matrix, S=A^{-1}*B */
42: BV Q; /* orthogonal basis of converged eigenvectors */
43: } EPS_LYAPII_MATSHELL;
45: typedef struct {
46: Mat S; /* the matrix from which the implicit operator is built */
47: PetscInt n; /* the size of matrix S, the operator is nxn */
48: LME lme; /* dummy LME object */
49: #if defined(PETSC_USE_COMPLEX)
50: Mat A,B,F;
51: Vec w;
52: #endif
53: } EPS_EIG_MATSHELL;
55: PetscErrorCode EPSSetUp_LyapII(EPS eps)
56: {
58: PetscRandom rand;
59: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
62: EPSCheckSinvert(eps);
63: if (eps->ncv!=PETSC_DEFAULT) {
64: if (eps->ncv<eps->nev+1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nev+1");
65: } else eps->ncv = eps->nev+1;
66: if (eps->mpd!=PETSC_DEFAULT) { PetscInfo(eps,"Warning: parameter mpd ignored\n"); }
67: if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(1000*eps->nev,100*eps->n);
68: if (!eps->which) eps->which=EPS_LARGEST_REAL;
69: if (eps->which!=EPS_LARGEST_REAL) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver supports only largest real eigenvalues");
70: EPSCheckUnsupported(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION | EPS_FEATURE_TWOSIDED);
72: if (!ctx->rkc) ctx->rkc = 10;
73: if (!ctx->rkl) ctx->rkl = 3*ctx->rkc;
74: if (!ctx->lme) { EPSLyapIIGetLME(eps,&ctx->lme); }
75: LMESetProblemType(ctx->lme,LME_LYAPUNOV);
76: LMESetErrorIfNotConverged(ctx->lme,PETSC_TRUE);
78: if (!ctx->ds) {
79: DSCreate(PetscObjectComm((PetscObject)eps),&ctx->ds);
80: PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->ds);
81: DSSetType(ctx->ds,DSSVD);
82: }
83: DSAllocate(ctx->ds,ctx->rkl);
85: DSSetType(eps->ds,DSNHEP);
86: DSAllocate(eps->ds,eps->ncv);
88: EPSAllocateSolution(eps,0);
89: BVGetRandomContext(eps->V,&rand); /* make sure the random context is available when duplicating */
90: EPSSetWorkVecs(eps,3);
91: return(0);
92: }
94: static PetscErrorCode MatMult_EPSLyapIIOperator(Mat M,Vec x,Vec r)
95: {
96: PetscErrorCode ierr;
97: EPS_LYAPII_MATSHELL *matctx;
100: MatShellGetContext(M,&matctx);
101: MatMult(matctx->S,x,r);
102: BVOrthogonalizeVec(matctx->Q,r,NULL,NULL,NULL);
103: return(0);
104: }
106: static PetscErrorCode MatDestroy_EPSLyapIIOperator(Mat M)
107: {
108: PetscErrorCode ierr;
109: EPS_LYAPII_MATSHELL *matctx;
112: MatShellGetContext(M,&matctx);
113: MatDestroy(&matctx->S);
114: PetscFree(matctx);
115: return(0);
116: }
118: static PetscErrorCode MatMult_EigOperator(Mat M,Vec x,Vec y)
119: {
120: PetscErrorCode ierr;
121: EPS_EIG_MATSHELL *matctx;
122: #if !defined(PETSC_USE_COMPLEX)
123: PetscInt n;
124: PetscScalar *Y,*C,zero=0.0,done=1.0,dtwo=2.0;
125: const PetscScalar *S,*X;
126: PetscBLASInt n_;
127: #endif
130: MatShellGetContext(M,&matctx);
132: #if defined(PETSC_USE_COMPLEX)
133: MatMult(matctx->B,x,matctx->w);
134: MatSolve(matctx->F,matctx->w,y);
135: #else
136: VecGetArrayRead(x,&X);
137: VecGetArray(y,&Y);
138: MatDenseGetArrayRead(matctx->S,&S);
140: n = matctx->n;
141: PetscCalloc1(n*n,&C);
142: PetscBLASIntCast(n,&n_);
144: /* C = 2*S*X*S.' */
145: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&dtwo,S,&n_,X,&n_,&zero,Y,&n_));
146: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&n_,&n_,&n_,&done,Y,&n_,S,&n_,&zero,C,&n_));
148: /* Solve S*Y + Y*S' = -C */
149: LMEDenseLyapunov(matctx->lme,n,(PetscScalar*)S,n,C,n,Y,n);
151: PetscFree(C);
152: VecRestoreArrayRead(x,&X);
153: VecRestoreArray(y,&Y);
154: MatDenseRestoreArrayRead(matctx->S,&S);
155: #endif
156: return(0);
157: }
159: static PetscErrorCode MatDestroy_EigOperator(Mat M)
160: {
161: PetscErrorCode ierr;
162: EPS_EIG_MATSHELL *matctx;
165: MatShellGetContext(M,&matctx);
166: #if defined(PETSC_USE_COMPLEX)
167: MatDestroy(&matctx->A);
168: MatDestroy(&matctx->B);
169: MatDestroy(&matctx->F);
170: VecDestroy(&matctx->w);
171: #endif
172: PetscFree(matctx);
173: return(0);
174: }
176: /*
177: EV2x2: solve the eigenproblem for a 2x2 matrix M
178: */
179: static PetscErrorCode EV2x2(PetscScalar *M,PetscInt ld,PetscScalar *wr,PetscScalar *wi,PetscScalar *vec)
180: {
182: PetscBLASInt lwork=10,ld_;
183: PetscScalar work[10];
184: PetscBLASInt two=2,info;
185: #if defined(PETSC_USE_COMPLEX)
186: PetscReal rwork[6];
187: #endif
190: PetscBLASIntCast(ld,&ld_);
191: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
192: #if !defined(PETSC_USE_COMPLEX)
193: PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","V",&two,M,&ld_,wr,wi,NULL,&ld_,vec,&ld_,work,&lwork,&info));
194: #else
195: PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","V",&two,M,&ld_,wr,NULL,&ld_,vec,&ld_,work,&lwork,rwork,&info));
196: #endif
197: SlepcCheckLapackInfo("geev",info);
198: PetscFPTrapPop();
199: return(0);
200: }
202: /*
203: LyapIIBuildRHS: prepare the right-hand side of the Lyapunov equation SY + YS' = -2*S*Z*S'
204: in factored form:
205: if (V) U=sqrt(2)*S*V (uses 1 work vector)
206: else U=sqrt(2)*S*U (uses 2 work vectors)
207: where U,V are assumed to have rk columns.
208: */
209: static PetscErrorCode LyapIIBuildRHS(Mat S,PetscInt rk,Mat U,BV V,Vec *work)
210: {
212: PetscScalar *array,*uu;
213: PetscInt i,nloc;
214: Vec v,u=work[0];
217: MatGetLocalSize(U,&nloc,NULL);
218: for (i=0;i<rk;i++) {
219: MatDenseGetColumn(U,i,&array);
220: if (V) {
221: BVGetColumn(V,i,&v);
222: } else {
223: v = work[1];
224: VecPlaceArray(v,array);
225: }
226: MatMult(S,v,u);
227: if (V) {
228: BVRestoreColumn(V,i,&v);
229: } else {
230: VecResetArray(v);
231: }
232: VecScale(u,PETSC_SQRT2);
233: VecGetArray(u,&uu);
234: PetscArraycpy(array,uu,nloc);
235: VecRestoreArray(u,&uu);
236: MatDenseRestoreColumn(U,&array);
237: }
238: return(0);
239: }
241: /*
242: LyapIIBuildEigenMat: create shell matrix Op=A\B with A = kron(I,S)+kron(S,I), B = -2*kron(S,S)
243: where S is a sequential square dense matrix of order n.
244: v0 is the initial vector, should have the form v0 = w*w' (for instance 1*1')
245: */
246: static PetscErrorCode LyapIIBuildEigenMat(LME lme,Mat S,Mat *Op,Vec *v0)
247: {
248: PetscErrorCode ierr;
249: PetscInt n,m;
250: PetscBool create=PETSC_FALSE;
251: EPS_EIG_MATSHELL *matctx;
252: #if defined(PETSC_USE_COMPLEX)
253: PetscScalar theta,*aa,*bb;
254: const PetscScalar *ss;
255: PetscInt i,j,f,c,off,ld;
256: IS perm;
257: #endif
260: MatGetSize(S,&n,NULL);
261: if (!*Op) create=PETSC_TRUE;
262: else {
263: MatGetSize(*Op,&m,NULL);
264: if (m!=n*n) create=PETSC_TRUE;
265: }
266: if (create) {
267: MatDestroy(Op);
268: VecDestroy(v0);
269: PetscNew(&matctx);
270: #if defined(PETSC_USE_COMPLEX)
271: MatCreateSeqDense(PETSC_COMM_SELF,n*n,n*n,NULL,&matctx->A);
272: MatCreateSeqDense(PETSC_COMM_SELF,n*n,n*n,NULL,&matctx->B);
273: MatCreateVecs(matctx->A,NULL,&matctx->w);
274: #endif
275: MatCreateShell(PETSC_COMM_SELF,n*n,n*n,PETSC_DETERMINE,PETSC_DETERMINE,matctx,Op);
276: MatShellSetOperation(*Op,MATOP_MULT,(void(*)(void))MatMult_EigOperator);
277: MatShellSetOperation(*Op,MATOP_DESTROY,(void(*)(void))MatDestroy_EigOperator);
278: MatCreateVecs(*Op,NULL,v0);
279: } else {
280: MatShellGetContext(*Op,&matctx);
281: #if defined(PETSC_USE_COMPLEX)
282: MatZeroEntries(matctx->A);
283: #endif
284: }
285: #if defined(PETSC_USE_COMPLEX)
286: MatDenseGetArray(matctx->A,&aa);
287: MatDenseGetArray(matctx->B,&bb);
288: MatDenseGetArrayRead(S,&ss);
289: ld = n*n;
290: for (f=0;f<n;f++) {
291: off = f*n+f*n*ld;
292: for (i=0;i<n;i++) for (j=0;j<n;j++) aa[off+i+j*ld] = ss[i+j*n];
293: for (c=0;c<n;c++) {
294: off = f*n+c*n*ld;
295: theta = ss[f+c*n];
296: for (i=0;i<n;i++) aa[off+i+i*ld] += theta;
297: for (i=0;i<n;i++) for (j=0;j<n;j++) bb[off+i+j*ld] = -2*theta*ss[i+j*n];
298: }
299: }
300: MatDenseRestoreArray(matctx->A,&aa);
301: MatDenseRestoreArray(matctx->B,&bb);
302: MatDenseRestoreArrayRead(S,&ss);
303: ISCreateStride(PETSC_COMM_SELF,n*n,0,1,&perm);
304: MatDestroy(&matctx->F);
305: MatDuplicate(matctx->A,MAT_COPY_VALUES,&matctx->F);
306: MatLUFactor(matctx->F,perm,perm,0);
307: ISDestroy(&perm);
308: #endif
309: matctx->lme = lme;
310: matctx->S = S;
311: matctx->n = n;
312: VecSet(*v0,1.0);
313: return(0);
314: }
316: PetscErrorCode EPSSolve_LyapII(EPS eps)
317: {
318: PetscErrorCode ierr;
319: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
320: PetscInt i,ldds,rk,nloc,mloc,nv,idx,k;
321: Vec v,w,z=eps->work[0],v0=NULL;
322: Mat S,C,Ux[2],Y,Y1,R,U,W,X,Op=NULL;
323: BV V;
324: BVOrthogType type;
325: BVOrthogRefineType refine;
326: PetscScalar eigr[2],eigi[2],*array,er,ei,*uu,*s,*xx,*aa,pM[4],vec[4];
327: PetscReal eta;
328: EPS epsrr;
329: PetscReal norm;
330: EPS_LYAPII_MATSHELL *matctx;
333: DSGetLeadingDimension(ctx->ds,&ldds);
335: /* Operator for the Lyapunov equation */
336: PetscNew(&matctx);
337: STGetOperator(eps->st,&matctx->S);
338: MatGetLocalSize(matctx->S,&mloc,&nloc);
339: MatCreateShell(PetscObjectComm((PetscObject)eps),mloc,nloc,PETSC_DETERMINE,PETSC_DETERMINE,matctx,&S);
340: matctx->Q = eps->V;
341: MatShellSetOperation(S,MATOP_MULT,(void(*)(void))MatMult_EPSLyapIIOperator);
342: MatShellSetOperation(S,MATOP_DESTROY,(void(*)(void))MatDestroy_EPSLyapIIOperator);
343: LMESetCoefficients(ctx->lme,S,NULL,NULL,NULL);
345: /* Right-hand side */
346: BVDuplicateResize(eps->V,ctx->rkl,&V);
347: BVGetOrthogonalization(V,&type,&refine,&eta,NULL);
348: BVSetOrthogonalization(V,type,refine,eta,BV_ORTHOG_BLOCK_TSQR);
349: MatCreateDense(PetscObjectComm((PetscObject)eps),eps->nloc,PETSC_DECIDE,PETSC_DECIDE,1,NULL,&Ux[0]);
350: MatCreateDense(PetscObjectComm((PetscObject)eps),eps->nloc,PETSC_DECIDE,PETSC_DECIDE,2,NULL,&Ux[1]);
351: nv = ctx->rkl;
352: PetscMalloc1(nv,&s);
354: /* Initialize first column */
355: EPSGetStartVector(eps,0,NULL);
356: BVGetColumn(eps->V,0,&v);
357: BVInsertVec(V,0,v);
358: BVRestoreColumn(eps->V,0,&v);
359: BVSetActiveColumns(eps->V,0,0); /* no deflation at the beginning */
360: LyapIIBuildRHS(S,1,Ux[0],V,eps->work);
361: idx = 0;
363: /* EPS for rank reduction */
364: EPSCreate(PETSC_COMM_SELF,&epsrr);
365: EPSSetOptionsPrefix(epsrr,((PetscObject)eps)->prefix);
366: EPSAppendOptionsPrefix(epsrr,"eps_lyapii_");
367: EPSSetDimensions(epsrr,1,PETSC_DEFAULT,PETSC_DEFAULT);
368: EPSSetTolerances(epsrr,PETSC_MACHINE_EPSILON*100,PETSC_DEFAULT);
370: while (eps->reason == EPS_CONVERGED_ITERATING) {
371: eps->its++;
373: /* Matrix for placing the solution of the Lyapunov equation (an alias of V) */
374: BVSetActiveColumns(V,0,nv);
375: BVGetMat(V,&Y1);
376: MatZeroEntries(Y1);
377: MatCreateLRC(NULL,Y1,NULL,NULL,&Y);
378: LMESetSolution(ctx->lme,Y);
380: /* Solve the Lyapunov equation SY + YS' = -2*S*Z*S' */
381: MatCreateLRC(NULL,Ux[idx],NULL,NULL,&C);
382: LMESetRHS(ctx->lme,C);
383: MatDestroy(&C);
384: LMESolve(ctx->lme);
385: BVRestoreMat(V,&Y1);
386: MatDestroy(&Y);
388: /* SVD of the solution: [Q,R]=qr(V); [U,Sigma,~]=svd(R) */
389: DSSetDimensions(ctx->ds,nv,0,0);
390: DSSVDSetDimensions(ctx->ds,nv);
391: DSGetMat(ctx->ds,DS_MAT_A,&R);
392: BVOrthogonalize(V,R);
393: DSRestoreMat(ctx->ds,DS_MAT_A,&R);
394: DSSetState(ctx->ds,DS_STATE_RAW);
395: DSSolve(ctx->ds,s,NULL);
397: /* Determine rank */
398: rk = nv;
399: for (i=1;i<nv;i++) if (PetscAbsScalar(s[i]/s[0])<PETSC_SQRT_MACHINE_EPSILON) {rk=i; break;}
400: PetscInfo1(eps,"The computed solution of the Lyapunov equation has rank %D\n",rk);
401: rk = PetscMin(rk,ctx->rkc);
402: DSGetMat(ctx->ds,DS_MAT_U,&U);
403: BVMultInPlace(V,U,0,rk);
404: BVSetActiveColumns(V,0,rk);
405: MatDestroy(&U);
407: /* Rank reduction */
408: DSSetDimensions(ctx->ds,rk,0,0);
409: DSSVDSetDimensions(ctx->ds,rk);
410: DSGetMat(ctx->ds,DS_MAT_A,&W);
411: BVMatProject(V,S,V,W);
412: LyapIIBuildEigenMat(ctx->lme,W,&Op,&v0); /* Op=A\B, A=kron(I,S)+kron(S,I), B=-2*kron(S,S) */
413: EPSSetOperators(epsrr,Op,NULL);
414: EPSSetInitialSpace(epsrr,1,&v0);
415: EPSSolve(epsrr);
416: MatDestroy(&W);
417: EPSComputeVectors(epsrr);
418: /* Copy first eigenvector, vec(A)=x */
419: BVGetArray(epsrr->V,&xx);
420: DSGetArray(ctx->ds,DS_MAT_A,&aa);
421: for (i=0;i<rk;i++) {
422: PetscArraycpy(aa+i*ldds,xx+i*rk,rk);
423: }
424: DSRestoreArray(ctx->ds,DS_MAT_A,&aa);
425: BVRestoreArray(epsrr->V,&xx);
426: DSSetState(ctx->ds,DS_STATE_RAW);
427: /* Compute [U,Sigma,~] = svd(A), its rank should be 1 or 2 */
428: DSSolve(ctx->ds,s,NULL);
429: if (PetscAbsScalar(s[1]/s[0])<PETSC_SQRT_MACHINE_EPSILON) rk=1;
430: else rk = 2;
431: PetscInfo1(eps,"The eigenvector has rank %D\n",rk);
432: DSGetMat(ctx->ds,DS_MAT_U,&U);
433: BVMultInPlace(V,U,0,rk);
434: MatDestroy(&U);
436: /* Save V in Ux */
437: idx = (rk==2)?1:0;
438: for (i=0;i<rk;i++) {
439: BVGetColumn(V,i,&v);
440: VecGetArray(v,&uu);
441: MatDenseGetColumn(Ux[idx],i,&array);
442: PetscArraycpy(array,uu,eps->nloc);
443: MatDenseRestoreColumn(Ux[idx],&array);
444: VecRestoreArray(v,&uu);
445: BVRestoreColumn(V,i,&v);
446: }
448: /* Eigenpair approximation */
449: BVGetColumn(V,0,&v);
450: MatMult(S,v,z);
451: VecDot(z,v,pM);
452: BVRestoreColumn(V,0,&v);
453: if (rk>1) {
454: BVGetColumn(V,1,&w);
455: VecDot(z,w,pM+1);
456: MatMult(S,w,z);
457: VecDot(z,w,pM+3);
458: BVGetColumn(V,0,&v);
459: VecDot(z,v,pM+2);
460: BVRestoreColumn(V,0,&v);
461: BVRestoreColumn(V,1,&w);
462: EV2x2(pM,2,eigr,eigi,vec);
463: MatCreateSeqDense(PETSC_COMM_SELF,2,2,vec,&X);
464: BVSetActiveColumns(V,0,rk);
465: BVMultInPlace(V,X,0,rk);
466: MatDestroy(&X);
467: #if !defined(PETSC_USE_COMPLEX)
468: norm = eigr[0]*eigr[0]+eigi[0]*eigi[0];
469: er = eigr[0]/norm; ei = -eigi[0]/norm;
470: #else
471: er =1.0/eigr[0]; ei = 0.0;
472: #endif
473: } else {
474: eigr[0] = pM[0]; eigi[0] = 0.0;
475: er = 1.0/eigr[0]; ei = 0.0;
476: }
477: BVGetColumn(V,0,&v);
478: if (eigi[0]!=0.0) {
479: BVGetColumn(V,1,&w);
480: } else w = NULL;
481: eps->eigr[eps->nconv] = eigr[0]; eps->eigi[eps->nconv] = eigi[0];
482: EPSComputeResidualNorm_Private(eps,PETSC_FALSE,er,ei,v,w,eps->work,&norm);
483: BVRestoreColumn(V,0,&v);
484: if (w) {
485: BVRestoreColumn(V,1,&w);
486: }
487: (*eps->converged)(eps,er,ei,norm,&eps->errest[eps->nconv],eps->convergedctx);
488: k = 0;
489: if (eps->errest[eps->nconv]<eps->tol) {
490: k++;
491: if (rk==2) {
492: #if !defined (PETSC_USE_COMPLEX)
493: eps->eigr[eps->nconv+k] = eigr[0]; eps->eigi[eps->nconv+k] = -eigi[0];
494: #else
495: eps->eigr[eps->nconv+k] = PetscConj(eps->eigr[eps->nconv]);
496: #endif
497: k++;
498: }
499: /* Store converged eigenpairs and vectors for deflation */
500: for (i=0;i<k;i++) {
501: BVGetColumn(V,i,&v);
502: BVInsertVec(eps->V,eps->nconv+i,v);
503: BVRestoreColumn(V,i,&v);
504: }
505: eps->nconv += k;
506: BVSetActiveColumns(eps->V,eps->nconv-rk,eps->nconv);
507: BVOrthogonalize(eps->V,NULL);
508: DSSetDimensions(eps->ds,eps->nconv,0,0);
509: DSGetMat(eps->ds,DS_MAT_A,&W);
510: BVMatProject(eps->V,matctx->S,eps->V,W);
511: DSRestoreMat(eps->ds,DS_MAT_A,&W);
512: if (eps->nconv<eps->nev) {
513: idx = 0;
514: BVSetRandomColumn(V,0);
515: BVNormColumn(V,0,NORM_2,&norm);
516: BVScaleColumn(V,0,1.0/norm);
517: LyapIIBuildRHS(S,1,Ux[idx],V,eps->work);
518: }
519: } else {
520: /* Prepare right-hand side */
521: LyapIIBuildRHS(S,rk,Ux[idx],NULL,eps->work);
522: }
523: (*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx);
524: EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,eps->nconv+1);
525: }
526: STRestoreOperator(eps->st,&matctx->S);
527: MatDestroy(&S);
528: MatDestroy(&Ux[0]);
529: MatDestroy(&Ux[1]);
530: MatDestroy(&Op);
531: VecDestroy(&v0);
532: BVDestroy(&V);
533: EPSDestroy(&epsrr);
534: PetscFree(s);
535: return(0);
536: }
538: PetscErrorCode EPSSetFromOptions_LyapII(PetscOptionItems *PetscOptionsObject,EPS eps)
539: {
541: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
542: PetscInt k,array[2]={PETSC_DEFAULT,PETSC_DEFAULT};
543: PetscBool flg;
546: PetscOptionsHead(PetscOptionsObject,"EPS Lyapunov Inverse Iteration Options");
548: k = 2;
549: PetscOptionsIntArray("-eps_lyapii_ranks","Ranks for Lyapunov equation (one or two comma-separated integers)","EPSLyapIISetRanks",array,&k,&flg);
550: if (flg) {
551: EPSLyapIISetRanks(eps,array[0],array[1]);
552: }
554: PetscOptionsTail();
556: if (!ctx->lme) { EPSLyapIIGetLME(eps,&ctx->lme); }
557: LMESetFromOptions(ctx->lme);
558: return(0);
559: }
561: static PetscErrorCode EPSLyapIISetRanks_LyapII(EPS eps,PetscInt rkc,PetscInt rkl)
562: {
563: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
566: if (rkc==PETSC_DEFAULT) rkc = 10;
567: if (rkc<2) SETERRQ1(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The compressed rank %D must be larger than 1",rkc);
568: if (rkl==PETSC_DEFAULT) rkl = 3*rkc;
569: if (rkl<rkc) SETERRQ2(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The Lyapunov rank %D cannot be smaller than the compressed rank %D",rkl,rkc);
570: if (rkc != ctx->rkc) {
571: ctx->rkc = rkc;
572: eps->state = EPS_STATE_INITIAL;
573: }
574: if (rkl != ctx->rkl) {
575: ctx->rkl = rkl;
576: eps->state = EPS_STATE_INITIAL;
577: }
578: return(0);
579: }
581: /*@
582: EPSLyapIISetRanks - Set the ranks used in the solution of the Lyapunov equation.
584: Logically Collective on EPS
586: Input Parameters:
587: + eps - the eigenproblem solver context
588: . rkc - the compressed rank
589: - rkl - the Lyapunov rank
591: Options Database Key:
592: . -eps_lyapii_ranks <rkc,rkl> - Sets the rank parameters
594: Notes:
595: Lyapunov inverse iteration needs to solve a large-scale Lyapunov equation
596: at each iteration of the eigensolver. For this, an iterative solver (LME)
597: is used, which requires to prescribe the rank of the solution matrix X. This
598: is the meaning of parameter rkl. Later, this matrix is compressed into
599: another matrix of rank rkc. If not provided, rkl is a small multiple of rkc.
601: Level: intermediate
603: .seealso: EPSLyapIIGetRanks()
604: @*/
605: PetscErrorCode EPSLyapIISetRanks(EPS eps,PetscInt rkc,PetscInt rkl)
606: {
613: PetscTryMethod(eps,"EPSLyapIISetRanks_C",(EPS,PetscInt,PetscInt),(eps,rkc,rkl));
614: return(0);
615: }
617: static PetscErrorCode EPSLyapIIGetRanks_LyapII(EPS eps,PetscInt *rkc,PetscInt *rkl)
618: {
619: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
622: if (rkc) *rkc = ctx->rkc;
623: if (rkl) *rkl = ctx->rkl;
624: return(0);
625: }
627: /*@
628: EPSLyapIIGetRanks - Return the rank values used for the Lyapunov step.
630: Not Collective
632: Input Parameter:
633: . eps - the eigenproblem solver context
635: Output Parameters:
636: + rkc - the compressed rank
637: - rkl - the Lyapunov rank
639: Level: intermediate
641: .seealso: EPSLyapIISetRanks()
642: @*/
643: PetscErrorCode EPSLyapIIGetRanks(EPS eps,PetscInt *rkc,PetscInt *rkl)
644: {
649: PetscUseMethod(eps,"EPSLyapIIGetRanks_C",(EPS,PetscInt*,PetscInt*),(eps,rkc,rkl));
650: return(0);
651: }
653: static PetscErrorCode EPSLyapIISetLME_LyapII(EPS eps,LME lme)
654: {
656: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
659: PetscObjectReference((PetscObject)lme);
660: LMEDestroy(&ctx->lme);
661: ctx->lme = lme;
662: PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->lme);
663: eps->state = EPS_STATE_INITIAL;
664: return(0);
665: }
667: /*@
668: EPSLyapIISetLME - Associate a linear matrix equation solver object (LME) to the
669: eigenvalue solver.
671: Collective on EPS
673: Input Parameters:
674: + eps - the eigenproblem solver context
675: - lme - the linear matrix equation solver object
677: Level: advanced
679: .seealso: EPSLyapIIGetLME()
680: @*/
681: PetscErrorCode EPSLyapIISetLME(EPS eps,LME lme)
682: {
689: PetscTryMethod(eps,"EPSLyapIISetLME_C",(EPS,LME),(eps,lme));
690: return(0);
691: }
693: static PetscErrorCode EPSLyapIIGetLME_LyapII(EPS eps,LME *lme)
694: {
696: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
699: if (!ctx->lme) {
700: LMECreate(PetscObjectComm((PetscObject)eps),&ctx->lme);
701: LMESetOptionsPrefix(ctx->lme,((PetscObject)eps)->prefix);
702: LMEAppendOptionsPrefix(ctx->lme,"eps_lyapii_");
703: PetscObjectIncrementTabLevel((PetscObject)ctx->lme,(PetscObject)eps,1);
704: PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->lme);
705: }
706: *lme = ctx->lme;
707: return(0);
708: }
710: /*@
711: EPSLyapIIGetLME - Retrieve the linear matrix equation solver object (LME)
712: associated with the eigenvalue solver.
714: Not Collective
716: Input Parameter:
717: . eps - the eigenproblem solver context
719: Output Parameter:
720: . lme - the linear matrix equation solver object
722: Level: advanced
724: .seealso: EPSLyapIISetLME()
725: @*/
726: PetscErrorCode EPSLyapIIGetLME(EPS eps,LME *lme)
727: {
733: PetscUseMethod(eps,"EPSLyapIIGetLME_C",(EPS,LME*),(eps,lme));
734: return(0);
735: }
737: PetscErrorCode EPSView_LyapII(EPS eps,PetscViewer viewer)
738: {
740: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
741: PetscBool isascii;
744: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
745: if (isascii) {
746: PetscViewerASCIIPrintf(viewer," ranks: for Lyapunov solver=%D, after compression=%D\n",ctx->rkl,ctx->rkc);
747: if (!ctx->lme) { EPSLyapIIGetLME(eps,&ctx->lme); }
748: PetscViewerASCIIPushTab(viewer);
749: LMEView(ctx->lme,viewer);
750: PetscViewerASCIIPopTab(viewer);
751: }
752: return(0);
753: }
755: PetscErrorCode EPSReset_LyapII(EPS eps)
756: {
758: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
761: if (!ctx->lme) { LMEReset(ctx->lme); }
762: return(0);
763: }
765: PetscErrorCode EPSDestroy_LyapII(EPS eps)
766: {
768: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
771: LMEDestroy(&ctx->lme);
772: DSDestroy(&ctx->ds);
773: PetscFree(eps->data);
774: PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetLME_C",NULL);
775: PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetLME_C",NULL);
776: PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetRanks_C",NULL);
777: PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetRanks_C",NULL);
778: return(0);
779: }
781: PetscErrorCode EPSSetDefaultST_LyapII(EPS eps)
782: {
786: if (!((PetscObject)eps->st)->type_name) {
787: STSetType(eps->st,STSINVERT);
788: }
789: return(0);
790: }
792: SLEPC_EXTERN PetscErrorCode EPSCreate_LyapII(EPS eps)
793: {
794: EPS_LYAPII *ctx;
798: PetscNewLog(eps,&ctx);
799: eps->data = (void*)ctx;
801: eps->useds = PETSC_TRUE;
803: eps->ops->solve = EPSSolve_LyapII;
804: eps->ops->setup = EPSSetUp_LyapII;
805: eps->ops->setupsort = EPSSetUpSort_Default;
806: eps->ops->setfromoptions = EPSSetFromOptions_LyapII;
807: eps->ops->reset = EPSReset_LyapII;
808: eps->ops->destroy = EPSDestroy_LyapII;
809: eps->ops->view = EPSView_LyapII;
810: eps->ops->setdefaultst = EPSSetDefaultST_LyapII;
811: eps->ops->backtransform = EPSBackTransform_Default;
812: eps->ops->computevectors = EPSComputeVectors_Schur;
814: PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetLME_C",EPSLyapIISetLME_LyapII);
815: PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetLME_C",EPSLyapIIGetLME_LyapII);
816: PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetRanks_C",EPSLyapIISetRanks_LyapII);
817: PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetRanks_C",EPSLyapIIGetRanks_LyapII);
818: return(0);
819: }