Actual source code: lyapii.c

slepc-3.18.2 2023-01-26
Report Typos and Errors
  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: PetscErrorCode EPSSetUp_LyapII(EPS eps)
 56: {
 57:   PetscRandom    rand;
 58:   EPS_LYAPII     *ctx = (EPS_LYAPII*)eps->data;

 60:   EPSCheckSinvert(eps);
 61:   if (eps->ncv!=PETSC_DEFAULT) {
 63:   } else eps->ncv = eps->nev+1;
 64:   if (eps->mpd!=PETSC_DEFAULT) PetscInfo(eps,"Warning: parameter mpd ignored\n");
 65:   if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(1000*eps->nev,100*eps->n);
 66:   if (!eps->which) eps->which=EPS_LARGEST_REAL;
 68:   EPSCheckUnsupported(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION | EPS_FEATURE_TWOSIDED);

 70:   if (!ctx->rkc) ctx->rkc = 10;
 71:   if (!ctx->rkl) ctx->rkl = 3*ctx->rkc;
 72:   if (!ctx->lme) EPSLyapIIGetLME(eps,&ctx->lme);
 73:   LMESetProblemType(ctx->lme,LME_LYAPUNOV);
 74:   LMESetErrorIfNotConverged(ctx->lme,PETSC_TRUE);

 76:   if (!ctx->ds) {
 77:     DSCreate(PetscObjectComm((PetscObject)eps),&ctx->ds);
 78:     DSSetType(ctx->ds,DSSVD);
 79:   }
 80:   DSAllocate(ctx->ds,ctx->rkl);

 82:   DSSetType(eps->ds,DSNHEP);
 83:   DSAllocate(eps->ds,eps->ncv);

 85:   EPSAllocateSolution(eps,0);
 86:   BVGetRandomContext(eps->V,&rand);  /* make sure the random context is available when duplicating */
 87:   EPSSetWorkVecs(eps,3);
 88:   return 0;
 89: }

 91: static PetscErrorCode MatMult_EPSLyapIIOperator(Mat M,Vec x,Vec r)
 92: {
 93:   EPS_LYAPII_MATSHELL *matctx;

 95:   MatShellGetContext(M,&matctx);
 96:   MatMult(matctx->S,x,r);
 97:   BVOrthogonalizeVec(matctx->Q,r,NULL,NULL,NULL);
 98:   return 0;
 99: }

101: static PetscErrorCode MatDestroy_EPSLyapIIOperator(Mat M)
102: {
103:   EPS_LYAPII_MATSHELL *matctx;

105:   MatShellGetContext(M,&matctx);
106:   MatDestroy(&matctx->S);
107:   PetscFree(matctx);
108:   return 0;
109: }

111: static PetscErrorCode MatMult_EigOperator(Mat M,Vec x,Vec y)
112: {
113:   EPS_EIG_MATSHELL  *matctx;
114: #if !defined(PETSC_USE_COMPLEX)
115:   PetscInt          n,lds;
116:   PetscScalar       *Y,*C,zero=0.0,done=1.0,dtwo=2.0;
117:   const PetscScalar *S,*X;
118:   PetscBLASInt      n_,lds_;
119: #endif

121:   MatShellGetContext(M,&matctx);

123: #if defined(PETSC_USE_COMPLEX)
124:   MatMult(matctx->B,x,matctx->w);
125:   MatSolve(matctx->F,matctx->w,y);
126: #else
127:   VecGetArrayRead(x,&X);
128:   VecGetArray(y,&Y);
129:   MatDenseGetArrayRead(matctx->S,&S);
130:   MatDenseGetLDA(matctx->S,&lds);

132:   n = matctx->n;
133:   PetscCalloc1(n*n,&C);
134:   PetscBLASIntCast(n,&n_);
135:   PetscBLASIntCast(lds,&lds_);

137:   /* C = 2*S*X*S.' */
138:   PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&dtwo,S,&lds_,X,&n_,&zero,Y,&n_));
139:   PetscCallBLAS("BLASgemm",BLASgemm_("N","T",&n_,&n_,&n_,&done,Y,&n_,S,&lds_,&zero,C,&n_));

