Actual source code: lyapii.c
slepc-3.22.1 2024-10-28
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->ncv!=PETSC_DETERMINE) {
64: PetscCheck(eps->ncv>=eps->nev+1,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nev+1");
65: } else eps->ncv = eps->nev+1;
66: if (eps->mpd!=PETSC_DETERMINE) PetscCall(PetscInfo(eps,"Warning: parameter mpd ignored\n"));
67: if (eps->max_it==PETSC_DETERMINE) eps->max_it = PetscMax(1000*eps->nev,100*eps->n);
68: if (!eps->which) eps->which=EPS_LARGEST_REAL;
69: PetscCheck(eps->which==EPS_LARGEST_REAL,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver supports only largest real eigenvalues");
70: EPSCheckUnsupported(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION | EPS_FEATURE_TWOSIDED);
72: if (!ctx->rkc) ctx->rkc = 10;
73: if (!ctx->rkl) ctx->rkl = 3*ctx->rkc;
74: if (!ctx->lme) PetscCall(EPSLyapIIGetLME(eps,&ctx->lme));
75: PetscCall(LMESetProblemType(ctx->lme,LME_LYAPUNOV));
76: PetscCall(LMESetErrorIfNotConverged(ctx->lme,PETSC_TRUE));
78: if (!ctx->ds) {
79: PetscCall(DSCreate(PetscObjectComm((PetscObject)eps),&ctx->ds));
80: PetscCall(DSSetType(ctx->ds,DSSVD));
81: }
82: PetscCall(DSAllocate(ctx->ds,ctx->rkl));
84: PetscCall(DSSetType(eps->ds,DSNHEP));
85: PetscCall(DSAllocate(eps->ds,eps->ncv));
87: PetscCall(EPSAllocateSolution(eps,0));
88: PetscCall(BVGetRandomContext(eps->V,&rand)); /* make sure the random context is available when duplicating */
89: PetscCall(EPSSetWorkVecs(eps,3));
90: PetscFunctionReturn(PETSC_SUCCESS);
91: }
93: static PetscErrorCode MatMult_EPSLyapIIOperator(Mat M,Vec x,Vec r)
94: {
95: EPS_LYAPII_MATSHELL *matctx;
97: PetscFunctionBegin;
98: PetscCall(MatShellGetContext(M,&matctx));
99: PetscCall(MatMult(matctx->S,x,r));
100: PetscCall(BVOrthogonalizeVec(matctx->Q,r,NULL,NULL,NULL));
101: PetscFunctionReturn(PETSC_SUCCESS);
102: }
104: static PetscErrorCode MatDestroy_EPSLyapIIOperator(Mat M)
105: {
106: EPS_LYAPII_MATSHELL *matctx;
108: PetscFunctionBegin;
109: PetscCall(MatShellGetContext(M,&matctx));
110: PetscCall(MatDestroy(&matctx->S));
111: PetscCall(PetscFree(matctx));
112: PetscFunctionReturn(PETSC_SUCCESS);
113: }
115: static PetscErrorCode MatMult_EigOperator(Mat M,Vec x,Vec y)
116: {
117: EPS_EIG_MATSHELL *matctx;
118: #if !defined(PETSC_USE_COMPLEX)
119: PetscInt n,lds;
120: PetscScalar *Y,*C,zero=0.0,done=1.0,dtwo=2.0;
121: const PetscScalar *S,*X;
122: PetscBLASInt n_,lds_;
123: #endif
125: PetscFunctionBegin;
126: PetscCall(MatShellGetContext(M,&matctx));
128: #if defined(PETSC_USE_COMPLEX)
129: PetscCall(MatMult(matctx->B,x,matctx->w));
130: PetscCall(MatSolve(matctx->F,matctx->w,y));
131: #else
132: PetscCall(VecGetArrayRead(x,&X));
133: PetscCall(VecGetArray(y,&Y));
134: PetscCall(MatDenseGetArrayRead(matctx->S,&S));
135: PetscCall(MatDenseGetLDA(matctx->S,&lds));
137: n = matctx->n;
138: PetscCall(PetscCalloc1(n*n,&C));
139: PetscCall(PetscBLASIntCast(n,&n_));
140: PetscCall(PetscBLASIntCast(lds,&lds_));
142: /* C = 2*S*X*S.' */
143: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&dtwo,S,&lds_,X,&n_,&zero,Y,&n_));
144: PetscCallBLAS("BLASgemm",BLASgemm_("N","T",&n_,&n_,&n_,&done,Y,&n_,S,&lds_,&zero,C,&n_));
146: /* Solve S*Y + Y*S' = -C */
147: PetscCall(LMEDenseLyapunov(matctx->lme,n,(PetscScalar*)S,lds,C,n,Y,n));
149: PetscCall(PetscFree(C));
150: PetscCall(VecRestoreArrayRead(x,&X));
151: PetscCall(VecRestoreArray(y,&Y));
152: PetscCall(MatDenseRestoreArrayRead(matctx->S,&S));
153: #endif
154: PetscFunctionReturn(PETSC_SUCCESS);
155: }
157: static PetscErrorCode MatDestroy_EigOperator(Mat M)
158: {
159: EPS_EIG_MATSHELL *matctx;
161: PetscFunctionBegin;
162: PetscCall(MatShellGetContext(M,&matctx));
163: #if defined(PETSC_USE_COMPLEX)
164: PetscCall(MatDestroy(&matctx->A));
165: PetscCall(MatDestroy(&matctx->B));
166: PetscCall(MatDestroy(&matctx->F));
167: PetscCall(VecDestroy(&matctx->w));
168: #else
169: PetscCall(MatDestroy(&matctx->S));
170: #endif
171: PetscCall(PetscFree(matctx));
172: PetscFunctionReturn(PETSC_SUCCESS);
173: }
175: /*
176: EV2x2: solve the eigenproblem for a 2x2 matrix M
177: */
178: static PetscErrorCode EV2x2(PetscScalar *M,PetscInt ld,PetscScalar *wr,PetscScalar *wi,PetscScalar *vec)
179: {
180: PetscBLASInt lwork=10,ld_;
181: PetscScalar work[10];
182: PetscBLASInt two=2,info;
183: #if defined(PETSC_USE_COMPLEX)
184: PetscReal rwork[6];
185: #endif
187: PetscFunctionBegin;
188: PetscCall(PetscBLASIntCast(ld,&ld_));
189: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
190: #if !defined(PETSC_USE_COMPLEX)
191: PetscCallBLAS("LAPACKgeev",LAPACKgeev_("N","V",&two,M,&ld_,wr,wi,NULL,&ld_,vec,&ld_,work,&lwork,&info));
192: #else
193: PetscCallBLAS("LAPACKgeev",LAPACKgeev_("N","V",&two,M,&ld_,wr,NULL,&ld_,vec,&ld_,work,&lwork,rwork,&info));
194: #endif
195: SlepcCheckLapackInfo("geev",info);
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 defined(PETSC_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 defined(PETSC_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,(void(*)(void))MatMult_EigOperator));
271: PetscCall(MatShellSetOperation(*Op,MATOP_DESTROY,(void(*)(void))MatDestroy_EigOperator));
272: PetscCall(MatCreateVecs(*Op,NULL,v0));
273: } else {
274: PetscCall(MatShellGetContext(*Op,&matctx));
275: #if defined(PETSC_USE_COMPLEX)
276: PetscCall(MatZeroEntries(matctx->A));
277: #endif
278: }
279: #if defined(PETSC_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: Mat S,C,Ux[2],Y,Y1,R,U,W,X,Op=NULL;
318: BV V;
319: BVOrthogType type;
320: BVOrthogRefineType refine;
321: PetscScalar eigr[2],eigi[2],*array,er,ei,*uu,*s,*xx,*aa,pM[4],vec[4];
322: PetscReal eta;
323: EPS epsrr;
324: PetscReal norm;
325: EPS_LYAPII_MATSHELL *matctx;
327: PetscFunctionBegin;
328: PetscCall(DSGetLeadingDimension(ctx->ds,&ldds));
330: /* Operator for the Lyapunov equation */
331: PetscCall(PetscNew(&matctx));
332: PetscCall(STGetOperator(eps->st,&matctx->S));
333: PetscCall(MatGetLocalSize(matctx->S,&mloc,&nloc));
334: PetscCall(MatCreateShell(PetscObjectComm((PetscObject)eps),mloc,nloc,PETSC_DETERMINE,PETSC_DETERMINE,matctx,&S));
335: matctx->Q = eps->V;
336: PetscCall(MatShellSetOperation(S,MATOP_MULT,(void(*)(void))MatMult_EPSLyapIIOperator));
337: PetscCall(MatShellSetOperation(S,MATOP_DESTROY,(void(*)(void))MatDestroy_EPSLyapIIOperator));
338: PetscCall(LMESetCoefficients(ctx->lme,S,NULL,NULL,NULL));
340: /* Right-hand side */
341: PetscCall(BVDuplicateResize(eps->V,ctx->rkl,&V));
342: PetscCall(BVGetOrthogonalization(V,&type,&refine,&eta,NULL));
343: PetscCall(BVSetOrthogonalization(V,type,refine,eta,BV_ORTHOG_BLOCK_TSQR));
344: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)eps),eps->nloc,PETSC_DECIDE,PETSC_DECIDE,1,NULL,&Ux[0]));
345: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)eps),eps->nloc,PETSC_DECIDE,PETSC_DECIDE,2,NULL,&Ux[1]));
346: nv = ctx->rkl;
347: PetscCall(PetscMalloc1(nv,&s));
349: /* Initialize first column */
350: PetscCall(EPSGetStartVector(eps,0,NULL));
351: PetscCall(BVGetColumn(eps->V,0,&v));
352: PetscCall(BVInsertVec(V,0,v));
353: PetscCall(BVRestoreColumn(eps->V,0,&v));
354: PetscCall(BVSetActiveColumns(eps->V,0,0)); /* no deflation at the beginning */
355: PetscCall(LyapIIBuildRHS(S,1,Ux[0],V,eps->work));
356: idx = 0;
358: /* EPS for rank reduction */
359: PetscCall(EPSCreate(PETSC_COMM_SELF,&epsrr));
360: PetscCall(EPSSetOptionsPrefix(epsrr,((PetscObject)eps)->prefix));
361: PetscCall(EPSAppendOptionsPrefix(epsrr,"eps_lyapii_"));
362: PetscCall(EPSSetDimensions(epsrr,1,PETSC_CURRENT,PETSC_CURRENT));
363: PetscCall(EPSSetTolerances(epsrr,PETSC_MACHINE_EPSILON*100,PETSC_CURRENT));
365: while (eps->reason == EPS_CONVERGED_ITERATING) {
366: eps->its++;
368: /* Matrix for placing the solution of the Lyapunov equation (an alias of V) */
369: PetscCall(BVSetActiveColumns(V,0,nv));
370: PetscCall(BVGetMat(V,&Y1));
371: PetscCall(MatZeroEntries(Y1));
372: PetscCall(MatCreateLRC(NULL,Y1,NULL,NULL,&Y));
373: PetscCall(LMESetSolution(ctx->lme,Y));
375: /* Solve the Lyapunov equation SY + YS' = -2*S*Z*S' */
376: PetscCall(MatCreateLRC(NULL,Ux[idx],NULL,NULL,&C));
377: PetscCall(LMESetRHS(ctx->lme,C));
378: PetscCall(MatDestroy(&C));
379: PetscCall(LMESolve(ctx->lme));
380: PetscCall(BVRestoreMat(V,&Y1));
381: PetscCall(MatDestroy(&Y));
383: /* SVD of the solution: [Q,R]=qr(V); [U,Sigma,~]=svd(R) */
384: PetscCall(DSSetDimensions(ctx->ds,nv,0,0));
385: PetscCall(DSSVDSetDimensions(ctx->ds,nv));
386: PetscCall(DSGetMat(ctx->ds,DS_MAT_A,&R));
387: PetscCall(BVOrthogonalize(V,R));
388: PetscCall(DSRestoreMat(ctx->ds,DS_MAT_A,&R));
389: PetscCall(DSSetState(ctx->ds,DS_STATE_RAW));
390: PetscCall(DSSolve(ctx->ds,s,NULL));
392: /* Determine rank */
393: rk = nv;
394: for (i=1;i<nv;i++) if (PetscAbsScalar(s[i]/s[0])<PETSC_SQRT_MACHINE_EPSILON) {rk=i; break;}
395: PetscCall(PetscInfo(eps,"The computed solution of the Lyapunov equation has rank %" PetscInt_FMT "\n",rk));
396: rk = PetscMin(rk,ctx->rkc);
397: PetscCall(DSGetMat(ctx->ds,DS_MAT_U,&U));
398: PetscCall(BVMultInPlace(V,U,0,rk));
399: PetscCall(DSRestoreMat(ctx->ds,DS_MAT_U,&U));
400: PetscCall(BVSetActiveColumns(V,0,rk));
402: /* Rank reduction */
403: PetscCall(DSSetDimensions(ctx->ds,rk,0,0));
404: PetscCall(DSSVDSetDimensions(ctx->ds,rk));
405: PetscCall(DSGetMat(ctx->ds,DS_MAT_A,&W));
406: PetscCall(BVMatProject(V,S,V,W));
407: PetscCall(LyapIIBuildEigenMat(ctx->lme,W,&Op,&v0)); /* Op=A\B, A=kron(I,S)+kron(S,I), B=-2*kron(S,S) */
408: PetscCall(DSRestoreMat(ctx->ds,DS_MAT_A,&W));
409: PetscCall(EPSSetOperators(epsrr,Op,NULL));
410: PetscCall(EPSSetInitialSpace(epsrr,1,&v0));
411: PetscCall(EPSSolve(epsrr));
412: PetscCall(EPSComputeVectors(epsrr));
413: /* Copy first eigenvector, vec(A)=x */
414: PetscCall(BVGetArray(epsrr->V,&xx));
415: PetscCall(DSGetArray(ctx->ds,DS_MAT_A,&aa));
416: for (i=0;i<rk;i++) PetscCall(PetscArraycpy(aa+i*ldds,xx+i*rk,rk));
417: PetscCall(DSRestoreArray(ctx->ds,DS_MAT_A,&aa));
418: PetscCall(BVRestoreArray(epsrr->V,&xx));
419: PetscCall(DSSetState(ctx->ds,DS_STATE_RAW));
420: /* Compute [U,Sigma,~] = svd(A), its rank should be 1 or 2 */
421: PetscCall(DSSolve(ctx->ds,s,NULL));
422: if (PetscAbsScalar(s[1]/s[0])<PETSC_SQRT_MACHINE_EPSILON) rk=1;
423: else rk = 2;
424: PetscCall(PetscInfo(eps,"The eigenvector has rank %" PetscInt_FMT "\n",rk));
425: PetscCall(DSGetMat(ctx->ds,DS_MAT_U,&U));
426: PetscCall(BVMultInPlace(V,U,0,rk));
427: PetscCall(DSRestoreMat(ctx->ds,DS_MAT_U,&U));
429: /* Save V in Ux */
430: idx = (rk==2)?1:0;
431: for (i=0;i<rk;i++) {
432: PetscCall(BVGetColumn(V,i,&v));
433: PetscCall(VecGetArray(v,&uu));
434: PetscCall(MatDenseGetColumn(Ux[idx],i,&array));
435: PetscCall(PetscArraycpy(array,uu,eps->nloc));
436: PetscCall(MatDenseRestoreColumn(Ux[idx],&array));
437: PetscCall(VecRestoreArray(v,&uu));
438: PetscCall(BVRestoreColumn(V,i,&v));
439: }
441: /* Eigenpair approximation */
442: PetscCall(BVGetColumn(V,0,&v));
443: PetscCall(MatMult(S,v,z));
444: PetscCall(VecDot(z,v,pM));
445: PetscCall(BVRestoreColumn(V,0,&v));
446: if (rk>1) {
447: PetscCall(BVGetColumn(V,1,&w));
448: PetscCall(VecDot(z,w,pM+1));
449: PetscCall(MatMult(S,w,z));
450: PetscCall(VecDot(z,w,pM+3));
451: PetscCall(BVGetColumn(V,0,&v));
452: PetscCall(VecDot(z,v,pM+2));
453: PetscCall(BVRestoreColumn(V,0,&v));
454: PetscCall(BVRestoreColumn(V,1,&w));
455: PetscCall(EV2x2(pM,2,eigr,eigi,vec));
456: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,2,2,vec,&X));
457: PetscCall(BVSetActiveColumns(V,0,rk));
458: PetscCall(BVMultInPlace(V,X,0,rk));
459: PetscCall(MatDestroy(&X));
460: #if !defined(PETSC_USE_COMPLEX)
461: norm = eigr[0]*eigr[0]+eigi[0]*eigi[0];
462: er = eigr[0]/norm; ei = -eigi[0]/norm;
463: #else
464: er =1.0/eigr[0]; ei = 0.0;
465: #endif
466: } else {
467: eigr[0] = pM[0]; eigi[0] = 0.0;
468: er = 1.0/eigr[0]; ei = 0.0;
469: }
470: PetscCall(BVGetColumn(V,0,&v));
471: if (eigi[0]!=0.0) PetscCall(BVGetColumn(V,1,&w));
472: else w = NULL;
473: eps->eigr[eps->nconv] = eigr[0]; eps->eigi[eps->nconv] = eigi[0];
474: PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_FALSE,er,ei,v,w,eps->work,&norm));
475: PetscCall(BVRestoreColumn(V,0,&v));
476: if (w) PetscCall(BVRestoreColumn(V,1,&w));
477: PetscCall((*eps->converged)(eps,er,ei,norm,&eps->errest[eps->nconv],eps->convergedctx));
478: k = 0;
479: if (eps->errest[eps->nconv]<eps->tol) {
480: k++;
481: if (rk==2) {
482: #if !defined (PETSC_USE_COMPLEX)
483: eps->eigr[eps->nconv+k] = eigr[0]; eps->eigi[eps->nconv+k] = -eigi[0];
484: #else
485: eps->eigr[eps->nconv+k] = PetscConj(eps->eigr[eps->nconv]);
486: #endif
487: k++;
488: }
489: /* Store converged eigenpairs and vectors for deflation */
490: for (i=0;i<k;i++) {
491: PetscCall(BVGetColumn(V,i,&v));
492: PetscCall(BVInsertVec(eps->V,eps->nconv+i,v));
493: PetscCall(BVRestoreColumn(V,i,&v));
494: }
495: eps->nconv += k;
496: PetscCall(BVSetActiveColumns(eps->V,eps->nconv-rk,eps->nconv));
497: PetscCall(BVOrthogonalize(eps->V,NULL));
498: PetscCall(DSSetDimensions(eps->ds,eps->nconv,0,0));
499: PetscCall(DSGetMat(eps->ds,DS_MAT_A,&W));
500: PetscCall(BVMatProject(eps->V,matctx->S,eps->V,W));
501: PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&W));
502: if (eps->nconv<eps->nev) {
503: idx = 0;
504: PetscCall(BVSetRandomColumn(V,0));
505: PetscCall(BVNormColumn(V,0,NORM_2,&norm));
506: PetscCall(BVScaleColumn(V,0,1.0/norm));
507: PetscCall(LyapIIBuildRHS(S,1,Ux[idx],V,eps->work));
508: }
509: } else {
510: /* Prepare right-hand side */
511: PetscCall(LyapIIBuildRHS(S,rk,Ux[idx],NULL,eps->work));
512: }
513: PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx));
514: PetscCall(EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,eps->nconv+1));
515: }
516: PetscCall(STRestoreOperator(eps->st,&matctx->S));
517: PetscCall(MatDestroy(&S));
518: PetscCall(MatDestroy(&Ux[0]));
519: PetscCall(MatDestroy(&Ux[1]));
520: PetscCall(MatDestroy(&Op));
521: PetscCall(VecDestroy(&v0));
522: PetscCall(BVDestroy(&V));
523: PetscCall(EPSDestroy(&epsrr));
524: PetscCall(PetscFree(s));
525: PetscFunctionReturn(PETSC_SUCCESS);
526: }
528: static PetscErrorCode EPSSetFromOptions_LyapII(EPS eps,PetscOptionItems *PetscOptionsObject)
529: {
530: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
531: PetscInt k,array[2]={PETSC_DETERMINE,PETSC_DETERMINE};
532: PetscBool flg;
534: PetscFunctionBegin;
535: PetscOptionsHeadBegin(PetscOptionsObject,"EPS Lyapunov Inverse Iteration Options");
537: k = 2;
538: PetscCall(PetscOptionsIntArray("-eps_lyapii_ranks","Ranks for Lyapunov equation (one or two comma-separated integers)","EPSLyapIISetRanks",array,&k,&flg));
539: if (flg) PetscCall(EPSLyapIISetRanks(eps,array[0],array[1]));
541: PetscOptionsHeadEnd();
543: if (!ctx->lme) PetscCall(EPSLyapIIGetLME(eps,&ctx->lme));
544: PetscCall(LMESetFromOptions(ctx->lme));
545: PetscFunctionReturn(PETSC_SUCCESS);
546: }
548: static PetscErrorCode EPSLyapIISetRanks_LyapII(EPS eps,PetscInt rkc,PetscInt rkl)
549: {
550: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
552: PetscFunctionBegin;
553: if (rkc==PETSC_DETERMINE) {
554: if (ctx->rkc != 10) eps->state = EPS_STATE_INITIAL;
555: ctx->rkc = 10;
556: } else if (rkc!=PETSC_CURRENT) {
557: PetscCheck(rkc>1,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The compressed rank %" PetscInt_FMT " must be larger than 1",rkc);
558: if (ctx->rkc != rkc) eps->state = EPS_STATE_INITIAL;
559: ctx->rkc = rkc;
560: }
561: if (rkl==PETSC_DETERMINE) {
562: if (ctx->rkl != 3*rkc) eps->state = EPS_STATE_INITIAL;
563: ctx->rkl = 3*rkc;
564: } else if (rkl!=PETSC_CURRENT) {
565: 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);
566: if (ctx->rkl != rkl) eps->state = EPS_STATE_INITIAL;
567: ctx->rkl = rkl;
568: }
569: PetscFunctionReturn(PETSC_SUCCESS);
570: }
572: /*@
573: EPSLyapIISetRanks - Set the ranks used in the solution of the Lyapunov equation.
575: Logically Collective
577: Input Parameters:
578: + eps - the eigenproblem solver context
579: . rkc - the compressed rank
580: - rkl - the Lyapunov rank
582: Options Database Key:
583: . -eps_lyapii_ranks <rkc,rkl> - Sets the rank parameters
585: Notes:
586: PETSC_CURRENT can be used to preserve the current value of any of the
587: arguments, and PETSC_DETERMINE to set them to a default value.
589: Lyapunov inverse iteration needs to solve a large-scale Lyapunov equation
590: at each iteration of the eigensolver. For this, an iterative solver (LME)
591: is used, which requires to prescribe the rank of the solution matrix X. This
592: is the meaning of parameter rkl. Later, this matrix is compressed into
593: another matrix of rank rkc. If not provided, rkl is a small multiple of rkc.
595: Level: intermediate
597: .seealso: EPSLyapIIGetRanks()
598: @*/
599: PetscErrorCode EPSLyapIISetRanks(EPS eps,PetscInt rkc,PetscInt rkl)
600: {
601: PetscFunctionBegin;
605: PetscTryMethod(eps,"EPSLyapIISetRanks_C",(EPS,PetscInt,PetscInt),(eps,rkc,rkl));
606: PetscFunctionReturn(PETSC_SUCCESS);
607: }
609: static PetscErrorCode EPSLyapIIGetRanks_LyapII(EPS eps,PetscInt *rkc,PetscInt *rkl)
610: {
611: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
613: PetscFunctionBegin;
614: if (rkc) *rkc = ctx->rkc;
615: if (rkl) *rkl = ctx->rkl;
616: PetscFunctionReturn(PETSC_SUCCESS);
617: }
619: /*@
620: EPSLyapIIGetRanks - Return the rank values used for the Lyapunov step.
622: Not Collective
624: Input Parameter:
625: . eps - the eigenproblem solver context
627: Output Parameters:
628: + rkc - the compressed rank
629: - rkl - the Lyapunov rank
631: Level: intermediate
633: .seealso: EPSLyapIISetRanks()
634: @*/
635: PetscErrorCode EPSLyapIIGetRanks(EPS eps,PetscInt *rkc,PetscInt *rkl)
636: {
637: PetscFunctionBegin;
639: PetscUseMethod(eps,"EPSLyapIIGetRanks_C",(EPS,PetscInt*,PetscInt*),(eps,rkc,rkl));
640: PetscFunctionReturn(PETSC_SUCCESS);
641: }
643: static PetscErrorCode EPSLyapIISetLME_LyapII(EPS eps,LME lme)
644: {
645: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
647: PetscFunctionBegin;
648: PetscCall(PetscObjectReference((PetscObject)lme));
649: PetscCall(LMEDestroy(&ctx->lme));
650: ctx->lme = lme;
651: eps->state = EPS_STATE_INITIAL;
652: PetscFunctionReturn(PETSC_SUCCESS);
653: }
655: /*@
656: EPSLyapIISetLME - Associate a linear matrix equation solver object (LME) to the
657: eigenvalue solver.
659: Collective
661: Input Parameters:
662: + eps - the eigenproblem solver context
663: - lme - the linear matrix equation solver object
665: Level: advanced
667: .seealso: EPSLyapIIGetLME()
668: @*/
669: PetscErrorCode EPSLyapIISetLME(EPS eps,LME lme)
670: {
671: PetscFunctionBegin;
674: PetscCheckSameComm(eps,1,lme,2);
675: PetscTryMethod(eps,"EPSLyapIISetLME_C",(EPS,LME),(eps,lme));
676: PetscFunctionReturn(PETSC_SUCCESS);
677: }
679: static PetscErrorCode EPSLyapIIGetLME_LyapII(EPS eps,LME *lme)
680: {
681: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
683: PetscFunctionBegin;
684: if (!ctx->lme) {
685: PetscCall(LMECreate(PetscObjectComm((PetscObject)eps),&ctx->lme));
686: PetscCall(LMESetOptionsPrefix(ctx->lme,((PetscObject)eps)->prefix));
687: PetscCall(LMEAppendOptionsPrefix(ctx->lme,"eps_lyapii_"));
688: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->lme,(PetscObject)eps,1));
689: }
690: *lme = ctx->lme;
691: PetscFunctionReturn(PETSC_SUCCESS);
692: }
694: /*@
695: EPSLyapIIGetLME - Retrieve the linear matrix equation solver object (LME)
696: associated with the eigenvalue solver.
698: Not Collective
700: Input Parameter:
701: . eps - the eigenproblem solver context
703: Output Parameter:
704: . lme - the linear matrix equation solver object
706: Level: advanced
708: .seealso: EPSLyapIISetLME()
709: @*/
710: PetscErrorCode EPSLyapIIGetLME(EPS eps,LME *lme)
711: {
712: PetscFunctionBegin;
714: PetscAssertPointer(lme,2);
715: PetscUseMethod(eps,"EPSLyapIIGetLME_C",(EPS,LME*),(eps,lme));
716: PetscFunctionReturn(PETSC_SUCCESS);
717: }
719: static PetscErrorCode EPSView_LyapII(EPS eps,PetscViewer viewer)
720: {
721: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
722: PetscBool isascii;
724: PetscFunctionBegin;
725: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
726: if (isascii) {
727: PetscCall(PetscViewerASCIIPrintf(viewer," ranks: for Lyapunov solver=%" PetscInt_FMT ", after compression=%" PetscInt_FMT "\n",ctx->rkl,ctx->rkc));
728: if (!ctx->lme) PetscCall(EPSLyapIIGetLME(eps,&ctx->lme));
729: PetscCall(PetscViewerASCIIPushTab(viewer));
730: PetscCall(LMEView(ctx->lme,viewer));
731: PetscCall(PetscViewerASCIIPopTab(viewer));
732: }
733: PetscFunctionReturn(PETSC_SUCCESS);
734: }
736: static PetscErrorCode EPSReset_LyapII(EPS eps)
737: {
738: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
740: PetscFunctionBegin;
741: if (!ctx->lme) PetscCall(LMEReset(ctx->lme));
742: PetscFunctionReturn(PETSC_SUCCESS);
743: }
745: static PetscErrorCode EPSDestroy_LyapII(EPS eps)
746: {
747: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
749: PetscFunctionBegin;
750: PetscCall(LMEDestroy(&ctx->lme));
751: PetscCall(DSDestroy(&ctx->ds));
752: PetscCall(PetscFree(eps->data));
753: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetLME_C",NULL));
754: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetLME_C",NULL));
755: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetRanks_C",NULL));
756: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetRanks_C",NULL));
757: PetscFunctionReturn(PETSC_SUCCESS);
758: }
760: static PetscErrorCode EPSSetDefaultST_LyapII(EPS eps)
761: {
762: PetscFunctionBegin;
763: if (!((PetscObject)eps->st)->type_name) PetscCall(STSetType(eps->st,STSINVERT));
764: PetscFunctionReturn(PETSC_SUCCESS);
765: }
767: SLEPC_EXTERN PetscErrorCode EPSCreate_LyapII(EPS eps)
768: {
769: EPS_LYAPII *ctx;
771: PetscFunctionBegin;
772: PetscCall(PetscNew(&ctx));
773: eps->data = (void*)ctx;
775: eps->useds = PETSC_TRUE;
777: eps->ops->solve = EPSSolve_LyapII;
778: eps->ops->setup = EPSSetUp_LyapII;
779: eps->ops->setupsort = EPSSetUpSort_Default;
780: eps->ops->setfromoptions = EPSSetFromOptions_LyapII;
781: eps->ops->reset = EPSReset_LyapII;
782: eps->ops->destroy = EPSDestroy_LyapII;
783: eps->ops->view = EPSView_LyapII;
784: eps->ops->setdefaultst = EPSSetDefaultST_LyapII;
785: eps->ops->backtransform = EPSBackTransform_Default;
786: eps->ops->computevectors = EPSComputeVectors_Schur;
788: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetLME_C",EPSLyapIISetLME_LyapII));
789: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetLME_C",EPSLyapIIGetLME_LyapII));
790: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetRanks_C",EPSLyapIISetRanks_LyapII));
791: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetRanks_C",EPSLyapIIGetRanks_LyapII));
792: PetscFunctionReturn(PETSC_SUCCESS);
793: }