Actual source code: lyapii.c
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 PetscDefined(USE_COMPLEX)
50: Mat A,B,F;
51: Vec w;
52: #endif
53: } EPS_EIG_MATSHELL;
55: static PetscErrorCode EPSSetUp_LyapII(EPS eps)
56: {
57: PetscRandom rand;
58: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
60: PetscFunctionBegin;
61: EPSCheckSinvert(eps);
62: EPSCheckNotStructured(eps);
63: if (eps->nev==0) eps->nev = 1;
64: if (eps->ncv!=PETSC_DETERMINE) {
65: PetscCheck(eps->ncv>=eps->nev+1,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nev+1");
66: } else eps->ncv = eps->nev+1;
67: if (eps->mpd!=PETSC_DETERMINE) PetscCall(PetscInfo(eps,"Warning: parameter mpd ignored\n"));
68: if (eps->max_it==PETSC_DETERMINE) eps->max_it = PetscMax(1000*eps->nev,100*eps->n);
69: if (!eps->which) eps->which=EPS_LARGEST_REAL;
70: PetscCheck(eps->which==EPS_LARGEST_REAL,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver supports only largest real eigenvalues");
71: EPSCheckUnsupported(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION | EPS_FEATURE_THRESHOLD | EPS_FEATURE_TWOSIDED);
73: if (!ctx->rkc) ctx->rkc = 10;
74: if (!ctx->rkl) ctx->rkl = 3*ctx->rkc;
75: if (!ctx->lme) PetscCall(EPSLyapIIGetLME(eps,&ctx->lme));
76: PetscCall(LMESetProblemType(ctx->lme,LME_LYAPUNOV));
77: PetscCall(LMESetErrorIfNotConverged(ctx->lme,PETSC_TRUE));
79: if (!ctx->ds) {
80: PetscCall(DSCreate(PetscObjectComm((PetscObject)eps),&ctx->ds));
81: PetscCall(DSSetType(ctx->ds,DSSVD));
82: }
83: PetscCall(DSAllocate(ctx->ds,ctx->rkl));
85: PetscCall(DSSetType(eps->ds,DSNHEP));
86: PetscCall(DSAllocate(eps->ds,eps->ncv));
88: PetscCall(EPSAllocateSolution(eps,0));
89: PetscCall(BVGetRandomContext(eps->V,&rand)); /* make sure the random context is available when duplicating */
90: PetscCall(EPSSetWorkVecs(eps,3));
91: PetscFunctionReturn(PETSC_SUCCESS);
92: }
94: static PetscErrorCode MatMult_EPSLyapIIOperator(Mat M,Vec x,Vec r)
95: {
96: EPS_LYAPII_MATSHELL *matctx;
98: PetscFunctionBegin;
99: PetscCall(MatShellGetContext(M,&matctx));
100: PetscCall(MatMult(matctx->S,x,r));
101: PetscCall(BVOrthogonalizeVec(matctx->Q,r,NULL,NULL,NULL));
102: PetscFunctionReturn(PETSC_SUCCESS);
103: }
105: static PetscErrorCode MatDestroy_EPSLyapIIOperator(Mat M)
106: {
107: EPS_LYAPII_MATSHELL *matctx;
109: PetscFunctionBegin;
110: PetscCall(MatShellGetContext(M,&matctx));
111: PetscCall(MatDestroy(&matctx->S));
112: PetscCall(PetscFree(matctx));
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: static PetscErrorCode MatMult_EigOperator(Mat M,Vec x,Vec y)
117: {
118: EPS_EIG_MATSHELL *matctx;
119: #if !PetscDefined(USE_COMPLEX)
120: PetscInt n,lds;
121: PetscScalar *Y,*C,zero=0.0,done=1.0,dtwo=2.0;
122: const PetscScalar *S,*X;
123: PetscBLASInt n_,lds_;
124: #endif
126: PetscFunctionBegin;
127: PetscCall(MatShellGetContext(M,&matctx));
129: #if PetscDefined(USE_COMPLEX)
130: PetscCall(MatMult(matctx->B,x,matctx->w));
131: PetscCall(MatSolve(matctx->F,matctx->w,y));
132: #else
133: PetscCall(VecGetArrayRead(x,&X));
134: PetscCall(VecGetArray(y,&Y));
135: PetscCall(MatDenseGetArrayRead(matctx->S,&S));
136: PetscCall(MatDenseGetLDA(matctx->S,&lds));
138: n = matctx->n;
139: PetscCall(PetscCalloc1(n*n,&C));
140: PetscCall(PetscBLASIntCast(n,&n_));
141: PetscCall(PetscBLASIntCast(lds,&lds_));
143: /* C = 2*S*X*S.' */
144: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&dtwo,S,&lds_,X,&n_,&zero,Y,&n_));
145: PetscCallBLAS("BLASgemm",BLASgemm_("N","T",&n_,&n_,&n_,&done,Y,&n_,S,&lds_,&zero,C,&n_));
147: /* Solve S*Y + Y*S' = -C */
148: PetscCall(LMEDenseLyapunov(matctx->lme,n,(PetscScalar*)S,lds,C,n,Y,n));
150: PetscCall(PetscFree(C));
151: PetscCall(VecRestoreArrayRead(x,&X));
152: PetscCall(VecRestoreArray(y,&Y));
153: PetscCall(MatDenseRestoreArrayRead(matctx->S,&S));
154: #endif
155: PetscFunctionReturn(PETSC_SUCCESS);
156: }
158: static PetscErrorCode MatDestroy_EigOperator(Mat M)
159: {
160: EPS_EIG_MATSHELL *matctx;
162: PetscFunctionBegin;
163: PetscCall(MatShellGetContext(M,&matctx));
164: #if PetscDefined(USE_COMPLEX)
165: PetscCall(MatDestroy(&matctx->A));
166: PetscCall(MatDestroy(&matctx->B));
167: PetscCall(MatDestroy(&matctx->F));
168: PetscCall(VecDestroy(&matctx->w));
169: #else
170: PetscCall(MatDestroy(&matctx->S));
171: #endif
172: PetscCall(PetscFree(matctx));
173: PetscFunctionReturn(PETSC_SUCCESS);
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: {
181: PetscBLASInt lwork=10,ld_;
182: PetscScalar work[10];
183: PetscBLASInt two=2;
184: #if PetscDefined(USE_COMPLEX)
185: PetscReal rwork[6];
186: #endif
188: PetscFunctionBegin;
189: PetscCall(PetscBLASIntCast(ld,&ld_));
190: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
191: #if !PetscDefined(USE_COMPLEX)
192: PetscCallLAPACKInfo("LAPACKgeev",LAPACKgeev_("N","V",&two,M,&ld_,wr,wi,NULL,&ld_,vec,&ld_,work,&lwork,&info));
193: #else
194: PetscCallLAPACKInfo("LAPACKgeev",LAPACKgeev_("N","V",&two,M,&ld_,wr,NULL,&ld_,vec,&ld_,work,&lwork,rwork,&info));
195: #endif
196: PetscCall(PetscFPTrapPop());
197: PetscFunctionReturn(PETSC_SUCCESS);
198: }
200: /*
201: LyapIIBuildRHS: prepare the right-hand side of the Lyapunov equation SY + YS' = -2*S*Z*S'
202: in factored form:
203: if (V) U=sqrt(2)*S*V (uses 1 work vector)
204: else U=sqrt(2)*S*U (uses 2 work vectors)
205: where U,V are assumed to have rk columns.
206: */
207: static PetscErrorCode LyapIIBuildRHS(Mat S,PetscInt rk,Mat U,BV V,Vec *work)
208: {
209: PetscScalar *array,*uu;
210: PetscInt i,nloc;
211: Vec v,u=work[0];
213: PetscFunctionBegin;
214: PetscCall(MatGetLocalSize(U,&nloc,NULL));
215: for (i=0;i<rk;i++) {
216: PetscCall(MatDenseGetColumn(U,i,&array));
217: if (V) PetscCall(BVGetColumn(V,i,&v));
218: else {
219: v = work[1];
220: PetscCall(VecPlaceArray(v,array));
221: }
222: PetscCall(MatMult(S,v,u));
223: if (V) PetscCall(BVRestoreColumn(V,i,&v));
224: else PetscCall(VecResetArray(v));
225: PetscCall(VecScale(u,PETSC_SQRT2));
226: PetscCall(VecGetArray(u,&uu));
227: PetscCall(PetscArraycpy(array,uu,nloc));
228: PetscCall(VecRestoreArray(u,&uu));
229: PetscCall(MatDenseRestoreColumn(U,&array));
230: }
231: PetscFunctionReturn(PETSC_SUCCESS);
232: }
234: /*
235: LyapIIBuildEigenMat: create shell matrix Op=A\B with A = kron(I,S)+kron(S,I), B = -2*kron(S,S)
236: where S is a sequential square dense matrix of order n.
237: v0 is the initial vector, should have the form v0 = w*w' (for instance 1*1')
238: */
239: static PetscErrorCode LyapIIBuildEigenMat(LME lme,Mat S,Mat *Op,Vec *v0)
240: {
241: PetscInt n,m;
242: PetscBool create=PETSC_FALSE;
243: EPS_EIG_MATSHELL *matctx;
244: #if PetscDefined(USE_COMPLEX)
245: PetscScalar theta,*aa,*bb;
246: const PetscScalar *ss;
247: PetscInt i,j,f,c,off,ld,lds;
248: IS perm;
249: #endif
251: PetscFunctionBegin;
252: PetscCall(MatGetSize(S,&n,NULL));
253: if (!*Op) create=PETSC_TRUE;
254: else {
255: PetscCall(MatGetSize(*Op,&m,NULL));
256: if (m!=n*n) create=PETSC_TRUE;
257: }
258: if (create) {
259: PetscCall(MatDestroy(Op));
260: PetscCall(VecDestroy(v0));
261: PetscCall(PetscNew(&matctx));
262: #if PetscDefined(USE_COMPLEX)
263: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,n*n,n*n,NULL,&matctx->A));
264: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,n*n,n*n,NULL,&matctx->B));
265: PetscCall(MatCreateVecs(matctx->A,NULL,&matctx->w));
266: #else
267: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&matctx->S));
268: #endif
269: PetscCall(MatCreateShell(PETSC_COMM_SELF,n*n,n*n,PETSC_DETERMINE,PETSC_DETERMINE,matctx,Op));
270: PetscCall(MatShellSetOperation(*Op,MATOP_MULT,(PetscErrorCodeFn*)MatMult_EigOperator));
271: PetscCall(MatShellSetOperation(*Op,MATOP_DESTROY,(PetscErrorCodeFn*)MatDestroy_EigOperator));
272: PetscCall(MatCreateVecs(*Op,NULL,v0));
273: } else {
274: PetscCall(MatShellGetContext(*Op,&matctx));
275: #if PetscDefined(USE_COMPLEX)
276: PetscCall(MatZeroEntries(matctx->A));
277: #endif
278: }
279: #if PetscDefined(USE_COMPLEX)
280: PetscCall(MatDenseGetArray(matctx->A,&aa));
281: PetscCall(MatDenseGetArray(matctx->B,&bb));
282: PetscCall(MatDenseGetArrayRead(S,&ss));
283: PetscCall(MatDenseGetLDA(S,&lds));
284: ld = n*n;
285: for (f=0;f<n;f++) {
286: off = f*n+f*n*ld;
287: for (i=0;i<n;i++) for (j=0;j<n;j++) aa[off+i+j*ld] = ss[i+j*lds];
288: for (c=0;c<n;c++) {
289: off = f*n+c*n*ld;
290: theta = ss[f+c*lds];
291: for (i=0;i<n;i++) aa[off+i+i*ld] += theta;
292: for (i=0;i<n;i++) for (j=0;j<n;j++) bb[off+i+j*ld] = -2*theta*ss[i+j*lds];
293: }
294: }
295: PetscCall(MatDenseRestoreArray(matctx->A,&aa));
296: PetscCall(MatDenseRestoreArray(matctx->B,&bb));
297: PetscCall(MatDenseRestoreArrayRead(S,&ss));
298: PetscCall(ISCreateStride(PETSC_COMM_SELF,n*n,0,1,&perm));
299: PetscCall(MatDestroy(&matctx->F));
300: PetscCall(MatDuplicate(matctx->A,MAT_COPY_VALUES,&matctx->F));
301: PetscCall(MatLUFactor(matctx->F,perm,perm,NULL));
302: PetscCall(ISDestroy(&perm));
303: #else
304: PetscCall(MatCopy(S,matctx->S,SAME_NONZERO_PATTERN));
305: #endif
306: matctx->lme = lme;
307: matctx->n = n;
308: PetscCall(VecSet(*v0,1.0));
309: PetscFunctionReturn(PETSC_SUCCESS);
310: }
312: static PetscErrorCode EPSSolve_LyapII(EPS eps)
313: {
314: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
315: PetscInt i,ldds,rk,nloc,mloc,nv,idx,k;
316: Vec v,w,z=eps->work[0],v0=NULL;
317: VecType vtype;
318: Mat S,C,Ux[2],Y,Y1,R,U,W,X,Op=NULL;
319: BV V;
320: BVOrthogType type;
321: BVOrthogRefineType refine;
322: PetscScalar eigr[2],eigi[2],*array,er,ei,*uu,*s,*xx,*aa,pM[4],vec[4];
323: PetscReal eta;
324: EPS epsrr;
325: PetscReal norm;
326: EPS_LYAPII_MATSHELL *matctx;
328: PetscFunctionBegin;
329: PetscCall(DSGetLeadingDimension(ctx->ds,&ldds));
331: /* Operator for the Lyapunov equation */
332: PetscCall(PetscNew(&matctx));
333: PetscCall(STGetOperator(eps->st,&matctx->S));
334: PetscCall(MatGetLocalSize(matctx->S,&mloc,&nloc));
335: PetscCall(MatCreateShell(PetscObjectComm((PetscObject)eps),mloc,nloc,PETSC_DETERMINE,PETSC_DETERMINE,matctx,&S));
336: matctx->Q = eps->V;
337: PetscCall(MatShellSetOperation(S,MATOP_MULT,(PetscErrorCodeFn*)MatMult_EPSLyapIIOperator));
338: PetscCall(MatShellSetOperation(S,MATOP_DESTROY,(PetscErrorCodeFn*)MatDestroy_EPSLyapIIOperator));
339: /* make sure the shell matrix generates a vector of the same type as the problem matrices */
340: PetscCall(MatGetVecType(matctx->S,&vtype));
341: PetscCall(MatShellSetVecType(S,vtype));
342: PetscCall(LMESetCoefficients(ctx->lme,S,NULL,NULL,NULL));
344: /* Right-hand side */
345: PetscCall(BVDuplicateResize(eps->V,ctx->rkl,&V));
346: PetscCall(BVGetOrthogonalization(V,&type,&refine,&eta,NULL));
347: PetscCall(BVSetOrthogonalization(V,type,refine,eta,BV_ORTHOG_BLOCK_TSQR));
348: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)eps),eps->nloc,PETSC_DECIDE,PETSC_DECIDE,1,NULL,&Ux[0]));
349: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)eps),eps->nloc,PETSC_DECIDE,PETSC_DECIDE,2,NULL,&Ux[1]));
350: nv = ctx->rkl;
351: PetscCall(PetscMalloc1(nv,&s));
353: /* Initialize first column */
354: PetscCall(EPSGetStartVector(eps,0,NULL));
355: PetscCall(BVGetColumn(eps->V,0,&v));
356: PetscCall(BVInsertVec(V,0,v));
357: PetscCall(BVRestoreColumn(eps->V,0,&v));
358: PetscCall(BVSetActiveColumns(eps->V,0,0)); /* no deflation at the beginning */
359: PetscCall(LyapIIBuildRHS(S,1,Ux[0],V,eps->work));
360: idx = 0;
362: /* EPS for rank reduction */
363: PetscCall(EPSCreate(PETSC_COMM_SELF,&epsrr));
364: PetscCall(EPSSetOptionsPrefix(epsrr,((PetscObject)eps)->prefix));
365: PetscCall(EPSAppendOptionsPrefix(epsrr,"eps_lyapii_"));
366: PetscCall(EPSSetDimensions(epsrr,1,PETSC_CURRENT,PETSC_CURRENT));
367: PetscCall(EPSSetTolerances(epsrr,PETSC_MACHINE_EPSILON*100,PETSC_CURRENT));
369: while (eps->reason == EPS_CONVERGED_ITERATING) {
370: eps->its++;
372: /* Matrix for placing the solution of the Lyapunov equation (an alias of V) */
373: PetscCall(BVSetActiveColumns(V,0,nv));
374: PetscCall(BVGetMat(V,&Y1));
375: PetscCall(MatZeroEntries(Y1));
376: PetscCall(MatCreateLRC(NULL,Y1,NULL,NULL,&Y));
377: PetscCall(LMESetSolution(ctx->lme,Y));
379: /* Solve the Lyapunov equation SY + YS' = -2*S*Z*S' */
380: PetscCall(MatCreateLRC(NULL,Ux[idx],NULL,NULL,&C));
381: PetscCall(LMESetRHS(ctx->lme,C));
382: PetscCall(MatDestroy(&C));
383: PetscCall(LMESolve(ctx->lme));
384: PetscCall(BVRestoreMat(V,&Y1));
385: PetscCall(MatDestroy(&Y));
387: /* SVD of the solution: [Q,R]=qr(V); [U,Sigma,~]=svd(R) */
388: PetscCall(DSSetDimensions(ctx->ds,nv,0,0));
389: PetscCall(DSSVDSetDimensions(ctx->ds,nv));
390: PetscCall(DSGetMat(ctx->ds,DS_MAT_A,&R));
391: PetscCall(BVOrthogonalize(V,R));
392: PetscCall(DSRestoreMat(ctx->ds,DS_MAT_A,&R));
393: PetscCall(DSSetState(ctx->ds,DS_STATE_RAW));
394: PetscCall(DSSolve(ctx->ds,s,NULL));
396: /* Determine rank */
397: rk = nv;
398: for (i=1;i<nv;i++) if (PetscAbsScalar(s[i]/s[0])<PETSC_SQRT_MACHINE_EPSILON) {rk=i; break;}
399: PetscCall(PetscInfo(eps,"The computed solution of the Lyapunov equation has rank %" PetscInt_FMT "\n",rk));
400: rk = PetscMin(rk,ctx->rkc);
401: PetscCall(DSGetMat(ctx->ds,DS_MAT_U,&U));
402: PetscCall(BVMultInPlace(V,U,0,rk));
403: PetscCall(DSRestoreMat(ctx->ds,DS_MAT_U,&U));
404: PetscCall(BVSetActiveColumns(V,0,rk));
406: /* Rank reduction */
407: PetscCall(DSSetDimensions(ctx->ds,rk,0,0));
408: PetscCall(DSSVDSetDimensions(ctx->ds,rk));
409: PetscCall(DSGetMat(ctx->ds,DS_MAT_A,&W));
410: PetscCall(BVMatProject(V,S,V,W));
411: PetscCall(LyapIIBuildEigenMat(ctx->lme,W,&Op,&v0)); /* Op=A\B, A=kron(I,S)+kron(S,I), B=-2*kron(S,S) */
412: PetscCall(DSRestoreMat(ctx->ds,DS_MAT_A,&W));
413: PetscCall(EPSSetOperators(epsrr,Op,NULL));
414: PetscCall(EPSSetInitialSpace(epsrr,1,&v0));
415: PetscCall(EPSSolve(epsrr));
416: PetscCall(EPSComputeVectors(epsrr));
417: /* Copy first eigenvector, vec(A)=x */
418: PetscCall(BVGetArray(epsrr->V,&xx));
419: PetscCall(DSGetArray(ctx->ds,DS_MAT_A,&aa));
420: for (i=0;i<rk;i++) PetscCall(PetscArraycpy(aa+i*ldds,xx+i*rk,rk));
421: PetscCall(DSRestoreArray(ctx->ds,DS_MAT_A,&aa));
422: PetscCall(BVRestoreArray(epsrr->V,&xx));
423: PetscCall(DSSetState(ctx->ds,DS_STATE_RAW));
424: /* Compute [U,Sigma,~] = svd(A), its rank should be 1 or 2 */
425: PetscCall(DSSolve(ctx->ds,s,NULL));
426: if (PetscAbsScalar(s[1]/s[0])<PETSC_SQRT_MACHINE_EPSILON) rk=1;
427: else rk = 2;
428: PetscCall(PetscInfo(eps,"The eigenvector has rank %" PetscInt_FMT "\n",rk));
429: PetscCall(DSGetMat(ctx->ds,DS_MAT_U,&U));
430: PetscCall(BVMultInPlace(V,U,0,rk));
431: PetscCall(DSRestoreMat(ctx->ds,DS_MAT_U,&U));
433: /* Save V in Ux */
434: idx = (rk==2)?1:0;
435: for (i=0;i<rk;i++) {
436: PetscCall(BVGetColumn(V,i,&v));
437: PetscCall(VecGetArray(v,&uu));
438: PetscCall(MatDenseGetColumn(Ux[idx],i,&array));
439: PetscCall(PetscArraycpy(array,uu,eps->nloc));
440: PetscCall(MatDenseRestoreColumn(Ux[idx],&array));
441: PetscCall(VecRestoreArray(v,&uu));
442: PetscCall(BVRestoreColumn(V,i,&v));
443: }
445: /* Eigenpair approximation */
446: PetscCall(BVGetColumn(V,0,&v));
447: PetscCall(MatMult(S,v,z));
448: PetscCall(VecDot(z,v,pM));
449: PetscCall(BVRestoreColumn(V,0,&v));
450: if (rk>1) {
451: PetscCall(BVGetColumn(V,1,&w));
452: PetscCall(VecDot(z,w,pM+1));
453: PetscCall(MatMult(S,w,z));
454: PetscCall(VecDot(z,w,pM+3));
455: PetscCall(BVGetColumn(V,0,&v));
456: PetscCall(VecDot(z,v,pM+2));
457: PetscCall(BVRestoreColumn(V,0,&v));
458: PetscCall(BVRestoreColumn(V,1,&w));
459: PetscCall(EV2x2(pM,2,eigr,eigi,vec));
460: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,2,2,vec,&X));
461: PetscCall(BVSetActiveColumns(V,0,rk));
462: PetscCall(BVMultInPlace(V,X,0,rk));
463: PetscCall(MatDestroy(&X));
464: #if !PetscDefined(USE_COMPLEX)
465: norm = eigr[0]*eigr[0]+eigi[0]*eigi[0];
466: er = eigr[0]/norm; ei = -eigi[0]/norm;
467: #else
468: er =1.0/eigr[0]; ei = 0.0;
469: #endif
470: } else {
471: eigr[0] = pM[0]; eigi[0] = 0.0;
472: er = 1.0/eigr[0]; ei = 0.0;
473: }
474: PetscCall(BVGetColumn(V,0,&v));
475: if (eigi[0]!=0.0) PetscCall(BVGetColumn(V,1,&w));
476: else w = NULL;
477: eps->eigr[eps->nconv] = eigr[0]; eps->eigi[eps->nconv] = eigi[0];
478: PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_FALSE,er,ei,v,w,eps->work,&norm));
479: PetscCall(BVRestoreColumn(V,0,&v));
480: if (w) PetscCall(BVRestoreColumn(V,1,&w));
481: PetscCall((*eps->converged)(eps,er,ei,norm,&eps->errest[eps->nconv],eps->convergedctx));
482: k = 0;
483: if (eps->errest[eps->nconv]<eps->tol) {
484: k++;
485: if (rk==2) {
486: #if !PetscDefined(USE_COMPLEX)
487: eps->eigr[eps->nconv+k] = eigr[0]; eps->eigi[eps->nconv+k] = -eigi[0];
488: #else
489: eps->eigr[eps->nconv+k] = PetscConj(eps->eigr[eps->nconv]);
490: #endif
491: k++;
492: }
493: /* Store converged eigenpairs and vectors for deflation */
494: for (i=0;i<k;i++) {
495: PetscCall(BVGetColumn(V,i,&v));
496: PetscCall(BVInsertVec(eps->V,eps->nconv+i,v));
497: PetscCall(BVRestoreColumn(V,i,&v));
498: }
499: eps->nconv += k;
500: PetscCall(BVSetActiveColumns(eps->V,eps->nconv-rk,eps->nconv));
501: PetscCall(BVOrthogonalize(eps->V,NULL));
502: PetscCall(DSSetDimensions(eps->ds,eps->nconv,0,0));
503: PetscCall(DSGetMat(eps->ds,DS_MAT_A,&W));
504: PetscCall(BVMatProject(eps->V,matctx->S,eps->V,W));
505: PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&W));
506: if (eps->nconv<eps->nev) {
507: idx = 0;
508: PetscCall(BVSetRandomColumn(V,0));
509: PetscCall(BVNormColumn(V,0,NORM_2,&norm));
510: PetscCall(BVScaleColumn(V,0,1.0/norm));
511: PetscCall(LyapIIBuildRHS(S,1,Ux[idx],V,eps->work));
512: }
513: } else {
514: /* Prepare right-hand side */
515: PetscCall(LyapIIBuildRHS(S,rk,Ux[idx],NULL,eps->work));
516: }
517: PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx));
518: PetscCall(EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,eps->nconv+1));
519: }
520: PetscCall(STRestoreOperator(eps->st,&matctx->S));
521: PetscCall(MatDestroy(&S));
522: PetscCall(MatDestroy(&Ux[0]));
523: PetscCall(MatDestroy(&Ux[1]));
524: PetscCall(MatDestroy(&Op));
525: PetscCall(VecDestroy(&v0));
526: PetscCall(BVDestroy(&V));
527: PetscCall(EPSDestroy(&epsrr));
528: PetscCall(PetscFree(s));
529: PetscFunctionReturn(PETSC_SUCCESS);
530: }
532: static PetscErrorCode EPSSetFromOptions_LyapII(EPS eps,PetscOptionItems PetscOptionsObject)
533: {
534: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
535: PetscInt k,array[2]={PETSC_DETERMINE,PETSC_DETERMINE};
536: PetscBool flg;
538: PetscFunctionBegin;
539: PetscOptionsHeadBegin(PetscOptionsObject,"EPS Lyapunov Inverse Iteration Options");
541: k = 2;
542: PetscCall(PetscOptionsIntArray("-eps_lyapii_ranks","Ranks for Lyapunov equation (one or two comma-separated integers)","EPSLyapIISetRanks",array,&k,&flg));
543: if (flg) PetscCall(EPSLyapIISetRanks(eps,array[0],array[1]));
545: PetscOptionsHeadEnd();
547: if (!ctx->lme) PetscCall(EPSLyapIIGetLME(eps,&ctx->lme));
548: PetscCall(LMESetFromOptions(ctx->lme));
549: PetscFunctionReturn(PETSC_SUCCESS);
550: }
552: static PetscErrorCode EPSLyapIISetRanks_LyapII(EPS eps,PetscInt rkc,PetscInt rkl)
553: {
554: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
556: PetscFunctionBegin;
557: if (rkc==PETSC_DETERMINE) {
558: if (ctx->rkc != 10) eps->state = EPS_STATE_INITIAL;
559: ctx->rkc = 10;
560: } else if (rkc!=PETSC_CURRENT) {
561: PetscCheck(rkc>1,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The compressed rank %" PetscInt_FMT " must be larger than 1",rkc);
562: if (ctx->rkc != rkc) eps->state = EPS_STATE_INITIAL;
563: ctx->rkc = rkc;
564: }
565: if (rkl==PETSC_DETERMINE) {
566: if (ctx->rkl != 3*rkc) eps->state = EPS_STATE_INITIAL;
567: ctx->rkl = 3*rkc;
568: } else if (rkl!=PETSC_CURRENT) {
569: PetscCheck(rkl>=rkc,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The Lyapunov rank %" PetscInt_FMT " cannot be smaller than the compressed rank %" PetscInt_FMT,rkl,rkc);
570: if (ctx->rkl != rkl) eps->state = EPS_STATE_INITIAL;
571: ctx->rkl = rkl;
572: }
573: PetscFunctionReturn(PETSC_SUCCESS);
574: }
576: /*@
577: EPSLyapIISetRanks - Set the ranks used in the solution of the Lyapunov equation.
579: Logically Collective
581: Input Parameters:
582: + eps - the linear eigensolver context
583: . rkc - the compressed rank
584: - rkl - the Lyapunov rank
586: Options Database Key:
587: . -eps_lyapii_ranks rkc,rkl - sets the rank parameters
589: Notes:
590: `PETSC_CURRENT` can be used to preserve the current value of any of the
591: arguments, and `PETSC_DETERMINE` to set them to a default value.
593: Lyapunov inverse iteration needs to solve a large-scale Lyapunov equation
594: at each iteration of the eigensolver. For this, an iterative solver (`LME`)
595: is used, which requires to prescribe the rank of the solution matrix $X$. This
596: is the meaning of parameter `rkl`. Later, this matrix is compressed into
597: another matrix of rank `rkc`. If not provided, `rkl` is a small multiple of `rkc`.
599: Level: intermediate
601: .seealso: [](ch:eps), `EPSLYAPII`, `EPSLyapIIGetRanks()`
602: @*/
603: PetscErrorCode EPSLyapIISetRanks(EPS eps,PetscInt rkc,PetscInt rkl)
604: {
605: PetscFunctionBegin;
609: PetscTryMethod(eps,"EPSLyapIISetRanks_C",(EPS,PetscInt,PetscInt),(eps,rkc,rkl));
610: PetscFunctionReturn(PETSC_SUCCESS);
611: }
613: static PetscErrorCode EPSLyapIIGetRanks_LyapII(EPS eps,PetscInt *rkc,PetscInt *rkl)
614: {
615: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
617: PetscFunctionBegin;
618: if (rkc) *rkc = ctx->rkc;
619: if (rkl) *rkl = ctx->rkl;
620: PetscFunctionReturn(PETSC_SUCCESS);
621: }
623: /*@
624: EPSLyapIIGetRanks - Return the rank values used for the Lyapunov step.
626: Not Collective
628: Input Parameter:
629: . eps - the linear eigensolver context
631: Output Parameters:
632: + rkc - the compressed rank
633: - rkl - the Lyapunov rank
635: Level: intermediate
637: .seealso: [](ch:eps), `EPSLYAPII`, `EPSLyapIISetRanks()`
638: @*/
639: PetscErrorCode EPSLyapIIGetRanks(EPS eps,PetscInt *rkc,PetscInt *rkl)
640: {
641: PetscFunctionBegin;
643: PetscUseMethod(eps,"EPSLyapIIGetRanks_C",(EPS,PetscInt*,PetscInt*),(eps,rkc,rkl));
644: PetscFunctionReturn(PETSC_SUCCESS);
645: }
647: static PetscErrorCode EPSLyapIISetLME_LyapII(EPS eps,LME lme)
648: {
649: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
651: PetscFunctionBegin;
652: PetscCall(PetscObjectReference((PetscObject)lme));
653: PetscCall(LMEDestroy(&ctx->lme));
654: ctx->lme = lme;
655: eps->state = EPS_STATE_INITIAL;
656: PetscFunctionReturn(PETSC_SUCCESS);
657: }
659: /*@
660: EPSLyapIISetLME - Associate a linear matrix equation solver object (`LME`) to the
661: eigenvalue solver.
663: Collective
665: Input Parameters:
666: + eps - the linear eigensolver context
667: - lme - the linear matrix equation solver context
669: Level: advanced
671: .seealso: [](ch:eps), `EPSLYAPII`, `EPSLyapIIGetLME()`
672: @*/
673: PetscErrorCode EPSLyapIISetLME(EPS eps,LME lme)
674: {
675: PetscFunctionBegin;
678: PetscCheckSameComm(eps,1,lme,2);
679: PetscTryMethod(eps,"EPSLyapIISetLME_C",(EPS,LME),(eps,lme));
680: PetscFunctionReturn(PETSC_SUCCESS);
681: }
683: static PetscErrorCode EPSLyapIIGetLME_LyapII(EPS eps,LME *lme)
684: {
685: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
687: PetscFunctionBegin;
688: if (!ctx->lme) {
689: PetscCall(LMECreate(PetscObjectComm((PetscObject)eps),&ctx->lme));
690: PetscCall(LMESetOptionsPrefix(ctx->lme,((PetscObject)eps)->prefix));
691: PetscCall(LMEAppendOptionsPrefix(ctx->lme,"eps_lyapii_"));
692: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->lme,(PetscObject)eps,1));
693: }
694: *lme = ctx->lme;
695: PetscFunctionReturn(PETSC_SUCCESS);
696: }
698: /*@
699: EPSLyapIIGetLME - Retrieve the linear matrix equation solver object (`LME`)
700: associated with the eigenvalue solver.
702: Not Collective
704: Input Parameter:
705: . eps - the linear eigensolver context
707: Output Parameter:
708: . lme - the linear matrix equation solver context
710: Level: advanced
712: .seealso: [](ch:eps), `EPSLYAPII`, `EPSLyapIISetLME()`
713: @*/
714: PetscErrorCode EPSLyapIIGetLME(EPS eps,LME *lme)
715: {
716: PetscFunctionBegin;
718: PetscAssertPointer(lme,2);
719: PetscUseMethod(eps,"EPSLyapIIGetLME_C",(EPS,LME*),(eps,lme));
720: PetscFunctionReturn(PETSC_SUCCESS);
721: }
723: static PetscErrorCode EPSView_LyapII(EPS eps,PetscViewer viewer)
724: {
725: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
726: PetscBool isascii;
728: PetscFunctionBegin;
729: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
730: if (isascii) {
731: PetscCall(PetscViewerASCIIPrintf(viewer," ranks: for Lyapunov solver=%" PetscInt_FMT ", after compression=%" PetscInt_FMT "\n",ctx->rkl,ctx->rkc));
732: if (!ctx->lme) PetscCall(EPSLyapIIGetLME(eps,&ctx->lme));
733: PetscCall(PetscViewerASCIIPushTab(viewer));
734: PetscCall(LMEView(ctx->lme,viewer));
735: PetscCall(PetscViewerASCIIPopTab(viewer));
736: }
737: PetscFunctionReturn(PETSC_SUCCESS);
738: }
740: static PetscErrorCode EPSReset_LyapII(EPS eps)
741: {
742: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
744: PetscFunctionBegin;
745: if (!ctx->lme) PetscCall(LMEReset(ctx->lme));
746: PetscFunctionReturn(PETSC_SUCCESS);
747: }
749: static PetscErrorCode EPSDestroy_LyapII(EPS eps)
750: {
751: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
753: PetscFunctionBegin;
754: PetscCall(LMEDestroy(&ctx->lme));
755: PetscCall(DSDestroy(&ctx->ds));
756: PetscCall(PetscFree(eps->data));
757: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetLME_C",NULL));
758: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetLME_C",NULL));
759: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetRanks_C",NULL));
760: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetRanks_C",NULL));
761: PetscFunctionReturn(PETSC_SUCCESS);
762: }
764: static PetscErrorCode EPSSetDefaultST_LyapII(EPS eps)
765: {
766: PetscFunctionBegin;
767: if (!((PetscObject)eps->st)->type_name) PetscCall(STSetType(eps->st,STSINVERT));
768: PetscFunctionReturn(PETSC_SUCCESS);
769: }
771: /*MC
772: EPSLYAPII - EPSLYAPII = "lyapii" - The Lyapunov inverse iteration.
774: Notes:
775: This solver implements the method of Lyapunov inverse iteration
776: {cite:p}`Mee10,Elm13` to compute rightmost eigenvalues of
777: non-Hermitian matrices (or matrix pencils).
779: At each step of the eigensolver, a Lyapunov equation must be solved,
780: and this is done with an `LME` object, see `EPSLyapIIGetLME()`.
782: This solver may be useful for analyzing the stability of PDE's.
783: Note that the method requires the input matrix to be stable, so
784: it generally requires to shift the matrix before passing it in
785: `EPSSetOperators()`.
787: Level: beginner
789: .seealso: [](ch:eps), `EPS`, `EPSType`, `EPSSetType()`, `EPSLyapIIGetLME()`, `EPSSetOperators()`
790: M*/
791: SLEPC_EXTERN PetscErrorCode EPSCreate_LyapII(EPS eps)
792: {
793: EPS_LYAPII *ctx;
795: PetscFunctionBegin;
796: PetscCall(PetscNew(&ctx));
797: eps->data = (void*)ctx;
799: eps->useds = PETSC_TRUE;
801: eps->ops->solve = EPSSolve_LyapII;
802: eps->ops->setup = EPSSetUp_LyapII;
803: eps->ops->setupsort = EPSSetUpSort_Default;
804: eps->ops->setfromoptions = EPSSetFromOptions_LyapII;
805: eps->ops->reset = EPSReset_LyapII;
806: eps->ops->destroy = EPSDestroy_LyapII;
807: eps->ops->view = EPSView_LyapII;
808: eps->ops->setdefaultst = EPSSetDefaultST_LyapII;
809: eps->ops->backtransform = EPSBackTransform_Default;
810: eps->ops->computevectors = EPSComputeVectors_Schur;
812: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetLME_C",EPSLyapIISetLME_LyapII));
813: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetLME_C",EPSLyapIIGetLME_LyapII));
814: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetRanks_C",EPSLyapIISetRanks_LyapII));
815: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetRanks_C",EPSLyapIIGetRanks_LyapII));
816: PetscFunctionReturn(PETSC_SUCCESS);
817: }