141:   /* Solve S*Y + Y*S' = -C */
142:   LMEDenseLyapunov(matctx->lme,n,(PetscScalar*)S,lds,C,n,Y,n);

144:   PetscFree(C);
145:   VecRestoreArrayRead(x,&X);
146:   VecRestoreArray(y,&Y);
147:   MatDenseRestoreArrayRead(matctx->S,&S);
148: #endif
149:   return 0;
150: }

152: static PetscErrorCode MatDestroy_EigOperator(Mat M)
153: {
154:   EPS_EIG_MATSHELL *matctx;

156:   MatShellGetContext(M,&matctx);
157: #if defined(PETSC_USE_COMPLEX)
158:   MatDestroy(&matctx->A);
159:   MatDestroy(&matctx->B);
160:   MatDestroy(&matctx->F);
161:   VecDestroy(&matctx->w);
162: #else
163:   MatDestroy(&matctx->S);
164: #endif
165:   PetscFree(matctx);
166:   return 0;
167: }

169: /*
170:    EV2x2: solve the eigenproblem for a 2x2 matrix M
171:  */
172: static PetscErrorCode EV2x2(PetscScalar *M,PetscInt ld,PetscScalar *wr,PetscScalar *wi,PetscScalar *vec)
173: {
174:   PetscBLASInt   lwork=10,ld_;
175:   PetscScalar    work[10];
176:   PetscBLASInt   two=2,info;
177: #if defined(PETSC_USE_COMPLEX)
178:   PetscReal      rwork[6];
179: #endif

181:   PetscBLASIntCast(ld,&ld_);
182:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
183: #if !defined(PETSC_USE_COMPLEX)
184:   PetscCallBLAS("LAPACKgeev",LAPACKgeev_("N","V",&two,M,&ld_,wr,wi,NULL,&ld_,vec,&ld_,work,&lwork,&info));
185: #else
186:   PetscCallBLAS("LAPACKgeev",LAPACKgeev_("N","V",&two,M,&ld_,wr,NULL,&ld_,vec,&ld_,work,&lwork,rwork,&info));
187: #endif
188:   SlepcCheckLapackInfo("geev",info);
189:   PetscFPTrapPop();
190:   return 0;
191: }

193: /*
194:    LyapIIBuildRHS: prepare the right-hand side of the Lyapunov equation SY + YS' = -2*S*Z*S'
195:    in factored form:
196:       if (V)  U=sqrt(2)*S*V    (uses 1 work vector)
197:       else    U=sqrt(2)*S*U    (uses 2 work vectors)
198:    where U,V are assumed to have rk columns.
199:  */
200: static PetscErrorCode LyapIIBuildRHS(Mat S,PetscInt rk,Mat U,BV V,Vec *work)
201: {
202:   PetscScalar    *array,*uu;
203:   PetscInt       i,nloc;
204:   Vec            v,u=work[0];

206:   MatGetLocalSize(U,&nloc,NULL);
207:   for (i=0;i<rk;i++) {
208:     MatDenseGetColumn(U,i,&array);
209:     if (V) BVGetColumn(V,i,&v);
210:     else {
211:       v = work[1];
212:       VecPlaceArray(v,array);
213:     }
214:     MatMult(S,v,u);
215:     if (V) BVRestoreColumn(V,i,&v);
216:     else VecResetArray(v);
217:     VecScale(u,PETSC_SQRT2);
218:     VecGetArray(u,&uu);
219:     PetscArraycpy(array,uu,nloc);
220:     VecRestoreArray(u,&uu);
221:     MatDenseRestoreColumn(U,&array);
222:   }
223:   return 0;
224: }

