Actual source code: lyapii.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 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: {
57: PetscRandom rand;
58: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
60: EPSCheckSinvert(eps);
61: if (eps->ncv!=PETSC_DEFAULT) {
63: } else eps->ncv = eps->nev+1;
64: if (eps->mpd!=PETSC_DEFAULT) PetscInfo(eps,"Warning: parameter mpd ignored\n");
65: if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(1000*eps->nev,100*eps->n);
66: if (!eps->which) eps->which=EPS_LARGEST_REAL;
68: EPSCheckUnsupported(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION | EPS_FEATURE_TWOSIDED);
70: if (!ctx->rkc) ctx->rkc = 10;
71: if (!ctx->rkl) ctx->rkl = 3*ctx->rkc;
72: if (!ctx->lme) EPSLyapIIGetLME(eps,&ctx->lme);
73: LMESetProblemType(ctx->lme,LME_LYAPUNOV);
74: LMESetErrorIfNotConverged(ctx->lme,PETSC_TRUE);
76: if (!ctx->ds) {
77: DSCreate(PetscObjectComm((PetscObject)eps),&ctx->ds);
78: DSSetType(ctx->ds,DSSVD);
79: }
80: DSAllocate(ctx->ds,ctx->rkl);
82: DSSetType(eps->ds,DSNHEP);
83: DSAllocate(eps->ds,eps->ncv);
85: EPSAllocateSolution(eps,0);
86: BVGetRandomContext(eps->V,&rand); /* make sure the random context is available when duplicating */
87: EPSSetWorkVecs(eps,3);
88: return 0;
89: }
91: static PetscErrorCode MatMult_EPSLyapIIOperator(Mat M,Vec x,Vec r)
92: {
93: EPS_LYAPII_MATSHELL *matctx;
95: MatShellGetContext(M,&matctx);
96: MatMult(matctx->S,x,r);
97: BVOrthogonalizeVec(matctx->Q,r,NULL,NULL,NULL);
98: return 0;
99: }
101: static PetscErrorCode MatDestroy_EPSLyapIIOperator(Mat M)
102: {
103: EPS_LYAPII_MATSHELL *matctx;
105: MatShellGetContext(M,&matctx);
106: MatDestroy(&matctx->S);
107: PetscFree(matctx);
108: return 0;
109: }
111: static PetscErrorCode MatMult_EigOperator(Mat M,Vec x,Vec y)
112: {
113: EPS_EIG_MATSHELL *matctx;
114: #if !defined(PETSC_USE_COMPLEX)
115: PetscInt n,lds;
116: PetscScalar *Y,*C,zero=0.0,done=1.0,dtwo=2.0;
117: const PetscScalar *S,*X;
118: PetscBLASInt n_,lds_;
119: #endif
121: MatShellGetContext(M,&matctx);
123: #if defined(PETSC_USE_COMPLEX)
124: MatMult(matctx->B,x,matctx->w);
125: MatSolve(matctx->F,matctx->w,y);
126: #else
127: VecGetArrayRead(x,&X);
128: VecGetArray(y,&Y);
129: MatDenseGetArrayRead(matctx->S,&S);
130: MatDenseGetLDA(matctx->S,&lds);
132: n = matctx->n;
133: PetscCalloc1(n*n,&C);
134: PetscBLASIntCast(n,&n_);
135: PetscBLASIntCast(lds,&lds_);
137: /* C = 2*S*X*S.' */
138: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&dtwo,S,&lds_,X,&n_,&zero,Y,&n_));
139: PetscCallBLAS("BLASgemm",BLASgemm_("N","T",&n_,&n_,&n_,&done,Y,&n_,S,&lds_,&zero,C,&n_));
141: /* Solve S*Y + Y*S' = -C */
142: LMEDenseLyapunov(matctx->lme,n,(PetscScalar*)S,lds,C,n,Y,n);
144: PetscFree(C);
145: VecRestoreArrayRead(x,&X);
146: VecRestoreArray(y,&Y);
147: MatDenseRestoreArrayRead(matctx->S,&S);
148: #endif
149: return 0;
150: }
152: static PetscErrorCode MatDestroy_EigOperator(Mat M)
153: {
154: EPS_EIG_MATSHELL *matctx;
156: MatShellGetContext(M,&matctx);
157: #if defined(PETSC_USE_COMPLEX)
158: MatDestroy(&matctx->A);
159: MatDestroy(&matctx->B);
160: MatDestroy(&matctx->F);
161: VecDestroy(&matctx->w);
162: #else
163: MatDestroy(&matctx->S);
164: #endif
165: PetscFree(matctx);
166: return 0;
167: }
169: /*
170: EV2x2: solve the eigenproblem for a 2x2 matrix M
171: */
172: static PetscErrorCode EV2x2(PetscScalar *M,PetscInt ld,PetscScalar *wr,PetscScalar *wi,PetscScalar *vec)
173: {
174: PetscBLASInt lwork=10,ld_;
175: PetscScalar work[10];
176: PetscBLASInt two=2,info;
177: #if defined(PETSC_USE_COMPLEX)
178: PetscReal rwork[6];
179: #endif
181: PetscBLASIntCast(ld,&ld_);
182: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
183: #if !defined(PETSC_USE_COMPLEX)
184: PetscCallBLAS("LAPACKgeev",LAPACKgeev_("N","V",&two,M,&ld_,wr,wi,NULL,&ld_,vec,&ld_,work,&lwork,&info));
185: #else
186: PetscCallBLAS("LAPACKgeev",LAPACKgeev_("N","V",&two,M,&ld_,wr,NULL,&ld_,vec,&ld_,work,&lwork,rwork,&info));
187: #endif
188: SlepcCheckLapackInfo("geev",info);
189: PetscFPTrapPop();
190: return 0;
191: }
193: /*
194: LyapIIBuildRHS: prepare the right-hand side of the Lyapunov equation SY + YS' = -2*S*Z*S'
195: in factored form:
196: if (V) U=sqrt(2)*S*V (uses 1 work vector)
197: else U=sqrt(2)*S*U (uses 2 work vectors)
198: where U,V are assumed to have rk columns.
199: */
200: static PetscErrorCode LyapIIBuildRHS(Mat S,PetscInt rk,Mat U,BV V,Vec *work)
201: {
202: PetscScalar *array,*uu;
203: PetscInt i,nloc;
204: Vec v,u=work[0];
206: MatGetLocalSize(U,&nloc,NULL);
207: for (i=0;i<rk;i++) {
208: MatDenseGetColumn(U,i,&array);
209: if (V) BVGetColumn(V,i,&v);
210: else {
211: v = work[1];
212: VecPlaceArray(v,array);
213: }
214: MatMult(S,v,u);
215: if (V) BVRestoreColumn(V,i,&v);
216: else VecResetArray(v);
217: VecScale(u,PETSC_SQRT2);
218: VecGetArray(u,&uu);
219: PetscArraycpy(array,uu,nloc);
220: VecRestoreArray(u,&uu);
221: MatDenseRestoreColumn(U,&array);
222: }
223: return 0;
224: }
226: /*
227: LyapIIBuildEigenMat: create shell matrix Op=A\B with A = kron(I,S)+kron(S,I), B = -2*kron(S,S)
228: where S is a sequential square dense matrix of order n.
229: v0 is the initial vector, should have the form v0 = w*w' (for instance 1*1')
230: */
231: static PetscErrorCode LyapIIBuildEigenMat(LME lme,Mat S,Mat *Op,Vec *v0)
232: {
233: PetscInt n,m;
234: PetscBool create=PETSC_FALSE;
235: EPS_EIG_MATSHELL *matctx;
236: #if defined(PETSC_USE_COMPLEX)
237: PetscScalar theta,*aa,*bb;
238: const PetscScalar *ss;
239: PetscInt i,j,f,c,off,ld,lds;
240: IS perm;
241: #endif
243: MatGetSize(S,&n,NULL);
244: if (!*Op) create=PETSC_TRUE;
245: else {
246: MatGetSize(*Op,&m,NULL);
247: if (m!=n*n) create=PETSC_TRUE;
248: }
249: if (create) {
250: MatDestroy(Op);
251: VecDestroy(v0);
252: PetscNew(&matctx);
253: #if defined(PETSC_USE_COMPLEX)
254: MatCreateSeqDense(PETSC_COMM_SELF,n*n,n*n,NULL,&matctx->A);
255: MatCreateSeqDense(PETSC_COMM_SELF,n*n,n*n,NULL,&matctx->B);
256: MatCreateVecs(matctx->A,NULL,&matctx->w);
257: #else
258: MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&matctx->S);
259: #endif
260: MatCreateShell(PETSC_COMM_SELF,n*n,n*n,PETSC_DETERMINE,PETSC_DETERMINE,matctx,Op);
261: MatShellSetOperation(*Op,MATOP_MULT,(void(*)(void))MatMult_EigOperator);
262: MatShellSetOperation(*Op,MATOP_DESTROY,(void(*)(void))MatDestroy_EigOperator);
263: MatCreateVecs(*Op,NULL,v0);
264: } else {
265: MatShellGetContext(*Op,&matctx);
266: #if defined(PETSC_USE_COMPLEX)
267: MatZeroEntries(matctx->A);
268: #endif
269: }
270: #if defined(PETSC_USE_COMPLEX)
271: MatDenseGetArray(matctx->A,&aa);
272: MatDenseGetArray(matctx->B,&bb);
273: MatDenseGetArrayRead(S,&ss);
274: MatDenseGetLDA(S,&lds);
275: ld = n*n;
276: for (f=0;f<n;f++) {
277: off = f*n+f*n*ld;
278: for (i=0;i<n;i++) for (j=0;j<n;j++) aa[off+i+j*ld] = ss[i+j*lds];
279: for (c=0;c<n;c++) {
280: off = f*n+c*n*ld;
281: theta = ss[f+c*lds];
282: for (i=0;i<n;i++) aa[off+i+i*ld] += theta;
283: for (i=0;i<n;i++) for (j=0;j<n;j++) bb[off+i+j*ld] = -2*theta*ss[i+j*lds];
284: }
285: }
286: MatDenseRestoreArray(matctx->A,&aa);
287: MatDenseRestoreArray(matctx->B,&bb);
288: MatDenseRestoreArrayRead(S,&ss);
289: ISCreateStride(PETSC_COMM_SELF,n*n,0,1,&perm);
290: MatDestroy(&matctx->F);
291: MatDuplicate(matctx->A,MAT_COPY_VALUES,&matctx->F);
292: MatLUFactor(matctx->F,perm,perm,NULL);
293: ISDestroy(&perm);
294: #else
295: MatCopy(S,matctx->S,SAME_NONZERO_PATTERN);
296: #endif
297: matctx->lme = lme;
298: matctx->n = n;
299: VecSet(*v0,1.0);
300: return 0;
301: }
303: PetscErrorCode EPSSolve_LyapII(EPS eps)
304: {
305: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
306: PetscInt i,ldds,rk,nloc,mloc,nv,idx,k;
307: Vec v,w,z=eps->work[0],v0=NULL;
308: Mat S,C,Ux[2],Y,Y1,R,U,W,X,Op=NULL;
309: BV V;
310: BVOrthogType type;
311: BVOrthogRefineType refine;
312: PetscScalar eigr[2],eigi[2],*array,er,ei,*uu,*s,*xx,*aa,pM[4],vec[4];
313: PetscReal eta;
314: EPS epsrr;
315: PetscReal norm;
316: EPS_LYAPII_MATSHELL *matctx;
318: DSGetLeadingDimension(ctx->ds,&ldds);
320: /* Operator for the Lyapunov equation */
321: PetscNew(&matctx);
322: STGetOperator(eps->st,&matctx->S);
323: MatGetLocalSize(matctx->S,&mloc,&nloc);
324: MatCreateShell(PetscObjectComm((PetscObject)eps),mloc,nloc,PETSC_DETERMINE,PETSC_DETERMINE,matctx,&S);
325: matctx->Q = eps->V;
326: MatShellSetOperation(S,MATOP_MULT,(void(*)(void))MatMult_EPSLyapIIOperator);
327: MatShellSetOperation(S,MATOP_DESTROY,(void(*)(void))MatDestroy_EPSLyapIIOperator);
328: LMESetCoefficients(ctx->lme,S,NULL,NULL,NULL);
330: /* Right-hand side */
331: BVDuplicateResize(eps->V,ctx->rkl,&V);
332: BVGetOrthogonalization(V,&type,&refine,&eta,NULL);
333: BVSetOrthogonalization(V,type,refine,eta,BV_ORTHOG_BLOCK_TSQR);
334: MatCreateDense(PetscObjectComm((PetscObject)eps),eps->nloc,PETSC_DECIDE,PETSC_DECIDE,1,NULL,&Ux[0]);
335: MatCreateDense(PetscObjectComm((PetscObject)eps),eps->nloc,PETSC_DECIDE,PETSC_DECIDE,2,NULL,&Ux[1]);
336: nv = ctx->rkl;
337: PetscMalloc1(nv,&s);
339: /* Initialize first column */
340: EPSGetStartVector(eps,0,NULL);
341: BVGetColumn(eps->V,0,&v);
342: BVInsertVec(V,0,v);
343: BVRestoreColumn(eps->V,0,&v);
344: BVSetActiveColumns(eps->V,0,0); /* no deflation at the beginning */
345: LyapIIBuildRHS(S,1,Ux[0],V,eps->work);
346: idx = 0;
348: /* EPS for rank reduction */
349: EPSCreate(PETSC_COMM_SELF,&epsrr);
350: EPSSetOptionsPrefix(epsrr,((PetscObject)eps)->prefix);
351: EPSAppendOptionsPrefix(epsrr,"eps_lyapii_");
352: EPSSetDimensions(epsrr,1,PETSC_DEFAULT,PETSC_DEFAULT);
353: EPSSetTolerances(epsrr,PETSC_MACHINE_EPSILON*100,PETSC_DEFAULT);
355: while (eps->reason == EPS_CONVERGED_ITERATING) {
356: eps->its++;
358: /* Matrix for placing the solution of the Lyapunov equation (an alias of V) */
359: BVSetActiveColumns(V,0,nv);
360: BVGetMat(V,&Y1);
361: MatZeroEntries(Y1);
362: MatCreateLRC(NULL,Y1,NULL,NULL,&Y);
363: LMESetSolution(ctx->lme,Y);
365: /* Solve the Lyapunov equation SY + YS' = -2*S*Z*S' */
366: MatCreateLRC(NULL,Ux[idx],NULL,NULL,&C);
367: LMESetRHS(ctx->lme,C);
368: MatDestroy(&C);
369: LMESolve(ctx->lme);
370: BVRestoreMat(V,&Y1);
371: MatDestroy(&Y);
373: /* SVD of the solution: [Q,R]=qr(V); [U,Sigma,~]=svd(R) */
374: DSSetDimensions(ctx->ds,nv,0,0);
375: DSSVDSetDimensions(ctx->ds,nv);
376: DSGetMat(ctx->ds,DS_MAT_A,&R);
377: BVOrthogonalize(V,R);
378: DSRestoreMat(ctx->ds,DS_MAT_A,&R);
379: DSSetState(ctx->ds,DS_STATE_RAW);
380: DSSolve(ctx->ds,s,NULL);
382: /* Determine rank */
383: rk = nv;
384: for (i=1;i<nv;i++) if (PetscAbsScalar(s[i]/s[0])<PETSC_SQRT_MACHINE_EPSILON) {rk=i; break;}
385: PetscInfo(eps,"The computed solution of the Lyapunov equation has rank %" PetscInt_FMT "\n",rk);
386: rk = PetscMin(rk,ctx->rkc);
387: DSGetMat(ctx->ds,DS_MAT_U,&U);
388: BVMultInPlace(V,U,0,rk);
389: DSRestoreMat(ctx->ds,DS_MAT_U,&U);
390: BVSetActiveColumns(V,0,rk);
392: /* Rank reduction */
393: DSSetDimensions(ctx->ds,rk,0,0);
394: DSSVDSetDimensions(ctx->ds,rk);
395: DSGetMat(ctx->ds,DS_MAT_A,&W);
396: BVMatProject(V,S,V,W);
397: LyapIIBuildEigenMat(ctx->lme,W,&Op,&v0); /* Op=A\B, A=kron(I,S)+kron(S,I), B=-2*kron(S,S) */
398: DSRestoreMat(ctx->ds,DS_MAT_A,&W);
399: EPSSetOperators(epsrr,Op,NULL);
400: EPSSetInitialSpace(epsrr,1,&v0);
401: EPSSolve(epsrr);
402: EPSComputeVectors(epsrr);
403: /* Copy first eigenvector, vec(A)=x */
404: BVGetArray(epsrr->V,&xx);
405: DSGetArray(ctx->ds,DS_MAT_A,&aa);
406: for (i=0;i<rk;i++) PetscArraycpy(aa+i*ldds,xx+i*rk,rk);
407: DSRestoreArray(ctx->ds,DS_MAT_A,&aa);
408: BVRestoreArray(epsrr->V,&xx);
409: DSSetState(ctx->ds,DS_STATE_RAW);
410: /* Compute [U,Sigma,~] = svd(A), its rank should be 1 or 2 */
411: DSSolve(ctx->ds,s,NULL);
412: if (PetscAbsScalar(s[1]/s[0])<PETSC_SQRT_MACHINE_EPSILON) rk=1;
413: else rk = 2;
414: PetscInfo(eps,"The eigenvector has rank %" PetscInt_FMT "\n",rk);
415: DSGetMat(ctx->ds,DS_MAT_U,&U);
416: BVMultInPlace(V,U,0,rk);
417: DSRestoreMat(ctx->ds,DS_MAT_U,&U);
419: /* Save V in Ux */
420: idx = (rk==2)?1:0;
421: for (i=0;i<rk;i++) {
422: BVGetColumn(V,i,&v);
423: VecGetArray(v,&uu);
424: MatDenseGetColumn(Ux[idx],i,&array);
425: PetscArraycpy(array,uu,eps->nloc);
426: MatDenseRestoreColumn(Ux[idx],&array);
427: VecRestoreArray(v,&uu);
428: BVRestoreColumn(V,i,&v);
429: }
431: /* Eigenpair approximation */
432: BVGetColumn(V,0,&v);
433: MatMult(S,v,z);
434: VecDot(z,v,pM);
435: BVRestoreColumn(V,0,&v);
436: if (rk>1) {
437: BVGetColumn(V,1,&w);
438: VecDot(z,w,pM+1);
439: MatMult(S,w,z);
440: VecDot(z,w,pM+3);
441: BVGetColumn(V,0,&v);
442: VecDot(z,v,pM+2);
443: BVRestoreColumn(V,0,&v);
444: BVRestoreColumn(V,1,&w);
445: EV2x2(pM,2,eigr,eigi,vec);
446: MatCreateSeqDense(PETSC_COMM_SELF,2,2,vec,&X);
447: BVSetActiveColumns(V,0,rk);
448: BVMultInPlace(V,X,0,rk);
449: MatDestroy(&X);
450: #if !defined(PETSC_USE_COMPLEX)
451: norm = eigr[0]*eigr[0]+eigi[0]*eigi[0];
452: er = eigr[0]/norm; ei = -eigi[0]/norm;
453: #else
454: er =1.0/eigr[0]; ei = 0.0;
455: #endif
456: } else {
457: eigr[0] = pM[0]; eigi[0] = 0.0;
458: er = 1.0/eigr[0]; ei = 0.0;
459: }
460: BVGetColumn(V,0,&v);
461: if (eigi[0]!=0.0) BVGetColumn(V,1,&w);
462: else w = NULL;
463: eps->eigr[eps->nconv] = eigr[0]; eps->eigi[eps->nconv] = eigi[0];
464: EPSComputeResidualNorm_Private(eps,PETSC_FALSE,er,ei,v,w,eps->work,&norm);
465: BVRestoreColumn(V,0,&v);
466: if (w) BVRestoreColumn(V,1,&w);
467: (*eps->converged)(eps,er,ei,norm,&eps->errest[eps->nconv],eps->convergedctx);
468: k = 0;
469: if (eps->errest[eps->nconv]<eps->tol) {
470: k++;
471: if (rk==2) {
472: #if !defined (PETSC_USE_COMPLEX)
473: eps->eigr[eps->nconv+k] = eigr[0]; eps->eigi[eps->nconv+k] = -eigi[0];
474: #else
475: eps->eigr[eps->nconv+k] = PetscConj(eps->eigr[eps->nconv]);
476: #endif
477: k++;
478: }
479: /* Store converged eigenpairs and vectors for deflation */
480: for (i=0;i<k;i++) {
481: BVGetColumn(V,i,&v);
482: BVInsertVec(eps->V,eps->nconv+i,v);
483: BVRestoreColumn(V,i,&v);
484: }
485: eps->nconv += k;
486: BVSetActiveColumns(eps->V,eps->nconv-rk,eps->nconv);
487: BVOrthogonalize(eps->V,NULL);
488: DSSetDimensions(eps->ds,eps->nconv,0,0);
489: DSGetMat(eps->ds,DS_MAT_A,&W);
490: BVMatProject(eps->V,matctx->S,eps->V,W);
491: DSRestoreMat(eps->ds,DS_MAT_A,&W);
492: if (eps->nconv<eps->nev) {
493: idx = 0;
494: BVSetRandomColumn(V,0);
495: BVNormColumn(V,0,NORM_2,&norm);
496: BVScaleColumn(V,0,1.0/norm);
497: LyapIIBuildRHS(S,1,Ux[idx],V,eps->work);
498: }
499: } else {
500: /* Prepare right-hand side */
501: LyapIIBuildRHS(S,rk,Ux[idx],NULL,eps->work);
502: }
503: (*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx);
504: EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,eps->nconv+1);
505: }
506: STRestoreOperator(eps->st,&matctx->S);
507: MatDestroy(&S);
508: MatDestroy(&Ux[0]);
509: MatDestroy(&Ux[1]);
510: MatDestroy(&Op);
511: VecDestroy(&v0);
512: BVDestroy(&V);
513: EPSDestroy(&epsrr);
514: PetscFree(s);
515: return 0;
516: }
518: PetscErrorCode EPSSetFromOptions_LyapII(EPS eps,PetscOptionItems *PetscOptionsObject)
519: {
520: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
521: PetscInt k,array[2]={PETSC_DEFAULT,PETSC_DEFAULT};
522: PetscBool flg;
524: PetscOptionsHeadBegin(PetscOptionsObject,"EPS Lyapunov Inverse Iteration Options");
526: k = 2;
527: PetscOptionsIntArray("-eps_lyapii_ranks","Ranks for Lyapunov equation (one or two comma-separated integers)","EPSLyapIISetRanks",array,&k,&flg);
528: if (flg) EPSLyapIISetRanks(eps,array[0],array[1]);
530: PetscOptionsHeadEnd();
532: if (!ctx->lme) EPSLyapIIGetLME(eps,&ctx->lme);
533: LMESetFromOptions(ctx->lme);
534: return 0;
535: }
537: static PetscErrorCode EPSLyapIISetRanks_LyapII(EPS eps,PetscInt rkc,PetscInt rkl)
538: {
539: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
541: if (rkc==PETSC_DEFAULT) rkc = 10;
543: if (rkl==PETSC_DEFAULT) rkl = 3*rkc;
545: if (rkc != ctx->rkc) {
546: ctx->rkc = rkc;
547: eps->state = EPS_STATE_INITIAL;
548: }
549: if (rkl != ctx->rkl) {
550: ctx->rkl = rkl;
551: eps->state = EPS_STATE_INITIAL;
552: }
553: return 0;
554: }
556: /*@
557: EPSLyapIISetRanks - Set the ranks used in the solution of the Lyapunov equation.
559: Logically Collective on EPS
561: Input Parameters:
562: + eps - the eigenproblem solver context
563: . rkc - the compressed rank
564: - rkl - the Lyapunov rank
566: Options Database Key:
567: . -eps_lyapii_ranks <rkc,rkl> - Sets the rank parameters
569: Notes:
570: Lyapunov inverse iteration needs to solve a large-scale Lyapunov equation
571: at each iteration of the eigensolver. For this, an iterative solver (LME)
572: is used, which requires to prescribe the rank of the solution matrix X. This
573: is the meaning of parameter rkl. Later, this matrix is compressed into
574: another matrix of rank rkc. If not provided, rkl is a small multiple of rkc.
576: Level: intermediate
578: .seealso: EPSLyapIIGetRanks()
579: @*/
580: PetscErrorCode EPSLyapIISetRanks(EPS eps,PetscInt rkc,PetscInt rkl)
581: {
585: PetscTryMethod(eps,"EPSLyapIISetRanks_C",(EPS,PetscInt,PetscInt),(eps,rkc,rkl));
586: return 0;
587: }
589: static PetscErrorCode EPSLyapIIGetRanks_LyapII(EPS eps,PetscInt *rkc,PetscInt *rkl)
590: {
591: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
593: if (rkc) *rkc = ctx->rkc;
594: if (rkl) *rkl = ctx->rkl;
595: return 0;
596: }
598: /*@
599: EPSLyapIIGetRanks - Return the rank values used for the Lyapunov step.
601: Not Collective
603: Input Parameter:
604: . eps - the eigenproblem solver context
606: Output Parameters:
607: + rkc - the compressed rank
608: - rkl - the Lyapunov rank
610: Level: intermediate
612: .seealso: EPSLyapIISetRanks()
613: @*/
614: PetscErrorCode EPSLyapIIGetRanks(EPS eps,PetscInt *rkc,PetscInt *rkl)
615: {
617: PetscUseMethod(eps,"EPSLyapIIGetRanks_C",(EPS,PetscInt*,PetscInt*),(eps,rkc,rkl));
618: return 0;
619: }
621: static PetscErrorCode EPSLyapIISetLME_LyapII(EPS eps,LME lme)
622: {
623: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
625: PetscObjectReference((PetscObject)lme);
626: LMEDestroy(&ctx->lme);
627: ctx->lme = lme;
628: eps->state = EPS_STATE_INITIAL;
629: return 0;
630: }
632: /*@
633: EPSLyapIISetLME - Associate a linear matrix equation solver object (LME) to the
634: eigenvalue solver.
636: Collective on EPS
638: Input Parameters:
639: + eps - the eigenproblem solver context
640: - lme - the linear matrix equation solver object
642: Level: advanced
644: .seealso: EPSLyapIIGetLME()
645: @*/
646: PetscErrorCode EPSLyapIISetLME(EPS eps,LME lme)
647: {
651: PetscTryMethod(eps,"EPSLyapIISetLME_C",(EPS,LME),(eps,lme));
652: return 0;
653: }
655: static PetscErrorCode EPSLyapIIGetLME_LyapII(EPS eps,LME *lme)
656: {
657: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
659: if (!ctx->lme) {
660: LMECreate(PetscObjectComm((PetscObject)eps),&ctx->lme);
661: LMESetOptionsPrefix(ctx->lme,((PetscObject)eps)->prefix);
662: LMEAppendOptionsPrefix(ctx->lme,"eps_lyapii_");
663: PetscObjectIncrementTabLevel((PetscObject)ctx->lme,(PetscObject)eps,1);
664: }
665: *lme = ctx->lme;
666: return 0;
667: }
669: /*@
670: EPSLyapIIGetLME - Retrieve the linear matrix equation solver object (LME)
671: associated with the eigenvalue solver.
673: Not Collective
675: Input Parameter:
676: . eps - the eigenproblem solver context
678: Output Parameter:
679: . lme - the linear matrix equation solver object
681: Level: advanced
683: .seealso: EPSLyapIISetLME()
684: @*/
685: PetscErrorCode EPSLyapIIGetLME(EPS eps,LME *lme)
686: {
689: PetscUseMethod(eps,"EPSLyapIIGetLME_C",(EPS,LME*),(eps,lme));
690: return 0;
691: }
693: PetscErrorCode EPSView_LyapII(EPS eps,PetscViewer viewer)
694: {
695: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
696: PetscBool isascii;
698: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
699: if (isascii) {
700: PetscViewerASCIIPrintf(viewer," ranks: for Lyapunov solver=%" PetscInt_FMT ", after compression=%" PetscInt_FMT "\n",ctx->rkl,ctx->rkc);
701: if (!ctx->lme) EPSLyapIIGetLME(eps,&ctx->lme);
702: PetscViewerASCIIPushTab(viewer);
703: LMEView(ctx->lme,viewer);
704: PetscViewerASCIIPopTab(viewer);
705: }
706: return 0;
707: }
709: PetscErrorCode EPSReset_LyapII(EPS eps)
710: {
711: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
713: if (!ctx->lme) LMEReset(ctx->lme);
714: return 0;
715: }
717: PetscErrorCode EPSDestroy_LyapII(EPS eps)
718: {
719: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
721: LMEDestroy(&ctx->lme);
722: DSDestroy(&ctx->ds);
723: PetscFree(eps->data);
724: PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetLME_C",NULL);
725: PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetLME_C",NULL);
726: PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetRanks_C",NULL);
727: PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetRanks_C",NULL);
728: return 0;
729: }
731: PetscErrorCode EPSSetDefaultST_LyapII(EPS eps)
732: {
733: if (!((PetscObject)eps->st)->type_name) STSetType(eps->st,STSINVERT);
734: return 0;
735: }
737: SLEPC_EXTERN PetscErrorCode EPSCreate_LyapII(EPS eps)
738: {
739: EPS_LYAPII *ctx;
741: PetscNew(&ctx);
742: eps->data = (void*)ctx;
744: eps->useds = PETSC_TRUE;
746: eps->ops->solve = EPSSolve_LyapII;
747: eps->ops->setup = EPSSetUp_LyapII;
748: eps->ops->setupsort = EPSSetUpSort_Default;
749: eps->ops->setfromoptions = EPSSetFromOptions_LyapII;
750: eps->ops->reset = EPSReset_LyapII;
751: eps->ops->destroy = EPSDestroy_LyapII;
752: eps->ops->view = EPSView_LyapII;
753: eps->ops->setdefaultst = EPSSetDefaultST_LyapII;
754: eps->ops->backtransform = EPSBackTransform_Default;
755: eps->ops->computevectors = EPSComputeVectors_Schur;
757: PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetLME_C",EPSLyapIISetLME_LyapII);
758: PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetLME_C",EPSLyapIIGetLME_LyapII);
759: PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetRanks_C",EPSLyapIISetRanks_LyapII);
760: PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetRanks_C",EPSLyapIIGetRanks_LyapII);
761: return 0;
762: }