Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
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"
12 :
13 : Method: Lyapunov inverse iteration
14 :
15 : Algorithm:
16 :
17 : Lyapunov inverse iteration using LME solvers
18 :
19 : References:
20 :
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.
24 :
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 : */
29 :
30 : #include <slepc/private/epsimpl.h> /*I "slepceps.h" I*/
31 : #include <slepcblaslapack.h>
32 :
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;
39 :
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;
44 :
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;
54 :
55 3 : static PetscErrorCode EPSSetUp_LyapII(EPS eps)
56 : {
57 3 : PetscRandom rand;
58 3 : EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
59 :
60 3 : PetscFunctionBegin;
61 3 : EPSCheckSinvert(eps);
62 3 : EPSCheckNotStructured(eps);
63 3 : if (eps->nev==0) eps->nev = 1;
64 3 : if (eps->ncv!=PETSC_DETERMINE) {
65 0 : PetscCheck(eps->ncv>=eps->nev+1,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nev+1");
66 3 : } else eps->ncv = eps->nev+1;
67 3 : if (eps->mpd!=PETSC_DETERMINE) PetscCall(PetscInfo(eps,"Warning: parameter mpd ignored\n"));
68 3 : if (eps->max_it==PETSC_DETERMINE) eps->max_it = PetscMax(1000*eps->nev,100*eps->n);
69 3 : if (!eps->which) eps->which=EPS_LARGEST_REAL;
70 3 : PetscCheck(eps->which==EPS_LARGEST_REAL,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver supports only largest real eigenvalues");
71 3 : EPSCheckUnsupported(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION | EPS_FEATURE_THRESHOLD | EPS_FEATURE_TWOSIDED);
72 :
73 3 : if (!ctx->rkc) ctx->rkc = 10;
74 3 : if (!ctx->rkl) ctx->rkl = 3*ctx->rkc;
75 3 : if (!ctx->lme) PetscCall(EPSLyapIIGetLME(eps,&ctx->lme));
76 3 : PetscCall(LMESetProblemType(ctx->lme,LME_LYAPUNOV));
77 3 : PetscCall(LMESetErrorIfNotConverged(ctx->lme,PETSC_TRUE));
78 :
79 3 : if (!ctx->ds) {
80 3 : PetscCall(DSCreate(PetscObjectComm((PetscObject)eps),&ctx->ds));
81 3 : PetscCall(DSSetType(ctx->ds,DSSVD));
82 : }
83 3 : PetscCall(DSAllocate(ctx->ds,ctx->rkl));
84 :
85 3 : PetscCall(DSSetType(eps->ds,DSNHEP));
86 3 : PetscCall(DSAllocate(eps->ds,eps->ncv));
87 :
88 3 : PetscCall(EPSAllocateSolution(eps,0));
89 3 : PetscCall(BVGetRandomContext(eps->V,&rand)); /* make sure the random context is available when duplicating */
90 3 : PetscCall(EPSSetWorkVecs(eps,3));
91 3 : PetscFunctionReturn(PETSC_SUCCESS);
92 : }
93 :
94 3174 : static PetscErrorCode MatMult_EPSLyapIIOperator(Mat M,Vec x,Vec r)
95 : {
96 3174 : EPS_LYAPII_MATSHELL *matctx;
97 :
98 3174 : PetscFunctionBegin;
99 3174 : PetscCall(MatShellGetContext(M,&matctx));
100 3174 : PetscCall(MatMult(matctx->S,x,r));
101 3174 : PetscCall(BVOrthogonalizeVec(matctx->Q,r,NULL,NULL,NULL));
102 3174 : PetscFunctionReturn(PETSC_SUCCESS);
103 : }
104 :
105 3 : static PetscErrorCode MatDestroy_EPSLyapIIOperator(Mat M)
106 : {
107 3 : EPS_LYAPII_MATSHELL *matctx;
108 :
109 3 : PetscFunctionBegin;
110 3 : PetscCall(MatShellGetContext(M,&matctx));
111 3 : PetscCall(MatDestroy(&matctx->S));
112 3 : PetscCall(PetscFree(matctx));
113 3 : PetscFunctionReturn(PETSC_SUCCESS);
114 : }
115 :
116 704 : static PetscErrorCode MatMult_EigOperator(Mat M,Vec x,Vec y)
117 : {
118 704 : EPS_EIG_MATSHELL *matctx;
119 : #if !defined(PETSC_USE_COMPLEX)
120 : PetscInt n,lds;
121 : PetscScalar *Y,*C,zero=0.0,done=1.0,dtwo=2.0;
122 : const PetscScalar *S,*X;
123 : PetscBLASInt n_,lds_;
124 : #endif
125 :
126 704 : PetscFunctionBegin;
127 704 : PetscCall(MatShellGetContext(M,&matctx));
128 :
129 : #if defined(PETSC_USE_COMPLEX)
130 704 : PetscCall(MatMult(matctx->B,x,matctx->w));
131 704 : PetscCall(MatSolve(matctx->F,matctx->w,y));
132 : #else
133 : PetscCall(VecGetArrayRead(x,&X));
134 : PetscCall(VecGetArray(y,&Y));
135 : PetscCall(MatDenseGetArrayRead(matctx->S,&S));
136 : PetscCall(MatDenseGetLDA(matctx->S,&lds));
137 :
138 : n = matctx->n;
139 : PetscCall(PetscCalloc1(n*n,&C));
140 : PetscCall(PetscBLASIntCast(n,&n_));
141 : PetscCall(PetscBLASIntCast(lds,&lds_));
142 :
143 : /* C = 2*S*X*S.' */
144 : PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&dtwo,S,&lds_,X,&n_,&zero,Y,&n_));
145 : PetscCallBLAS("BLASgemm",BLASgemm_("N","T",&n_,&n_,&n_,&done,Y,&n_,S,&lds_,&zero,C,&n_));
146 :
147 : /* Solve S*Y + Y*S' = -C */
148 : PetscCall(LMEDenseLyapunov(matctx->lme,n,(PetscScalar*)S,lds,C,n,Y,n));
149 :
150 : PetscCall(PetscFree(C));
151 : PetscCall(VecRestoreArrayRead(x,&X));
152 : PetscCall(VecRestoreArray(y,&Y));
153 : PetscCall(MatDenseRestoreArrayRead(matctx->S,&S));
154 : #endif
155 704 : PetscFunctionReturn(PETSC_SUCCESS);
156 : }
157 :
158 3 : static PetscErrorCode MatDestroy_EigOperator(Mat M)
159 : {
160 3 : EPS_EIG_MATSHELL *matctx;
161 :
162 3 : PetscFunctionBegin;
163 3 : PetscCall(MatShellGetContext(M,&matctx));
164 : #if defined(PETSC_USE_COMPLEX)
165 3 : PetscCall(MatDestroy(&matctx->A));
166 3 : PetscCall(MatDestroy(&matctx->B));
167 3 : PetscCall(MatDestroy(&matctx->F));
168 3 : PetscCall(VecDestroy(&matctx->w));
169 : #else
170 : PetscCall(MatDestroy(&matctx->S));
171 : #endif
172 3 : PetscCall(PetscFree(matctx));
173 3 : PetscFunctionReturn(PETSC_SUCCESS);
174 : }
175 :
176 : /*
177 : EV2x2: solve the eigenproblem for a 2x2 matrix M
178 : */
179 25 : static PetscErrorCode EV2x2(PetscScalar *M,PetscInt ld,PetscScalar *wr,PetscScalar *wi,PetscScalar *vec)
180 : {
181 25 : PetscBLASInt lwork=10,ld_;
182 25 : PetscScalar work[10];
183 25 : PetscBLASInt two=2,info;
184 : #if defined(PETSC_USE_COMPLEX)
185 25 : PetscReal rwork[6];
186 : #endif
187 :
188 25 : PetscFunctionBegin;
189 25 : PetscCall(PetscBLASIntCast(ld,&ld_));
190 25 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
191 : #if !defined(PETSC_USE_COMPLEX)
192 : PetscCallBLAS("LAPACKgeev",LAPACKgeev_("N","V",&two,M,&ld_,wr,wi,NULL,&ld_,vec,&ld_,work,&lwork,&info));
193 : #else
194 25 : PetscCallBLAS("LAPACKgeev",LAPACKgeev_("N","V",&two,M,&ld_,wr,NULL,&ld_,vec,&ld_,work,&lwork,rwork,&info));
195 : #endif
196 25 : SlepcCheckLapackInfo("geev",info);
197 25 : PetscCall(PetscFPTrapPop());
198 25 : PetscFunctionReturn(PETSC_SUCCESS);
199 : }
200 :
201 : /*
202 : LyapIIBuildRHS: prepare the right-hand side of the Lyapunov equation SY + YS' = -2*S*Z*S'
203 : in factored form:
204 : if (V) U=sqrt(2)*S*V (uses 1 work vector)
205 : else U=sqrt(2)*S*U (uses 2 work vectors)
206 : where U,V are assumed to have rk columns.
207 : */
208 28 : static PetscErrorCode LyapIIBuildRHS(Mat S,PetscInt rk,Mat U,BV V,Vec *work)
209 : {
210 28 : PetscScalar *array,*uu;
211 28 : PetscInt i,nloc;
212 28 : Vec v,u=work[0];
213 :
214 28 : PetscFunctionBegin;
215 28 : PetscCall(MatGetLocalSize(U,&nloc,NULL));
216 75 : for (i=0;i<rk;i++) {
217 47 : PetscCall(MatDenseGetColumn(U,i,&array));
218 47 : if (V) PetscCall(BVGetColumn(V,i,&v));
219 : else {
220 40 : v = work[1];
221 40 : PetscCall(VecPlaceArray(v,array));
222 : }
223 47 : PetscCall(MatMult(S,v,u));
224 47 : if (V) PetscCall(BVRestoreColumn(V,i,&v));
225 40 : else PetscCall(VecResetArray(v));
226 47 : PetscCall(VecScale(u,PETSC_SQRT2));
227 47 : PetscCall(VecGetArray(u,&uu));
228 47 : PetscCall(PetscArraycpy(array,uu,nloc));
229 47 : PetscCall(VecRestoreArray(u,&uu));
230 47 : PetscCall(MatDenseRestoreColumn(U,&array));
231 : }
232 28 : PetscFunctionReturn(PETSC_SUCCESS);
233 : }
234 :
235 : /*
236 : LyapIIBuildEigenMat: create shell matrix Op=A\B with A = kron(I,S)+kron(S,I), B = -2*kron(S,S)
237 : where S is a sequential square dense matrix of order n.
238 : v0 is the initial vector, should have the form v0 = w*w' (for instance 1*1')
239 : */
240 28 : static PetscErrorCode LyapIIBuildEigenMat(LME lme,Mat S,Mat *Op,Vec *v0)
241 : {
242 28 : PetscInt n,m;
243 28 : PetscBool create=PETSC_FALSE;
244 28 : EPS_EIG_MATSHELL *matctx;
245 : #if defined(PETSC_USE_COMPLEX)
246 28 : PetscScalar theta,*aa,*bb;
247 28 : const PetscScalar *ss;
248 28 : PetscInt i,j,f,c,off,ld,lds;
249 28 : IS perm;
250 : #endif
251 :
252 28 : PetscFunctionBegin;
253 28 : PetscCall(MatGetSize(S,&n,NULL));
254 28 : if (!*Op) create=PETSC_TRUE;
255 : else {
256 25 : PetscCall(MatGetSize(*Op,&m,NULL));
257 25 : if (m!=n*n) create=PETSC_TRUE;
258 : }
259 25 : if (create) {
260 3 : PetscCall(MatDestroy(Op));
261 3 : PetscCall(VecDestroy(v0));
262 3 : PetscCall(PetscNew(&matctx));
263 : #if defined(PETSC_USE_COMPLEX)
264 3 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,n*n,n*n,NULL,&matctx->A));
265 3 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,n*n,n*n,NULL,&matctx->B));
266 3 : PetscCall(MatCreateVecs(matctx->A,NULL,&matctx->w));
267 : #else
268 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&matctx->S));
269 : #endif
270 3 : PetscCall(MatCreateShell(PETSC_COMM_SELF,n*n,n*n,PETSC_DETERMINE,PETSC_DETERMINE,matctx,Op));
271 3 : PetscCall(MatShellSetOperation(*Op,MATOP_MULT,(void(*)(void))MatMult_EigOperator));
272 3 : PetscCall(MatShellSetOperation(*Op,MATOP_DESTROY,(void(*)(void))MatDestroy_EigOperator));
273 3 : PetscCall(MatCreateVecs(*Op,NULL,v0));
274 : } else {
275 25 : PetscCall(MatShellGetContext(*Op,&matctx));
276 : #if defined(PETSC_USE_COMPLEX)
277 25 : PetscCall(MatZeroEntries(matctx->A));
278 : #endif
279 : }
280 : #if defined(PETSC_USE_COMPLEX)
281 28 : PetscCall(MatDenseGetArray(matctx->A,&aa));
282 28 : PetscCall(MatDenseGetArray(matctx->B,&bb));
283 28 : PetscCall(MatDenseGetArrayRead(S,&ss));
284 28 : PetscCall(MatDenseGetLDA(S,&lds));
285 28 : ld = n*n;
286 282 : for (f=0;f<n;f++) {
287 254 : off = f*n+f*n*ld;
288 24242 : for (i=0;i<n;i++) for (j=0;j<n;j++) aa[off+i+j*ld] = ss[i+j*lds];
289 2586 : for (c=0;c<n;c++) {
290 2332 : off = f*n+c*n*ld;
291 2332 : theta = ss[f+c*lds];
292 23988 : for (i=0;i<n;i++) aa[off+i+i*ld] += theta;
293 227236 : for (i=0;i<n;i++) for (j=0;j<n;j++) bb[off+i+j*ld] = -2*theta*ss[i+j*lds];
294 : }
295 : }
296 28 : PetscCall(MatDenseRestoreArray(matctx->A,&aa));
297 28 : PetscCall(MatDenseRestoreArray(matctx->B,&bb));
298 28 : PetscCall(MatDenseRestoreArrayRead(S,&ss));
299 28 : PetscCall(ISCreateStride(PETSC_COMM_SELF,n*n,0,1,&perm));
300 28 : PetscCall(MatDestroy(&matctx->F));
301 28 : PetscCall(MatDuplicate(matctx->A,MAT_COPY_VALUES,&matctx->F));
302 28 : PetscCall(MatLUFactor(matctx->F,perm,perm,NULL));
303 28 : PetscCall(ISDestroy(&perm));
304 : #else
305 : PetscCall(MatCopy(S,matctx->S,SAME_NONZERO_PATTERN));
306 : #endif
307 28 : matctx->lme = lme;
308 28 : matctx->n = n;
309 28 : PetscCall(VecSet(*v0,1.0));
310 28 : PetscFunctionReturn(PETSC_SUCCESS);
311 : }
312 :
313 3 : static PetscErrorCode EPSSolve_LyapII(EPS eps)
314 : {
315 3 : EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
316 3 : PetscInt i,ldds,rk,nloc,mloc,nv,idx,k;
317 3 : Vec v,w,z=eps->work[0],v0=NULL;
318 3 : Mat S,C,Ux[2],Y,Y1,R,U,W,X,Op=NULL;
319 3 : BV V;
320 3 : BVOrthogType type;
321 3 : BVOrthogRefineType refine;
322 3 : PetscScalar eigr[2],eigi[2],*array,er,ei,*uu,*s,*xx,*aa,pM[4],vec[4];
323 3 : PetscReal eta;
324 3 : EPS epsrr;
325 3 : PetscReal norm;
326 3 : EPS_LYAPII_MATSHELL *matctx;
327 :
328 3 : PetscFunctionBegin;
329 3 : PetscCall(DSGetLeadingDimension(ctx->ds,&ldds));
330 :
331 : /* Operator for the Lyapunov equation */
332 3 : PetscCall(PetscNew(&matctx));
333 3 : PetscCall(STGetOperator(eps->st,&matctx->S));
334 3 : PetscCall(MatGetLocalSize(matctx->S,&mloc,&nloc));
335 3 : PetscCall(MatCreateShell(PetscObjectComm((PetscObject)eps),mloc,nloc,PETSC_DETERMINE,PETSC_DETERMINE,matctx,&S));
336 3 : matctx->Q = eps->V;
337 3 : PetscCall(MatShellSetOperation(S,MATOP_MULT,(void(*)(void))MatMult_EPSLyapIIOperator));
338 3 : PetscCall(MatShellSetOperation(S,MATOP_DESTROY,(void(*)(void))MatDestroy_EPSLyapIIOperator));
339 3 : PetscCall(LMESetCoefficients(ctx->lme,S,NULL,NULL,NULL));
340 :
341 : /* Right-hand side */
342 3 : PetscCall(BVDuplicateResize(eps->V,ctx->rkl,&V));
343 3 : PetscCall(BVGetOrthogonalization(V,&type,&refine,&eta,NULL));
344 3 : PetscCall(BVSetOrthogonalization(V,type,refine,eta,BV_ORTHOG_BLOCK_TSQR));
345 3 : PetscCall(MatCreateDense(PetscObjectComm((PetscObject)eps),eps->nloc,PETSC_DECIDE,PETSC_DECIDE,1,NULL,&Ux[0]));
346 3 : PetscCall(MatCreateDense(PetscObjectComm((PetscObject)eps),eps->nloc,PETSC_DECIDE,PETSC_DECIDE,2,NULL,&Ux[1]));
347 3 : nv = ctx->rkl;
348 3 : PetscCall(PetscMalloc1(nv,&s));
349 :
350 : /* Initialize first column */
351 3 : PetscCall(EPSGetStartVector(eps,0,NULL));
352 3 : PetscCall(BVGetColumn(eps->V,0,&v));
353 3 : PetscCall(BVInsertVec(V,0,v));
354 3 : PetscCall(BVRestoreColumn(eps->V,0,&v));
355 3 : PetscCall(BVSetActiveColumns(eps->V,0,0)); /* no deflation at the beginning */
356 3 : PetscCall(LyapIIBuildRHS(S,1,Ux[0],V,eps->work));
357 3 : idx = 0;
358 :
359 : /* EPS for rank reduction */
360 3 : PetscCall(EPSCreate(PETSC_COMM_SELF,&epsrr));
361 3 : PetscCall(EPSSetOptionsPrefix(epsrr,((PetscObject)eps)->prefix));
362 3 : PetscCall(EPSAppendOptionsPrefix(epsrr,"eps_lyapii_"));
363 3 : PetscCall(EPSSetDimensions(epsrr,1,PETSC_CURRENT,PETSC_CURRENT));
364 3 : PetscCall(EPSSetTolerances(epsrr,PETSC_MACHINE_EPSILON*100,PETSC_CURRENT));
365 :
366 31 : while (eps->reason == EPS_CONVERGED_ITERATING) {
367 28 : eps->its++;
368 :
369 : /* Matrix for placing the solution of the Lyapunov equation (an alias of V) */
370 28 : PetscCall(BVSetActiveColumns(V,0,nv));
371 28 : PetscCall(BVGetMat(V,&Y1));
372 28 : PetscCall(MatZeroEntries(Y1));
373 28 : PetscCall(MatCreateLRC(NULL,Y1,NULL,NULL,&Y));
374 28 : PetscCall(LMESetSolution(ctx->lme,Y));
375 :
376 : /* Solve the Lyapunov equation SY + YS' = -2*S*Z*S' */
377 28 : PetscCall(MatCreateLRC(NULL,Ux[idx],NULL,NULL,&C));
378 28 : PetscCall(LMESetRHS(ctx->lme,C));
379 28 : PetscCall(MatDestroy(&C));
380 28 : PetscCall(LMESolve(ctx->lme));
381 28 : PetscCall(BVRestoreMat(V,&Y1));
382 28 : PetscCall(MatDestroy(&Y));
383 :
384 : /* SVD of the solution: [Q,R]=qr(V); [U,Sigma,~]=svd(R) */
385 28 : PetscCall(DSSetDimensions(ctx->ds,nv,0,0));
386 28 : PetscCall(DSSVDSetDimensions(ctx->ds,nv));
387 28 : PetscCall(DSGetMat(ctx->ds,DS_MAT_A,&R));
388 28 : PetscCall(BVOrthogonalize(V,R));
389 28 : PetscCall(DSRestoreMat(ctx->ds,DS_MAT_A,&R));
390 28 : PetscCall(DSSetState(ctx->ds,DS_STATE_RAW));
391 28 : PetscCall(DSSolve(ctx->ds,s,NULL));
392 :
393 : /* Determine rank */
394 519 : rk = nv;
395 519 : for (i=1;i<nv;i++) if (PetscAbsScalar(s[i]/s[0])<PETSC_SQRT_MACHINE_EPSILON) {rk=i; break;}
396 28 : PetscCall(PetscInfo(eps,"The computed solution of the Lyapunov equation has rank %" PetscInt_FMT "\n",rk));
397 28 : rk = PetscMin(rk,ctx->rkc);
398 28 : PetscCall(DSGetMat(ctx->ds,DS_MAT_U,&U));
399 28 : PetscCall(BVMultInPlace(V,U,0,rk));
400 28 : PetscCall(DSRestoreMat(ctx->ds,DS_MAT_U,&U));
401 28 : PetscCall(BVSetActiveColumns(V,0,rk));
402 :
403 : /* Rank reduction */
404 28 : PetscCall(DSSetDimensions(ctx->ds,rk,0,0));
405 28 : PetscCall(DSSVDSetDimensions(ctx->ds,rk));
406 28 : PetscCall(DSGetMat(ctx->ds,DS_MAT_A,&W));
407 28 : PetscCall(BVMatProject(V,S,V,W));
408 28 : PetscCall(LyapIIBuildEigenMat(ctx->lme,W,&Op,&v0)); /* Op=A\B, A=kron(I,S)+kron(S,I), B=-2*kron(S,S) */
409 28 : PetscCall(DSRestoreMat(ctx->ds,DS_MAT_A,&W));
410 28 : PetscCall(EPSSetOperators(epsrr,Op,NULL));
411 28 : PetscCall(EPSSetInitialSpace(epsrr,1,&v0));
412 28 : PetscCall(EPSSolve(epsrr));
413 28 : PetscCall(EPSComputeVectors(epsrr));
414 : /* Copy first eigenvector, vec(A)=x */
415 28 : PetscCall(BVGetArray(epsrr->V,&xx));
416 28 : PetscCall(DSGetArray(ctx->ds,DS_MAT_A,&aa));
417 282 : for (i=0;i<rk;i++) PetscCall(PetscArraycpy(aa+i*ldds,xx+i*rk,rk));
418 28 : PetscCall(DSRestoreArray(ctx->ds,DS_MAT_A,&aa));
419 28 : PetscCall(BVRestoreArray(epsrr->V,&xx));
420 28 : PetscCall(DSSetState(ctx->ds,DS_STATE_RAW));
421 : /* Compute [U,Sigma,~] = svd(A), its rank should be 1 or 2 */
422 28 : PetscCall(DSSolve(ctx->ds,s,NULL));
423 28 : if (PetscAbsScalar(s[1]/s[0])<PETSC_SQRT_MACHINE_EPSILON) rk=1;
424 25 : else rk = 2;
425 28 : PetscCall(PetscInfo(eps,"The eigenvector has rank %" PetscInt_FMT "\n",rk));
426 28 : PetscCall(DSGetMat(ctx->ds,DS_MAT_U,&U));
427 28 : PetscCall(BVMultInPlace(V,U,0,rk));
428 28 : PetscCall(DSRestoreMat(ctx->ds,DS_MAT_U,&U));
429 :
430 : /* Save V in Ux */
431 28 : idx = (rk==2)?1:0;
432 81 : for (i=0;i<rk;i++) {
433 53 : PetscCall(BVGetColumn(V,i,&v));
434 53 : PetscCall(VecGetArray(v,&uu));
435 53 : PetscCall(MatDenseGetColumn(Ux[idx],i,&array));
436 53 : PetscCall(PetscArraycpy(array,uu,eps->nloc));
437 53 : PetscCall(MatDenseRestoreColumn(Ux[idx],&array));
438 53 : PetscCall(VecRestoreArray(v,&uu));
439 53 : PetscCall(BVRestoreColumn(V,i,&v));
440 : }
441 :
442 : /* Eigenpair approximation */
443 28 : PetscCall(BVGetColumn(V,0,&v));
444 28 : PetscCall(MatMult(S,v,z));
445 28 : PetscCall(VecDot(z,v,pM));
446 28 : PetscCall(BVRestoreColumn(V,0,&v));
447 28 : if (rk>1) {
448 25 : PetscCall(BVGetColumn(V,1,&w));
449 25 : PetscCall(VecDot(z,w,pM+1));
450 25 : PetscCall(MatMult(S,w,z));
451 25 : PetscCall(VecDot(z,w,pM+3));
452 25 : PetscCall(BVGetColumn(V,0,&v));
453 25 : PetscCall(VecDot(z,v,pM+2));
454 25 : PetscCall(BVRestoreColumn(V,0,&v));
455 25 : PetscCall(BVRestoreColumn(V,1,&w));
456 25 : PetscCall(EV2x2(pM,2,eigr,eigi,vec));
457 25 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,2,2,vec,&X));
458 25 : PetscCall(BVSetActiveColumns(V,0,rk));
459 25 : PetscCall(BVMultInPlace(V,X,0,rk));
460 25 : PetscCall(MatDestroy(&X));
461 : #if !defined(PETSC_USE_COMPLEX)
462 : norm = eigr[0]*eigr[0]+eigi[0]*eigi[0];
463 : er = eigr[0]/norm; ei = -eigi[0]/norm;
464 : #else
465 25 : er =1.0/eigr[0]; ei = 0.0;
466 : #endif
467 : } else {
468 3 : eigr[0] = pM[0]; eigi[0] = 0.0;
469 3 : er = 1.0/eigr[0]; ei = 0.0;
470 : }
471 28 : PetscCall(BVGetColumn(V,0,&v));
472 28 : if (eigi[0]!=0.0) PetscCall(BVGetColumn(V,1,&w));
473 28 : else w = NULL;
474 28 : eps->eigr[eps->nconv] = eigr[0]; eps->eigi[eps->nconv] = eigi[0];
475 28 : PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_FALSE,er,ei,v,w,eps->work,&norm));
476 28 : PetscCall(BVRestoreColumn(V,0,&v));
477 28 : if (w) PetscCall(BVRestoreColumn(V,1,&w));
478 28 : PetscCall((*eps->converged)(eps,er,ei,norm,&eps->errest[eps->nconv],eps->convergedctx));
479 28 : k = 0;
480 28 : if (eps->errest[eps->nconv]<eps->tol) {
481 7 : k++;
482 7 : if (rk==2) {
483 : #if !defined (PETSC_USE_COMPLEX)
484 : eps->eigr[eps->nconv+k] = eigr[0]; eps->eigi[eps->nconv+k] = -eigi[0];
485 : #else
486 6 : eps->eigr[eps->nconv+k] = PetscConj(eps->eigr[eps->nconv]);
487 : #endif
488 6 : k++;
489 : }
490 : /* Store converged eigenpairs and vectors for deflation */
491 20 : for (i=0;i<k;i++) {
492 13 : PetscCall(BVGetColumn(V,i,&v));
493 13 : PetscCall(BVInsertVec(eps->V,eps->nconv+i,v));
494 13 : PetscCall(BVRestoreColumn(V,i,&v));
495 : }
496 7 : eps->nconv += k;
497 7 : PetscCall(BVSetActiveColumns(eps->V,eps->nconv-rk,eps->nconv));
498 7 : PetscCall(BVOrthogonalize(eps->V,NULL));
499 7 : PetscCall(DSSetDimensions(eps->ds,eps->nconv,0,0));
500 7 : PetscCall(DSGetMat(eps->ds,DS_MAT_A,&W));
501 7 : PetscCall(BVMatProject(eps->V,matctx->S,eps->V,W));
502 7 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&W));
503 7 : if (eps->nconv<eps->nev) {
504 4 : idx = 0;
505 4 : PetscCall(BVSetRandomColumn(V,0));
506 4 : PetscCall(BVNormColumn(V,0,NORM_2,&norm));
507 4 : PetscCall(BVScaleColumn(V,0,1.0/norm));
508 4 : PetscCall(LyapIIBuildRHS(S,1,Ux[idx],V,eps->work));
509 : }
510 : } else {
511 : /* Prepare right-hand side */
512 21 : PetscCall(LyapIIBuildRHS(S,rk,Ux[idx],NULL,eps->work));
513 : }
514 28 : PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx));
515 31 : PetscCall(EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,eps->nconv+1));
516 : }
517 3 : PetscCall(STRestoreOperator(eps->st,&matctx->S));
518 3 : PetscCall(MatDestroy(&S));
519 3 : PetscCall(MatDestroy(&Ux[0]));
520 3 : PetscCall(MatDestroy(&Ux[1]));
521 3 : PetscCall(MatDestroy(&Op));
522 3 : PetscCall(VecDestroy(&v0));
523 3 : PetscCall(BVDestroy(&V));
524 3 : PetscCall(EPSDestroy(&epsrr));
525 3 : PetscCall(PetscFree(s));
526 3 : PetscFunctionReturn(PETSC_SUCCESS);
527 : }
528 :
529 3 : static PetscErrorCode EPSSetFromOptions_LyapII(EPS eps,PetscOptionItems *PetscOptionsObject)
530 : {
531 3 : EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
532 3 : PetscInt k,array[2]={PETSC_DETERMINE,PETSC_DETERMINE};
533 3 : PetscBool flg;
534 :
535 3 : PetscFunctionBegin;
536 3 : PetscOptionsHeadBegin(PetscOptionsObject,"EPS Lyapunov Inverse Iteration Options");
537 :
538 3 : k = 2;
539 3 : PetscCall(PetscOptionsIntArray("-eps_lyapii_ranks","Ranks for Lyapunov equation (one or two comma-separated integers)","EPSLyapIISetRanks",array,&k,&flg));
540 3 : if (flg) PetscCall(EPSLyapIISetRanks(eps,array[0],array[1]));
541 :
542 3 : PetscOptionsHeadEnd();
543 :
544 3 : if (!ctx->lme) PetscCall(EPSLyapIIGetLME(eps,&ctx->lme));
545 3 : PetscCall(LMESetFromOptions(ctx->lme));
546 3 : PetscFunctionReturn(PETSC_SUCCESS);
547 : }
548 :
549 1 : static PetscErrorCode EPSLyapIISetRanks_LyapII(EPS eps,PetscInt rkc,PetscInt rkl)
550 : {
551 1 : EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
552 :
553 1 : PetscFunctionBegin;
554 1 : if (rkc==PETSC_DETERMINE) {
555 0 : if (ctx->rkc != 10) eps->state = EPS_STATE_INITIAL;
556 0 : ctx->rkc = 10;
557 1 : } else if (rkc!=PETSC_CURRENT) {
558 1 : PetscCheck(rkc>1,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The compressed rank %" PetscInt_FMT " must be larger than 1",rkc);
559 1 : if (ctx->rkc != rkc) eps->state = EPS_STATE_INITIAL;
560 1 : ctx->rkc = rkc;
561 : }
562 1 : if (rkl==PETSC_DETERMINE) {
563 0 : if (ctx->rkl != 3*rkc) eps->state = EPS_STATE_INITIAL;
564 0 : ctx->rkl = 3*rkc;
565 1 : } else if (rkl!=PETSC_CURRENT) {
566 1 : 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);
567 1 : if (ctx->rkl != rkl) eps->state = EPS_STATE_INITIAL;
568 1 : ctx->rkl = rkl;
569 : }
570 1 : PetscFunctionReturn(PETSC_SUCCESS);
571 : }
572 :
573 : /*@
574 : EPSLyapIISetRanks - Set the ranks used in the solution of the Lyapunov equation.
575 :
576 : Logically Collective
577 :
578 : Input Parameters:
579 : + eps - the eigenproblem solver context
580 : . rkc - the compressed rank
581 : - rkl - the Lyapunov rank
582 :
583 : Options Database Key:
584 : . -eps_lyapii_ranks <rkc,rkl> - Sets the rank parameters
585 :
586 : Notes:
587 : PETSC_CURRENT can be used to preserve the current value of any of the
588 : arguments, and PETSC_DETERMINE to set them to a default value.
589 :
590 : Lyapunov inverse iteration needs to solve a large-scale Lyapunov equation
591 : at each iteration of the eigensolver. For this, an iterative solver (LME)
592 : is used, which requires to prescribe the rank of the solution matrix X. This
593 : is the meaning of parameter rkl. Later, this matrix is compressed into
594 : another matrix of rank rkc. If not provided, rkl is a small multiple of rkc.
595 :
596 : Level: intermediate
597 :
598 : .seealso: EPSLyapIIGetRanks()
599 : @*/
600 1 : PetscErrorCode EPSLyapIISetRanks(EPS eps,PetscInt rkc,PetscInt rkl)
601 : {
602 1 : PetscFunctionBegin;
603 1 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
604 3 : PetscValidLogicalCollectiveInt(eps,rkc,2);
605 3 : PetscValidLogicalCollectiveInt(eps,rkl,3);
606 1 : PetscTryMethod(eps,"EPSLyapIISetRanks_C",(EPS,PetscInt,PetscInt),(eps,rkc,rkl));
607 1 : PetscFunctionReturn(PETSC_SUCCESS);
608 : }
609 :
610 1 : static PetscErrorCode EPSLyapIIGetRanks_LyapII(EPS eps,PetscInt *rkc,PetscInt *rkl)
611 : {
612 1 : EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
613 :
614 1 : PetscFunctionBegin;
615 1 : if (rkc) *rkc = ctx->rkc;
616 1 : if (rkl) *rkl = ctx->rkl;
617 1 : PetscFunctionReturn(PETSC_SUCCESS);
618 : }
619 :
620 : /*@
621 : EPSLyapIIGetRanks - Return the rank values used for the Lyapunov step.
622 :
623 : Not Collective
624 :
625 : Input Parameter:
626 : . eps - the eigenproblem solver context
627 :
628 : Output Parameters:
629 : + rkc - the compressed rank
630 : - rkl - the Lyapunov rank
631 :
632 : Level: intermediate
633 :
634 : .seealso: EPSLyapIISetRanks()
635 : @*/
636 1 : PetscErrorCode EPSLyapIIGetRanks(EPS eps,PetscInt *rkc,PetscInt *rkl)
637 : {
638 1 : PetscFunctionBegin;
639 1 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
640 1 : PetscUseMethod(eps,"EPSLyapIIGetRanks_C",(EPS,PetscInt*,PetscInt*),(eps,rkc,rkl));
641 1 : PetscFunctionReturn(PETSC_SUCCESS);
642 : }
643 :
644 0 : static PetscErrorCode EPSLyapIISetLME_LyapII(EPS eps,LME lme)
645 : {
646 0 : EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
647 :
648 0 : PetscFunctionBegin;
649 0 : PetscCall(PetscObjectReference((PetscObject)lme));
650 0 : PetscCall(LMEDestroy(&ctx->lme));
651 0 : ctx->lme = lme;
652 0 : eps->state = EPS_STATE_INITIAL;
653 0 : PetscFunctionReturn(PETSC_SUCCESS);
654 : }
655 :
656 : /*@
657 : EPSLyapIISetLME - Associate a linear matrix equation solver object (LME) to the
658 : eigenvalue solver.
659 :
660 : Collective
661 :
662 : Input Parameters:
663 : + eps - the eigenproblem solver context
664 : - lme - the linear matrix equation solver object
665 :
666 : Level: advanced
667 :
668 : .seealso: EPSLyapIIGetLME()
669 : @*/
670 0 : PetscErrorCode EPSLyapIISetLME(EPS eps,LME lme)
671 : {
672 0 : PetscFunctionBegin;
673 0 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
674 0 : PetscValidHeaderSpecific(lme,LME_CLASSID,2);
675 0 : PetscCheckSameComm(eps,1,lme,2);
676 0 : PetscTryMethod(eps,"EPSLyapIISetLME_C",(EPS,LME),(eps,lme));
677 0 : PetscFunctionReturn(PETSC_SUCCESS);
678 : }
679 :
680 3 : static PetscErrorCode EPSLyapIIGetLME_LyapII(EPS eps,LME *lme)
681 : {
682 3 : EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
683 :
684 3 : PetscFunctionBegin;
685 3 : if (!ctx->lme) {
686 3 : PetscCall(LMECreate(PetscObjectComm((PetscObject)eps),&ctx->lme));
687 3 : PetscCall(LMESetOptionsPrefix(ctx->lme,((PetscObject)eps)->prefix));
688 3 : PetscCall(LMEAppendOptionsPrefix(ctx->lme,"eps_lyapii_"));
689 3 : PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->lme,(PetscObject)eps,1));
690 : }
691 3 : *lme = ctx->lme;
692 3 : PetscFunctionReturn(PETSC_SUCCESS);
693 : }
694 :
695 : /*@
696 : EPSLyapIIGetLME - Retrieve the linear matrix equation solver object (LME)
697 : associated with the eigenvalue solver.
698 :
699 : Not Collective
700 :
701 : Input Parameter:
702 : . eps - the eigenproblem solver context
703 :
704 : Output Parameter:
705 : . lme - the linear matrix equation solver object
706 :
707 : Level: advanced
708 :
709 : .seealso: EPSLyapIISetLME()
710 : @*/
711 3 : PetscErrorCode EPSLyapIIGetLME(EPS eps,LME *lme)
712 : {
713 3 : PetscFunctionBegin;
714 3 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
715 3 : PetscAssertPointer(lme,2);
716 3 : PetscUseMethod(eps,"EPSLyapIIGetLME_C",(EPS,LME*),(eps,lme));
717 3 : PetscFunctionReturn(PETSC_SUCCESS);
718 : }
719 :
720 1 : static PetscErrorCode EPSView_LyapII(EPS eps,PetscViewer viewer)
721 : {
722 1 : EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
723 1 : PetscBool isascii;
724 :
725 1 : PetscFunctionBegin;
726 1 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
727 1 : if (isascii) {
728 1 : PetscCall(PetscViewerASCIIPrintf(viewer," ranks: for Lyapunov solver=%" PetscInt_FMT ", after compression=%" PetscInt_FMT "\n",ctx->rkl,ctx->rkc));
729 1 : if (!ctx->lme) PetscCall(EPSLyapIIGetLME(eps,&ctx->lme));
730 1 : PetscCall(PetscViewerASCIIPushTab(viewer));
731 1 : PetscCall(LMEView(ctx->lme,viewer));
732 1 : PetscCall(PetscViewerASCIIPopTab(viewer));
733 : }
734 1 : PetscFunctionReturn(PETSC_SUCCESS);
735 : }
736 :
737 3 : static PetscErrorCode EPSReset_LyapII(EPS eps)
738 : {
739 3 : EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
740 :
741 3 : PetscFunctionBegin;
742 3 : if (!ctx->lme) PetscCall(LMEReset(ctx->lme));
743 3 : PetscFunctionReturn(PETSC_SUCCESS);
744 : }
745 :
746 3 : static PetscErrorCode EPSDestroy_LyapII(EPS eps)
747 : {
748 3 : EPS_LYAPII *ctx = (EPS_LYAPII*)eps->data;
749 :
750 3 : PetscFunctionBegin;
751 3 : PetscCall(LMEDestroy(&ctx->lme));
752 3 : PetscCall(DSDestroy(&ctx->ds));
753 3 : PetscCall(PetscFree(eps->data));
754 3 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetLME_C",NULL));
755 3 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetLME_C",NULL));
756 3 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetRanks_C",NULL));
757 3 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetRanks_C",NULL));
758 3 : PetscFunctionReturn(PETSC_SUCCESS);
759 : }
760 :
761 6 : static PetscErrorCode EPSSetDefaultST_LyapII(EPS eps)
762 : {
763 6 : PetscFunctionBegin;
764 6 : if (!((PetscObject)eps->st)->type_name) PetscCall(STSetType(eps->st,STSINVERT));
765 6 : PetscFunctionReturn(PETSC_SUCCESS);
766 : }
767 :
768 3 : SLEPC_EXTERN PetscErrorCode EPSCreate_LyapII(EPS eps)
769 : {
770 3 : EPS_LYAPII *ctx;
771 :
772 3 : PetscFunctionBegin;
773 3 : PetscCall(PetscNew(&ctx));
774 3 : eps->data = (void*)ctx;
775 :
776 3 : eps->useds = PETSC_TRUE;
777 :
778 3 : eps->ops->solve = EPSSolve_LyapII;
779 3 : eps->ops->setup = EPSSetUp_LyapII;
780 3 : eps->ops->setupsort = EPSSetUpSort_Default;
781 3 : eps->ops->setfromoptions = EPSSetFromOptions_LyapII;
782 3 : eps->ops->reset = EPSReset_LyapII;
783 3 : eps->ops->destroy = EPSDestroy_LyapII;
784 3 : eps->ops->view = EPSView_LyapII;
785 3 : eps->ops->setdefaultst = EPSSetDefaultST_LyapII;
786 3 : eps->ops->backtransform = EPSBackTransform_Default;
787 3 : eps->ops->computevectors = EPSComputeVectors_Schur;
788 :
789 3 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetLME_C",EPSLyapIISetLME_LyapII));
790 3 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetLME_C",EPSLyapIIGetLME_LyapII));
791 3 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIISetRanks_C",EPSLyapIISetRanks_LyapII));
792 3 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLyapIIGetRanks_C",EPSLyapIIGetRanks_LyapII));
793 3 : PetscFunctionReturn(PETSC_SUCCESS);
794 : }
|