226: /*
227:    LyapIIBuildEigenMat: create shell matrix Op=A\B with A = kron(I,S)+kron(S,I), B = -2*kron(S,S)
228:    where S is a sequential square dense matrix of order n.
229:    v0 is the initial vector, should have the form v0 = w*w' (for instance 1*1')
230:  */
231: static PetscErrorCode LyapIIBuildEigenMat(LME lme,Mat S,Mat *Op,Vec *v0)
232: {
233:   PetscInt          n,m;
234:   PetscBool         create=PETSC_FALSE;
235:   EPS_EIG_MATSHELL  *matctx;
236: #if defined(PETSC_USE_COMPLEX)
237:   PetscScalar       theta,*aa,*bb;
238:   const PetscScalar *ss;
239:   PetscInt          i,j,f,c,off,ld,lds;
240:   IS                perm;
241: #endif

243:   MatGetSize(S,&n,NULL);
244:   if (!*Op) create=PETSC_TRUE;
245:   else {
246:     MatGetSize(*Op,&m,NULL);
247:     if (m!=n*n) create=PETSC_TRUE;
248:   }
249:   if (create) {
250:     MatDestroy(Op);
251:     VecDestroy(v0);
252:     PetscNew(&matctx);
253: #if defined(PETSC_USE_COMPLEX)
254:     MatCreateSeqDense(PETSC_COMM_SELF,n*n,n*n,NULL,&matctx->A);
255:     MatCreateSeqDense(PETSC_COMM_SELF,n*n,n*n,NULL,&matctx->B);
256:     MatCreateVecs(matctx->A,NULL,&matctx->w);
257: #else
258:     MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&matctx->S);
259: #endif
260:     MatCreateShell(PETSC_COMM_SELF,n*n,n*n,PETSC_DETERMINE,PETSC_DETERMINE,matctx,Op);
261:     MatShellSetOperation(*Op,MATOP_MULT,(void(*)(void))MatMult_EigOperator);
262:     MatShellSetOperation(*Op,MATOP_DESTROY,(void(*)(void))MatDestroy_EigOperator);
263:     MatCreateVecs(*Op,NULL,v0);
264:   } else {
265:     MatShellGetContext(*Op,&matctx);
266: #if defined(PETSC_USE_COMPLEX)
267:     MatZeroEntries(matctx->A);
268: #endif
269:   }
270: #if defined(PETSC_USE_COMPLEX)
271:   MatDenseGetArray(matctx->A,&aa);
272:   MatDenseGetArray(matctx->B,&bb);
273:   MatDenseGetArrayRead(S,&ss);
274:   MatDenseGetLDA(S,&lds);
275:   ld = n*n;
276:   for (f=0;f<n;f++) {
277:     off = f*n+f*n*ld;
278:     for (i=0;i<n;i++) for (j=0;j<n;j++) aa[off+i+j*ld] = ss[i+j*lds];
279:     for (c=0;c<n;c++) {
280:       off = f*n+c*n*ld;
281:       theta = ss[f+c*lds];
282:       for (i=0;i<n;i++) aa[off+i+i*ld] += theta;
283:       for (i=0;i<n;i++) for (j=0;j<n;j++) bb[off+i+j*ld] = -2*theta*ss[i+j*lds];
284:     }
285:   }
286:   MatDenseRestoreArray(matctx->A,&aa);
287:   MatDenseRestoreArray(matctx->B,&bb);
288:   MatDenseRestoreArrayRead(S,&ss);
289:   ISCreateStride(PETSC_COMM_SELF,n*n,0,1,&perm);
290:   MatDestroy(&matctx->F);
291:   MatDuplicate(matctx->A,MAT_COPY_VALUES,&matctx->F);
292:   MatLUFactor(matctx->F,perm,perm,NULL);
293:   ISDestroy(&perm);
294: #else
295:   MatCopy(S,matctx->S,SAME_NONZERO_PATTERN);
296: #endif
297:   matctx->lme = lme;
298:   matctx->n = n;
299:   VecSet(*v0,1.0);
300:   return 0;
301: }

