Actual source code: lmedense.c
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_,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 !PetscDefined(USE_COMPLEX)
36: PetscCallLAPACKInfo("LAPACKgesvd",LAPACKgesvd_("S","O",&n_,&n_,A,&lda_,sg,U,&ldu_,NULL,&n_,work,&lw_,&info));
37: #else
38: PetscCallLAPACKInfo("LAPACKgesvd",LAPACKgesvd_("S","O",&n_,&n_,A,&lda_,sg,U,&ldu_,NULL,&n_,work,&lw_,rwork,&info));
39: #endif
40: PetscCall(PetscFPTrapPop());
41: tol = 10*PETSC_MACHINE_EPSILON*n*sg[0];
42: for (j=0;j<n;j++) {
43: if (sg[j]>tol) {
44: for (i=0;i<n;i++) U[i+j*n] *= sg[j];
45: rk++;
46: } else break;
47: }
48: *rank = rk;
49: PetscCall(PetscFree3(sg,work,rwork));
50: PetscFunctionReturn(PETSC_SUCCESS);
51: }
53: #if PetscDefined(USE_INFO)
54: /*
55: LyapunovCholResidual - compute the residual norm ||A*U'*U+U'*U*A'+B*B'||
56: */
57: static PetscErrorCode LyapunovCholResidual(PetscInt m,PetscScalar *A,PetscInt lda,PetscInt k,PetscScalar *B,PetscInt ldb,PetscScalar *U,PetscInt ldu,PetscReal *res)
58: {
59: PetscBLASInt n,kk,la,lb,lu;
60: PetscScalar *M,*R,zero=0.0,done=1.0;
62: PetscFunctionBegin;
63: *res = 0;
64: PetscCall(PetscBLASIntCast(lda,&la));
65: PetscCall(PetscBLASIntCast(ldb,&lb));
66: PetscCall(PetscBLASIntCast(ldu,&lu));
67: PetscCall(PetscBLASIntCast(m,&n));
68: PetscCall(PetscBLASIntCast(k,&kk));
69: PetscCall(PetscMalloc2(m*m,&M,m*m,&R));
71: /* R = B*B' */
72: PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&n,&n,&kk,&done,B,&lb,B,&lb,&zero,R,&n));
73: /* M = A*U' */
74: PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&n,&n,&n,&done,A,&la,U,&lu,&zero,M,&n));
75: /* R = R+M*U */
76: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&done,M,&n,U,&lu,&done,R,&n));
77: /* R = R+U'*M' */
78: PetscCallBLAS("BLASgemm",BLASgemm_("C","C",&n,&n,&n,&done,U,&lu,M,&n,&done,R,&n));
80: *res = LAPACKlange_("F",&n,&n,R,&n,NULL);
81: PetscCall(PetscFree2(M,R));
82: PetscFunctionReturn(PETSC_SUCCESS);
83: }
85: /*
86: LyapunovResidual - compute the residual norm ||A*X+X*A'+B||
87: */
88: static PetscErrorCode LyapunovResidual(PetscInt m,PetscScalar *A,PetscInt lda,PetscScalar *B,PetscInt ldb,PetscScalar *X,PetscInt ldx,PetscReal *res)
89: {
90: PetscInt i;
91: PetscBLASInt n,la,lb,lx;
92: PetscScalar *R,done=1.0;
94: PetscFunctionBegin;
95: *res = 0;
96: PetscCall(PetscBLASIntCast(lda,&la));
97: PetscCall(PetscBLASIntCast(ldb,&lb));
98: PetscCall(PetscBLASIntCast(ldx,&lx));
99: PetscCall(PetscBLASIntCast(m,&n));
100: PetscCall(PetscMalloc1(m*m,&R));
102: /* R = B+A*X */
103: for (i=0;i<m;i++) PetscCall(PetscArraycpy(R+i*m,B+i*ldb,m));
104: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&done,A,&la,X,&lx,&done,R,&n));
105: /* R = R+X*A' */
106: PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&n,&n,&n,&done,X,&lx,A,&la,&done,R,&n));
108: *res = LAPACKlange_("F",&n,&n,R,&n,NULL);
109: PetscCall(PetscFree(R));
110: PetscFunctionReturn(PETSC_SUCCESS);
111: }
112: #endif
114: #if defined(SLEPC_HAVE_SLICOT)
115: /*
116: HessLyapunovChol_SLICOT - implementation used when SLICOT is available
117: */
118: static PetscErrorCode HessLyapunovChol_SLICOT(PetscInt m,PetscScalar *H,PetscInt ldh,PetscInt k,PetscScalar *B,PetscInt ldb,PetscScalar *U,PetscInt ldu,PetscReal *res)
119: {
120: PetscBLASInt lwork,info,n,kk,lu,ione=1,sdim;
121: PetscInt i,j;
122: PetscReal scal;
123: PetscScalar *Q,*W,*wr,*wi,*work;
125: PetscFunctionBegin;
126: PetscCall(PetscBLASIntCast(ldu,&lu));
127: PetscCall(PetscBLASIntCast(m,&n));
128: PetscCall(PetscBLASIntCast(k,&kk));
129: PetscCall(PetscBLASIntCast(6*m,&lwork));
130: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
132: /* transpose W = H' */
133: PetscCall(PetscMalloc5(m*m,&W,m*m,&Q,m,&wr,m,&wi,lwork,&work));
134: for (j=0;j<m;j++) {
135: for (i=0;i<m;i++) W[i+j*m] = H[j+i*ldh];
136: }
138: /* compute the real Schur form of W */
139: PetscCallLAPACKInfo("LAPACKgees",LAPACKgees_("V","N",NULL,&n,W,&n,&sdim,wr,wi,Q,&n,work,&lwork,NULL,&info));
140: #if PetscDefined(USE_DEBUG)
141: 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");
142: #endif
144: /* copy B' into first rows of U */
145: for (i=0;i<k;i++) {
146: for (j=0;j<m;j++) U[i+j*ldu] = B[j+i*ldb];
147: }
149: /* solve Lyapunov equation (Hammarling) */
150: PetscCallBLAS("SLICOTsb03od",SLICOTsb03od_("C","F","N",&n,&kk,W,&n,Q,&n,U,&lu,&scal,wr,wi,work,&lwork,&info));
151: PetscCheck(!info,PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in SLICOT subroutine SB03OD: info=%" PetscBLASInt_FMT,info);
152: PetscCheck(scal==1.0,PETSC_COMM_SELF,PETSC_ERR_SUP,"Current implementation cannot handle scale factor %g",scal);
154: /* resnorm = norm(H(m+1,:)*U'*U), use Q(:,1) = U'*U(:,m) */
155: if (res) {
156: for (j=0;j<m;j++) Q[j] = U[j+(m-1)*ldu];
157: PetscCallBLAS("BLAStrmv",BLAStrmv_("U","C","N",&n,U,&lu,Q,&ione));
158: PetscCheck(k==1,PETSC_COMM_SELF,PETSC_ERR_LIB,"Residual error is intended for k=1 only, but you set k=%" PetscInt_FMT,k);
159: *res *= BLASnrm2_(&n,Q,&ione);
160: }
162: PetscCall(PetscFPTrapPop());
163: PetscCall(PetscFree5(W,Q,wr,wi,work));
164: PetscFunctionReturn(PETSC_SUCCESS);
165: }
167: #else
169: /*
170: Compute the upper Cholesky factor of A
171: */
172: static PetscErrorCode CholeskyFactor(PetscInt m,PetscScalar *A,PetscInt lda)
173: {
174: PetscInt i;
175: PetscScalar *S;
176: PetscBLASInt info,n,ld;
178: PetscFunctionBegin;
179: PetscCall(PetscBLASIntCast(m,&n));
180: PetscCall(PetscBLASIntCast(lda,&ld));
181: PetscCall(PetscMalloc1(m*m,&S));
182: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
184: /* save a copy of matrix in S */
185: for (i=0;i<m;i++) PetscCall(PetscArraycpy(S+i*m,A+i*lda,m));
187: /* compute upper Cholesky factor in R */
188: PetscCallBLAS("LAPACKpotrf",LAPACKpotrf_("U",&n,A,&ld,&info));
189: PetscCall(PetscLogFlops((1.0*n*n*n)/3.0));
191: if (info) {
192: for (i=0;i<m;i++) {
193: PetscCall(PetscArraycpy(A+i*lda,S+i*m,m));
194: A[i+i*lda] += 50.0*PETSC_MACHINE_EPSILON;
195: }
196: PetscCallLAPACKInfo("LAPACKpotrf",LAPACKpotrf_("U",&n,A,&ld,&info));
197: PetscCall(PetscLogFlops((1.0*n*n*n)/3.0));
198: }
200: /* Zero out entries below the diagonal */
201: for (i=0;i<m-1;i++) PetscCall(PetscArrayzero(A+i*lda+i+1,m-i-1));
202: PetscCall(PetscFPTrapPop());
203: PetscCall(PetscFree(S));
204: PetscFunctionReturn(PETSC_SUCCESS);
205: }
207: /*
208: HessLyapunovChol_LAPACK - alternative implementation when SLICOT is not available
209: */
210: static PetscErrorCode HessLyapunovChol_LAPACK(PetscInt m,PetscScalar *H,PetscInt ldh,PetscInt k,PetscScalar *B,PetscInt ldb,PetscScalar *U,PetscInt ldu,PetscReal *res)
211: {
212: PetscBLASInt ilo=1,lwork,n,kk,lu,lb,ione=1;
213: PetscInt i,j;
214: PetscReal scal;
215: PetscScalar *Q,*C,*W,*Z,*wr,*work,zero=0.0,done=1.0,dmone=-1.0;
216: #if !PetscDefined(USE_COMPLEX)
217: PetscScalar *wi;
218: #endif
220: PetscFunctionBegin;
221: PetscCall(PetscBLASIntCast(ldb,&lb));
222: PetscCall(PetscBLASIntCast(ldu,&lu));
223: PetscCall(PetscBLASIntCast(m,&n));
224: PetscCall(PetscBLASIntCast(k,&kk));
225: PetscCall(PetscBLASIntCast(6*m,&lwork));
226: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
227: C = U;
229: #if !PetscDefined(USE_COMPLEX)
230: PetscCall(PetscMalloc6(m*m,&Q,m*m,&W,m*k,&Z,m,&wr,m,&wi,lwork,&work));
231: #else
232: PetscCall(PetscMalloc5(m*m,&Q,m*m,&W,m*k,&Z,m,&wr,lwork,&work));
233: #endif
235: /* save a copy W = H */
236: for (j=0;j<m;j++) {
237: for (i=0;i<m;i++) W[i+j*m] = H[i+j*ldh];
238: }
240: /* compute the (real) Schur form of W */
241: #if !PetscDefined(USE_COMPLEX)
242: PetscCallLAPACKInfo("LAPACKhseqr",LAPACKhseqr_("S","I",&n,&ilo,&n,W,&n,wr,wi,Q,&n,work,&lwork,&info));
243: #else
244: PetscCallLAPACKInfo("LAPACKhseqr",LAPACKhseqr_("S","I",&n,&ilo,&n,W,&n,wr,Q,&n,work,&lwork,&info));
245: #endif
246: #if PetscDefined(USE_DEBUG)
247: 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]));
248: #endif
250: /* C = -Z*Z', Z = Q'*B */
251: PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&n,&kk,&n,&done,Q,&n,B,&lb,&zero,Z,&n));
252: PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&n,&n,&kk,&dmone,Z,&n,Z,&n,&zero,C,&lu));
254: /* solve triangular Sylvester equation */
255: PetscCallLAPACKInfo("LAPACKtrsyl",LAPACKtrsyl_("N","C",&ione,&n,&n,W,&n,W,&n,C,&lu,&scal,&info));
256: PetscCheck(scal==1.0,PETSC_COMM_SELF,PETSC_ERR_SUP,"Current implementation cannot handle scale factor %g",(double)scal);
258: /* back-transform C = Q*C*Q' */
259: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&done,Q,&n,C,&n,&zero,W,&n));
260: PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&n,&n,&n,&done,W,&n,Q,&n,&zero,C,&lu));
262: /* resnorm = norm(H(m+1,:)*Y) */
263: if (res) {
264: PetscCheck(k==1,PETSC_COMM_SELF,PETSC_ERR_LIB,"Residual error is intended for k=1 only, but you set k=%" PetscInt_FMT,k);
265: *res *= BLASnrm2_(&n,C+m-1,&n);
266: }
268: /* U = chol(C) */
269: PetscCall(CholeskyFactor(m,C,ldu));
271: PetscCall(PetscFPTrapPop());
272: #if !PetscDefined(USE_COMPLEX)
273: PetscCall(PetscFree6(Q,W,Z,wr,wi,work));
274: #else
275: PetscCall(PetscFree5(Q,W,Z,wr,work));
276: #endif
277: PetscFunctionReturn(PETSC_SUCCESS);
278: }
280: #endif /* SLEPC_HAVE_SLICOT */
282: /*@
283: LMEDenseHessLyapunovChol - Computes the Cholesky factor of the solution of a
284: dense Lyapunov equation with an upper Hessenberg coefficient matrix.
286: Logically Collective
288: Input Parameters:
289: + lme - the linear matrix equation solver context
290: . m - number of rows and columns of `H`
291: . H - coefficient matrix
292: . ldh - leading dimension of `H`
293: . k - number of columns of `B`
294: . B - right-hand side matrix
295: . ldb - leading dimension of `B`
296: - ldu - leading dimension of `U`
298: Output Parameters:
299: + U - Cholesky factor of the solution
300: - res - (optional) residual norm, on input it should contain $h_{m+1,m}$
302: Note:
303: The Lyapunov equation has the form $HX + XH^* = -BB^*$, where $H$ is an $m\times m$
304: upper Hessenberg matrix, $B$ is an $m\times k$ matrix and the solution is expressed
305: as $X = U^*U$, where $U$ is upper triangular. $H$ is supposed to be stable.
307: When `k=1` and the `res` argument is provided, the last row of `X` is used to
308: compute the residual norm of a Lyapunov equation projected via Arnoldi.
310: Level: developer
312: .seealso: [](ch:lme), `LMEDenseLyapunov()`, `LMESolve()`
313: @*/
314: PetscErrorCode LMEDenseHessLyapunovChol(LME lme,PetscInt m,PetscScalar H[],PetscInt ldh,PetscInt k,PetscScalar B[],PetscInt ldb,PetscScalar U[],PetscInt ldu,PetscReal *res)
315: {
316: #if PetscDefined(USE_INFO)
317: PetscReal error;
318: #endif
320: PetscFunctionBegin;
323: PetscAssertPointer(H,3);
326: PetscAssertPointer(B,6);
328: PetscAssertPointer(U,8);
332: #if defined(SLEPC_HAVE_SLICOT)
333: PetscCall(HessLyapunovChol_SLICOT(m,H,ldh,k,B,ldb,U,ldu,res));
334: #else
335: PetscCall(HessLyapunovChol_LAPACK(m,H,ldh,k,B,ldb,U,ldu,res));
336: #endif
338: #if PetscDefined(USE_INFO)
339: if (PetscLogPrintInfo) {
340: PetscCall(LyapunovCholResidual(m,H,ldh,k,B,ldb,U,ldu,&error));
341: PetscCall(PetscInfo(lme,"Residual norm of dense Lyapunov equation = %g\n",(double)error));
342: }
343: #endif
344: PetscFunctionReturn(PETSC_SUCCESS);
345: }
347: #if defined(SLEPC_HAVE_SLICOT)
348: /*
349: Lyapunov_SLICOT - implementation used when SLICOT is available
350: */
351: static PetscErrorCode Lyapunov_SLICOT(PetscInt m,PetscScalar *H,PetscInt ldh,PetscScalar *B,PetscInt ldb,PetscScalar *X,PetscInt ldx)
352: {
353: PetscBLASInt sdim,lwork,info,n,lx,*iwork;
354: PetscInt i,j;
355: PetscReal scal,sep,ferr,*work;
356: PetscScalar *Q,*W,*wr,*wi;
358: PetscFunctionBegin;
359: PetscCall(PetscBLASIntCast(ldx,&lx));
360: PetscCall(PetscBLASIntCast(m,&n));
361: PetscCall(PetscBLASIntCast(PetscMax(20,m*m),&lwork));
362: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
364: /* transpose W = H' */
365: PetscCall(PetscMalloc6(m*m,&W,m*m,&Q,m,&wr,m,&wi,m*m,&iwork,lwork,&work));
366: for (j=0;j<m;j++) {
367: for (i=0;i<m;i++) W[i+j*m] = H[j+i*ldh];
368: }
370: /* compute the real Schur form of W */
371: PetscCallLAPACKInfo("LAPACKgees",LAPACKgees_("V","N",NULL,&n,W,&n,&sdim,wr,wi,Q,&n,work,&lwork,NULL,&info));
373: /* copy -B into X */
374: for (i=0;i<m;i++) {
375: for (j=0;j<m;j++) X[i+j*ldx] = -B[i+j*ldb];
376: }
378: /* solve Lyapunov equation (Hammarling) */
379: PetscCallBLAS("SLICOTsb03md",SLICOTsb03md_("C","X","F","N",&n,W,&n,Q,&n,X,&lx,&scal,&sep,&ferr,wr,wi,iwork,work,&lwork,&info));
380: PetscCheck(!info,PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in SLICOT subroutine SB03OD: info=%" PetscBLASInt_FMT,info);
381: PetscCheck(scal==1.0,PETSC_COMM_SELF,PETSC_ERR_SUP,"Current implementation cannot handle scale factor %g",scal);
383: PetscCall(PetscFPTrapPop());
384: PetscCall(PetscFree6(W,Q,wr,wi,iwork,work));
385: PetscFunctionReturn(PETSC_SUCCESS);
386: }
388: #else
390: /*
391: Lyapunov_LAPACK - alternative implementation when SLICOT is not available
392: */
393: static PetscErrorCode Lyapunov_LAPACK(PetscInt m,PetscScalar *A,PetscInt lda,PetscScalar *B,PetscInt ldb,PetscScalar *X,PetscInt ldx)
394: {
395: PetscBLASInt sdim,lwork,n,lx,lb,ione=1;
396: PetscInt i,j;
397: PetscReal scal;
398: PetscScalar *Q,*W,*Z,*wr,*work,zero=0.0,done=1.0,dmone=-1.0;
399: #if PetscDefined(USE_COMPLEX)
400: PetscReal *rwork;
401: #else
402: PetscScalar *wi;
403: #endif
405: PetscFunctionBegin;
406: PetscCall(PetscBLASIntCast(ldb,&lb));
407: PetscCall(PetscBLASIntCast(ldx,&lx));
408: PetscCall(PetscBLASIntCast(m,&n));
409: PetscCall(PetscBLASIntCast(6*m,&lwork));
410: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
412: #if !PetscDefined(USE_COMPLEX)
413: PetscCall(PetscMalloc6(m*m,&Q,m*m,&W,m*m,&Z,m,&wr,m,&wi,lwork,&work));
414: #else
415: PetscCall(PetscMalloc6(m*m,&Q,m*m,&W,m*m,&Z,m,&wr,lwork,&work,m,&rwork));
416: #endif
418: /* save a copy W = A */
419: for (j=0;j<m;j++) {
420: for (i=0;i<m;i++) W[i+j*m] = A[i+j*lda];
421: }
423: /* compute the (real) Schur form of W */
424: #if !PetscDefined(USE_COMPLEX)
425: PetscCallLAPACKInfo("LAPACKgees",LAPACKgees_("V","N",NULL,&n,W,&n,&sdim,wr,wi,Q,&n,work,&lwork,NULL,&info));
426: #else
427: PetscCallLAPACKInfo("LAPACKgees",LAPACKgees_("V","N",NULL,&n,W,&n,&sdim,wr,Q,&n,work,&lwork,rwork,NULL,&info));
428: #endif
430: /* X = -Q'*B*Q */
431: PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&n,&n,&n,&done,Q,&n,B,&lb,&zero,Z,&n));
432: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&dmone,Z,&n,Q,&n,&zero,X,&lx));
434: /* solve triangular Sylvester equation */
435: PetscCallLAPACKInfo("LAPACKtrsyl",LAPACKtrsyl_("N","C",&ione,&n,&n,W,&n,W,&n,X,&lx,&scal,&info));
436: PetscCheck(scal==1.0,PETSC_COMM_SELF,PETSC_ERR_SUP,"Current implementation cannot handle scale factor %g",(double)scal);
438: /* back-transform X = Q*X*Q' */
439: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&n,&n,&done,Q,&n,X,&n,&zero,W,&n));
440: PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&n,&n,&n,&done,W,&n,Q,&n,&zero,X,&lx));
442: PetscCall(PetscFPTrapPop());
443: #if !PetscDefined(USE_COMPLEX)
444: PetscCall(PetscFree6(Q,W,Z,wr,wi,work));
445: #else
446: PetscCall(PetscFree6(Q,W,Z,wr,work,rwork));
447: #endif
448: PetscFunctionReturn(PETSC_SUCCESS);
449: }
451: #endif /* SLEPC_HAVE_SLICOT */
453: /*@
454: LMEDenseLyapunov - Computes the solution of a dense continuous-time Lyapunov
455: equation.
457: Logically Collective
459: Input Parameters:
460: + lme - the linear matrix equation solver context
461: . m - number of rows and columns of `A`
462: . A - coefficient matrix
463: . lda - leading dimension of `A`
464: . B - right-hand side matrix
465: . ldb - leading dimension of `B`
466: - ldx - leading dimension of `X`
468: Output Parameter:
469: . X - the solution
471: Note:
472: The Lyapunov equation has the form $AX + XA^* = -B$, where all are $m\times m$
473: matrices, and $B$ is symmetric.
475: Level: developer
477: .seealso: [](ch:lme), `LMEDenseHessLyapunovChol()`, `LMESolve()`
478: @*/
479: PetscErrorCode LMEDenseLyapunov(LME lme,PetscInt m,PetscScalar A[],PetscInt lda,PetscScalar B[],PetscInt ldb,PetscScalar X[],PetscInt ldx)
480: {
481: #if PetscDefined(USE_INFO)
482: PetscReal error;
483: #endif
485: PetscFunctionBegin;
488: PetscAssertPointer(A,3);
490: PetscAssertPointer(B,5);
492: PetscAssertPointer(X,7);
495: #if defined(SLEPC_HAVE_SLICOT)
496: PetscCall(Lyapunov_SLICOT(m,A,lda,B,ldb,X,ldx));
497: #else
498: PetscCall(Lyapunov_LAPACK(m,A,lda,B,ldb,X,ldx));
499: #endif
501: #if PetscDefined(USE_INFO)
502: if (PetscLogPrintInfo) {
503: PetscCall(LyapunovResidual(m,A,lda,B,ldb,X,ldx,&error));
504: PetscCall(PetscInfo(lme,"Residual norm of dense Lyapunov equation = %g\n",(double)error));
505: }
506: #endif
507: PetscFunctionReturn(PETSC_SUCCESS);
508: }