Actual source code: lyapii.c

slepc-3.15.2 2021-09-20
Report Typos and Errors
  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: }