303: PetscErrorCode EPSSolve_LyapII(EPS eps)
304: {
305:   EPS_LYAPII          *ctx = (EPS_LYAPII*)eps->data;
306:   PetscInt            i,ldds,rk,nloc,mloc,nv,idx,k;
307:   Vec                 v,w,z=eps->work[0],v0=NULL;
308:   Mat                 S,C,Ux[2],Y,Y1,R,U,W,X,Op=NULL;
309:   BV                  V;
310:   BVOrthogType        type;
311:   BVOrthogRefineType  refine;
312:   PetscScalar         eigr[2],eigi[2],*array,er,ei,*uu,*s,*xx,*aa,pM[4],vec[4];
313:   PetscReal           eta;
314:   EPS                 epsrr;
315:   PetscReal           norm;
316:   EPS_LYAPII_MATSHELL *matctx;

318:   DSGetLeadingDimension(ctx->ds,&ldds);

320:   /* Operator for the Lyapunov equation */
321:   PetscNew(&matctx);
322:   STGetOperator(eps->st,&matctx->S);
323:   MatGetLocalSize(matctx->S,&mloc,&nloc);
324:   MatCreateShell(PetscObjectComm((PetscObject)eps),mloc,nloc,PETSC_DETERMINE,PETSC_DETERMINE,matctx,&S);
325:   matctx->Q = eps->V;
326:   MatShellSetOperation(S,MATOP_MULT,(void(*)(void))MatMult_EPSLyapIIOperator);
327:   MatShellSetOperation(S,MATOP_DESTROY,(void(*)(void))MatDestroy_EPSLyapIIOperator);
328:   LMESetCoefficients(ctx->lme,S,NULL,NULL,NULL);

330:   /* Right-hand side */
331:   BVDuplicateResize(eps->V,ctx->rkl,&V);
332:   BVGetOrthogonalization(V,&type,&refine,&eta,NULL);
333:   BVSetOrthogonalization(V,type,refine,eta,BV_ORTHOG_BLOCK_TSQR);
334:   MatCreateDense(PetscObjectComm((PetscObject)eps),eps->nloc,PETSC_DECIDE,PETSC_DECIDE,1,NULL,&Ux[0]);
335:   MatCreateDense(PetscObjectComm((PetscObject)eps),eps->nloc,PETSC_DECIDE,PETSC_DECIDE,2,NULL,&Ux[1]);
336:   nv = ctx->rkl;
337:   PetscMalloc1(nv,&s);

339:   /* Initialize first column */
340:   EPSGetStartVector(eps,0,NULL);
341:   BVGetColumn(eps->V,0,&v);
342:   BVInsertVec(V,0,v);
343:   BVRestoreColumn(eps->V,0,&v);
344:   BVSetActiveColumns(eps->V,0,0);  /* no deflation at the beginning */
345:   LyapIIBuildRHS(S,1,Ux[0],V,eps->work);
346:   idx = 0;

348:   /* EPS for rank reduction */
349:   EPSCreate(PETSC_COMM_SELF,&epsrr);
350:   EPSSetOptionsPrefix(epsrr,((PetscObject)eps)->prefix);
351:   EPSAppendOptionsPrefix(epsrr,"eps_lyapii_");
352:   EPSSetDimensions(epsrr,1,PETSC_DEFAULT,PETSC_DEFAULT);
353:   EPSSetTolerances(epsrr,PETSC_MACHINE_EPSILON*100,PETSC_DEFAULT);

355:   while (eps->reason == EPS_CONVERGED_ITERATING) {
356:     eps->its++;

358:     /* Matrix for placing the solution of the Lyapunov equation (an alias of V) */
359:     BVSetActiveColumns(V,0,nv);
360:     BVGetMat(V,&Y1);
361:     MatZeroEntries(Y1);
362:     MatCreateLRC(NULL,Y1,NULL,NULL,&Y);
363:     LMESetSolution(ctx->lme,Y);

365:     /* Solve the Lyapunov equation SY + YS' = -2*S*Z*S' */
366:     MatCreateLRC(NULL,Ux[idx],NULL,NULL,&C);
367:     LMESetRHS(ctx->lme,C);
368:     MatDestroy(&C);
369:     LMESolve(ctx->lme);
370:     BVRestoreMat(V,&Y1);
371:     MatDestroy(&Y);

373:     /* SVD of the solution: [Q,R]=qr(V); [U,Sigma,~]=svd(R) */
374:     DSSetDimensions(ctx->ds,nv,0,0);
375:     DSSVDSetDimensions(ctx->ds,nv);
376:     DSGetMat(ctx->ds,DS_MAT_A,&R);
377:     BVOrthogonalize(V,R);
378:     DSRestoreMat(ctx->ds,DS_MAT_A,&R);
379:     DSSetState(ctx->ds,DS_STATE_RAW);
380:     DSSolve(ctx->ds,s,NULL);

382:     /* Determine rank */
383:     rk = nv;
384:     for (i=1;i<nv;i++) if (PetscAbsScalar(s[i]/s[0])<PETSC_SQRT_MACHINE_EPSILON) {rk=i; break;}
385:     PetscInfo(eps,"The computed solution of the Lyapunov equation has rank %" PetscInt_FMT "\n",rk);
386:     rk = PetscMin(rk,ctx->rkc);
387:     DSGetMat(ctx->ds,DS_MAT_U,&U);
388:     BVMultInPlace(V,U,0,rk);
389:     DSRestoreMat(ctx->ds,DS_MAT_U,&U);
390:     BVSetActiveColumns(V,0,rk);

392:     /* Rank reduction */
393:     DSSetDimensions(ctx->ds,rk,0,0);
394:     DSSVDSetDimensions(ctx->ds,rk);
395:     DSGetMat(ctx->ds,DS_MAT_A,&W);
396:     BVMatProject(V,S,V,W);
397:     LyapIIBuildEigenMat(ctx->lme,W,&Op,&v0); /* Op=A\B, A=kron(I,S)+kron(S,I), B=-2*kron(S,S) */
398:     DSRestoreMat(ctx->ds,DS_MAT_A,&W);
399:     EPSSetOperators(epsrr,Op,NULL);
400:     EPSSetInitialSpace(epsrr,1,&v0);
401:     EPSSolve(epsrr);
402:     EPSComputeVectors(epsrr);
403:     /* Copy first eigenvector, vec(A)=x */
404:     BVGetArray(epsrr->V,&xx);
405:     DSGetArray(ctx->ds,DS_MAT_A,&aa);
406:     for (i=0;i<rk;i++) PetscArraycpy(aa+i*ldds,xx+i*rk,rk);
407:     DSRestoreArray(ctx->ds,DS_MAT_A,&aa);
408:     BVRestoreArray(epsrr->V,&xx);
409:     DSSetState(ctx->ds,DS_STATE_RAW);
410:     /* Compute [U,Sigma,~] = svd(A), its rank should be 1 or 2 */
411:     DSSolve(ctx->ds,s,NULL);
412:     if (PetscAbsScalar(s[1]/s[0])<PETSC_SQRT_MACHINE_EPSILON) rk=1;
413:     else rk = 2;
414:     PetscInfo(eps,"The eigenvector has rank %" PetscInt_FMT "\n",rk);
415:     DSGetMat(ctx->ds,DS_MAT_U,&U);
416:     BVMultInPlace(V,U,0,rk);
417:     DSRestoreMat(ctx->ds,DS_MAT_U,&U);

419:     /* Save V in Ux */
420:     idx = (rk==2)?1:0;
421:     for (i=0;i<rk;i++) {
422:       BVGetColumn(V,i,&v);
423:       VecGetArray(v,&uu);
424:       MatDenseGetColumn(Ux[idx],i,&array);
425:       PetscArraycpy(array,uu,eps->nloc);
426:       MatDenseRestoreColumn(Ux[idx],&array);
427:       VecRestoreArray(v,&uu);
428:       BVRestoreColumn(V,i,&v);
429:     }

431:     /* Eigenpair approximation */
432:     BVGetColumn(V,0,&v);
433:     MatMult(S,v,z);
434:     VecDot(z,v,pM);
435:     BVRestoreColumn(V,0,&v);
436:     if (rk>1) {
437:       BVGetColumn(V,1,&w);
438:       VecDot(z,w,pM+1);
439:       MatMult(S,w,z);
440:       VecDot(z,w,pM+3);
441:       BVGetColumn(V,0,&v);
442:       VecDot(z,v,pM+2);
443:       BVRestoreColumn(V,0,&v);
444:       BVRestoreColumn(V,1,&w);
445:       EV2x2(pM,2,eigr,eigi,vec);
446:       MatCreateSeqDense(PETSC_COMM_SELF,2,2,vec,&X);
447:       BVSetActiveColumns(V,0,rk);
448:       BVMultInPlace(V,X,0,rk);
449:       MatDestroy(&X);
450: #if !defined(PETSC_USE_COMPLEX)
451:       norm = eigr[0]*eigr[0]+eigi[0]*eigi[0];
452:       er = eigr[0]/norm; ei = -eigi[0]/norm;
453: #else
454:       er =1.0/eigr[0]; ei = 0.0;
455: #endif
456:     } else {
457:       eigr[0] = pM[0]; eigi[0] = 0.0;
458:       er = 1.0/eigr[0]; ei = 0.0;
459:     }
460:     BVGetColumn(V,0,&v);
461:     if (eigi[0]!=0.0) BVGetColumn(V,1,&w);
462:     else w = NULL;
463:     eps->eigr[eps->nconv] = eigr[0]; eps->eigi[eps->nconv] = eigi[0];
464:     EPSComputeResidualNorm_Private(eps,PETSC_FALSE,er,ei,v,w,eps->work,&norm);
465:     BVRestoreColumn(V,0,&v);
466:     if (w) BVRestoreColumn(V,1,&w);
467:     (*eps->converged)(eps,er,ei,norm,&eps->errest[eps->nconv],eps->convergedctx);
468:     k = 0;
469:     if (eps->errest[eps->nconv]<eps->tol) {
470:       k++;
471:       if (rk==2) {
472: #if !defined (PETSC_USE_COMPLEX)
473:         eps->eigr[eps->nconv+k] = eigr[0]; eps->eigi[eps->nconv+k] = -eigi[0];
474: #else
475:         eps->eigr[eps->nconv+k] = PetscConj(eps->eigr[eps->nconv]);
476: #endif
477:         k++;
478:       }
479:       /* Store converged eigenpairs and vectors for deflation */
480:       for (i=0;i<k;i++) {
481:         BVGetColumn(V,i,&v);
482:         BVInsertVec(eps->V,eps->nconv+i,v);
483:         BVRestoreColumn(V,i,&v);
484:       }
485:       eps->nconv += k;
486:       BVSetActiveColumns(eps->V,eps->nconv-rk,eps->nconv);
487:       BVOrthogonalize(eps->V,NULL);
488:       DSSetDimensions(eps->ds,eps->nconv,0,0);
489:       DSGetMat(eps->ds,DS_MAT_A,&W);
490:       BVMatProject(eps->V,matctx->S,eps->V,W);
491:       DSRestoreMat(eps->ds,DS_MAT_A,&W);
492:       if (eps->nconv<eps->nev) {
493:         idx = 0;
494:         BVSetRandomColumn(V,0);
495:         BVNormColumn(V,0,NORM_2,&norm);
496:         BVScaleColumn(V,0,1.0/norm);
497:         LyapIIBuildRHS(S,1,Ux[idx],V,eps->work);
498:       }
499:     } else {
500:       /* Prepare right-hand side */
501:       LyapIIBuildRHS(S,rk,Ux[idx],NULL,eps->work);
502:     }
503:     (*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx);
504:     EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,eps->nconv+1);
505:   }
506:   STRestoreOperator(eps->st,&matctx->S);
507:   MatDestroy(&S);
508:   MatDestroy(&Ux[0]);
509:   MatDestroy(&Ux[1]);
510:   MatDestroy(&Op);
511:   VecDestroy(&v0);
512:   BVDestroy(&V);
513:   EPSDestroy(&epsrr);
514:   PetscFree(s);
515:   return 0;
516: }

