From 7a47225de4ca135247790707bf01f71d9514eaac Mon Sep 17 00:00:00 2001 From: mschauer <moritzschauer@web.de> Date: Fri, 27 Jun 2014 18:16:27 +0200 Subject: [PATCH 1/6] Add lyap() and sylvester(). --- base/exports.jl | 2 ++ base/linalg.jl | 2 ++ base/linalg/dense.jl | 26 ++++++++++++++++++++++++++ base/linalg/lapack.jl | 37 +++++++++++++++++++++++++++++++++++++ 4 files changed, 67 insertions(+) diff --git a/base/exports.jl b/base/exports.jl index 816aaf07c6704..ef56cc6dd23f6 100644 --- a/base/exports.jl +++ b/base/exports.jl @@ -654,6 +654,7 @@ export lu, lufact!, lufact, + lyap, norm, null, peakflops, @@ -677,6 +678,7 @@ export svdfact, svdvals!, svdvals, + sylvester, trace, transpose, tril!, diff --git a/base/linalg.jl b/base/linalg.jl index 3ccc3ff508521..e463f50f17b33 100644 --- a/base/linalg.jl +++ b/base/linalg.jl @@ -85,6 +85,7 @@ export lu, lufact, lufact!, + lyap, norm, null, peakflops, @@ -108,6 +109,7 @@ export svdfact, svdvals!, svdvals, + sylvester, trace, transpose, tril, diff --git a/base/linalg/dense.jl b/base/linalg/dense.jl index f28dc5a24a10e..4e8de64f2f027 100644 --- a/base/linalg/dense.jl +++ b/base/linalg/dense.jl @@ -457,3 +457,29 @@ function cond(A::StridedMatrix, p::Real=2) end throw(ArgumentError("invalid p-norm p=$p. Valid: 1, 2 or Inf")) end + +## Lyapunov and Sylvester equation + +# AX + XB + C = 0 +function sylvester{T<:BlasFloat}(A::StridedMatrix{T},B::StridedMatrix{T},C::StridedMatrix{T}) + RA, QA = schur(A) + RB, QB = schur(B') + + D = -QA'*C*QB + Y, scale = LAPACK.trsyl!('N', T <: Complex ? 'C' : 'T', RA, RB, D) + scale!(QA*Y*QB', inv(scale)) +end +sylvester{T<:Integer}(A::StridedMatrix{T},B::StridedMatrix{T},C::StridedMatrix{T}) = sylvester(float(A), float(B), float(C)) + +# AX + XA' + C = 0 +function lyap{T<:BlasFloat}(A::StridedMatrix{T},C::StridedMatrix{T}) + chksquare(A, C) + R, Q = schur(A) + + D = -Q'*C*Q + Y, scale = LAPACK.trsyl!('N', T <: Complex ? 'C' : 'T', R, R, D) + scale!(Q*Y*Q', inv(scale)) +end +lyap{T<:Integer}(A::StridedMatrix{T},C::StridedMatrix{T}) = lyap(float(A), float(C)) +lyap{T<:Number}(a::T, c::T) = -c/(2a) + diff --git a/base/linalg/lapack.jl b/base/linalg/lapack.jl index c9a8f71ca1d38..7f79660a16871 100644 --- a/base/linalg/lapack.jl +++ b/base/linalg/lapack.jl @@ -3665,4 +3665,41 @@ for (fn, elty) in ((:dtrttf_, :Float64), end end +# Solves the real Sylvester matrix equation: op(A)*X +- X*op(B) = scale*C and A and B are both upper quasi triangular. +for (fn, elty) in ((:dtrsyl_, :Float64), + (:strsyl_, :Float32), + (:ztrsyl_, :Complex128), + (:ctrsyl_, :Complex64)) + @eval begin + function trsyl!(transa::BlasChar, transb::BlasChar, A::StridedMatrix{$elty}, B::StridedMatrix{$elty}, C::StridedMatrix{$elty}, isgn::BlasInt=1) + chkstride1(A) + chkstride1(B) + chkstride1(C) + m = size(A, 1) + n = size(B, 1) + lda = max(1, stride(A, 2)) + ldb = max(1, stride(B, 2)) + if lda < m throw(DimensionMismatch("")) end + if ldb < n throw(DimensionMismatch("")) end + m1, n1 = size(C) + if m != m1 || n != n1 throw(DimensionMismatch("")) end + ldc = max(1, stride(C, 2)) + + scale = Array($elty, 1) + info = Array(BlasInt, 1) + + ccall(($(string(fn)), liblapack), Void, + (Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, + Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, + Ptr{$elty}, Ptr{BlasInt}), + &transa, &transb, &isgn, &m, &n, + A, &lda, B, &ldb, C, &ldc, + scale, info) + @lapackerror + C, scale[1] + end + end +end + + end # module From 8231bdeda24a7c63063ebffc55ef216c4e2f1110 Mon Sep 17 00:00:00 2001 From: mschauer <moritzschauer@web.de> Date: Fri, 27 Jun 2014 19:14:48 +0200 Subject: [PATCH 2/6] Test for lyap() and sylvester() and taking some tests out of the inner loop in linalg1.jl --- test/linalg1.jl | 37 +++++++++++++++++++++++++++---------- 1 file changed, 27 insertions(+), 10 deletions(-) diff --git a/test/linalg1.jl b/test/linalg1.jl index bb6b2a67990bd..2b4b76113128c 100644 --- a/test/linalg1.jl +++ b/test/linalg1.jl @@ -2,7 +2,7 @@ debug = false import Base.LinAlg: BlasComplex, BlasFloat, BlasReal, QRPivoted -n = 10 +n = 10 srand(1234321) a = rand(n,n) @@ -17,16 +17,20 @@ end areal = randn(n,n)/2 aimg = randn(n,n)/2 +a2real = randn(n,n)/2 +a2img = randn(n,n)/2 breal = randn(n,2)/2 bimg = randn(n,2)/2 + for eltya in (Float32, Float64, Complex64, Complex128, BigFloat, Int) + a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(areal, aimg) : areal) + a2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(a2real, a2img) : a2real) + asym = a'+a # symmetric indefinite + apd = a'*a # symmetric positive-definite + ε = εa = eps(abs(float(one(eltya)))) + for eltyb in (Float32, Float64, Complex64, Complex128, Int) - a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(areal, aimg) : areal) b = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex(breal, bimg) : breal) - asym = a'+a # symmetric indefinite - apd = a'*a # symmetric positive-definite - - εa = eps(abs(float(one(eltya)))) εb = eps(abs(float(one(eltyb)))) ε = max(εa,εb) @@ -249,8 +253,12 @@ debug && println("Test null") @test size(null(b), 2) == 0 end + end # for eltyb + +debug && println("\ntype of a: ", eltya, "\n") + debug && println("Test pinv") - if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia + if eltya != BigFloat # Revisit when implemented in julia pinva15 = pinv(a[:,1:5]) @test_approx_eq a[:,1:5]*pinva15*a[:,1:5] a[:,1:5] @test_approx_eq pinva15*a[:,1:5]*pinva15 pinva15 @@ -258,14 +266,23 @@ debug && println("Test pinv") # if isreal(a) debug && println("Matrix square root") - if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia + if eltya != BigFloat # Revisit when implemented in julia asq = sqrtm(a) @test_approx_eq asq*asq a asymsq = sqrtm(asym) @test_approx_eq asymsq*asymsq asym end -end -end + +debug && println("Lyapunov/Sylvester") + if eltya != BigFloat + let + x = lyap(a, a2) + @test_approx_eq -a2 a*x + x*a' + x2 = sylvester(a[1:3, 1:3], a[4:n, 4:n], a2[1:3,4:n]) + @test_approx_eq -a2[1:3, 4:n] a[1:3, 1:3]*x2 + x2*a[4:n, 4:n] + end + end +end # for eltya #6941 #@test (ones(10^7,4)*ones(4))[3] == 4.0 From ea8f889bc6682ff1cd61d8016ce9e510200fc481 Mon Sep 17 00:00:00 2001 From: mschauer <moritzschauer@web.de> Date: Fri, 27 Jun 2014 20:43:58 +0200 Subject: [PATCH 3/6] Fix: Scale is real --- base/linalg/lapack.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/base/linalg/lapack.jl b/base/linalg/lapack.jl index 7f79660a16871..1eb6d99e38186 100644 --- a/base/linalg/lapack.jl +++ b/base/linalg/lapack.jl @@ -3666,12 +3666,12 @@ for (fn, elty) in ((:dtrttf_, :Float64), end # Solves the real Sylvester matrix equation: op(A)*X +- X*op(B) = scale*C and A and B are both upper quasi triangular. -for (fn, elty) in ((:dtrsyl_, :Float64), - (:strsyl_, :Float32), - (:ztrsyl_, :Complex128), - (:ctrsyl_, :Complex64)) +for (fn, elty, relty) in ((:dtrsyl_, :Float64, :Float64), + (:strsyl_, :Float32, :Float32), + (:ztrsyl_, :Complex128, :Float64), + (:ctrsyl_, :Complex64, :Float32)) @eval begin - function trsyl!(transa::BlasChar, transb::BlasChar, A::StridedMatrix{$elty}, B::StridedMatrix{$elty}, C::StridedMatrix{$elty}, isgn::BlasInt=1) + function trsyl!(transa::BlasChar, transb::BlasChar, A::StridedMatrix{$elty}, B::StridedMatrix{$elty}, C::StridedMatrix{$elty}, isgn::Int=1) chkstride1(A) chkstride1(B) chkstride1(C) @@ -3685,13 +3685,13 @@ for (fn, elty) in ((:dtrsyl_, :Float64), if m != m1 || n != n1 throw(DimensionMismatch("")) end ldc = max(1, stride(C, 2)) - scale = Array($elty, 1) + scale = Array($relty, 1) info = Array(BlasInt, 1) ccall(($(string(fn)), liblapack), Void, (Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, - Ptr{$elty}, Ptr{BlasInt}), + Ptr{$relty}, Ptr{BlasInt}), &transa, &transb, &isgn, &m, &n, A, &lda, B, &ldb, C, &ldc, scale, info) From cceeeda4c917c540eca5b8eef7206e97e6326ec5 Mon Sep 17 00:00:00 2001 From: mschauer <moritzschauer@web.de> Date: Fri, 27 Jun 2014 21:27:08 +0200 Subject: [PATCH 4/6] Doc for lyap() and sylvester(). --- doc/stdlib/linalg.rst | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/doc/stdlib/linalg.rst b/doc/stdlib/linalg.rst index 92be96b0d4795..ed4b2cd523a3c 100644 --- a/doc/stdlib/linalg.rst +++ b/doc/stdlib/linalg.rst @@ -463,6 +463,14 @@ Linear algebra functions in Julia are largely implemented by calling functions f Matrix exponential. +.. function:: lyap(A, C) + + Computes the solution ``X`` to the continuous Lyapunov equation ``AX + XA' + C = 0``, where no eigenvalue of ``A`` has a zero real part and no two eigenvalues are negative complex conjugates of each other. + +.. function:: sylvester(A, B, C) + + Computes the solution ``X`` to the Sylvester equation ``AX + XB + C = 0``, where ``A``, ``B`` and ``C`` have compatible dimensions and ``A`` and ``-B`` have no eigenvalues with equal real part. + .. function:: issym(A) -> Bool Test whether a matrix is symmetric. From 26e7b87709a0c589a6d32cedc63198e9fa34e912 Mon Sep 17 00:00:00 2001 From: mschauer <moritzschauer@web.de> Date: Sat, 28 Jun 2014 02:09:33 +0200 Subject: [PATCH 5/6] Avoid explicit transposes --- base/linalg/dense.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/base/linalg/dense.jl b/base/linalg/dense.jl index 4e8de64f2f027..3da6522702b77 100644 --- a/base/linalg/dense.jl +++ b/base/linalg/dense.jl @@ -463,11 +463,11 @@ end # AX + XB + C = 0 function sylvester{T<:BlasFloat}(A::StridedMatrix{T},B::StridedMatrix{T},C::StridedMatrix{T}) RA, QA = schur(A) - RB, QB = schur(B') + RB, QB = schur(B) - D = -QA'*C*QB - Y, scale = LAPACK.trsyl!('N', T <: Complex ? 'C' : 'T', RA, RB, D) - scale!(QA*Y*QB', inv(scale)) + D = -Ac_mul_B(QA,C*QB) + Y, scale = LAPACK.trsyl!('N','N', RA, RB, D) + scale!(QA*A_mul_Bc(Y,QB), inv(scale)) end sylvester{T<:Integer}(A::StridedMatrix{T},B::StridedMatrix{T},C::StridedMatrix{T}) = sylvester(float(A), float(B), float(C)) @@ -476,9 +476,9 @@ function lyap{T<:BlasFloat}(A::StridedMatrix{T},C::StridedMatrix{T}) chksquare(A, C) R, Q = schur(A) - D = -Q'*C*Q + D = -Ac_mul_B(Q,C*Q) Y, scale = LAPACK.trsyl!('N', T <: Complex ? 'C' : 'T', R, R, D) - scale!(Q*Y*Q', inv(scale)) + scale!(Q*A_mul_Bc(Y,Q), inv(scale)) end lyap{T<:Integer}(A::StridedMatrix{T},C::StridedMatrix{T}) = lyap(float(A), float(C)) lyap{T<:Number}(a::T, c::T) = -c/(2a) From 80abe17f9062ff17a727c1028284e85c2a105876 Mon Sep 17 00:00:00 2001 From: mschauer <moritzschauer@web.de> Date: Mon, 30 Jun 2014 15:59:25 +0200 Subject: [PATCH 6/6] Use chksquare in trsyl! --- base/linalg/dense.jl | 1 - base/linalg/lapack.jl | 9 ++------- 2 files changed, 2 insertions(+), 8 deletions(-) diff --git a/base/linalg/dense.jl b/base/linalg/dense.jl index 3da6522702b77..135813eda5f1a 100644 --- a/base/linalg/dense.jl +++ b/base/linalg/dense.jl @@ -473,7 +473,6 @@ sylvester{T<:Integer}(A::StridedMatrix{T},B::StridedMatrix{T},C::StridedMatrix{T # AX + XA' + C = 0 function lyap{T<:BlasFloat}(A::StridedMatrix{T},C::StridedMatrix{T}) - chksquare(A, C) R, Q = schur(A) D = -Ac_mul_B(Q,C*Q) diff --git a/base/linalg/lapack.jl b/base/linalg/lapack.jl index 1eb6d99e38186..f271b1d91adb8 100644 --- a/base/linalg/lapack.jl +++ b/base/linalg/lapack.jl @@ -3672,15 +3672,10 @@ for (fn, elty, relty) in ((:dtrsyl_, :Float64, :Float64), (:ctrsyl_, :Complex64, :Float32)) @eval begin function trsyl!(transa::BlasChar, transb::BlasChar, A::StridedMatrix{$elty}, B::StridedMatrix{$elty}, C::StridedMatrix{$elty}, isgn::Int=1) - chkstride1(A) - chkstride1(B) - chkstride1(C) - m = size(A, 1) - n = size(B, 1) + chkstride1(A, B, C) + m, n = chksquare(A, B) lda = max(1, stride(A, 2)) ldb = max(1, stride(B, 2)) - if lda < m throw(DimensionMismatch("")) end - if ldb < n throw(DimensionMismatch("")) end m1, n1 = size(C) if m != m1 || n != n1 throw(DimensionMismatch("")) end ldc = max(1, stride(C, 2))