Actual source code: lyapii.c
slepc-main 2024-12-17
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: 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 !defined(PETSC_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 defined(PETSC_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 defined(PETSC_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,info;
184: #if defined(PETSC_USE_COMPLEX)
185: PetscReal rwork[6];
186: #endif
188: PetscFunctionBegin;
189: PetscCall(PetscBLASIntCast(ld,&ld_));
190: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
191: #if !defined(PETSC_USE_COMPLEX)
192: PetscCallBLAS("LAPACKgeev",LAPACKgeev_("N","V",&two,M,&ld_,wr,wi,NULL,&ld_,vec,&ld_,work,&lwork,&info));
193: #else
194: PetscCallBLAS("LAPACKgeev",LAPACKgeev_("N","V",&two,M,&ld_,wr,NULL,&ld_,vec,&ld_,work,&lwork,rwork,&info));
195: #endif
196: SlepcCheckLapackInfo("geev",info);
197: PetscCall(PetscFPTrapPop());
198: PetscFunctionReturn(PETSC_SUCCESS);
199: }
201: /*
202: LyapIIBuildRHS: prepare the right-hand side of the Lyapunov equation SY + YS' = -2*S*Z*S'
203: in factored form:
204: if (V) U=sqrt(2)*S*V (uses 1 work vector)
205: else U=sqrt(2)*S*U (uses 2 work vectors)
206: where U,V are assumed to have rk columns.
207: */
208: static PetscErrorCode LyapIIBuildRHS(Mat S,PetscInt rk,Mat U,BV V,Vec *work)
209: {
210: PetscScalar *array,*uu;
211: PetscInt i,nloc;
212: Vec v,u=work[0];
214: PetscFunctionBegin;
215: PetscCall(MatGetLocalSize(U,&nloc,NULL));
216: for (i=0;i<rk;i++) {
217: PetscCall(MatDenseGetColumn(U,i,&array));
218: if (V) PetscCall(BVGetColumn(V,i,&v));
219: else {
220: v = work[1];
221: PetscCall(VecPlaceArray(v,array));
222: }
223: PetscCall(MatMult(S,v,u));
224: if (V) PetscCall(BVRestoreColumn(V,i,&v));
225: else PetscCall(VecResetArray(v));
226: PetscCall(VecScale(u,PETSC_SQRT2));
227: PetscCall(VecGetArray(u,&uu));
228: PetscCall(PetscArraycpy(array,uu,nloc));
229: PetscCall(VecRestoreArray(u,&uu));
230: PetscCall(MatDenseRestoreColumn(U,&array));
231: }
232: PetscFunctionReturn(PETSC_SUCCESS);
233: }
235: /*
236: LyapIIBuildEigenMat: create shell matrix Op=A\B with A = kron(I,S)+kron(S,I), B = -2*kron(S,S)
237: where S is a sequential square dense matrix of order n.
238: v0 is the initial vector, should have the form v0 = w*w' (for instance 1*1')
239: */
240: static PetscErrorCode LyapIIBuildEigenMat(LME lme,Mat S,Mat *Op,Vec *v0)
241: {
242: PetscInt n,m;
243: PetscBool create=PETSC_FALSE;
244: EPS_EIG_MATSHELL *matctx;
245: #if defined(PETSC_USE_COMPLEX)
246: PetscScalar theta,*aa,*bb;
247: const PetscScalar *ss;
248: PetscInt i,j,f,c,off,ld,lds;
249: IS perm;
250: #endif
252: PetscFunctionBegin;
253: PetscCall(MatGetSize(S,&n,NULL));
254: if (!*Op) create=PETSC_TRUE;
255: else {
256: PetscCall(MatGetSize(*Op,&m,NULL));
257: if (m!=n*n) create=PETSC_TRUE;
258: }
259: if (create) {
260: PetscCall(MatDestroy(Op));
261: PetscCall(VecDestroy(v0));
262: PetscCall(PetscNew(&matctx));
263: #if defined(PETSC_USE_COMPLEX)
264: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,n*n,n*n,NULL,&matctx->A));
265: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,n*n,n*n,NULL,&matctx->B));
266: PetscCall(MatCreateVecs(matctx->A,NULL,&matctx->w));
267: #else
268: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&matctx->S));
269: #endif
270: PetscCall(MatCreateShell(PETSC_COMM_SELF,n*n,n*n,PETSC_DETERMINE,PETSC_DETERMINE,matctx,Op));
271: PetscCall(MatShellSetOperation(*Op,MATOP_MULT,(void(*)(void))MatMult_EigOperator));
272: PetscCall(MatShellSetOperation(*Op,MATOP_DESTROY,(void(*)(void))MatDestroy_EigOperator));
273: PetscCall(MatCreateVecs(*Op,NULL,v0));
274: } else {
275: PetscCall(MatShellGetContext(*Op,&matctx));
276: #if defined(PETSC_USE_COMPLEX)
277: PetscCall(MatZeroEntries(matctx->A));
278: #endif
279: }
280: #if defined(PETSC_USE_COMPLEX)
281: PetscCall(MatDenseGetArray(matctx->A,&aa));
282: PetscCall(MatDenseGetArray(matctx->B,&bb));
283: PetscCall(MatDenseGetArrayRead(S,&ss));
284: PetscCall(MatDenseGetLDA(S,&lds));
285: ld = n*n;
286: for (f=0;f<n;f++) {
287: off = f*n+f*n*ld;
288: for (i=0;i<n;i++) for (j=0;j<n;j++) aa[off+i+j*ld] = ss[i+j*lds];
289: for (c=0;c<n;c++) {
290: off = f*n+c*n*ld;
291: theta = ss[f+c*lds];
292: for (i=0;i<n;i++) aa[off+i+i*ld] += theta;
293: for (i=0;i<n;i++) for (j=0;j<n;j++) bb[off+i+j*ld] = -2*theta*ss[i+j*lds];
294: }
295: }
296: PetscCall(MatDenseRestoreArray(matctx->A,&aa));
297: PetscCall(MatDenseRestoreArray(matctx->B,&bb));
298: PetscCall(MatDenseRestoreArrayRead(S,&ss));
299: PetscCall(ISCreateStride(PETSC_COMM_SELF,n*n,0,1,&perm));
300: PetscCall(MatDestroy(&matctx->F));
301: PetscCall(MatDuplicate(matctx->A,MAT_COPY_VALUES,&matctx->F));
302: PetscCall(MatLUFactor(matctx->F,perm,perm,NULL));
303: PetscCall(ISDestroy(&perm));
304: #else
305: PetscCall(MatCopy(S,matctx->S,SAME_NONZERO_PATTERN));
306: #endif
307: matctx->lme = lme;
308: matctx->n = n;
309: PetscCall(VecSet(*v0,1.0));
310: PetscFunctionReturn(PETSC_SUCCESS);
311: }
313: static PetscErrorCode EPSSolve_LyapII(EPS eps)
314: {
315: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
316: PetscInt i,ldds,rk,nloc,mloc,nv,idx,k;
317: Vec v,w,z=eps->work[0],v0=NULL;
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,(void(*)(void))MatMult_EPSLyapIIOperator));
338: PetscCall(MatShellSetOperation(S,MATOP_DESTROY,(void(*)(void))MatDestroy_EPSLyapIIOperator));
339: PetscCall(LMESetCoefficients(ctx->lme,S,NULL,NULL,NULL));
341: /* Right-hand side */
342: PetscCall(BVDuplicateResize(eps->V,ctx->rkl,&V));
343: PetscCall(BVGetOrthogonalization(V,&type,&refine,&eta,NULL));
344: PetscCall(BVSetOrthogonalization(V,type,refine,eta,BV_ORTHOG_BLOCK_TSQR));
345: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)eps),eps->nloc,PETSC_DECIDE,PETSC_DECIDE,1,NULL,&Ux[0]));
346: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)eps),eps->nloc,PETSC_DECIDE,PETSC_DECIDE,2,NULL,&Ux[1]));
347: nv = ctx->rkl;
348: PetscCall(PetscMalloc1(nv,&s));
350: /* Initialize first column */
351: PetscCall(EPSGetStartVector(eps,0,NULL));
352: PetscCall(BVGetColumn(eps->V,0,&v));
353: PetscCall(BVInsertVec(V,0,v));
354: PetscCall(BVRestoreColumn(eps->V,0,&v));
355: PetscCall(BVSetActiveColumns(eps->V,0,0)); /* no deflation at the beginning */
356: PetscCall(LyapIIBuildRHS(S,1,Ux[0],V,eps->work));
357: idx = 0;
359: /* EPS for rank reduction */
360: PetscCall(EPSCreate(PETSC_COMM_SELF,&epsrr));
361: PetscCall(EPSSetOptionsPrefix(epsrr,((PetscObject)eps)->prefix));
362: PetscCall(EPSAppendOptionsPrefix(epsrr,"eps_lyapii_"));
363: PetscCall(EPSSetDimensions(epsrr,1,PETSC_CURRENT,PETSC_CURRENT));
364: PetscCall(EPSSetTolerances(epsrr,PETSC_MACHINE_EPSILON*100,PETSC_CURRENT));
366: while (eps->reason == EPS_CONVERGED_ITERATING) {
367: eps->its++;
369: /* Matrix for placing the solution of the Lyapunov equation (an alias of V) */
370: PetscCall(BVSetActiveColumns(V,0,nv));
371: PetscCall(BVGetMat(V,&Y1));
372: PetscCall(MatZeroEntries(Y1));
373: PetscCall(MatCreateLRC(NULL,Y1,NULL,NULL,&Y));
374: PetscCall(LMESetSolution(ctx->lme,Y));
376: /* Solve the Lyapunov equation SY + YS' = -2*S*Z*S' */
377: PetscCall(MatCreateLRC(NULL,Ux[idx],NULL,NULL,&C));
378: PetscCall(LMESetRHS(ctx->lme,C));
379: PetscCall(MatDestroy(&C));
380: PetscCall(LMESolve(ctx->lme));
381: PetscCall(BVRestoreMat(V,&Y1));
382: PetscCall(MatDestroy(&Y));
384: /* SVD of the solution: [Q,R]=qr(V); [U,Sigma,~]=svd(R) */
385: PetscCall(DSSetDimensions(ctx->ds,nv,0,0));
386: PetscCall(DSSVDSetDimensions(ctx->ds,nv));
387: PetscCall(DSGetMat(ctx->ds,DS_MAT_A,&R));
388: PetscCall(BVOrthogonalize(V,R));
389: PetscCall(DSRestoreMat(ctx->ds,DS_MAT_A,&R));
390: PetscCall(DSSetState(ctx->ds,DS_STATE_RAW));
391: PetscCall(DSSolve(ctx->ds,s,NULL));
393: /* Determine rank */
394: rk = nv;
395: for (i=1;i<nv;i++) if (PetscAbsScalar(s[i]/s[0])<PETSC_SQRT_MACHINE_EPSILON) {rk=i; break;}
396: PetscCall(PetscInfo(eps,"The computed solution of the Lyapunov equation has rank %" PetscInt_FMT "\n",rk));
397: rk = PetscMin(rk,ctx->rkc);
398: PetscCall(DSGetMat(ctx->ds,DS_MAT_U,&U));
399: PetscCall(BVMultInPlace(V,U,0,rk));
400: PetscCall(DSRestoreMat(ctx->ds,DS_MAT_U,&U));
401: PetscCall(BVSetActiveColumns(V,0,rk));
403: /* Rank reduction */
404: PetscCall(DSSetDimensions(ctx->ds,rk,0,0));
405: PetscCall(DSSVDSetDimensions(ctx->ds,rk));
406: PetscCall(DSGetMat(ctx->ds,DS_MAT_A,&W));
407: PetscCall(BVMatProject(V,S,V,W));
408: PetscCall(LyapIIBuildEigenMat(ctx->lme,W,&Op,&v0)); /* Op=A\B, A=kron(I,S)+kron(S,I), B=-2*kron(S,S) */
409: PetscCall(DSRestoreMat(ctx->ds,DS_MAT_A,&W));
410: PetscCall(EPSSetOperators(epsrr,Op,NULL));
411: PetscCall(EPSSetInitialSpace(epsrr,1,&v0));
412: PetscCall(EPSSolve(epsrr));
413: PetscCall(EPSComputeVectors(epsrr));
414: /* Copy first eigenvector, vec(A)=x */
415: PetscCall(BVGetArray(epsrr->V,&xx));
416: PetscCall(DSGetArray(ctx->ds,DS_MAT_A,&aa));
417: for (i=0;i<rk;i++) PetscCall(PetscArraycpy(aa+i*ldds,xx+i*rk,rk));
418: PetscCall(DSRestoreArray(ctx->ds,DS_MAT_A,&aa));
419: PetscCall(BVRestoreArray(epsrr->V,&xx));
420: PetscCall(DSSetState(ctx->ds,DS_STATE_RAW));
421: /* Compute [U,Sigma,~] = svd(A), its rank should be 1 or 2 */
422: PetscCall(DSSolve(ctx->ds,s,NULL));
423: if (PetscAbsScalar(s[1]/s[0])<PETSC_SQRT_MACHINE_EPSILON) rk=1;
424: else rk = 2;
425: PetscCall(PetscInfo(eps,"The eigenvector has rank %" PetscInt_FMT "\n",rk));
426: PetscCall(DSGetMat(ctx->ds,DS_MAT_U,&U));
427: PetscCall(BVMultInPlace(V,U,0,rk));
428: PetscCall(DSRestoreMat(ctx->ds,DS_MAT_U,&U));
430: /* Save V in Ux */
431: idx = (rk==2)?1:0;
432: for (i=0;i<rk;i++) {
433: PetscCall(BVGetColumn(V,i,&v));
434: PetscCall(VecGetArray(v,&uu));
435: PetscCall(MatDenseGetColumn(Ux[idx],i,&array));
436: PetscCall(PetscArraycpy(array,uu,eps->nloc));
437: PetscCall(MatDenseRestoreColumn(Ux[idx],&array));
438: PetscCall(VecRestoreArray(v,&uu));
439: PetscCall(BVRestoreColumn(V,i,&v));
440: }
442: /* Eigenpair approximation */
443: PetscCall(BVGetColumn(V,0,&v));
444: PetscCall(MatMult(S,v,z));
445: PetscCall(VecDot(z,v,pM));
446: PetscCall(BVRestoreColumn(V,0,&v));
447: if (rk>1) {
448: PetscCall(BVGetColumn(V,1,&w));
449: PetscCall(VecDot(z,w,pM+1));
450: PetscCall(MatMult(S,w,z));
451: PetscCall(VecDot(z,w,pM+3));
452: PetscCall(BVGetColumn(V,0,&v));
453: PetscCall(VecDot(z,v,pM+2));
454: PetscCall(BVRestoreColumn(V,0,&v));
455: PetscCall(BVRestoreColumn(V,1,&w));
456: PetscCall(EV2x2(pM,2,eigr,eigi,vec));
457: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,2,2,vec,&X));
458: PetscCall(BVSetActiveColumns(V,0,rk));
459: PetscCall(BVMultInPlace(V,X,0,rk));
460: PetscCall(MatDestroy(&X));
461: #if !defined(PETSC_USE_COMPLEX)
462: norm = eigr[0]*eigr[0]+eigi[0]*eigi[0];
463: er = eigr[0]/norm; ei = -eigi[0]/norm;
464: #else
465: er =1.0/eigr[0]; ei = 0.0;
466: #endif
467: } else {
468: eigr[0] = pM[0]; eigi[0] = 0.0;
469: er = 1.0/eigr[0]; ei = 0.0;
470: }
471: PetscCall(BVGetColumn(V,0,&v));
472: if (eigi[0]!=0.0) PetscCall(BVGetColumn(V,1,&w));
473: else w = NULL;
474: eps->eigr[eps->nconv] = eigr[0]; eps->eigi[eps->nconv] = eigi[0];
475: PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_FALSE,er,ei,v,w,eps->work,&norm));
476: PetscCall(BVRestoreColumn(V,0,&v));
477: if (w) PetscCall(BVRestoreColumn(V,1,&w));
478: PetscCall((*eps->converged)(eps,er,ei,norm,&eps->errest[eps->nconv],eps->convergedctx));
479: k = 0;
480: if (eps->errest[eps->nconv]<eps->tol) {
481: k++;
482: if (rk==2) {
483: #if !defined (PETSC_USE_COMPLEX)
484: eps->eigr[eps->nconv+k] = eigr[0]; eps->eigi[eps->nconv+k] = -eigi[0];
485: #else
486: eps->eigr[eps->nconv+k] = PetscConj(eps->eigr[eps->nconv]);
487: #endif
488: k++;
489: }
490: /* Store converged eigenpairs and vectors for deflation */
491: for (i=0;i<k;i++) {
492: PetscCall(BVGetColumn(V,i,&v));
493: PetscCall(BVInsertVec(eps->V,eps->nconv+i,v));
494: PetscCall(BVRestoreColumn(V,i,&v));
495: }
496: eps->nconv += k;
497: PetscCall(BVSetActiveColumns(eps->V,eps->nconv-rk,eps->nconv));
498: PetscCall(BVOrthogonalize(eps->V,NULL));
499: PetscCall(DSSetDimensions(eps->ds,eps->nconv,0,0));
500: PetscCall(DSGetMat(eps->ds,DS_MAT_A,&W));
501: PetscCall(BVMatProject(eps->V,matctx->S,eps->V,W));
502: PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&W));
503: if (eps->nconv<eps->nev) {
504: idx = 0;
505: PetscCall(BVSetRandomColumn(V,0));
506: PetscCall(BVNormColumn(V,0,NORM_2,&norm));
507: PetscCall(BVScaleColumn(V,0,1.0/norm));
508: PetscCall(LyapIIBuildRHS(S,1,Ux[idx],V,eps->work));
509: }
510: } else {
511: /* Prepare right-hand side */
512: PetscCall(LyapIIBuildRHS(S,rk,Ux[idx],NULL,eps->work));
513: }
514: PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx));
515: PetscCall(EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,eps->nconv+1));
516: }
517: PetscCall(STRestoreOperator(eps->st,&matctx->S));
518: PetscCall(MatDestroy(&S));
519: PetscCall(MatDestroy(&Ux[0]));
520: PetscCall(MatDestroy(&Ux[1]));
521: PetscCall(MatDestroy(&Op));
522: PetscCall(VecDestroy(&v0));
523: PetscCall(BVDestroy(&V));
524: PetscCall(EPSDestroy(&epsrr));
525: PetscCall(PetscFree(s));
526: PetscFunctionReturn(PETSC_SUCCESS);
527: }
529: static PetscErrorCode EPSSetFromOptions_LyapII(EPS eps,PetscOptionItems *PetscOptionsObject)
530: {
531: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
532: PetscInt k,array[2]={PETSC_DETERMINE,PETSC_DETERMINE};
533: PetscBool flg;
535: PetscFunctionBegin;
536: PetscOptionsHeadBegin(PetscOptionsObject,"EPS Lyapunov Inverse Iteration Options");
538: k = 2;
539: PetscCall(PetscOptionsIntArray("-eps_lyapii_ranks","Ranks for Lyapunov equation (one or two comma-separated integers)","EPSLyapIISetRanks",array,&k,&flg));
540: if (flg) PetscCall(EPSLyapIISetRanks(eps,array[0],array[1]));
542: PetscOptionsHeadEnd();
544: if (!ctx->lme) PetscCall(EPSLyapIIGetLME(eps,&ctx->lme));
545: PetscCall(LMESetFromOptions(ctx->lme));
546: PetscFunctionReturn(PETSC_SUCCESS);
547: }
549: static PetscErrorCode EPSLyapIISetRanks_LyapII(EPS eps,PetscInt rkc,PetscInt rkl)
550: {
551: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
553: PetscFunctionBegin;
554: if (rkc==PETSC_DETERMINE) {
555: if (ctx->rkc != 10) eps->state = EPS_STATE_INITIAL;
556: ctx->rkc = 10;
557: } else if (rkc!=PETSC_CURRENT) {
558: PetscCheck(rkc>1,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The compressed rank %" PetscInt_FMT " must be larger than 1",rkc);
559: if (ctx->rkc != rkc) eps->state = EPS_STATE_INITIAL;
560: ctx->rkc = rkc;
561: }
562: if (rkl==PETSC_DETERMINE) {
563: if (ctx->rkl != 3*rkc) eps->state = EPS_STATE_INITIAL;
564: ctx->rkl = 3*rkc;
565: } else if (rkl!=PETSC_CURRENT) {
566: 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);
567: if (ctx->rkl != rkl) eps->state = EPS_STATE_INITIAL;
568: ctx->rkl = rkl;
569: }
570: PetscFunctionReturn(PETSC_SUCCESS);
571: }
573: /*@
574: EPSLyapIISetRanks - Set the ranks used in the solution of the Lyapunov equation.
576: Logically Collective
578: Input Parameters:
579: + eps - the eigenproblem solver context
580: . rkc - the compressed rank
581: - rkl - the Lyapunov rank
583: Options Database Key:
584: . -eps_lyapii_ranks <rkc,rkl> - Sets the rank parameters
586: Notes:
587: PETSC_CURRENT can be used to preserve the current value of any of the
588: arguments, and PETSC_DETERMINE to set them to a default value.
590: Lyapunov inverse iteration needs to solve a large-scale Lyapunov equation
591: at each iteration of the eigensolver. For this, an iterative solver (LME)
592: is used, which requires to prescribe the rank of the solution matrix X. This
593: is the meaning of parameter rkl. Later, this matrix is compressed into
594: another matrix of rank rkc. If not provided, rkl is a small multiple of rkc.
596: Level: intermediate
598: .seealso: EPSLyapIIGetRanks()
599: @*/
600: PetscErrorCode EPSLyapIISetRanks(EPS eps,PetscInt rkc,PetscInt rkl)
601: {
602: PetscFunctionBegin;
606: PetscTryMethod(eps,"EPSLyapIISetRanks_C",(EPS,PetscInt,PetscInt),(eps,rkc,rkl));
607: PetscFunctionReturn(PETSC_SUCCESS);
608: }
610: static PetscErrorCode EPSLyapIIGetRanks_LyapII(EPS eps,PetscInt *rkc,PetscInt *rkl)
611: {
612: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
614: PetscFunctionBegin;
615: if (rkc) *rkc = ctx->rkc;
616: if (rkl) *rkl = ctx->rkl;
617: PetscFunctionReturn(PETSC_SUCCESS);
618: }
620: /*@
621: EPSLyapIIGetRanks - Return the rank values used for the Lyapunov step.
623: Not Collective
625: Input Parameter:
626: . eps - the eigenproblem solver context
628: Output Parameters:
629: + rkc - the compressed rank
630: - rkl - the Lyapunov rank
632: Level: intermediate
634: .seealso: EPSLyapIISetRanks()
635: @*/
636: PetscErrorCode EPSLyapIIGetRanks(EPS eps,PetscInt *rkc,PetscInt *rkl)
637: {
638: PetscFunctionBegin;
640: PetscUseMethod(eps,"EPSLyapIIGetRanks_C",(EPS,PetscInt*,PetscInt*),(eps,rkc,rkl));
641: PetscFunctionReturn(PETSC_SUCCESS);
642: }
644: static PetscErrorCode EPSLyapIISetLME_LyapII(EPS eps,LME lme)
645: {
646: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
648: PetscFunctionBegin;
649: PetscCall(PetscObjectReference((PetscObject)lme));
650: PetscCall(LMEDestroy(&ctx->lme));
651: ctx->lme = lme;
652: eps->state = EPS_STATE_INITIAL;
653: PetscFunctionReturn(PETSC_SUCCESS);
654: }
656: /*@
657: EPSLyapIISetLME - Associate a linear matrix equation solver object (LME) to the
658: eigenvalue solver.
660: Collective
662: Input Parameters:
663: + eps - the eigenproblem solver context
664: - lme - the linear matrix equation solver object
666: Level: advanced
668: .seealso: EPSLyapIIGetLME()
669: @*/
670: PetscErrorCode EPSLyapIISetLME(EPS eps,LME lme)
671: {
672: PetscFunctionBegin;
675: PetscCheckSameComm(eps,1,lme,2);
676: PetscTryMethod(eps,"EPSLyapIISetLME_C",(EPS,LME),(eps,lme));
677: PetscFunctionReturn(PETSC_SUCCESS);
678: }
680: static PetscErrorCode EPSLyapIIGetLME_LyapII(EPS eps,LME *lme)
681: {
682: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
684: PetscFunctionBegin;
685: if (!ctx->lme) {
686: PetscCall(LMECreate(PetscObjectComm((PetscObject)eps),&ctx->lme));
687: PetscCall(LMESetOptionsPrefix(ctx->lme,((PetscObject)eps)->prefix));
688: PetscCall(LMEAppendOptionsPrefix(ctx->lme,"eps_lyapii_"));
689: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->lme,(PetscObject)eps,1));
690: }
691: *lme = ctx->lme;
692: PetscFunctionReturn(PETSC_SUCCESS);
693: }
695: /*@
696: EPSLyapIIGetLME - Retrieve the linear matrix equation solver object (LME)
697: associated with the eigenvalue solver.
699: Not Collective
701: Input Parameter:
702: . eps - the eigenproblem solver context
704: Output Parameter:
705: . lme - the linear matrix equation solver object
707: Level: advanced
709: .seealso: EPSLyapIISetLME()
710: @*/
711: PetscErrorCode EPSLyapIIGetLME(EPS eps,LME *lme)
712: {
713: PetscFunctionBegin;
715: PetscAssertPointer(lme,2);
716: PetscUseMethod(eps,"EPSLyapIIGetLME_C",(EPS,LME*),(eps,lme));
717: PetscFunctionReturn(PETSC_SUCCESS);
718: }
720: static PetscErrorCode EPSView_LyapII(EPS eps,PetscViewer viewer)
721: {
722: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
723: PetscBool isascii;
725: PetscFunctionBegin;
726: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
727: if (isascii) {
728: PetscCall(PetscViewerASCIIPrintf(viewer," ranks: for Lyapunov solver=%" PetscInt_FMT ", after compression=%" PetscInt_FMT "\n",ctx->rkl,ctx->rkc));
729: if (!ctx->lme) PetscCall(EPSLyapIIGetLME(eps,&ctx->lme));
730: PetscCall(PetscViewerASCIIPushTab(viewer));
731: PetscCall(LMEView(ctx->lme,viewer));
732: PetscCall(PetscViewerASCIIPopTab(viewer));
733: }
734: PetscFunctionReturn(PETSC_SUCCESS);
735: }
737: static PetscErrorCode EPSReset_LyapII(EPS eps)
738: {
739: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
741: PetscFunctionBegin;
742: if (!ctx->lme) PetscCall(LMEReset(ctx->lme));
743: PetscFunctionReturn(PETSC_SUCCESS);
744: }
746: static PetscErrorCode EPSDestroy_LyapII(EPS eps)
747: {
748: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
750: PetscFunctionBegin;
751: PetscCall(LMEDestroy(&ctx->lme));
752: PetscCall(DSDestroy(&ctx->ds));
753: PetscCall(PetscFree(eps->data));
754: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetLME_C",NULL));
755: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetLME_C",NULL));
756: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetRanks_C",NULL));
757: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetRanks_C",NULL));
758: PetscFunctionReturn(PETSC_SUCCESS);
759: }
761: static PetscErrorCode EPSSetDefaultST_LyapII(EPS eps)
762: {
763: PetscFunctionBegin;
764: if (!((PetscObject)eps->st)->type_name) PetscCall(STSetType(eps->st,STSINVERT));
765: PetscFunctionReturn(PETSC_SUCCESS);
766: }
768: SLEPC_EXTERN PetscErrorCode EPSCreate_LyapII(EPS eps)
769: {
770: EPS_LYAPII *ctx;
772: PetscFunctionBegin;
773: PetscCall(PetscNew(&ctx));
774: eps->data = (void*)ctx;
776: eps->useds = PETSC_TRUE;
778: eps->ops->solve = EPSSolve_LyapII;
779: eps->ops->setup = EPSSetUp_LyapII;
780: eps->ops->setupsort = EPSSetUpSort_Default;
781: eps->ops->setfromoptions = EPSSetFromOptions_LyapII;
782: eps->ops->reset = EPSReset_LyapII;
783: eps->ops->destroy = EPSDestroy_LyapII;
784: eps->ops->view = EPSView_LyapII;
785: eps->ops->setdefaultst = EPSSetDefaultST_LyapII;
786: eps->ops->backtransform = EPSBackTransform_Default;
787: eps->ops->computevectors = EPSComputeVectors_Schur;
789: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetLME_C",EPSLyapIISetLME_LyapII));
790: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetLME_C",EPSLyapIIGetLME_LyapII));
791: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetRanks_C",EPSLyapIISetRanks_LyapII));
792: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetRanks_C",EPSLyapIIGetRanks_LyapII));
793: PetscFunctionReturn(PETSC_SUCCESS);
794: }