518: PetscErrorCode EPSSetFromOptions_LyapII(EPS eps,PetscOptionItems *PetscOptionsObject)
519: {
520:   EPS_LYAPII     *ctx = (EPS_LYAPII*)eps->data;
521:   PetscInt       k,array[2]={PETSC_DEFAULT,PETSC_DEFAULT};
522:   PetscBool      flg;

524:   PetscOptionsHeadBegin(PetscOptionsObject,"EPS Lyapunov Inverse Iteration Options");

526:     k = 2;
527:     PetscOptionsIntArray("-eps_lyapii_ranks","Ranks for Lyapunov equation (one or two comma-separated integers)","EPSLyapIISetRanks",array,&k,&flg);
528:     if (flg) EPSLyapIISetRanks(eps,array[0],array[1]);

530:   PetscOptionsHeadEnd();

532:   if (!ctx->lme) EPSLyapIIGetLME(eps,&ctx->lme);
533:   LMESetFromOptions(ctx->lme);
534:   return 0;
535: }

537: static PetscErrorCode EPSLyapIISetRanks_LyapII(EPS eps,PetscInt rkc,PetscInt rkl)
538: {
539:   EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;

541:   if (rkc==PETSC_DEFAULT) rkc = 10;
543:   if (rkl==PETSC_DEFAULT) rkl = 3*rkc;
545:   if (rkc != ctx->rkc) {
546:     ctx->rkc   = rkc;
547:     eps->state = EPS_STATE_INITIAL;
548:   }
549:   if (rkl != ctx->rkl) {
550:     ctx->rkl   = rkl;
551:     eps->state = EPS_STATE_INITIAL;
552:   }
553:   return 0;
554: }

