Actual source code: lyapii.c
slepc-3.15.2 2021-09-20
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2021, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
10: /*
11: SLEPc eigensolver: "lyapii"
13: Method: Lyapunov inverse iteration
15: Algorithm:
17: Lyapunov inverse iteration using LME solvers
19: References:
21: [1] H.C. Elman and M. Wu, "Lyapunov inverse iteration for computing a
22: few rightmost eigenvalues of large generalized eigenvalue problems",
23: SIAM J. Matrix Anal. Appl. 34(4):1685-1707, 2013.
25: [2] K. Meerbergen and A. Spence, "Inverse iteration for purely imaginary
26: eigenvalues with application to the detection of Hopf bifurcations in
27: large-scale problems", SIAM J. Matrix Anal. Appl. 31:1982-1999, 2010.
28: */
30: #include <slepc/private/epsimpl.h>
31: #include <slepcblaslapack.h>
33: typedef struct {
34: LME lme; /* Lyapunov solver */
35: DS ds; /* used to compute the SVD for compression */
36: PetscInt rkl; /* prescribed rank for the Lyapunov solver */
37: PetscInt rkc; /* the compressed rank, cannot be larger than rkl */
38: } EPS_LYAPII;
40: typedef struct {
41: Mat S; /* the operator matrix, S=A^{-1}*B */
42: BV Q; /* orthogonal basis of converged eigenvectors */
43: } EPS_LYAPII_MATSHELL;
45: typedef struct {
46: Mat S; /* the matrix from which the implicit operator is built */
47: PetscInt n; /* the size of matrix S, the operator is nxn */
48: LME lme; /* dummy LME object */
49: #if defined(PETSC_USE_COMPLEX)
50: Mat A,B,F;
51: Vec w;
52: #endif
53: } EPS_EIG_MATSHELL;
55: PetscErrorCode EPSSetUp_LyapII(EPS eps)
56: {
58: PetscRandom rand;
59: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
62: EPSCheckSinvert(eps);
63: if (eps->ncv!=PETSC_DEFAULT) {
64: if (eps->ncv<eps->nev+1) SETERRQ(PetscObjectComm((PetscObject)eps),1,"The value of ncv must be at least nev+1");
65: } else eps->ncv = eps->nev+1;
66: if (eps->mpd!=PETSC_DEFAULT) { PetscInfo(eps,"Warning: parameter mpd ignored\n"); }
67: if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(1000*eps->nev,100*eps->n);
68: if (!eps->which) eps->which=EPS_LARGEST_REAL;
69: if (eps->which!=EPS_LARGEST_REAL) SETERRQ(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) { EPSLyapIIGetLME(eps,&ctx->lme); }
75: LMESetProblemType(ctx->lme,LME_LYAPUNOV);
76: LMESetErrorIfNotConverged(ctx->lme,PETSC_TRUE);
78: if (!ctx->ds) {
79: DSCreate(PetscObjectComm((PetscObject)eps),&ctx->ds);
80: PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->ds);
81: DSSetType(ctx->ds,DSSVD);
82: }
83: DSAllocate(ctx->ds,ctx->rkl);
85: DSSetType(eps->ds,DSNHEP);
86: DSAllocate(eps->ds,eps->ncv);
88: EPSAllocateSolution(eps,0);
89: BVGetRandomContext(eps->V,&rand); /* make sure the random context is available when duplicating */
90: EPSSetWorkVecs(eps,3);
91: return(0);
92: }
94: static PetscErrorCode MatMult_EPSLyapIIOperator(Mat M,Vec x,Vec r)
95: {
96: PetscErrorCode ierr;
97: EPS_LYAPII_MATSHELL *matctx;
100: MatShellGetContext(M,(void**)&matctx);
101: MatMult(matctx->S,x,r);
102: BVOrthogonalizeVec(matctx->Q,r,NULL,NULL,NULL);
103: return(0);
104: }
106: static PetscErrorCode MatDestroy_EPSLyapIIOperator(Mat M)
107: {
108: PetscErrorCode ierr;
109: EPS_LYAPII_MATSHELL *matctx;
112: MatShellGetContext(M,(void**)&matctx);
113: MatDestroy(&matctx->S);
114: PetscFree(matctx);
115: return(0);
116: }
118: static PetscErrorCode MatMult_EigOperator(Mat M,Vec x,Vec y)
119: {
120: PetscErrorCode ierr;
121: EPS_EIG_MATSHELL *matctx;
122: #if !defined(PETSC_USE_COMPLEX)
123: PetscInt n;
124: PetscScalar *Y,*C,zero=0.0,done=1.0,dtwo=2.0;
125: const PetscScalar *S,*X;
126: PetscBLASInt n_;
127: #endif
130: MatShellGetContext(M,(void**)&matctx);
132: #if defined(PETSC_USE_COMPLEX)
133: MatMult(matctx->B,x,matctx->w);
134: MatSolve(matctx->F,matctx->w,y);
135: #else
136: VecGetArrayRead(x,&X);
137: VecGetArray(y,&Y);
138: MatDenseGetArrayRead(matctx->S,&S);
140: n = matctx->n;
141: PetscCalloc1(n*n,&C);
142: PetscBLASIntCast(n,&n_);
144: /* C = 2*S*X*S.' */
145: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&dtwo,S,&n_,X,&n_,&zero,Y,&n_));
146: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&n_,&n_,&n_,&done,Y,&n_,S,&n_,&zero,C,&n_));
148: /* Solve S*Y + Y*S' = -C */
149: LMEDenseLyapunov(matctx->lme,n,(PetscScalar*)S,n,C,n,Y,n);
151: PetscFree(C);
152: VecRestoreArrayRead(x,&X);
153: VecRestoreArray(y,&Y);
154: MatDenseRestoreArrayRead(matctx->S,&S);
155: #endif
156: return(0);
157: }
159: static PetscErrorCode MatDestroy_EigOperator(Mat M)
160: {
161: PetscErrorCode ierr;
162: EPS_EIG_MATSHELL *matctx;
165: MatShellGetContext(M,(void**)&matctx);
166: #if defined(PETSC_USE_COMPLEX)
167: MatDestroy(&matctx->A);
168: MatDestroy(&matctx->B);
169: MatDestroy(&matctx->F);
170: VecDestroy(&matctx->w);
171: #endif
172: PetscFree(matctx);
173: return(0);
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: {
182: PetscBLASInt lwork=10,ld_;
183: #if !defined(PETSC_HAVE_ESSL)
184: PetscScalar work[10];
185: PetscBLASInt two=2,info;
186: #else
187: PetscInt i;
188: PetscBLASInt idummy,io=1;
189: PetscScalar wri[4];
190: #endif
191: #if defined(PETSC_HAVE_ESSL) || defined(PETSC_USE_COMPLEX)
192: PetscReal rwork[6];
193: #endif
196: PetscBLASIntCast(ld,&ld_);
197: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
198: #if !defined(PETSC_HAVE_ESSL)
199: #if !defined(PETSC_USE_COMPLEX)
200: PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","V",&two,M,&ld_,wr,wi,NULL,&ld_,vec,&ld_,work,&lwork,&info));
201: #else
202: PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","V",&two,M,&ld_,wr,NULL,&ld_,vec,&ld_,work,&lwork,rwork,&info));
203: #endif
204: SlepcCheckLapackInfo("geev",info);
205: #else /* defined(PETSC_HAVE_ESSL) */
206: PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_(&io,M,&ld_,wri,vec,&ld_,&idummy,&ld_,rwork,&lwork));
207: #if !defined(PETSC_USE_COMPLEX)
208: for (i=0;i<2;i++) {
209: wr[i] = wri[2*i];
210: wi[i] = wri[2*i+1];
211: }
212: #else
213: for (i=0;i<2;i++) wr[i] = wri[i];
214: #endif
215: #endif
216: PetscFPTrapPop();
217: return(0);
218: }
220: /*
221: LyapIIBuildRHS: prepare the right-hand side of the Lyapunov equation SY + YS' = -2*S*Z*S'
222: in factored form:
223: if (V) U=sqrt(2)*S*V (uses 1 work vector)
224: else U=sqrt(2)*S*U (uses 2 work vectors)
225: where U,V are assumed to have rk columns.
226: */
227: static PetscErrorCode LyapIIBuildRHS(Mat S,PetscInt rk,Mat U,BV V,Vec *work)
228: {
230: PetscScalar *array,*uu;
231: PetscInt i,nloc;
232: Vec v,u=work[0];
235: MatGetLocalSize(U,&nloc,NULL);
236: for (i=0;i<rk;i++) {
237: MatDenseGetColumn(U,i,&array);
238: if (V) {
239: BVGetColumn(V,i,&v);
240: } else {
241: v = work[1];
242: VecPlaceArray(v,array);
243: }
244: MatMult(S,v,u);
245: if (V) {
246: BVRestoreColumn(V,i,&v);
247: } else {
248: VecResetArray(v);
249: }
250: VecScale(u,PETSC_SQRT2);
251: VecGetArray(u,&uu);
252: PetscArraycpy(array,uu,nloc);
253: VecRestoreArray(u,&uu);
254: MatDenseRestoreColumn(U,&array);
255: }
256: return(0);
257: }
259: /*
260: LyapIIBuildEigenMat: create shell matrix Op=A\B with A = kron(I,S)+kron(S,I), B = -2*kron(S,S)
261: where S is a sequential square dense matrix of order n.
262: v0 is the initial vector, should have the form v0 = w*w' (for instance 1*1')
263: */
264: static PetscErrorCode LyapIIBuildEigenMat(LME lme,Mat S,Mat *Op,Vec *v0)
265: {
266: PetscErrorCode ierr;
267: PetscInt n,m;
268: PetscBool create=PETSC_FALSE;
269: EPS_EIG_MATSHELL *matctx;
270: #if defined(PETSC_USE_COMPLEX)
271: PetscScalar theta,*aa,*bb;
272: const PetscScalar *ss;
273: PetscInt i,j,f,c,off,ld;
274: IS perm;
275: #endif
278: MatGetSize(S,&n,NULL);
279: if (!*Op) create=PETSC_TRUE;
280: else {
281: MatGetSize(*Op,&m,NULL);
282: if (m!=n*n) create=PETSC_TRUE;
283: }
284: if (create) {
285: MatDestroy(Op);
286: VecDestroy(v0);
287: PetscNew(&matctx);
288: #if defined(PETSC_USE_COMPLEX)
289: MatCreateSeqDense(PETSC_COMM_SELF,n*n,n*n,NULL,&matctx->A);
290: MatCreateSeqDense(PETSC_COMM_SELF,n*n,n*n,NULL,&matctx->B);
291: MatCreateVecs(matctx->A,NULL,&matctx->w);
292: #endif
293: MatCreateShell(PETSC_COMM_SELF,n*n,n*n,PETSC_DETERMINE,PETSC_DETERMINE,matctx,Op);
294: MatShellSetOperation(*Op,MATOP_MULT,(void(*)(void))MatMult_EigOperator);
295: MatShellSetOperation(*Op,MATOP_DESTROY,(void(*)(void))MatDestroy_EigOperator);
296: MatCreateVecs(*Op,NULL,v0);
297: } else {
298: MatShellGetContext(*Op,(void**)&matctx);
299: #if defined(PETSC_USE_COMPLEX)
300: MatZeroEntries(matctx->A);
301: #endif
302: }
303: #if defined(PETSC_USE_COMPLEX)
304: MatDenseGetArray(matctx->A,&aa);
305: MatDenseGetArray(matctx->B,&bb);
306: MatDenseGetArrayRead(S,&ss);
307: ld = n*n;
308: for (f=0;f<n;f++) {
309: off = f*n+f*n*ld;
310: for (i=0;i<n;i++) for (j=0;j<n;j++) aa[off+i+j*ld] = ss[i+j*n];
311: for (c=0;c<n;c++) {
312: off = f*n+c*n*ld;
313: theta = ss[f+c*n];
314: for (i=0;i<n;i++) aa[off+i+i*ld] += theta;
315: for (i=0;i<n;i++) for (j=0;j<n;j++) bb[off+i+j*ld] = -2*theta*ss[i+j*n];
316: }
317: }
318: MatDenseRestoreArray(matctx->A,&aa);
319: MatDenseRestoreArray(matctx->B,&bb);
320: MatDenseRestoreArrayRead(S,&ss);
321: ISCreateStride(PETSC_COMM_SELF,n*n,0,1,&perm);
322: MatDestroy(&matctx->F);
323: MatDuplicate(matctx->A,MAT_COPY_VALUES,&matctx->F);
324: MatLUFactor(matctx->F,perm,perm,0);
325: ISDestroy(&perm);
326: #endif
327: matctx->lme = lme;
328: matctx->S = S;
329: matctx->n = n;
330: VecSet(*v0,1.0);
331: return(0);
332: }
334: PetscErrorCode EPSSolve_LyapII(EPS eps)
335: {
336: PetscErrorCode ierr;
337: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
338: PetscInt i,ldds,rk,nloc,mloc,nv,idx,k;
339: Vec v,w,z=eps->work[0],v0=NULL;
340: Mat S,C,Ux[2],Y,Y1,R,U,W,X,Op=NULL;
341: BV V;
342: BVOrthogType type;
343: BVOrthogRefineType refine;
344: PetscScalar eigr[2],eigi[2],*array,er,ei,*uu,*s,*xx,*aa,pM[4],vec[4];
345: PetscReal eta;
346: EPS epsrr;
347: PetscReal norm;
348: EPS_LYAPII_MATSHELL *matctx;
351: DSGetLeadingDimension(ctx->ds,&ldds);
353: /* Operator for the Lyapunov equation */
354: PetscNew(&matctx);
355: STGetOperator(eps->st,&matctx->S);
356: MatGetLocalSize(matctx->S,&mloc,&nloc);
357: MatCreateShell(PetscObjectComm((PetscObject)eps),mloc,nloc,PETSC_DETERMINE,PETSC_DETERMINE,matctx,&S);
358: matctx->Q = eps->V;
359: MatShellSetOperation(S,MATOP_MULT,(void(*)(void))MatMult_EPSLyapIIOperator);
360: MatShellSetOperation(S,MATOP_DESTROY,(void(*)(void))MatDestroy_EPSLyapIIOperator);
361: LMESetCoefficients(ctx->lme,S,NULL,NULL,NULL);
363: /* Right-hand side */
364: BVDuplicateResize(eps->V,ctx->rkl,&V);
365: BVGetOrthogonalization(V,&type,&refine,&eta,NULL);
366: BVSetOrthogonalization(V,type,refine,eta,BV_ORTHOG_BLOCK_TSQR);
367: MatCreateDense(PetscObjectComm((PetscObject)eps),eps->nloc,PETSC_DECIDE,PETSC_DECIDE,1,NULL,&Ux[0]);
368: MatCreateDense(PetscObjectComm((PetscObject)eps),eps->nloc,PETSC_DECIDE,PETSC_DECIDE,2,NULL,&Ux[1]);
369: nv = ctx->rkl;
370: PetscMalloc1(nv,&s);
372: /* Initialize first column */
373: EPSGetStartVector(eps,0,NULL);
374: BVGetColumn(eps->V,0,&v);
375: BVInsertVec(V,0,v);
376: BVRestoreColumn(eps->V,0,&v);
377: BVSetActiveColumns(eps->V,0,0); /* no deflation at the beginning */
378: LyapIIBuildRHS(S,1,Ux[0],V,eps->work);
379: idx = 0;
381: /* EPS for rank reduction */
382: EPSCreate(PETSC_COMM_SELF,&epsrr);
383: EPSSetOptionsPrefix(epsrr,((PetscObject)eps)->prefix);
384: EPSAppendOptionsPrefix(epsrr,"eps_lyapii_");
385: EPSSetDimensions(epsrr,1,PETSC_DEFAULT,PETSC_DEFAULT);
386: EPSSetTolerances(epsrr,PETSC_MACHINE_EPSILON*100,PETSC_DEFAULT);
388: while (eps->reason == EPS_CONVERGED_ITERATING) {
389: eps->its++;
391: /* Matrix for placing the solution of the Lyapunov equation (an alias of V) */
392: BVSetActiveColumns(V,0,nv);
393: BVGetMat(V,&Y1);
394: MatZeroEntries(Y1);
395: MatCreateLRC(NULL,Y1,NULL,NULL,&Y);
396: LMESetSolution(ctx->lme,Y);
398: /* Solve the Lyapunov equation SY + YS' = -2*S*Z*S' */
399: MatCreateLRC(NULL,Ux[idx],NULL,NULL,&C);
400: LMESetRHS(ctx->lme,C);
401: MatDestroy(&C);
402: LMESolve(ctx->lme);
403: BVRestoreMat(V,&Y1);
404: MatDestroy(&Y);
406: /* SVD of the solution: [Q,R]=qr(V); [U,Sigma,~]=svd(R) */
407: DSSetDimensions(ctx->ds,nv,nv,0,0);
408: DSGetMat(ctx->ds,DS_MAT_A,&R);
409: BVOrthogonalize(V,R);
410: DSRestoreMat(ctx->ds,DS_MAT_A,&R);
411: DSSetState(ctx->ds,DS_STATE_RAW);
412: DSSolve(ctx->ds,s,NULL);
414: /* Determine rank */
415: rk = nv;
416: for (i=1;i<nv;i++) if (PetscAbsScalar(s[i]/s[0])<PETSC_SQRT_MACHINE_EPSILON) {rk=i; break;}
417: PetscInfo1(eps,"The computed solution of the Lyapunov equation has rank %D\n",rk);
418: rk = PetscMin(rk,ctx->rkc);
419: DSGetMat(ctx->ds,DS_MAT_U,&U);
420: BVMultInPlace(V,U,0,rk);
421: BVSetActiveColumns(V,0,rk);
422: MatDestroy(&U);
424: /* Rank reduction */
425: DSSetDimensions(ctx->ds,rk,rk,0,0);
426: DSGetMat(ctx->ds,DS_MAT_A,&W);
427: BVMatProject(V,S,V,W);
428: LyapIIBuildEigenMat(ctx->lme,W,&Op,&v0); /* Op=A\B, A=kron(I,S)+kron(S,I), B=-2*kron(S,S) */
429: EPSSetOperators(epsrr,Op,NULL);
430: EPSSetInitialSpace(epsrr,1,&v0);
431: EPSSolve(epsrr);
432: MatDestroy(&W);
433: EPSComputeVectors(epsrr);
434: /* Copy first eigenvector, vec(A)=x */
435: BVGetArray(epsrr->V,&xx);
436: DSGetArray(ctx->ds,DS_MAT_A,&aa);
437: for (i=0;i<rk;i++) {
438: PetscArraycpy(aa+i*ldds,xx+i*rk,rk);
439: }
440: DSRestoreArray(ctx->ds,DS_MAT_A,&aa);
441: BVRestoreArray(epsrr->V,&xx);
442: DSSetState(ctx->ds,DS_STATE_RAW);
443: /* Compute [U,Sigma,~] = svd(A), its rank should be 1 or 2 */
444: DSSolve(ctx->ds,s,NULL);
445: if (PetscAbsScalar(s[1]/s[0])<PETSC_SQRT_MACHINE_EPSILON) rk=1;
446: else rk = 2;
447: PetscInfo1(eps,"The eigenvector has rank %D\n",rk);
448: DSGetMat(ctx->ds,DS_MAT_U,&U);
449: BVMultInPlace(V,U,0,rk);
450: MatDestroy(&U);
452: /* Save V in Ux */
453: idx = (rk==2)?1:0;
454: for (i=0;i<rk;i++) {
455: BVGetColumn(V,i,&v);
456: VecGetArray(v,&uu);
457: MatDenseGetColumn(Ux[idx],i,&array);
458: PetscArraycpy(array,uu,eps->nloc);
459: MatDenseRestoreColumn(Ux[idx],&array);
460: VecRestoreArray(v,&uu);
461: BVRestoreColumn(V,i,&v);
462: }
464: /* Eigenpair approximation */
465: BVGetColumn(V,0,&v);
466: MatMult(S,v,z);
467: VecDot(z,v,pM);
468: BVRestoreColumn(V,0,&v);
469: if (rk>1) {
470: BVGetColumn(V,1,&w);
471: VecDot(z,w,pM+1);
472: MatMult(S,w,z);
473: VecDot(z,w,pM+3);
474: BVGetColumn(V,0,&v);
475: VecDot(z,v,pM+2);
476: BVRestoreColumn(V,0,&v);
477: BVRestoreColumn(V,1,&w);
478: EV2x2(pM,2,eigr,eigi,vec);
479: MatCreateSeqDense(PETSC_COMM_SELF,2,2,vec,&X);
480: BVSetActiveColumns(V,0,rk);
481: BVMultInPlace(V,X,0,rk);
482: MatDestroy(&X);
483: #if !defined(PETSC_USE_COMPLEX)
484: norm = eigr[0]*eigr[0]+eigi[0]*eigi[0];
485: er = eigr[0]/norm; ei = -eigi[0]/norm;
486: #else
487: er =1.0/eigr[0]; ei = 0.0;
488: #endif
489: } else {
490: eigr[0] = pM[0]; eigi[0] = 0.0;
491: er = 1.0/eigr[0]; ei = 0.0;
492: }
493: BVGetColumn(V,0,&v);
494: if (eigi[0]!=0.0) {
495: BVGetColumn(V,1,&w);
496: } else w = NULL;
497: eps->eigr[eps->nconv] = eigr[0]; eps->eigi[eps->nconv] = eigi[0];
498: EPSComputeResidualNorm_Private(eps,PETSC_FALSE,er,ei,v,w,eps->work,&norm);
499: BVRestoreColumn(V,0,&v);
500: if (w) {
501: BVRestoreColumn(V,1,&w);
502: }
503: (*eps->converged)(eps,er,ei,norm,&eps->errest[eps->nconv],eps->convergedctx);
504: k = 0;
505: if (eps->errest[eps->nconv]<eps->tol) {
506: k++;
507: if (rk==2) {
508: #if !defined (PETSC_USE_COMPLEX)
509: eps->eigr[eps->nconv+k] = eigr[0]; eps->eigi[eps->nconv+k] = -eigi[0];
510: #else
511: eps->eigr[eps->nconv+k] = PetscConj(eps->eigr[eps->nconv]);
512: #endif
513: k++;
514: }
515: /* Store converged eigenpairs and vectors for deflation */
516: for (i=0;i<k;i++) {
517: BVGetColumn(V,i,&v);
518: BVInsertVec(eps->V,eps->nconv+i,v);
519: BVRestoreColumn(V,i,&v);
520: }
521: eps->nconv += k;
522: BVSetActiveColumns(eps->V,eps->nconv-rk,eps->nconv);
523: BVOrthogonalize(eps->V,NULL);
524: DSSetDimensions(eps->ds,eps->nconv,0,0,0);
525: DSGetMat(eps->ds,DS_MAT_A,&W);
526: BVMatProject(eps->V,matctx->S,eps->V,W);
527: DSRestoreMat(eps->ds,DS_MAT_A,&W);
528: if (eps->nconv<eps->nev) {
529: idx = 0;
530: BVSetRandomColumn(V,0);
531: BVNormColumn(V,0,NORM_2,&norm);
532: BVScaleColumn(V,0,1.0/norm);
533: LyapIIBuildRHS(S,1,Ux[idx],V,eps->work);
534: }
535: } else {
536: /* Prepare right-hand side */
537: LyapIIBuildRHS(S,rk,Ux[idx],NULL,eps->work);
538: }
539: (*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx);
540: EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,eps->nconv+1);
541: }
542: STRestoreOperator(eps->st,&matctx->S);
543: MatDestroy(&S);
544: MatDestroy(&Ux[0]);
545: MatDestroy(&Ux[1]);
546: MatDestroy(&Op);
547: VecDestroy(&v0);
548: BVDestroy(&V);
549: EPSDestroy(&epsrr);
550: PetscFree(s);
551: return(0);
552: }
554: PetscErrorCode EPSSetFromOptions_LyapII(PetscOptionItems *PetscOptionsObject,EPS eps)
555: {
557: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
558: PetscInt k,array[2]={PETSC_DEFAULT,PETSC_DEFAULT};
559: PetscBool flg;
562: PetscOptionsHead(PetscOptionsObject,"EPS Lyapunov Inverse Iteration Options");
564: k = 2;
565: PetscOptionsIntArray("-eps_lyapii_ranks","Ranks for Lyapunov equation (one or two comma-separated integers)","EPSLyapIISetRanks",array,&k,&flg);
566: if (flg) {
567: EPSLyapIISetRanks(eps,array[0],array[1]);
568: }
570: PetscOptionsTail();
572: if (!ctx->lme) { EPSLyapIIGetLME(eps,&ctx->lme); }
573: LMESetFromOptions(ctx->lme);
574: return(0);
575: }
577: static PetscErrorCode EPSLyapIISetRanks_LyapII(EPS eps,PetscInt rkc,PetscInt rkl)
578: {
579: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
582: if (rkc==PETSC_DEFAULT) rkc = 10;
583: if (rkc<2) SETERRQ1(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The compressed rank %D must be larger than 1",rkc);
584: if (rkl==PETSC_DEFAULT) rkl = 3*rkc;
585: if (rkl<rkc) SETERRQ2(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The Lyapunov rank %D cannot be smaller than the compressed rank %D",rkl,rkc);
586: if (rkc != ctx->rkc) {
587: ctx->rkc = rkc;
588: eps->state = EPS_STATE_INITIAL;
589: }
590: if (rkl != ctx->rkl) {
591: ctx->rkl = rkl;
592: eps->state = EPS_STATE_INITIAL;
593: }
594: return(0);
595: }
597: /*@
598: EPSLyapIISetRanks - Set the ranks used in the solution of the Lyapunov equation.
600: Logically Collective on EPS
602: Input Parameters:
603: + eps - the eigenproblem solver context
604: . rkc - the compressed rank
605: - rkl - the Lyapunov rank
607: Options Database Key:
608: . -eps_lyapii_ranks <rkc,rkl> - Sets the rank parameters
610: Notes:
611: Lyapunov inverse iteration needs to solve a large-scale Lyapunov equation
612: at each iteration of the eigensolver. For this, an iterative solver (LME)
613: is used, which requires to prescribe the rank of the solution matrix X. This
614: is the meaning of parameter rkl. Later, this matrix is compressed into
615: another matrix of rank rkc. If not provided, rkl is a small multiple of rkc.
617: Level: intermediate
619: .seealso: EPSLyapIIGetRanks()
620: @*/
621: PetscErrorCode EPSLyapIISetRanks(EPS eps,PetscInt rkc,PetscInt rkl)
622: {
629: PetscTryMethod(eps,"EPSLyapIISetRanks_C",(EPS,PetscInt,PetscInt),(eps,rkc,rkl));
630: return(0);
631: }
633: static PetscErrorCode EPSLyapIIGetRanks_LyapII(EPS eps,PetscInt *rkc,PetscInt *rkl)
634: {
635: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
638: if (rkc) *rkc = ctx->rkc;
639: if (rkl) *rkl = ctx->rkl;
640: return(0);
641: }
643: /*@
644: EPSLyapIIGetRanks - Return the rank values used for the Lyapunov step.
646: Not Collective
648: Input Parameter:
649: . eps - the eigenproblem solver context
651: Output Parameters:
652: + rkc - the compressed rank
653: - rkl - the Lyapunov rank
655: Level: intermediate
657: .seealso: EPSLyapIISetRanks()
658: @*/
659: PetscErrorCode EPSLyapIIGetRanks(EPS eps,PetscInt *rkc,PetscInt *rkl)
660: {
665: PetscUseMethod(eps,"EPSLyapIIGetRanks_C",(EPS,PetscInt*,PetscInt*),(eps,rkc,rkl));
666: return(0);
667: }
669: static PetscErrorCode EPSLyapIISetLME_LyapII(EPS eps,LME lme)
670: {
672: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
675: PetscObjectReference((PetscObject)lme);
676: LMEDestroy(&ctx->lme);
677: ctx->lme = lme;
678: PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->lme);
679: eps->state = EPS_STATE_INITIAL;
680: return(0);
681: }
683: /*@
684: EPSLyapIISetLME - Associate a linear matrix equation solver object (LME) to the
685: eigenvalue solver.
687: Collective on EPS
689: Input Parameters:
690: + eps - the eigenproblem solver context
691: - lme - the linear matrix equation solver object
693: Level: advanced
695: .seealso: EPSLyapIIGetLME()
696: @*/
697: PetscErrorCode EPSLyapIISetLME(EPS eps,LME lme)
698: {
705: PetscTryMethod(eps,"EPSLyapIISetLME_C",(EPS,LME),(eps,lme));
706: return(0);
707: }
709: static PetscErrorCode EPSLyapIIGetLME_LyapII(EPS eps,LME *lme)
710: {
712: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
715: if (!ctx->lme) {
716: LMECreate(PetscObjectComm((PetscObject)eps),&ctx->lme);
717: LMESetOptionsPrefix(ctx->lme,((PetscObject)eps)->prefix);
718: LMEAppendOptionsPrefix(ctx->lme,"eps_lyapii_");
719: PetscObjectIncrementTabLevel((PetscObject)ctx->lme,(PetscObject)eps,1);
720: PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->lme);
721: }
722: *lme = ctx->lme;
723: return(0);
724: }
726: /*@
727: EPSLyapIIGetLME - Retrieve the linear matrix equation solver object (LME)
728: associated with the eigenvalue solver.
730: Not Collective
732: Input Parameter:
733: . eps - the eigenproblem solver context
735: Output Parameter:
736: . lme - the linear matrix equation solver object
738: Level: advanced
740: .seealso: EPSLyapIISetLME()
741: @*/
742: PetscErrorCode EPSLyapIIGetLME(EPS eps,LME *lme)
743: {
749: PetscUseMethod(eps,"EPSLyapIIGetLME_C",(EPS,LME*),(eps,lme));
750: return(0);
751: }
753: PetscErrorCode EPSView_LyapII(EPS eps,PetscViewer viewer)
754: {
756: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
757: PetscBool isascii;
760: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
761: if (isascii) {
762: PetscViewerASCIIPrintf(viewer," ranks: for Lyapunov solver=%D, after compression=%D\n",ctx->rkl,ctx->rkc);
763: if (!ctx->lme) { EPSLyapIIGetLME(eps,&ctx->lme); }
764: PetscViewerASCIIPushTab(viewer);
765: LMEView(ctx->lme,viewer);
766: PetscViewerASCIIPopTab(viewer);
767: }
768: return(0);
769: }
771: PetscErrorCode EPSReset_LyapII(EPS eps)
772: {
774: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
777: if (!ctx->lme) { LMEReset(ctx->lme); }
778: return(0);
779: }
781: PetscErrorCode EPSDestroy_LyapII(EPS eps)
782: {
784: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
787: LMEDestroy(&ctx->lme);
788: DSDestroy(&ctx->ds);
789: PetscFree(eps->data);
790: PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetLME_C",NULL);
791: PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetLME_C",NULL);
792: PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetRanks_C",NULL);
793: PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetRanks_C",NULL);
794: return(0);
795: }
797: PetscErrorCode EPSSetDefaultST_LyapII(EPS eps)
798: {
802: if (!((PetscObject)eps->st)->type_name) {
803: STSetType(eps->st,STSINVERT);
804: }
805: return(0);
806: }
808: SLEPC_EXTERN PetscErrorCode EPSCreate_LyapII(EPS eps)
809: {
810: EPS_LYAPII *ctx;
814: PetscNewLog(eps,&ctx);
815: eps->data = (void*)ctx;
817: eps->useds = PETSC_TRUE;
819: eps->ops->solve = EPSSolve_LyapII;
820: eps->ops->setup = EPSSetUp_LyapII;
821: eps->ops->setupsort = EPSSetUpSort_Default;
822: eps->ops->setfromoptions = EPSSetFromOptions_LyapII;
823: eps->ops->reset = EPSReset_LyapII;
824: eps->ops->destroy = EPSDestroy_LyapII;
825: eps->ops->view = EPSView_LyapII;
826: eps->ops->setdefaultst = EPSSetDefaultST_LyapII;
827: eps->ops->backtransform = EPSBackTransform_Default;
828: eps->ops->computevectors = EPSComputeVectors_Schur;
830: PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetLME_C",EPSLyapIISetLME_LyapII);
831: PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetLME_C",EPSLyapIIGetLME_LyapII);
832: PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetRanks_C",EPSLyapIISetRanks_LyapII);
833: PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetRanks_C",EPSLyapIIGetRanks_LyapII);
834: return(0);
835: }