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 : Routines for solving dense matrix equations, in some cases calling SLICOT
12 : */
13 :
14 : #include <slepc/private/lmeimpl.h> /*I "slepclme.h" I*/
15 : #include <slepcblaslapack.h>
16 :
17 : /*
18 : LMEDenseRankSVD - given a square matrix A, compute its SVD U*S*V', and determine the
19 : numerical rank. On exit, U contains U*S and A is overwritten with V'
20 : */
21 51 : PetscErrorCode LMEDenseRankSVD(LME lme,PetscInt n,PetscScalar *A,PetscInt lda,PetscScalar *U,PetscInt ldu,PetscInt *rank)
22 : {
23 51 : PetscInt i,j,rk=0;
24 51 : PetscScalar *work;
25 51 : PetscReal tol,*sg,*rwork;
26 51 : PetscBLASInt n_,lda_,ldu_,info,lw_;
27 :
28 51 : PetscFunctionBegin;
29 51 : PetscCall(PetscCalloc3(n,&sg,10*n,&work,5*n,&rwork));
30 51 : PetscCall(PetscBLASIntCast(n,&n_));
31 51 : PetscCall(PetscBLASIntCast(lda,&lda_));
32 51 : PetscCall(PetscBLASIntCast(ldu,&ldu_));
33 51 : lw_ = 10*n_;
34 51 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
35 : #if !defined (PETSC_USE_COMPLEX)
36 51 : PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","O",&n_,&n_,A,&lda_,sg,U,&ldu_,NULL,&n_,work,&lw_,&info));
37 : #else
38 : PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","O",&n_,&n_,A,&lda_,sg,U,&ldu_,NULL,&n_,work,&lw_,rwork,&info));
39 : #endif
40 51 : PetscCall(PetscFPTrapPop());
41 51 : SlepcCheckLapackInfo("gesvd",info);
42 51 : tol = 10*PETSC_MACHINE_EPSILON*n*sg[0];
43 1629 : for (j=0;j<n;j++) {
44 1578 : if (sg[j]>tol) {
45 50934 : for (i=0;i<n;i++) U[i+j*n] *= sg[j];
46 1578 : rk++;
47 : } else break;
48 : }
49 51 : *rank = rk;
50 51 : PetscCall(PetscFree3(sg,work,rwork));
51 51 : PetscFunctionReturn(PETSC_SUCCESS);
52 : }
53 :
54 : #if defined(PETSC_USE_INFO)
55 : /*
56 : LyapunovCholResidual - compute the residual norm ||A*U'*U+U'*U*A'+B*B'||
57 : */
58 1 : static PetscErrorCode LyapunovCholResidual(PetscInt m,PetscScalar *A,PetscInt lda,PetscInt k,PetscScalar *B,PetscInt ldb,PetscScalar *U,PetscInt ldu,PetscReal *res)
59 : {
60 1 : PetscBLASInt n,kk,la,lb,lu;
61 1 : PetscScalar *M,*R,zero=0.0,done=1.0;
62 :
63 1 : PetscFunctionBegin;
64 1 : *res = 0;
65 1 : PetscCall(PetscBLASIntCast(lda,&la));
66 1 : PetscCall(PetscBLASIntCast(ldb,&lb));
67 1 : PetscCall(PetscBLASIntCast(ldu,&lu));
68 1 : PetscCall(PetscBLASIntCast(m,&n));
69 1 : PetscCall(PetscBLASIntCast(k,&kk));
70 1 : PetscCall(PetscMalloc2(m*m,&M,m*m,&R));
71 :
72 : /* R = B*B' */
73 1 : PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&n,&n,&kk,&done,B,&lb,B,&lb,&zero,R,&n));
74 : /* M = A*U' */
75 1 : PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&n,&n,&n,&done,A,&la,U,&lu,&zero,M,&n));
76 : /* R = R+M*U */
77 1 : PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&done,M,&n,U,&lu,&done,R,&n));
78 : /* R = R+U'*M' */
79 1 : PetscCallBLAS("BLASgemm",BLASgemm_("C","C",&n,&n,&n,&done,U,&lu,M,&n,&done,R,&n));
80 :
81 1 : *res = LAPACKlange_("F",&n,&n,R,&n,NULL);
82 1 : PetscCall(PetscFree2(M,R));
83 1 : PetscFunctionReturn(PETSC_SUCCESS);
84 : }
85 :
86 : /*
87 : LyapunovResidual - compute the residual norm ||A*X+X*A'+B||
88 : */
89 1 : static PetscErrorCode LyapunovResidual(PetscInt m,PetscScalar *A,PetscInt lda,PetscScalar *B,PetscInt ldb,PetscScalar *X,PetscInt ldx,PetscReal *res)
90 : {
91 1 : PetscInt i;
92 1 : PetscBLASInt n,la,lb,lx;
93 1 : PetscScalar *R,done=1.0;
94 :
95 1 : PetscFunctionBegin;
96 1 : *res = 0;
97 1 : PetscCall(PetscBLASIntCast(lda,&la));
98 1 : PetscCall(PetscBLASIntCast(ldb,&lb));
99 1 : PetscCall(PetscBLASIntCast(ldx,&lx));
100 1 : PetscCall(PetscBLASIntCast(m,&n));
101 1 : PetscCall(PetscMalloc1(m*m,&R));
102 :
103 : /* R = B+A*X */
104 11 : for (i=0;i<m;i++) PetscCall(PetscArraycpy(R+i*m,B+i*ldb,m));
105 1 : PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&done,A,&la,X,&lx,&done,R,&n));
106 : /* R = R+X*A' */
107 1 : PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&n,&n,&n,&done,X,&lx,A,&la,&done,R,&n));
108 :
109 1 : *res = LAPACKlange_("F",&n,&n,R,&n,NULL);
110 1 : PetscCall(PetscFree(R));
111 1 : PetscFunctionReturn(PETSC_SUCCESS);
112 : }
113 : #endif
114 :
115 : #if defined(SLEPC_HAVE_SLICOT)
116 : /*
117 : HessLyapunovChol_SLICOT - implementation used when SLICOT is available
118 : */
119 : static PetscErrorCode HessLyapunovChol_SLICOT(PetscInt m,PetscScalar *H,PetscInt ldh,PetscInt k,PetscScalar *B,PetscInt ldb,PetscScalar *U,PetscInt ldu,PetscReal *res)
120 : {
121 : PetscBLASInt lwork,info,n,kk,lu,ione=1,sdim;
122 : PetscInt i,j;
123 : PetscReal scal;
124 : PetscScalar *Q,*W,*wr,*wi,*work;
125 :
126 : PetscFunctionBegin;
127 : PetscCall(PetscBLASIntCast(ldu,&lu));
128 : PetscCall(PetscBLASIntCast(m,&n));
129 : PetscCall(PetscBLASIntCast(k,&kk));
130 : PetscCall(PetscBLASIntCast(6*m,&lwork));
131 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
132 :
133 : /* transpose W = H' */
134 : PetscCall(PetscMalloc5(m*m,&W,m*m,&Q,m,&wr,m,&wi,lwork,&work));
135 : for (j=0;j<m;j++) {
136 : for (i=0;i<m;i++) W[i+j*m] = H[j+i*ldh];
137 : }
138 :
139 : /* compute the real Schur form of W */
140 : PetscCallBLAS("LAPACKgees",LAPACKgees_("V","N",NULL,&n,W,&n,&sdim,wr,wi,Q,&n,work,&lwork,NULL,&info));
141 : SlepcCheckLapackInfo("gees",info);
142 : #if defined(PETSC_USE_DEBUG)
143 : for (i=0;i<m;i++) PetscCheck(PetscRealPart(wr[i])<0.0,PETSC_COMM_SELF,PETSC_ERR_USER_INPUT,"Eigenvalue with non-negative real part, the coefficient matrix is not stable");
144 : #endif
145 :
146 : /* copy B' into first rows of U */
147 : for (i=0;i<k;i++) {
148 : for (j=0;j<m;j++) U[i+j*ldu] = B[j+i*ldb];
149 : }
150 :
151 : /* solve Lyapunov equation (Hammarling) */
152 : PetscCallBLAS("SLICOTsb03od",SLICOTsb03od_("C","F","N",&n,&kk,W,&n,Q,&n,U,&lu,&scal,wr,wi,work,&lwork,&info));
153 : PetscCheck(!info,PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in SLICOT subroutine SB03OD: info=%" PetscBLASInt_FMT,info);
154 : PetscCheck(scal==1.0,PETSC_COMM_SELF,PETSC_ERR_SUP,"Current implementation cannot handle scale factor %g",scal);
155 :
156 : /* resnorm = norm(H(m+1,:)*U'*U), use Q(:,1) = U'*U(:,m) */
157 : if (res) {
158 : for (j=0;j<m;j++) Q[j] = U[j+(m-1)*ldu];
159 : PetscCallBLAS("BLAStrmv",BLAStrmv_("U","C","N",&n,U,&lu,Q,&ione));
160 : PetscCheck(k==1,PETSC_COMM_SELF,PETSC_ERR_LIB,"Residual error is intended for k=1 only, but you set k=%" PetscInt_FMT,k);
161 : *res *= BLASnrm2_(&n,Q,&ione);
162 : }
163 :
164 : PetscCall(PetscFPTrapPop());
165 : PetscCall(PetscFree5(W,Q,wr,wi,work));
166 : PetscFunctionReturn(PETSC_SUCCESS);
167 : }
168 :
169 : #else
170 :
171 : /*
172 : Compute the upper Cholesky factor of A
173 : */
174 56 : static PetscErrorCode CholeskyFactor(PetscInt m,PetscScalar *A,PetscInt lda)
175 : {
176 56 : PetscInt i;
177 56 : PetscScalar *S;
178 56 : PetscBLASInt info,n,ld;
179 :
180 56 : PetscFunctionBegin;
181 56 : PetscCall(PetscBLASIntCast(m,&n));
182 56 : PetscCall(PetscBLASIntCast(lda,&ld));
183 56 : PetscCall(PetscMalloc1(m*m,&S));
184 56 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
185 :
186 : /* save a copy of matrix in S */
187 1728 : for (i=0;i<m;i++) PetscCall(PetscArraycpy(S+i*m,A+i*lda,m));
188 :
189 : /* compute upper Cholesky factor in R */
190 56 : PetscCallBLAS("LAPACKpotrf",LAPACKpotrf_("U",&n,A,&ld,&info));
191 56 : PetscCall(PetscLogFlops((1.0*n*n*n)/3.0));
192 :
193 56 : if (info) {
194 1717 : for (i=0;i<m;i++) {
195 1662 : PetscCall(PetscArraycpy(A+i*lda,S+i*m,m));
196 1662 : A[i+i*lda] += 50.0*PETSC_MACHINE_EPSILON;
197 : }
198 55 : PetscCallBLAS("LAPACKpotrf",LAPACKpotrf_("U",&n,A,&ld,&info));
199 55 : SlepcCheckLapackInfo("potrf",info);
200 55 : PetscCall(PetscLogFlops((1.0*n*n*n)/3.0));
201 : }
202 :
203 : /* Zero out entries below the diagonal */
204 1672 : for (i=0;i<m-1;i++) PetscCall(PetscArrayzero(A+i*lda+i+1,m-i-1));
205 56 : PetscCall(PetscFPTrapPop());
206 56 : PetscCall(PetscFree(S));
207 56 : PetscFunctionReturn(PETSC_SUCCESS);
208 : }
209 :
210 : /*
211 : HessLyapunovChol_LAPACK - alternative implementation when SLICOT is not available
212 : */
213 56 : static PetscErrorCode HessLyapunovChol_LAPACK(PetscInt m,PetscScalar *H,PetscInt ldh,PetscInt k,PetscScalar *B,PetscInt ldb,PetscScalar *U,PetscInt ldu,PetscReal *res)
214 : {
215 56 : PetscBLASInt ilo=1,lwork,info,n,kk,lu,lb,ione=1;
216 56 : PetscInt i,j;
217 56 : PetscReal scal;
218 56 : PetscScalar *Q,*C,*W,*Z,*wr,*work,zero=0.0,done=1.0,dmone=-1.0;
219 : #if !defined(PETSC_USE_COMPLEX)
220 56 : PetscScalar *wi;
221 : #endif
222 :
223 56 : PetscFunctionBegin;
224 56 : PetscCall(PetscBLASIntCast(ldb,&lb));
225 56 : PetscCall(PetscBLASIntCast(ldu,&lu));
226 56 : PetscCall(PetscBLASIntCast(m,&n));
227 56 : PetscCall(PetscBLASIntCast(k,&kk));
228 56 : PetscCall(PetscBLASIntCast(6*m,&lwork));
229 56 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
230 56 : C = U;
231 :
232 : #if !defined(PETSC_USE_COMPLEX)
233 56 : PetscCall(PetscMalloc6(m*m,&Q,m*m,&W,m*k,&Z,m,&wr,m,&wi,lwork,&work));
234 : #else
235 : PetscCall(PetscMalloc5(m*m,&Q,m*m,&W,m*k,&Z,m,&wr,lwork,&work));
236 : #endif
237 :
238 : /* save a copy W = H */
239 1728 : for (j=0;j<m;j++) {
240 52892 : for (i=0;i<m;i++) W[i+j*m] = H[i+j*ldh];
241 : }
242 :
243 : /* compute the (real) Schur form of W */
244 : #if !defined(PETSC_USE_COMPLEX)
245 56 : PetscCallBLAS("LAPACKhseqr",LAPACKhseqr_("S","I",&n,&ilo,&n,W,&n,wr,wi,Q,&n,work,&lwork,&info));
246 : #else
247 : PetscCallBLAS("LAPACKhseqr",LAPACKhseqr_("S","I",&n,&ilo,&n,W,&n,wr,Q,&n,work,&lwork,&info));
248 : #endif
249 56 : SlepcCheckLapackInfo("hseqr",info);
250 : #if defined(PETSC_USE_DEBUG)
251 1728 : for (i=0;i<m;i++) PetscCheck(PetscRealPart(wr[i])<0.0,PETSC_COMM_SELF,PETSC_ERR_USER_INPUT,"Eigenvalue with non-negative real part %g, the coefficient matrix is not stable",(double)PetscRealPart(wr[i]));
252 : #endif
253 :
254 : /* C = -Z*Z', Z = Q'*B */
255 56 : PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&n,&kk,&n,&done,Q,&n,B,&lb,&zero,Z,&n));
256 56 : PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&n,&n,&kk,&dmone,Z,&n,Z,&n,&zero,C,&lu));
257 :
258 : /* solve triangular Sylvester equation */
259 56 : PetscCallBLAS("LAPACKtrsyl",LAPACKtrsyl_("N","C",&ione,&n,&n,W,&n,W,&n,C,&lu,&scal,&info));
260 56 : SlepcCheckLapackInfo("trsyl",info);
261 56 : PetscCheck(scal==1.0,PETSC_COMM_SELF,PETSC_ERR_SUP,"Current implementation cannot handle scale factor %g",(double)scal);
262 :
263 : /* back-transform C = Q*C*Q' */
264 56 : PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&done,Q,&n,C,&n,&zero,W,&n));
265 56 : PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&n,&n,&n,&done,W,&n,Q,&n,&zero,C,&lu));
266 :
267 : /* resnorm = norm(H(m+1,:)*Y) */
268 56 : if (res) {
269 55 : PetscCheck(k==1,PETSC_COMM_SELF,PETSC_ERR_LIB,"Residual error is intended for k=1 only, but you set k=%" PetscInt_FMT,k);
270 55 : *res *= BLASnrm2_(&n,C+m-1,&n);
271 : }
272 :
273 : /* U = chol(C) */
274 56 : PetscCall(CholeskyFactor(m,C,ldu));
275 :
276 56 : PetscCall(PetscFPTrapPop());
277 : #if !defined(PETSC_USE_COMPLEX)
278 56 : PetscCall(PetscFree6(Q,W,Z,wr,wi,work));
279 : #else
280 : PetscCall(PetscFree5(Q,W,Z,wr,work));
281 : #endif
282 56 : PetscFunctionReturn(PETSC_SUCCESS);
283 : }
284 :
285 : #endif /* SLEPC_HAVE_SLICOT */
286 :
287 : /*@C
288 : LMEDenseHessLyapunovChol - Computes the Cholesky factor of the solution of a
289 : dense Lyapunov equation with an upper Hessenberg coefficient matrix.
290 :
291 : Logically Collective
292 :
293 : Input Parameters:
294 : + lme - linear matrix equation solver context
295 : . m - number of rows and columns of H
296 : . H - coefficient matrix
297 : . ldh - leading dimension of H
298 : . k - number of columns of B
299 : . B - right-hand side matrix
300 : . ldb - leading dimension of B
301 : - ldu - leading dimension of U
302 :
303 : Output Parameters:
304 : + U - Cholesky factor of the solution
305 : - res - (optional) residual norm, on input it should contain H(m+1,m)
306 :
307 : Note:
308 : The Lyapunov equation has the form H*X + X*H' = -B*B', where H is an mxm
309 : upper Hessenberg matrix, B is an mxk matrix and the solution is expressed
310 : as X = U'*U, where U is upper triangular. H is supposed to be stable.
311 :
312 : When k=1 and the res argument is provided, the last row of X is used to
313 : compute the residual norm of a Lyapunov equation projected via Arnoldi.
314 :
315 : Level: developer
316 :
317 : .seealso: LMEDenseLyapunov(), LMESolve()
318 : @*/
319 56 : PetscErrorCode LMEDenseHessLyapunovChol(LME lme,PetscInt m,PetscScalar *H,PetscInt ldh,PetscInt k,PetscScalar *B,PetscInt ldb,PetscScalar *U,PetscInt ldu,PetscReal *res)
320 : {
321 : #if defined(PETSC_USE_INFO)
322 56 : PetscReal error;
323 : #endif
324 :
325 56 : PetscFunctionBegin;
326 56 : PetscValidHeaderSpecific(lme,LME_CLASSID,1);
327 224 : PetscValidLogicalCollectiveInt(lme,m,2);
328 56 : PetscAssertPointer(H,3);
329 224 : PetscValidLogicalCollectiveInt(lme,ldh,4);
330 224 : PetscValidLogicalCollectiveInt(lme,k,5);
331 56 : PetscAssertPointer(B,6);
332 224 : PetscValidLogicalCollectiveInt(lme,ldb,7);
333 56 : PetscAssertPointer(U,8);
334 224 : PetscValidLogicalCollectiveInt(lme,ldu,9);
335 221 : if (res) PetscValidLogicalCollectiveReal(lme,*res,10);
336 :
337 : #if defined(SLEPC_HAVE_SLICOT)
338 : PetscCall(HessLyapunovChol_SLICOT(m,H,ldh,k,B,ldb,U,ldu,res));
339 : #else
340 56 : PetscCall(HessLyapunovChol_LAPACK(m,H,ldh,k,B,ldb,U,ldu,res));
341 : #endif
342 :
343 : #if defined(PETSC_USE_INFO)
344 56 : if (PetscLogPrintInfo) {
345 1 : PetscCall(LyapunovCholResidual(m,H,ldh,k,B,ldb,U,ldu,&error));
346 1 : PetscCall(PetscInfo(lme,"Residual norm of dense Lyapunov equation = %g\n",(double)error));
347 : }
348 : #endif
349 56 : PetscFunctionReturn(PETSC_SUCCESS);
350 : }
351 :
352 : #if defined(SLEPC_HAVE_SLICOT)
353 : /*
354 : Lyapunov_SLICOT - implementation used when SLICOT is available
355 : */
356 : static PetscErrorCode Lyapunov_SLICOT(PetscInt m,PetscScalar *H,PetscInt ldh,PetscScalar *B,PetscInt ldb,PetscScalar *X,PetscInt ldx)
357 : {
358 : PetscBLASInt sdim,lwork,info,n,lx,*iwork;
359 : PetscInt i,j;
360 : PetscReal scal,sep,ferr,*work;
361 : PetscScalar *Q,*W,*wr,*wi;
362 :
363 : PetscFunctionBegin;
364 : PetscCall(PetscBLASIntCast(ldx,&lx));
365 : PetscCall(PetscBLASIntCast(m,&n));
366 : PetscCall(PetscBLASIntCast(PetscMax(20,m*m),&lwork));
367 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
368 :
369 : /* transpose W = H' */
370 : PetscCall(PetscMalloc6(m*m,&W,m*m,&Q,m,&wr,m,&wi,m*m,&iwork,lwork,&work));
371 : for (j=0;j<m;j++) {
372 : for (i=0;i<m;i++) W[i+j*m] = H[j+i*ldh];
373 : }
374 :
375 : /* compute the real Schur form of W */
376 : PetscCallBLAS("LAPACKgees",LAPACKgees_("V","N",NULL,&n,W,&n,&sdim,wr,wi,Q,&n,work,&lwork,NULL,&info));
377 : SlepcCheckLapackInfo("gees",info);
378 :
379 : /* copy -B into X */
380 : for (i=0;i<m;i++) {
381 : for (j=0;j<m;j++) X[i+j*ldx] = -B[i+j*ldb];
382 : }
383 :
384 : /* solve Lyapunov equation (Hammarling) */
385 : PetscCallBLAS("SLICOTsb03md",SLICOTsb03md_("C","X","F","N",&n,W,&n,Q,&n,X,&lx,&scal,&sep,&ferr,wr,wi,iwork,work,&lwork,&info));
386 : PetscCheck(!info,PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in SLICOT subroutine SB03OD: info=%" PetscBLASInt_FMT,info);
387 : PetscCheck(scal==1.0,PETSC_COMM_SELF,PETSC_ERR_SUP,"Current implementation cannot handle scale factor %g",scal);
388 :
389 : PetscCall(PetscFPTrapPop());
390 : PetscCall(PetscFree6(W,Q,wr,wi,iwork,work));
391 : PetscFunctionReturn(PETSC_SUCCESS);
392 : }
393 :
394 : #else
395 :
396 : /*
397 : Lyapunov_LAPACK - alternative implementation when SLICOT is not available
398 : */
399 620 : static PetscErrorCode Lyapunov_LAPACK(PetscInt m,PetscScalar *A,PetscInt lda,PetscScalar *B,PetscInt ldb,PetscScalar *X,PetscInt ldx)
400 : {
401 620 : PetscBLASInt sdim,lwork,info,n,lx,lb,ione=1;
402 620 : PetscInt i,j;
403 620 : PetscReal scal;
404 620 : PetscScalar *Q,*W,*Z,*wr,*work,zero=0.0,done=1.0,dmone=-1.0;
405 : #if defined(PETSC_USE_COMPLEX)
406 : PetscReal *rwork;
407 : #else
408 620 : PetscScalar *wi;
409 : #endif
410 :
411 620 : PetscFunctionBegin;
412 620 : PetscCall(PetscBLASIntCast(ldb,&lb));
413 620 : PetscCall(PetscBLASIntCast(ldx,&lx));
414 620 : PetscCall(PetscBLASIntCast(m,&n));
415 620 : PetscCall(PetscBLASIntCast(6*m,&lwork));
416 620 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
417 :
418 : #if !defined(PETSC_USE_COMPLEX)
419 620 : PetscCall(PetscMalloc6(m*m,&Q,m*m,&W,m*m,&Z,m,&wr,m,&wi,lwork,&work));
420 : #else
421 : PetscCall(PetscMalloc6(m*m,&Q,m*m,&W,m*m,&Z,m,&wr,lwork,&work,m,&rwork));
422 : #endif
423 :
424 : /* save a copy W = A */
425 6166 : for (j=0;j<m;j++) {
426 55822 : for (i=0;i<m;i++) W[i+j*m] = A[i+j*lda];
427 : }
428 :
429 : /* compute the (real) Schur form of W */
430 : #if !defined(PETSC_USE_COMPLEX)
431 620 : PetscCallBLAS("LAPACKgees",LAPACKgees_("V","N",NULL,&n,W,&n,&sdim,wr,wi,Q,&n,work,&lwork,NULL,&info));
432 : #else
433 : PetscCallBLAS("LAPACKgees",LAPACKgees_("V","N",NULL,&n,W,&n,&sdim,wr,Q,&n,work,&lwork,rwork,NULL,&info));
434 : #endif
435 620 : SlepcCheckLapackInfo("gees",info);
436 :
437 : /* X = -Q'*B*Q */
438 620 : PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&n,&n,&n,&done,Q,&n,B,&lb,&zero,Z,&n));
439 620 : PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&dmone,Z,&n,Q,&n,&zero,X,&lx));
440 :
441 : /* solve triangular Sylvester equation */
442 620 : PetscCallBLAS("LAPACKtrsyl",LAPACKtrsyl_("N","C",&ione,&n,&n,W,&n,W,&n,X,&lx,&scal,&info));
443 620 : SlepcCheckLapackInfo("trsyl",info);
444 620 : PetscCheck(scal==1.0,PETSC_COMM_SELF,PETSC_ERR_SUP,"Current implementation cannot handle scale factor %g",(double)scal);
445 :
446 : /* back-transform X = Q*X*Q' */
447 620 : PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&done,Q,&n,X,&n,&zero,W,&n));
448 620 : PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&n,&n,&n,&done,W,&n,Q,&n,&zero,X,&lx));
449 :
450 620 : PetscCall(PetscFPTrapPop());
451 : #if !defined(PETSC_USE_COMPLEX)
452 620 : PetscCall(PetscFree6(Q,W,Z,wr,wi,work));
453 : #else
454 : PetscCall(PetscFree6(Q,W,Z,wr,work,rwork));
455 : #endif
456 620 : PetscFunctionReturn(PETSC_SUCCESS);
457 : }
458 :
459 : #endif /* SLEPC_HAVE_SLICOT */
460 :
461 : /*@C
462 : LMEDenseLyapunov - Computes the solution of a dense continuous-time Lyapunov
463 : equation.
464 :
465 : Logically Collective
466 :
467 : Input Parameters:
468 : + lme - linear matrix equation solver context
469 : . m - number of rows and columns of A
470 : . A - coefficient matrix
471 : . lda - leading dimension of A
472 : . B - right-hand side matrix
473 : . ldb - leading dimension of B
474 : - ldx - leading dimension of X
475 :
476 : Output Parameter:
477 : . X - the solution
478 :
479 : Note:
480 : The Lyapunov equation has the form A*X + X*A' = -B, where all are mxm
481 : matrices, and B is symmetric.
482 :
483 : Level: developer
484 :
485 : .seealso: LMEDenseHessLyapunovChol(), LMESolve()
486 : @*/
487 620 : PetscErrorCode LMEDenseLyapunov(LME lme,PetscInt m,PetscScalar *A,PetscInt lda,PetscScalar *B,PetscInt ldb,PetscScalar *X,PetscInt ldx)
488 : {
489 : #if defined(PETSC_USE_INFO)
490 620 : PetscReal error;
491 : #endif
492 :
493 620 : PetscFunctionBegin;
494 620 : PetscValidHeaderSpecific(lme,LME_CLASSID,1);
495 2480 : PetscValidLogicalCollectiveInt(lme,m,2);
496 620 : PetscAssertPointer(A,3);
497 2480 : PetscValidLogicalCollectiveInt(lme,lda,4);
498 620 : PetscAssertPointer(B,5);
499 2480 : PetscValidLogicalCollectiveInt(lme,ldb,6);
500 620 : PetscAssertPointer(X,7);
501 2480 : PetscValidLogicalCollectiveInt(lme,ldx,8);
502 :
503 : #if defined(SLEPC_HAVE_SLICOT)
504 : PetscCall(Lyapunov_SLICOT(m,A,lda,B,ldb,X,ldx));
505 : #else
506 620 : PetscCall(Lyapunov_LAPACK(m,A,lda,B,ldb,X,ldx));
507 : #endif
508 :
509 : #if defined(PETSC_USE_INFO)
510 620 : if (PetscLogPrintInfo) {
511 1 : PetscCall(LyapunovResidual(m,A,lda,B,ldb,X,ldx,&error));
512 1 : PetscCall(PetscInfo(lme,"Residual norm of dense Lyapunov equation = %g\n",(double)error));
513 : }
514 : #endif
515 620 : PetscFunctionReturn(PETSC_SUCCESS);
516 : }
|
