Actual source code: lyapii.c
slepc-3.14.1 2020-12-08
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2020, 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: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
61: EPSCheckSinvert(eps);
62: if (eps->ncv!=PETSC_DEFAULT) {
63: if (eps->ncv<eps->nev+1) SETERRQ(PetscObjectComm((PetscObject)eps),1,"The value of ncv must be at least nev+1");
64: } else eps->ncv = eps->nev+1;
65: if (eps->mpd!=PETSC_DEFAULT) { 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: if (eps->which!=EPS_LARGEST_REAL) SETERRQ(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) { EPSLyapIIGetLME(eps,&ctx->lme); }
74: LMESetProblemType(ctx->lme,LME_LYAPUNOV);
75: LMESetErrorIfNotConverged(ctx->lme,PETSC_TRUE);
77: if (!ctx->ds) {
78: DSCreate(PetscObjectComm((PetscObject)eps),&ctx->ds);
79: PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->ds);
80: DSSetType(ctx->ds,DSSVD);
81: }
82: DSAllocate(ctx->ds,ctx->rkl);
84: DSSetType(eps->ds,DSNHEP);
85: DSAllocate(eps->ds,eps->ncv);
87: EPSAllocateSolution(eps,0);
88: EPSSetWorkVecs(eps,3);
89: return(0);
90: }
92: static PetscErrorCode MatMult_EPSLyapIIOperator(Mat M,Vec x,Vec r)
93: {
94: PetscErrorCode ierr;
95: EPS_LYAPII_MATSHELL *matctx;
98: MatShellGetContext(M,(void**)&matctx);
99: MatMult(matctx->S,x,r);
100: BVOrthogonalizeVec(matctx->Q,r,NULL,NULL,NULL);
101: return(0);
102: }
104: static PetscErrorCode MatDestroy_EPSLyapIIOperator(Mat M)
105: {
106: PetscErrorCode ierr;
107: EPS_LYAPII_MATSHELL *matctx;
110: MatShellGetContext(M,(void**)&matctx);
111: MatDestroy(&matctx->S);
112: PetscFree(matctx);
113: return(0);
114: }
116: static PetscErrorCode MatMult_EigOperator(Mat M,Vec x,Vec y)
117: {
118: PetscErrorCode ierr;
119: EPS_EIG_MATSHELL *matctx;
120: #if !defined(PETSC_USE_COMPLEX)
121: PetscInt n;
122: PetscScalar *S,*Y,*C,zero=0.0,done=1.0,dtwo=2.0;
123: const PetscScalar *X;
124: PetscBLASInt n_;
125: #endif
128: MatShellGetContext(M,(void**)&matctx);
130: #if defined(PETSC_USE_COMPLEX)
131: MatMult(matctx->B,x,matctx->w);
132: MatSolve(matctx->F,matctx->w,y);
133: #else
134: VecGetArrayRead(x,&X);
135: VecGetArray(y,&Y);
136: MatDenseGetArray(matctx->S,&S);
138: n = matctx->n;
139: PetscCalloc1(n*n,&C);
140: PetscBLASIntCast(n,&n_);
142: /* C = 2*S*X*S.' */
143: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&dtwo,S,&n_,X,&n_,&zero,Y,&n_));
144: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&n_,&n_,&n_,&done,Y,&n_,S,&n_,&zero,C,&n_));
146: /* Solve S*Y + Y*S' = -C */
147: LMEDenseLyapunov(matctx->lme,n,S,n,C,n,Y,n);
149: PetscFree(C);
150: VecRestoreArrayRead(x,&X);
151: VecRestoreArray(y,&Y);
152: MatDenseRestoreArray(matctx->S,&S);
153: #endif
154: return(0);
155: }
157: static PetscErrorCode MatDestroy_EigOperator(Mat M)
158: {
159: PetscErrorCode ierr;
160: EPS_EIG_MATSHELL *matctx;
163: MatShellGetContext(M,(void**)&matctx);
164: #if defined(PETSC_USE_COMPLEX)
165: MatDestroy(&matctx->A);
166: MatDestroy(&matctx->B);
167: MatDestroy(&matctx->F);
168: VecDestroy(&matctx->w);
169: #endif
170: PetscFree(matctx);
171: return(0);
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: {
180: PetscBLASInt lwork=10,ld_;
181: #if !defined(PETSC_HAVE_ESSL)
182: PetscScalar work[10];
183: PetscBLASInt two=2,info;
184: #else
185: PetscInt i;
186: PetscBLASInt idummy,io=1;
187: PetscScalar wri[4];
188: #endif
189: #if defined(PETSC_HAVE_ESSL) || defined(PETSC_USE_COMPLEX)
190: PetscReal rwork[6];
191: #endif
194: PetscBLASIntCast(ld,&ld_);
195: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
196: #if !defined(PETSC_HAVE_ESSL)
197: #if !defined(PETSC_USE_COMPLEX)
198: PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","V",&two,M,&ld_,wr,wi,NULL,&ld_,vec,&ld_,work,&lwork,&info));
199: #else
200: PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","V",&two,M,&ld_,wr,NULL,&ld_,vec,&ld_,work,&lwork,rwork,&info));
201: #endif
202: SlepcCheckLapackInfo("geev",info);
203: #else /* defined(PETSC_HAVE_ESSL) */
204: PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_(&io,M,&ld_,wri,vec,&ld_,&idummy,&ld_,rwork,&lwork));
205: #if !defined(PETSC_USE_COMPLEX)
206: for (i=0;i<2;i++) {
207: wr[i] = wri[2*i];
208: wi[i] = wri[2*i+1];
209: }
210: #else
211: for (i=0;i<2;i++) wr[i] = wri[i];
212: #endif
213: #endif
214: PetscFPTrapPop();
215: return(0);
216: }
218: /*
219: LyapIIBuildRHS: prepare the right-hand side of the Lyapunov equation SY + YS' = -2*S*Z*S'
220: in factored form:
221: if (V) U=sqrt(2)*S*V (uses 1 work vector)
222: else U=sqrt(2)*S*U (uses 2 work vectors)
223: where U,V are assumed to have rk columns.
224: */
225: static PetscErrorCode LyapIIBuildRHS(Mat S,PetscInt rk,Mat U,BV V,Vec *work)
226: {
228: PetscScalar *array,*uu;
229: PetscInt i,nloc;
230: Vec v,u=work[0];
233: MatGetLocalSize(U,&nloc,NULL);
234: for (i=0;i<rk;i++) {
235: MatDenseGetColumn(U,i,&array);
236: if (V) {
237: BVGetColumn(V,i,&v);
238: } else {
239: v = work[1];
240: VecPlaceArray(v,array);
241: }
242: MatMult(S,v,u);
243: if (V) {
244: BVRestoreColumn(V,i,&v);
245: } else {
246: VecResetArray(v);
247: }
248: VecScale(u,PetscSqrtReal(2.0));
249: VecGetArray(u,&uu);
250: PetscMemcpy(array,uu,nloc*sizeof(PetscScalar));
251: VecRestoreArray(u,&uu);
252: MatDenseRestoreColumn(U,&array);
253: }
254: return(0);
255: }
257: /*
258: LyapIIBuildEigenMat: create shell matrix Op=A\B with A = kron(I,S)+kron(S,I), B = -2*kron(S,S)
259: where S is a sequential square dense matrix of order n.
260: v0 is the initial vector, should have the form v0 = w*w' (for instance 1*1')
261: */
262: static PetscErrorCode LyapIIBuildEigenMat(LME lme,Mat S,Mat *Op,Vec *v0)
263: {
264: PetscErrorCode ierr;
265: PetscInt n,m;
266: PetscBool create=PETSC_FALSE;
267: EPS_EIG_MATSHELL *matctx;
268: #if defined(PETSC_USE_COMPLEX)
269: PetscScalar theta,*aa,*bb,*ss;
270: PetscInt i,j,f,c,off,ld;
271: IS perm;
272: #endif
275: MatGetSize(S,&n,NULL);
276: if (!*Op) create=PETSC_TRUE;
277: else {
278: MatGetSize(*Op,&m,NULL);
279: if (m!=n*n) create=PETSC_TRUE;
280: }
281: if (create) {
282: MatDestroy(Op);
283: VecDestroy(v0);
284: PetscNew(&matctx);
285: #if defined(PETSC_USE_COMPLEX)
286: MatCreateSeqDense(PETSC_COMM_SELF,n*n,n*n,NULL,&matctx->A);
287: MatCreateSeqDense(PETSC_COMM_SELF,n*n,n*n,NULL,&matctx->B);
288: MatCreateVecs(matctx->A,NULL,&matctx->w);
289: #endif
290: MatCreateShell(PETSC_COMM_SELF,n*n,n*n,PETSC_DETERMINE,PETSC_DETERMINE,matctx,Op);
291: MatShellSetOperation(*Op,MATOP_MULT,(void(*)(void))MatMult_EigOperator);
292: MatShellSetOperation(*Op,MATOP_DESTROY,(void(*)(void))MatDestroy_EigOperator);
293: MatCreateVecs(*Op,NULL,v0);
294: } else {
295: MatShellGetContext(*Op,(void**)&matctx);
296: #if defined(PETSC_USE_COMPLEX)
297: MatZeroEntries(matctx->A);
298: #endif
299: }
300: #if defined(PETSC_USE_COMPLEX)
301: MatDenseGetArray(matctx->A,&aa);
302: MatDenseGetArray(matctx->B,&bb);
303: MatDenseGetArray(S,&ss);
304: ld = n*n;
305: for (f=0;f<n;f++) {
306: off = f*n+f*n*ld;
307: for (i=0;i<n;i++) for (j=0;j<n;j++) aa[off+i+j*ld] = ss[i+j*n];
308: for (c=0;c<n;c++) {
309: off = f*n+c*n*ld;
310: theta = ss[f+c*n];
311: for (i=0;i<n;i++) aa[off+i+i*ld] += theta;
312: for (i=0;i<n;i++) for (j=0;j<n;j++) bb[off+i+j*ld] = -2*theta*ss[i+j*n];
313: }
314: }
315: MatDenseRestoreArray(matctx->A,&aa);
316: MatDenseRestoreArray(matctx->B,&bb);
317: MatDenseRestoreArray(S,&ss);
318: ISCreateStride(PETSC_COMM_SELF,n*n,0,1,&perm);
319: MatDestroy(&matctx->F);
320: MatDuplicate(matctx->A,MAT_COPY_VALUES,&matctx->F);
321: MatLUFactor(matctx->F,perm,perm,0);
322: ISDestroy(&perm);
323: #endif
324: matctx->lme = lme;
325: matctx->S = S;
326: matctx->n = n;
327: VecSet(*v0,1.0);
328: return(0);
329: }
331: PetscErrorCode EPSSolve_LyapII(EPS eps)
332: {
333: PetscErrorCode ierr;
334: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
335: PetscInt i,ldds,rk,nloc,mloc,nv,idx,k;
336: Vec v,w,z=eps->work[0],v0=NULL;
337: Mat S,C,Ux[2],Y,Y1,R,U,W,X,Op=NULL;
338: BV V;
339: BVOrthogType type;
340: BVOrthogRefineType refine;
341: PetscScalar eigr[2],eigi[2],*array,er,ei,*uu,*s,*xx,*aa,pM[4],vec[4];
342: PetscReal eta;
343: EPS epsrr;
344: PetscReal norm;
345: EPS_LYAPII_MATSHELL *matctx;
348: DSGetLeadingDimension(ctx->ds,&ldds);
350: /* Operator for the Lyapunov equation */
351: PetscNew(&matctx);
352: STGetOperator(eps->st,&matctx->S);
353: MatGetLocalSize(matctx->S,&mloc,&nloc);
354: MatCreateShell(PetscObjectComm((PetscObject)eps),mloc,nloc,PETSC_DETERMINE,PETSC_DETERMINE,matctx,&S);
355: matctx->Q = eps->V;
356: MatShellSetOperation(S,MATOP_MULT,(void(*)(void))MatMult_EPSLyapIIOperator);
357: MatShellSetOperation(S,MATOP_DESTROY,(void(*)(void))MatDestroy_EPSLyapIIOperator);
358: LMESetCoefficients(ctx->lme,S,NULL,NULL,NULL);
360: /* Right-hand side */
361: BVDuplicateResize(eps->V,ctx->rkl,&V);
362: BVGetOrthogonalization(V,&type,&refine,&eta,NULL);
363: BVSetOrthogonalization(V,type,refine,eta,BV_ORTHOG_BLOCK_TSQR);
364: MatCreateDense(PetscObjectComm((PetscObject)eps),eps->nloc,PETSC_DECIDE,PETSC_DECIDE,1,NULL,&Ux[0]);
365: MatCreateDense(PetscObjectComm((PetscObject)eps),eps->nloc,PETSC_DECIDE,PETSC_DECIDE,2,NULL,&Ux[1]);
366: nv = ctx->rkl;
367: PetscMalloc1(nv,&s);
369: /* Initialize first column */
370: EPSGetStartVector(eps,0,NULL);
371: BVGetColumn(eps->V,0,&v);
372: BVInsertVec(V,0,v);
373: BVRestoreColumn(eps->V,0,&v);
374: BVSetActiveColumns(eps->V,0,0); /* no deflation at the beginning */
375: LyapIIBuildRHS(S,1,Ux[0],V,eps->work);
376: idx = 0;
378: /* EPS for rank reduction */
379: EPSCreate(PETSC_COMM_SELF,&epsrr);
380: EPSSetOptionsPrefix(epsrr,((PetscObject)eps)->prefix);
381: EPSAppendOptionsPrefix(epsrr,"eps_lyapii_");
382: EPSSetDimensions(epsrr,1,PETSC_DEFAULT,PETSC_DEFAULT);
383: EPSSetTolerances(epsrr,PETSC_MACHINE_EPSILON*100,PETSC_DEFAULT);
385: while (eps->reason == EPS_CONVERGED_ITERATING) {
386: eps->its++;
388: /* Matrix for placing the solution of the Lyapunov equation (an alias of V) */
389: BVSetActiveColumns(V,0,nv);
390: BVGetMat(V,&Y1);
391: MatZeroEntries(Y1);
392: MatCreateLRC(NULL,Y1,NULL,NULL,&Y);
393: LMESetSolution(ctx->lme,Y);
395: /* Solve the Lyapunov equation SY + YS' = -2*S*Z*S' */
396: MatCreateLRC(NULL,Ux[idx],NULL,NULL,&C);
397: LMESetRHS(ctx->lme,C);
398: MatDestroy(&C);
399: LMESolve(ctx->lme);
400: BVRestoreMat(V,&Y1);
401: MatDestroy(&Y);
403: /* SVD of the solution: [Q,R]=qr(V); [U,Sigma,~]=svd(R) */
404: DSSetDimensions(ctx->ds,nv,nv,0,0);
405: DSGetMat(ctx->ds,DS_MAT_A,&R);
406: BVOrthogonalize(V,R);
407: DSRestoreMat(ctx->ds,DS_MAT_A,&R);
408: DSSetState(ctx->ds,DS_STATE_RAW);
409: DSSolve(ctx->ds,s,NULL);
411: /* Determine rank */
412: rk = nv;
413: for (i=1;i<nv;i++) if (PetscAbsScalar(s[i]/s[0])<PETSC_SQRT_MACHINE_EPSILON) {rk=i; break;}
414: PetscInfo1(eps,"The computed solution of the Lyapunov equation has rank %D\n",rk);
415: rk = PetscMin(rk,ctx->rkc);
416: DSGetMat(ctx->ds,DS_MAT_U,&U);
417: BVMultInPlace(V,U,0,rk);
418: BVSetActiveColumns(V,0,rk);
419: MatDestroy(&U);
421: /* Rank reduction */
422: DSSetDimensions(ctx->ds,rk,rk,0,0);
423: DSGetMat(ctx->ds,DS_MAT_A,&W);
424: BVMatProject(V,S,V,W);
425: LyapIIBuildEigenMat(ctx->lme,W,&Op,&v0); /* Op=A\B, A=kron(I,S)+kron(S,I), B=-2*kron(S,S) */
426: EPSSetOperators(epsrr,Op,NULL);
427: EPSSetInitialSpace(epsrr,1,&v0);
428: EPSSolve(epsrr);
429: MatDestroy(&W);
430: EPSComputeVectors(epsrr);
431: /* Copy first eigenvector, vec(A)=x */
432: BVGetArray(epsrr->V,&xx);
433: DSGetArray(ctx->ds,DS_MAT_A,&aa);
434: for (i=0;i<rk;i++) {
435: PetscMemcpy(aa+i*ldds,xx+i*rk,rk*sizeof(PetscScalar));
436: }
437: DSRestoreArray(ctx->ds,DS_MAT_A,&aa);
438: BVRestoreArray(epsrr->V,&xx);
439: DSSetState(ctx->ds,DS_STATE_RAW);
440: /* Compute [U,Sigma,~] = svd(A), its rank should be 1 or 2 */
441: DSSolve(ctx->ds,s,NULL);
442: if (PetscAbsScalar(s[1]/s[0])<PETSC_SQRT_MACHINE_EPSILON) rk=1;
443: else rk = 2;
444: PetscInfo1(eps,"The eigenvector has rank %D\n",rk);
445: DSGetMat(ctx->ds,DS_MAT_U,&U);
446: BVMultInPlace(V,U,0,rk);
447: MatDestroy(&U);
449: /* Save V in Ux */
450: idx = (rk==2)?1:0;
451: for (i=0;i<rk;i++) {
452: BVGetColumn(V,i,&v);
453: VecGetArray(v,&uu);
454: MatDenseGetColumn(Ux[idx],i,&array);
455: PetscMemcpy(array,uu,eps->nloc*sizeof(PetscScalar));
456: MatDenseRestoreColumn(Ux[idx],&array);
457: VecRestoreArray(v,&uu);
458: BVRestoreColumn(V,i,&v);
459: }
461: /* Eigenpair approximation */
462: BVGetColumn(V,0,&v);
463: MatMult(S,v,z);
464: VecDot(z,v,pM);
465: BVRestoreColumn(V,0,&v);
466: if (rk>1) {
467: BVGetColumn(V,1,&w);
468: VecDot(z,w,pM+1);
469: MatMult(S,w,z);
470: VecDot(z,w,pM+3);
471: BVGetColumn(V,0,&v);
472: VecDot(z,v,pM+2);
473: BVRestoreColumn(V,0,&v);
474: BVRestoreColumn(V,1,&w);
475: EV2x2(pM,2,eigr,eigi,vec);
476: MatCreateSeqDense(PETSC_COMM_SELF,2,2,vec,&X);
477: BVSetActiveColumns(V,0,rk);
478: BVMultInPlace(V,X,0,rk);
479: MatDestroy(&X);
480: #if !defined(PETSC_USE_COMPLEX)
481: norm = eigr[0]*eigr[0]+eigi[0]*eigi[0];
482: er = eigr[0]/norm; ei = -eigi[0]/norm;
483: #else
484: er =1.0/eigr[0]; ei = 0.0;
485: #endif
486: } else {
487: eigr[0] = pM[0]; eigi[0] = 0.0;
488: er = 1.0/eigr[0]; ei = 0.0;
489: }
490: BVGetColumn(V,0,&v);
491: if (eigi[0]!=0.0) {
492: BVGetColumn(V,1,&w);
493: } else w = NULL;
494: eps->eigr[eps->nconv] = eigr[0]; eps->eigi[eps->nconv] = eigi[0];
495: EPSComputeResidualNorm_Private(eps,PETSC_FALSE,er,ei,v,w,eps->work,&norm);
496: BVRestoreColumn(V,0,&v);
497: if (w) {
498: BVRestoreColumn(V,1,&w);
499: }
500: (*eps->converged)(eps,er,ei,norm,&eps->errest[eps->nconv],eps->convergedctx);
501: k = 0;
502: if (eps->errest[eps->nconv]<eps->tol) {
503: k++;
504: if (rk==2) {
505: #if !defined (PETSC_USE_COMPLEX)
506: eps->eigr[eps->nconv+k] = eigr[0]; eps->eigi[eps->nconv+k] = -eigi[0];
507: #else
508: eps->eigr[eps->nconv+k] = PetscConj(eps->eigr[eps->nconv]);
509: #endif
510: k++;
511: }
512: /* Store converged eigenpairs and vectors for deflation */
513: for (i=0;i<k;i++) {
514: BVGetColumn(V,i,&v);
515: BVInsertVec(eps->V,eps->nconv+i,v);
516: BVRestoreColumn(V,i,&v);
517: }
518: eps->nconv += k;
519: BVSetActiveColumns(eps->V,eps->nconv-rk,eps->nconv);
520: BVOrthogonalize(eps->V,NULL);
521: DSSetDimensions(eps->ds,eps->nconv,0,0,0);
522: DSGetMat(eps->ds,DS_MAT_A,&W);
523: BVMatProject(eps->V,matctx->S,eps->V,W);
524: DSRestoreMat(eps->ds,DS_MAT_A,&W);
525: if (eps->nconv<eps->nev) {
526: idx = 0;
527: BVSetRandomColumn(V,0);
528: BVNormColumn(V,0,NORM_2,&norm);
529: BVScaleColumn(V,0,1.0/norm);
530: LyapIIBuildRHS(S,1,Ux[idx],V,eps->work);
531: }
532: } else {
533: /* Prepare right-hand side */
534: LyapIIBuildRHS(S,rk,Ux[idx],NULL,eps->work);
535: }
536: (*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx);
537: EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,eps->nconv+1);
538: }
539: STRestoreOperator(eps->st,&matctx->S);
540: MatDestroy(&S);
541: MatDestroy(&Ux[0]);
542: MatDestroy(&Ux[1]);
543: MatDestroy(&Op);
544: VecDestroy(&v0);
545: BVDestroy(&V);
546: EPSDestroy(&epsrr);
547: PetscFree(s);
548: return(0);
549: }
551: PetscErrorCode EPSSetFromOptions_LyapII(PetscOptionItems *PetscOptionsObject,EPS eps)
552: {
554: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
555: PetscInt k,array[2]={PETSC_DEFAULT,PETSC_DEFAULT};
556: PetscBool flg;
559: PetscOptionsHead(PetscOptionsObject,"EPS Lyapunov Inverse Iteration Options");
561: k = 2;
562: PetscOptionsIntArray("-eps_lyapii_ranks","Ranks for Lyapunov equation (one or two comma-separated integers)","EPSLyapIISetRanks",array,&k,&flg);
563: if (flg) {
564: EPSLyapIISetRanks(eps,array[0],array[1]);
565: }
567: PetscOptionsTail();
569: if (!ctx->lme) { EPSLyapIIGetLME(eps,&ctx->lme); }
570: LMESetFromOptions(ctx->lme);
571: return(0);
572: }
574: static PetscErrorCode EPSLyapIISetRanks_LyapII(EPS eps,PetscInt rkc,PetscInt rkl)
575: {
576: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
579: if (rkc==PETSC_DEFAULT) rkc = 10;
580: if (rkc<2) SETERRQ1(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The compressed rank %D must be larger than 1",rkc);
581: if (rkl==PETSC_DEFAULT) rkl = 3*rkc;
582: 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);
583: if (rkc != ctx->rkc) {
584: ctx->rkc = rkc;
585: eps->state = EPS_STATE_INITIAL;
586: }
587: if (rkl != ctx->rkl) {
588: ctx->rkl = rkl;
589: eps->state = EPS_STATE_INITIAL;
590: }
591: return(0);
592: }
594: /*@
595: EPSLyapIISetRanks - Set the ranks used in the solution of the Lyapunov equation.
597: Collective on EPS
599: Input Parameters:
600: + eps - the eigenproblem solver context
601: . rkc - the compressed rank
602: - rkl - the Lyapunov rank
604: Options Database Key:
605: . -eps_lyapii_ranks <rkc,rkl> - Sets the rank parameters
607: Notes:
608: Lyapunov inverse iteration needs to solve a large-scale Lyapunov equation
609: at each iteration of the eigensolver. For this, an iterative solver (LME)
610: is used, which requires to prescribe the rank of the solution matrix X. This
611: is the meaning of parameter rkl. Later, this matrix is compressed into
612: another matrix of rank rkc. If not provided, rkl is a small multiple of rkc.
614: Level: intermediate
616: .seealso: EPSLyapIIGetRanks()
617: @*/
618: PetscErrorCode EPSLyapIISetRanks(EPS eps,PetscInt rkc,PetscInt rkl)
619: {
626: PetscTryMethod(eps,"EPSLyapIISetRanks_C",(EPS,PetscInt,PetscInt),(eps,rkc,rkl));
627: return(0);
628: }
630: static PetscErrorCode EPSLyapIIGetRanks_LyapII(EPS eps,PetscInt *rkc,PetscInt *rkl)
631: {
632: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
635: if (rkc) *rkc = ctx->rkc;
636: if (rkl) *rkl = ctx->rkl;
637: return(0);
638: }
640: /*@
641: EPSLyapIIGetRanks - Return the rank values used for the Lyapunov step.
643: Not Collective
645: Input Parameter:
646: . eps - the eigenproblem solver context
648: Output Parameters:
649: + rkc - the compressed rank
650: - rkl - the Lyapunov rank
652: Level: intermediate
654: .seealso: EPSLyapIISetRanks()
655: @*/
656: PetscErrorCode EPSLyapIIGetRanks(EPS eps,PetscInt *rkc,PetscInt *rkl)
657: {
662: PetscUseMethod(eps,"EPSLyapIIGetRanks_C",(EPS,PetscInt*,PetscInt*),(eps,rkc,rkl));
663: return(0);
664: }
666: static PetscErrorCode EPSLyapIISetLME_LyapII(EPS eps,LME lme)
667: {
669: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
672: PetscObjectReference((PetscObject)lme);
673: LMEDestroy(&ctx->lme);
674: ctx->lme = lme;
675: PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->lme);
676: eps->state = EPS_STATE_INITIAL;
677: return(0);
678: }
680: /*@
681: EPSLyapIISetLME - Associate a linear matrix equation solver object (LME) to the
682: eigenvalue solver.
684: Collective on EPS
686: Input Parameters:
687: + eps - the eigenproblem solver context
688: - lme - the linear matrix equation solver object
690: Level: advanced
692: .seealso: EPSLyapIIGetLME()
693: @*/
694: PetscErrorCode EPSLyapIISetLME(EPS eps,LME lme)
695: {
702: PetscTryMethod(eps,"EPSLyapIISetLME_C",(EPS,LME),(eps,lme));
703: return(0);
704: }
706: static PetscErrorCode EPSLyapIIGetLME_LyapII(EPS eps,LME *lme)
707: {
709: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
712: if (!ctx->lme) {
713: LMECreate(PetscObjectComm((PetscObject)eps),&ctx->lme);
714: LMESetOptionsPrefix(ctx->lme,((PetscObject)eps)->prefix);
715: LMEAppendOptionsPrefix(ctx->lme,"eps_lyapii_");
716: PetscObjectIncrementTabLevel((PetscObject)ctx->lme,(PetscObject)eps,1);
717: PetscLogObjectParent((PetscObject)eps,(PetscObject)ctx->lme);
718: }
719: *lme = ctx->lme;
720: return(0);
721: }
723: /*@
724: EPSLyapIIGetLME - Retrieve the linear matrix equation solver object (LME)
725: associated with the eigenvalue solver.
727: Not Collective
729: Input Parameter:
730: . eps - the eigenproblem solver context
732: Output Parameter:
733: . lme - the linear matrix equation solver object
735: Level: advanced
737: .seealso: EPSLyapIISetLME()
738: @*/
739: PetscErrorCode EPSLyapIIGetLME(EPS eps,LME *lme)
740: {
746: PetscUseMethod(eps,"EPSLyapIIGetLME_C",(EPS,LME*),(eps,lme));
747: return(0);
748: }
750: PetscErrorCode EPSView_LyapII(EPS eps,PetscViewer viewer)
751: {
753: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
754: PetscBool isascii;
757: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
758: if (isascii) {
759: PetscViewerASCIIPrintf(viewer," ranks: for Lyapunov solver=%D, after compression=%D\n",ctx->rkl,ctx->rkc);
760: if (!ctx->lme) { EPSLyapIIGetLME(eps,&ctx->lme); }
761: PetscViewerASCIIPushTab(viewer);
762: LMEView(ctx->lme,viewer);
763: PetscViewerASCIIPopTab(viewer);
764: }
765: return(0);
766: }
768: PetscErrorCode EPSReset_LyapII(EPS eps)
769: {
771: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
774: if (!ctx->lme) { LMEReset(ctx->lme); }
775: return(0);
776: }
778: PetscErrorCode EPSDestroy_LyapII(EPS eps)
779: {
781: EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
784: LMEDestroy(&ctx->lme);
785: DSDestroy(&ctx->ds);
786: PetscFree(eps->data);
787: PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetLME_C",NULL);
788: PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetLME_C",NULL);
789: PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetRanks_C",NULL);
790: PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetRanks_C",NULL);
791: return(0);
792: }
794: PetscErrorCode EPSSetDefaultST_LyapII(EPS eps)
795: {
799: if (!((PetscObject)eps->st)->type_name) {
800: STSetType(eps->st,STSINVERT);
801: }
802: return(0);
803: }
805: SLEPC_EXTERN PetscErrorCode EPSCreate_LyapII(EPS eps)
806: {
807: EPS_LYAPII *ctx;
811: PetscNewLog(eps,&ctx);
812: eps->data = (void*)ctx;
814: eps->useds = PETSC_TRUE;
816: eps->ops->solve = EPSSolve_LyapII;
817: eps->ops->setup = EPSSetUp_LyapII;
818: eps->ops->setupsort = EPSSetUpSort_Default;
819: eps->ops->setfromoptions = EPSSetFromOptions_LyapII;
820: eps->ops->reset = EPSReset_LyapII;
821: eps->ops->destroy = EPSDestroy_LyapII;
822: eps->ops->view = EPSView_LyapII;
823: eps->ops->setdefaultst = EPSSetDefaultST_LyapII;
824: eps->ops->backtransform = EPSBackTransform_Default;
825: eps->ops->computevectors = EPSComputeVectors_Schur;
827: PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetLME_C",EPSLyapIISetLME_LyapII);
828: PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetLME_C",EPSLyapIIGetLME_LyapII);
829: PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetRanks_C",EPSLyapIISetRanks_LyapII);
830: PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetRanks_C",EPSLyapIIGetRanks_LyapII);
831: return(0);
832: }