Added naive method for generic inputs.

Author: olof3

src/LyapTest.jl

@@ -80,7 +80,8 @@ end
 LinearAlgebra.schur(A::AbstractMatrix{T}) where T = schur!(LinearAlgebra.copy_oftype(A, LinearAlgebra.eigtype(T)))
 
 
-sylvcsoltype(A, B, C) = Base.promote_op(sylvc, eltype(A), eltype(B), eltype(C))
+sylvcsoltype(A, B, C) = Base.promote_op((a,b,c) -> c / (a + b), eltype(A), eltype(B), eltype(C))
+#sylvcsoltype(A, B, C) = Base.promote_op(sylvc, eltype(A), eltype(B), eltype(C))
 sylvdsoltype(A, B, C) = Base.promote_op(sylvd, eltype(A), eltype(B), eltype(C))
 
 """
@@ -236,6 +237,10 @@ function lyapc(A, Q)
 
      _check_lyap_inputs(A, Q)
 
+    if !hasmethod(schur!, (typeof(A),))
+        return lyapc(A, Q, Val(:naive))
+    end
+
     At2, U = schur(A')
 
     Q2 = U'*Q*U
@@ -244,6 +249,57 @@ function lyapc(A, Q)
 
     X = mul!(Y, U, Y*U')
 end
+
+
+function sylvc(A, B, C, ::Val{:naive})
+    Xv = kron(I(size(B,1)), A) + kron(transpose(B), I(size(A,1))) \ C[:]
+    return reshape(Xv, size(C))
+end
+
+# Mapping from Cartesian index into vector represenentation of Symmetric matrix
+@inline sub2triidx(i,j) = (j >= i ? (j*(j-1))>>>1 + i : (i*(i-1))>>>1 + j)
+
+function lyapc(A, Q, ::Val{:naive}) # Only works for real matrices A
+    # Sets up and solves a system of equations for the upper triangular part of X
+    # and solves that equation. This gives an n(n+1)/2 system instead of n^2
+    # ONLY WORKS FOR REAL A!
+    # Should be able to to base the discrete time version on this as well
+    _check_lyap_inputs(A, Q)
+
+    if !isreal(A); error("Only works for real A matrices"); end # Should call sylvc
+
+    n = size(Q, 1)
+    nR = (n*(n+1)) >>> 1 # Ssize of the system to be solved
+
+    R = zeros(eltype(A), nR, nR)
+    minusQv = Vector{eltype(Q)}(undef, nR)
+
+    for j=1:n
+        for i=1:j
+            m = sub2triidx(i,j) # Set up equation for X[i,j]
+            minusQv[m] = -Q[i,j]
+            # (AX + XA')[i,j] = sum(A[i,k]*X[k,j]) + sum(X[i,k]*A'[k,j])
+            # the r = trinum[i]+r gives the kth element in the upper traingle
+            # which correpsonds to X[i,j]
+            for k=1:n
+                R[m, sub2triidx(k,j)] += A[i,k]
+                R[m, sub2triidx(i,k)] += A[j,k]
+            end
+        end
+    end
+
+    Xv = R \ minusQv
+
+    # Fill the upper traingle of X from Xv
+    X = [Xv[sub2triidx(i,j)] for i=1:n, j=1:n]
+
+    return X
+end
+
+
+
+
+
 """
     X = lyapd(A, Q) -> X
 
@@ -316,6 +372,10 @@ function _sylvc_schur!(A::Matrix, B::Matrix, C::Matrix, alg::Union{Val{:sylv},Va
 
     @inbounds for j=1:nblocksb
         i0 = (alg === Val(:lyap) ? j : 1)
+
+        # if j > 1; mul!(C[i0:nblocksa,j], C[i0:nblocksa, 1:j-1], B[1:j-1, j], 1, -1); end # Could move this out?
+        # figure out the row indexing
+
         for i=i0:nblocksa
             if schurtype === Val(:complex)
                 if i > 1; C[i,j] -= sum(A[i, k] * C[k, j] for k=1:i-1); end

test/test_matrix_equations.jl

@@ -76,6 +76,9 @@ Cherm = [1.0 2+im; 2-im 1]
 @test LyapTest.lyapd(A, C) ≈ sylvd_diag(a, a, -C)
 @test LyapTest.lyapd(Ac, Cherm) ≈ sylvd_diag(ac, conj(ac), -Cherm)
 
+# A few tests for the naive version, should have more
+@test LyapTest.sylvc(A, B, C, Val(:naive)) ≈ sylvc_diag(a, a, -C)
+@test LyapTest.lyapc(A, C, Val(:naive)) ≈ sylvc_diag(a, a, -C)