Actual source code: lyapii.c

slepc-main 2024-11-09
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: 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: }