lyap function with 1x1 and 2x2 close form solutions (#246)

Author: mschauer

src/StaticArrays.jl

@@ -8,7 +8,7 @@ import Base: getindex, setindex!, size, similar, vec, show,
              length, convert, promote_op, promote_rule, map, map!, reduce, reducedim, mapreducedim,
              mapreduce, broadcast, broadcast!, conj, transpose, ctranspose,
              hcat, vcat, ones, zeros, eye, one, cross, vecdot, reshape, fill,
-             fill!, det, logdet, inv, eig, eigvals, expm, logm, sqrtm, trace, diag, vecnorm, norm, dot, diagm, diag,
+             fill!, det, logdet, inv, eig, eigvals, expm, logm, sqrtm, lyap, trace, diag, vecnorm, norm, dot, diagm, diag,
              lu, svd, svdvals, svdfact, factorize, ishermitian, issymmetric, isposdef,
              sum, diff, prod, count, any, all, minimum,
              maximum, extrema, mean, copy, rand, randn, randexp, rand!, randn!,
@@ -95,6 +95,7 @@ include("solve.jl")
 include("eigen.jl")
 include("expm.jl")
 include("sqrtm.jl")
+include("lyap.jl")
 include("triangular.jl")
 include("cholesky.jl")
 include("svd.jl")

src/lyap.jl

@@ -0,0 +1,12 @@
+lyap(a::StaticMatrix,  c::StaticMatrix) = _lyap(Size(a), Size(c), a, c)
+
+_lyap(::Size{(1,1)}, ::Size{(1,1)}, a::StaticMatrix,  c::StaticMatrix) = -c/(2a[1,1])
+
+@inline function _lyap(::Size{(2,2)}, ::Size{(2,2)}, a::StaticMatrix,  c::StaticMatrix)
+    d = det(a)
+    t = trace(a)
+     -(d*c  + (a - t*I)*c*(a-t*I)')/(2*d*t) # http://www.nber.org/papers/w8956.pdf
+end
+
+@inline _lyap(sa::Size, sc::Size, a::StaticMatrix, c::StaticMatrix) = sc(Base.lyap(Array(a),Array(c)))
+ 

test/lyap.jl

@@ -0,0 +1,9 @@
+@testset "Lyapunov equation" begin
+    @test lyap(@SMatrix([1]), @SMatrix [2]) === @SMatrix [-1.0]
+
+    @test lyap(@SMatrix([1 1; 0 1]), @SMatrix [2 1; 1 2]) == -eye(2)
+    @test isa(lyap(@SMatrix([1 1; 0 1]), @SMatrix [2 1; 1 2]), SArray{Tuple{2,2},Float64,2,4})
+    
+    @test lyap(@SMatrix([1 0 1; 0 1 0; 0 0 1]), @SMatrix [2 4 4; 0 4 20; 0 0 20]) == [-5.0 -2.0 3.0; 5.0 -2.0 -10.0; 5.0 0.0 -10.0]
+    @test isa(lyap(@SMatrix([1 0 1; 0 1 0; 0 0 1]), @SMatrix [2 4 4; 0 4 20; 0 0 20]), SizedArray{Tuple{3,3},Float64,2,2})
+end

test/runtests.jl

@@ -30,6 +30,7 @@ include("solve.jl")
 include("eigen.jl")
 include("expm.jl")
 include("sqrtm.jl")
+include("lyap.jl")
 include("chol.jl")
 include("deque.jl")
 include("io.jl")