556: /*@
557:    EPSLyapIISetRanks - Set the ranks used in the solution of the Lyapunov equation.

559:    Logically Collective on EPS

561:    Input Parameters:
562: +  eps - the eigenproblem solver context
563: .  rkc - the compressed rank
564: -  rkl - the Lyapunov rank

566:    Options Database Key:
567: .  -eps_lyapii_ranks <rkc,rkl> - Sets the rank parameters

569:    Notes:
570:    Lyapunov inverse iteration needs to solve a large-scale Lyapunov equation
571:    at each iteration of the eigensolver. For this, an iterative solver (LME)
572:    is used, which requires to prescribe the rank of the solution matrix X. This
573:    is the meaning of parameter rkl. Later, this matrix is compressed into
574:    another matrix of rank rkc. If not provided, rkl is a small multiple of rkc.

576:    Level: intermediate

578: .seealso: EPSLyapIIGetRanks()
579: @*/
580: PetscErrorCode EPSLyapIISetRanks(EPS eps,PetscInt rkc,PetscInt rkl)
581: {
585:   PetscTryMethod(eps,"EPSLyapIISetRanks_C",(EPS,PetscInt,PetscInt),(eps,rkc,rkl));
586:   return 0;
587: }

589: static PetscErrorCode EPSLyapIIGetRanks_LyapII(EPS eps,PetscInt *rkc,PetscInt *rkl)
590: {
591:   EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;

593:   if (rkc) *rkc = ctx->rkc;
594:   if (rkl) *rkl = ctx->rkl;
595:   return 0;
596: }

