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