Actual source code: lmedense.c
slepc-3.21.0 2024-03-30
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: Routines for solving dense matrix equations, in some cases calling SLICOT
12: */
14: #include <slepc/private/lmeimpl.h>
15: #include <slepcblaslapack.h>
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: PetscErrorCode LMEDenseRankSVD(LME lme,PetscInt n,PetscScalar *A,PetscInt lda,PetscScalar *U,PetscInt ldu,PetscInt *rank)
22: {
23: PetscInt i,j,rk=0;
24: PetscScalar *work;
25: PetscReal tol,*sg,*rwork;
26: PetscBLASInt n_,lda_,ldu_,info,lw_;
28: PetscFunctionBegin;
29: PetscCall(PetscCalloc3(n,&sg,10*n,&work,5*n,&rwork));
30: PetscCall(PetscBLASIntCast(n,&n_));
31: PetscCall(PetscBLASIntCast(lda,&lda_));
32: PetscCall(PetscBLASIntCast(ldu,&ldu_));
33: lw_ = 10*n_;
34: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
35: #if !defined (PETSC_USE_COMPLEX)
36: 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: PetscCall(PetscFPTrapPop());
41: SlepcCheckLapackInfo("gesvd",info);
42: tol = 10*PETSC_MACHINE_EPSILON*n*sg[0];
43: for (j=0;j<n;j++) {
44: if (sg[j]>tol) {
45: for (i=0;i<n;i++) U[i+j*n] *= sg[j];
46: rk++;
47: } else break;
48: }
49: *rank = rk;
50: PetscCall(PetscFree3(sg,work,rwork));
51: PetscFunctionReturn(PETSC_SUCCESS);
52: }
54: #if defined(PETSC_USE_INFO)
55: /*
56: LyapunovCholResidual - compute the residual norm ||A*U'*U+U'*U*A'+B*B'||
57: */
58: static PetscErrorCode LyapunovCholResidual(PetscInt m,PetscScalar *A,PetscInt lda,PetscInt k,PetscScalar *B,PetscInt ldb,PetscScalar *U,PetscInt ldu,PetscReal *res)
59: {
60: PetscBLASInt n,kk,la,lb,lu;
61: PetscScalar *M,*R,zero=0.0,done=1.0;
63: PetscFunctionBegin;
64: *res = 0;
65: PetscCall(PetscBLASIntCast(lda,&la));
66: PetscCall(PetscBLASIntCast(ldb,&lb));
67: PetscCall(PetscBLASIntCast(ldu,&lu));
68: PetscCall(PetscBLASIntCast(m,&n));
69: PetscCall(PetscBLASIntCast(k,&kk));
70: PetscCall(PetscMalloc2(m*m,&M,m*m,&R));
72: /* R = B*B' */
73: PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&n,&n,&kk,&done,B,&lb,B,&lb,&zero,R,&n));
74: /* M = A*U' */
75: PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&n,&n,&n,&done,A,&la,U,&lu,&zero,M,&n));
76: /* R = R+M*U */
77: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&done,M,&n,U,&lu,&done,R,&n));
78: /* R = R+U'*M' */
79: PetscCallBLAS("BLASgemm",BLASgemm_("C","C",&n,&n,&n,&done,U,&lu,M,&n,&done,R,&n));
81: *res = LAPACKlange_("F",&n,&n,R,&n,NULL);
82: PetscCall(PetscFree2(M,R));
83: PetscFunctionReturn(PETSC_SUCCESS);
84: }
86: /*
87: LyapunovResidual - compute the residual norm ||A*X+X*A'+B||
88: */
89: static PetscErrorCode LyapunovResidual(PetscInt m,PetscScalar *A,PetscInt lda,PetscScalar *B,PetscInt ldb,PetscScalar *X,PetscInt ldx,PetscReal *res)
90: {
91: PetscInt i;
92: PetscBLASInt n,la,lb,lx;
93: PetscScalar *R,done=1.0;
95: PetscFunctionBegin;
96: *res = 0;
97: PetscCall(PetscBLASIntCast(lda,&la));
98: PetscCall(PetscBLASIntCast(ldb,&lb));
99: PetscCall(PetscBLASIntCast(ldx,&lx));
100: PetscCall(PetscBLASIntCast(m,&n));
101: PetscCall(PetscMalloc1(m*m,&R));
103: /* R = B+A*X */
104: for (i=0;i<m;i++) PetscCall(PetscArraycpy(R+i*m,B+i*ldb,m));
105: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&done,A,&la,X,&lx,&done,R,&n));
106: /* R = R+X*A' */
107: PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&n,&n,&n,&done,X,&lx,A,&la,&done,R,&n));
109: *res = LAPACKlange_("F",&n,&n,R,&n,NULL);
110: PetscCall(PetscFree(R));
111: PetscFunctionReturn(PETSC_SUCCESS);
112: }
113: #endif
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;
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));
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: }
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
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: }
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);
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: }
164: PetscCall(PetscFPTrapPop());
165: PetscCall(PetscFree5(W,Q,wr,wi,work));
166: PetscFunctionReturn(PETSC_SUCCESS);
167: }
169: #else
171: /*
172: Compute the upper Cholesky factor of A
173: */
174: static PetscErrorCode CholeskyFactor(PetscInt m,PetscScalar *A,PetscInt lda)
175: {
176: PetscInt i;
177: PetscScalar *S;
178: PetscBLASInt info,n,ld;
180: PetscFunctionBegin;
181: PetscCall(PetscBLASIntCast(m,&n));
182: PetscCall(PetscBLASIntCast(lda,&ld));
183: PetscCall(PetscMalloc1(m*m,&S));
184: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
186: /* save a copy of matrix in S */
187: for (i=0;i<m;i++) PetscCall(PetscArraycpy(S+i*m,A+i*lda,m));
189: /* compute upper Cholesky factor in R */
190: PetscCallBLAS("LAPACKpotrf",LAPACKpotrf_("U",&n,A,&ld,&info));
191: PetscCall(PetscLogFlops((1.0*n*n*n)/3.0));
193: if (info) {
194: for (i=0;i<m;i++) {
195: PetscCall(PetscArraycpy(A+i*lda,S+i*m,m));
196: A[i+i*lda] += 50.0*PETSC_MACHINE_EPSILON;
197: }
198: PetscCallBLAS("LAPACKpotrf",LAPACKpotrf_("U",&n,A,&ld,&info));
199: SlepcCheckLapackInfo("potrf",info);
200: PetscCall(PetscLogFlops((1.0*n*n*n)/3.0));
201: }
203: /* Zero out entries below the diagonal */
204: for (i=0;i<m-1;i++) PetscCall(PetscArrayzero(A+i*lda+i+1,m-i-1));
205: PetscCall(PetscFPTrapPop());
206: PetscCall(PetscFree(S));
207: PetscFunctionReturn(PETSC_SUCCESS);
208: }
210: /*
211: HessLyapunovChol_LAPACK - alternative implementation when SLICOT is not available
212: */
213: static PetscErrorCode HessLyapunovChol_LAPACK(PetscInt m,PetscScalar *H,PetscInt ldh,PetscInt k,PetscScalar *B,PetscInt ldb,PetscScalar *U,PetscInt ldu,PetscReal *res)
214: {
215: PetscBLASInt ilo=1,lwork,info,n,kk,lu,lb,ione=1;
216: PetscInt i,j;
217: PetscReal scal;
218: PetscScalar *Q,*C,*W,*Z,*wr,*work,zero=0.0,done=1.0,dmone=-1.0;
219: #if !defined(PETSC_USE_COMPLEX)
220: PetscScalar *wi;
221: #endif
223: PetscFunctionBegin;
224: PetscCall(PetscBLASIntCast(ldb,&lb));
225: PetscCall(PetscBLASIntCast(ldu,&lu));
226: PetscCall(PetscBLASIntCast(m,&n));
227: PetscCall(PetscBLASIntCast(k,&kk));
228: PetscCall(PetscBLASIntCast(6*m,&lwork));
229: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
230: C = U;
232: #if !defined(PETSC_USE_COMPLEX)
233: 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
238: /* save a copy W = H */
239: for (j=0;j<m;j++) {
240: for (i=0;i<m;i++) W[i+j*m] = H[i+j*ldh];
241: }
243: /* compute the (real) Schur form of W */
244: #if !defined(PETSC_USE_COMPLEX)
245: 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: SlepcCheckLapackInfo("hseqr",info);
250: #if defined(PETSC_USE_DEBUG)
251: 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
254: /* C = -Z*Z', Z = Q'*B */
255: PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&n,&kk,&n,&done,Q,&n,B,&lb,&zero,Z,&n));
256: PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&n,&n,&kk,&dmone,Z,&n,Z,&n,&zero,C,&lu));
258: /* solve triangular Sylvester equation */
259: PetscCallBLAS("LAPACKtrsyl",LAPACKtrsyl_("N","C",&ione,&n,&n,W,&n,W,&n,C,&lu,&scal,&info));
260: SlepcCheckLapackInfo("trsyl",info);
261: PetscCheck(scal==1.0,PETSC_COMM_SELF,PETSC_ERR_SUP,"Current implementation cannot handle scale factor %g",(double)scal);
263: /* back-transform C = Q*C*Q' */
264: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&done,Q,&n,C,&n,&zero,W,&n));
265: PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&n,&n,&n,&done,W,&n,Q,&n,&zero,C,&lu));
267: /* resnorm = norm(H(m+1,:)*Y) */
268: if (res) {
269: 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: *res *= BLASnrm2_(&n,C+m-1,&n);
271: }
273: /* U = chol(C) */
274: PetscCall(CholeskyFactor(m,C,ldu));
276: PetscCall(PetscFPTrapPop());
277: #if !defined(PETSC_USE_COMPLEX)
278: PetscCall(PetscFree6(Q,W,Z,wr,wi,work));
279: #else
280: PetscCall(PetscFree5(Q,W,Z,wr,work));
281: #endif
282: PetscFunctionReturn(PETSC_SUCCESS);
283: }
285: #endif /* SLEPC_HAVE_SLICOT */
287: /*@C
288: LMEDenseHessLyapunovChol - Computes the Cholesky factor of the solution of a
289: dense Lyapunov equation with an upper Hessenberg coefficient matrix.
291: Logically Collective
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
303: Output Parameters:
304: + U - Cholesky factor of the solution
305: - res - (optional) residual norm, on input it should contain H(m+1,m)
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.
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.
315: Level: developer
317: .seealso: LMEDenseLyapunov(), LMESolve()
318: @*/
319: 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: PetscReal error;
323: #endif
325: PetscFunctionBegin;
328: PetscAssertPointer(H,3);
331: PetscAssertPointer(B,6);
333: PetscAssertPointer(U,8);
337: #if defined(SLEPC_HAVE_SLICOT)
338: PetscCall(HessLyapunovChol_SLICOT(m,H,ldh,k,B,ldb,U,ldu,res));
339: #else
340: PetscCall(HessLyapunovChol_LAPACK(m,H,ldh,k,B,ldb,U,ldu,res));
341: #endif
343: #if defined(PETSC_USE_INFO)
344: if (PetscLogPrintInfo) {
345: PetscCall(LyapunovCholResidual(m,H,ldh,k,B,ldb,U,ldu,&error));
346: PetscCall(PetscInfo(lme,"Residual norm of dense Lyapunov equation = %g\n",(double)error));
347: }
348: #endif
349: PetscFunctionReturn(PETSC_SUCCESS);
350: }
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;
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));
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: }
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);
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: }
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);
389: PetscCall(PetscFPTrapPop());
390: PetscCall(PetscFree6(W,Q,wr,wi,iwork,work));
391: PetscFunctionReturn(PETSC_SUCCESS);
392: }
394: #else
396: /*
397: Lyapunov_LAPACK - alternative implementation when SLICOT is not available
398: */
399: static PetscErrorCode Lyapunov_LAPACK(PetscInt m,PetscScalar *A,PetscInt lda,PetscScalar *B,PetscInt ldb,PetscScalar *X,PetscInt ldx)
400: {
401: PetscBLASInt sdim,lwork,info,n,lx,lb,ione=1;
402: PetscInt i,j;
403: PetscReal scal;
404: 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: PetscScalar *wi;
409: #endif
411: PetscFunctionBegin;
412: PetscCall(PetscBLASIntCast(ldb,&lb));
413: PetscCall(PetscBLASIntCast(ldx,&lx));
414: PetscCall(PetscBLASIntCast(m,&n));
415: PetscCall(PetscBLASIntCast(6*m,&lwork));
416: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
418: #if !defined(PETSC_USE_COMPLEX)
419: 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
424: /* save a copy W = A */
425: for (j=0;j<m;j++) {
426: for (i=0;i<m;i++) W[i+j*m] = A[i+j*lda];
427: }
429: /* compute the (real) Schur form of W */
430: #if !defined(PETSC_USE_COMPLEX)
431: 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: SlepcCheckLapackInfo("gees",info);
437: /* X = -Q'*B*Q */
438: PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&n,&n,&n,&done,Q,&n,B,&lb,&zero,Z,&n));
439: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&dmone,Z,&n,Q,&n,&zero,X,&lx));
441: /* solve triangular Sylvester equation */
442: PetscCallBLAS("LAPACKtrsyl",LAPACKtrsyl_("N","C",&ione,&n,&n,W,&n,W,&n,X,&lx,&scal,&info));
443: SlepcCheckLapackInfo("trsyl",info);
444: PetscCheck(scal==1.0,PETSC_COMM_SELF,PETSC_ERR_SUP,"Current implementation cannot handle scale factor %g",(double)scal);
446: /* back-transform X = Q*X*Q' */
447: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&done,Q,&n,X,&n,&zero,W,&n));
448: PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&n,&n,&n,&done,W,&n,Q,&n,&zero,X,&lx));
450: PetscCall(PetscFPTrapPop());
451: #if !defined(PETSC_USE_COMPLEX)
452: PetscCall(PetscFree6(Q,W,Z,wr,wi,work));
453: #else
454: PetscCall(PetscFree6(Q,W,Z,wr,work,rwork));
455: #endif
456: PetscFunctionReturn(PETSC_SUCCESS);
457: }
459: #endif /* SLEPC_HAVE_SLICOT */
461: /*@C
462: LMEDenseLyapunov - Computes the solution of a dense continuous-time Lyapunov
463: equation.
465: Logically Collective
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
476: Output Parameter:
477: . X - the solution
479: Note:
480: The Lyapunov equation has the form A*X + X*A' = -B, where all are mxm
481: matrices, and B is symmetric.
483: Level: developer
485: .seealso: LMEDenseHessLyapunovChol(), LMESolve()
486: @*/
487: 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: PetscReal error;
491: #endif
493: PetscFunctionBegin;
496: PetscAssertPointer(A,3);
498: PetscAssertPointer(B,5);
500: PetscAssertPointer(X,7);
503: #if defined(SLEPC_HAVE_SLICOT)
504: PetscCall(Lyapunov_SLICOT(m,A,lda,B,ldb,X,ldx));
505: #else
506: PetscCall(Lyapunov_LAPACK(m,A,lda,B,ldb,X,ldx));
507: #endif
509: #if defined(PETSC_USE_INFO)
510: if (PetscLogPrintInfo) {
511: PetscCall(LyapunovResidual(m,A,lda,B,ldb,X,ldx,&error));
512: PetscCall(PetscInfo(lme,"Residual norm of dense Lyapunov equation = %g\n",(double)error));
513: }
514: #endif
515: PetscFunctionReturn(PETSC_SUCCESS);
516: }