598: /*@
599:    EPSLyapIIGetRanks - Return the rank values used for the Lyapunov step.

601:    Not Collective

603:    Input Parameter:
604: .  eps - the eigenproblem solver context

606:    Output Parameters:
607: +  rkc - the compressed rank
608: -  rkl - the Lyapunov rank

610:    Level: intermediate

612: .seealso: EPSLyapIISetRanks()
613: @*/
614: PetscErrorCode EPSLyapIIGetRanks(EPS eps,PetscInt *rkc,PetscInt *rkl)
615: {
617:   PetscUseMethod(eps,"EPSLyapIIGetRanks_C",(EPS,PetscInt*,PetscInt*),(eps,rkc,rkl));
618:   return 0;
619: }

621: static PetscErrorCode EPSLyapIISetLME_LyapII(EPS eps,LME lme)
622: {
623:   EPS_LYAPII     *ctx = (EPS_LYAPII*)eps->data;

625:   PetscObjectReference((PetscObject)lme);
626:   LMEDestroy(&ctx->lme);
627:   ctx->lme = lme;
628:   eps->state = EPS_STATE_INITIAL;
629:   return 0;
630: }

632: /*@
633:    EPSLyapIISetLME - Associate a linear matrix equation solver object (LME) to the
634:    eigenvalue solver.

636:    Collective on EPS

638:    Input Parameters:
639: +  eps - the eigenproblem solver context
640: -  lme - the linear matrix equation solver object

642:    Level: advanced

644: .seealso: EPSLyapIIGetLME()
645: @*/
646: PetscErrorCode EPSLyapIISetLME(EPS eps,LME lme)
647: {
651:   PetscTryMethod(eps,"EPSLyapIISetLME_C",(EPS,LME),(eps,lme));
652:   return 0;
653: }

655: static PetscErrorCode EPSLyapIIGetLME_LyapII(EPS eps,LME *lme)
656: {
657:   EPS_LYAPII     *ctx = (EPS_LYAPII*)eps->data;

659:   if (!ctx->lme) {
660:     LMECreate(PetscObjectComm((PetscObject)eps),&ctx->lme);
661:     LMESetOptionsPrefix(ctx->lme,((PetscObject)eps)->prefix);
662:     LMEAppendOptionsPrefix(ctx->lme,"eps_lyapii_");
663:     PetscObjectIncrementTabLevel((PetscObject)ctx->lme,(PetscObject)eps,1);
664:   }
665:   *lme = ctx->lme;
666:   return 0;
667: }

669: /*@
670:    EPSLyapIIGetLME - Retrieve the linear matrix equation solver object (LME)
671:    associated with the eigenvalue solver.

673:    Not Collective

675:    Input Parameter:
676: .  eps - the eigenproblem solver context

678:    Output Parameter:
679: .  lme - the linear matrix equation solver object

681:    Level: advanced

683: .seealso: EPSLyapIISetLME()
684: @*/
685: PetscErrorCode EPSLyapIIGetLME(EPS eps,LME *lme)
686: {
689:   PetscUseMethod(eps,"EPSLyapIIGetLME_C",(EPS,LME*),(eps,lme));
690:   return 0;
691: }

693: PetscErrorCode EPSView_LyapII(EPS eps,PetscViewer viewer)
694: {
695:   EPS_LYAPII     *ctx = (EPS_LYAPII*)eps->data;
696:   PetscBool      isascii;

698:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
699:   if (isascii) {
700:     PetscViewerASCIIPrintf(viewer,"  ranks: for Lyapunov solver=%" PetscInt_FMT ", after compression=%" PetscInt_FMT "\n",ctx->rkl,ctx->rkc);
701:     if (!ctx->lme) EPSLyapIIGetLME(eps,&ctx->lme);
702:     PetscViewerASCIIPushTab(viewer);
703:     LMEView(ctx->lme,viewer);
704:     PetscViewerASCIIPopTab(viewer);
705:   }
706:   return 0;
707: }

709: PetscErrorCode EPSReset_LyapII(EPS eps)
710: {
711:   EPS_LYAPII     *ctx = (EPS_LYAPII*)eps->data;

713:   if (!ctx->lme) LMEReset(ctx->lme);
714:   return 0;
715: }

717: PetscErrorCode EPSDestroy_LyapII(EPS eps)
718: {
719:   EPS_LYAPII     *ctx = (EPS_LYAPII*)eps->data;

721:   LMEDestroy(&ctx->lme);
722:   DSDestroy(&ctx->ds);
723:   PetscFree(eps->data);
724:   PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetLME_C",NULL);
725:   PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetLME_C",NULL);
726:   PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetRanks_C",NULL);
727:   PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetRanks_C",NULL);
728:   return 0;
729: }

731: PetscErrorCode EPSSetDefaultST_LyapII(EPS eps)
732: {
733:   if (!((PetscObject)eps->st)->type_name) STSetType(eps->st,STSINVERT);
734:   return 0;
735: }

737: SLEPC_EXTERN PetscErrorCode EPSCreate_LyapII(EPS eps)
738: {
739:   EPS_LYAPII     *ctx;

741:   PetscNew(&ctx);
742:   eps->data = (void*)ctx;

744:   eps->useds = PETSC_TRUE;

746:   eps->ops->solve          = EPSSolve_LyapII;
747:   eps->ops->setup          = EPSSetUp_LyapII;
748:   eps->ops->setupsort      = EPSSetUpSort_Default;
749:   eps->ops->setfromoptions = EPSSetFromOptions_LyapII;
750:   eps->ops->reset          = EPSReset_LyapII;
751:   eps->ops->destroy        = EPSDestroy_LyapII;
752:   eps->ops->view           = EPSView_LyapII;
753:   eps->ops->setdefaultst   = EPSSetDefaultST_LyapII;
754:   eps->ops->backtransform  = EPSBackTransform_Default;
755:   eps->ops->computevectors = EPSComputeVectors_Schur;

757:   PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetLME_C",EPSLyapIISetLME_LyapII);
758:   PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetLME_C",EPSLyapIIGetLME_LyapII);
759:   PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetRanks_C",EPSLyapIISetRanks_LyapII);
760:   PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetRanks_C",EPSLyapIIGetRanks_LyapII);
761:   return 0;
762: }