fix for complex eltypes (#35)

* fix for complex eltypes

* version bump

Project.toml

@@ -1,7 +1,7 @@
 name = "FiniteHorizonGramians"
 uuid = "b59a298d-d283-4a37-9369-85a9f9a111a5"
 authors = ["Filip Tronarp <[email protected]>", "Tony Stillfjord <[email protected]>, and contributors"]
-version = "0.1.0"
+version = "0.1.1"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/exp_and_gram.jl

@@ -13,7 +13,7 @@ Current supported values of N are 3, 5, 7, 9, 13.
 struct ExpAndGram{T,N,A,B} <: AbstractExpAndGramAlgorithm where {T,N,A,B}
     pade_num::A
     gram_coeffs::B
-    normtol::T
+    normtol::Real
 end
 
 function exp_and_gram(
@@ -479,7 +479,7 @@ function ExpAndGram{T,3}() where {T}
         0.0 0.0 4.47213595499958 0.0
         0.0 0.0 0.0 0.37796447300922725
     ]
-    normtol = T(0.00067)
+    normtol = real(T)(0.00067)
     ExpAndGram{T,3,typeof(pade_num),typeof(gramcs)}(pade_num, gramcs, normtol)
 end
 
@@ -493,7 +493,7 @@ function ExpAndGram{T,5}() where {T}
         0.0 0.0 0.0 0.0 6.0 0.0
         0.0 0.0 0.0 0.0 0.0 0.30151134457776363
     ]
-    normtol = T(0.021)
+    normtol = real(T)(0.021)
     ExpAndGram{T,5,typeof(pade_num),typeof(gramcs)}(pade_num, gramcs, normtol)
 end
 
@@ -509,7 +509,7 @@ function ExpAndGram{T,7}() where {T}
         0.0 0.0 0.0 0.0 0.0 0.0 7.211102550927978 0.0
         0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.25819888974716115
     ]
-    normtol = T(0.13)
+    normtol = real(T)(0.13)
     ExpAndGram{T,7,typeof(pade_num),typeof(gramcs)}(pade_num, gramcs, normtol)
 end
 
@@ -538,7 +538,7 @@ function ExpAndGram{T,9}() where {T}
         0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 8.246211251235321 0.0
         0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.22941573387056177
     ]
-    normtol = T(0.41)
+    normtol = real(T)(0.41)
     ExpAndGram{T,9,typeof(pade_num),typeof(gramcs)}(pade_num, gramcs, normtol)
 end
 
@@ -575,6 +575,6 @@ function ExpAndGram{T,13}() where {T}
         0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 10.0 0.0
         0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.19245008972987526
     ]
-    normtol = T(1.57)
+    normtol = real(T)(1.57)
     ExpAndGram{T,13,typeof(pade_num),typeof(gramcs)}(pade_num, gramcs, normtol)
 end

src/utils.jl

@@ -1,9 +1,9 @@
 """
     _dims_if_compatible(A::AbstractMatrix, B::AbstractMatrix)
 
-Throws DimensionMismatch if A is not square or if the number of rows of B does not equal 
-the number of columns of A. 
-Is equivalent to size(B) if no error is thrown. 
+Throws DimensionMismatch if A is not square or if the number of rows of B does not equal
+the number of columns of A.
+Is equivalent to size(B) if no error is thrown.
 """
 function _dims_if_compatible(A::AbstractMatrix, B::AbstractMatrix)
     n, m = size(B)
@@ -19,9 +19,9 @@ end
 """
     triu2cholesky_factor!(A::AbstractMatrix{T})
 
-If A is an upper triangular matrix, it computes the product Q*A in-place, 
-where Q is a unitary transform such that Q*A is a valid Cholesky factor. 
-If A is not an upper triangular matrix, returns garbage. 
+If A is an upper triangular matrix, it computes the product Q*A in-place,
+where Q is a unitary transform such that Q*A is a valid Cholesky factor.
+If A is not an upper triangular matrix, returns garbage.
 """
 function triu2cholesky_factor!(A::AbstractMatrix{T}) where {T<:Number}
     LinearAlgebra.require_one_based_indexing(A)
@@ -39,7 +39,7 @@ end
 """
     _symmetrize!(A::AbstractMatrix{T}) where {T<:Number}
 
-Discards the skew-Hermitian part of A in-place. 
+Discards the skew-Hermitian part of A in-place.
 """
 function _symmetrize!(A::AbstractMatrix{T}) where {T<:Number}
     LinearAlgebra.require_one_based_indexing(A)
@@ -48,7 +48,7 @@ function _symmetrize!(A::AbstractMatrix{T}) where {T<:Number}
         for row = 1:col
             A[row, col] = (A[row, col] + A[col, row]') / 2
             if row != col
-                A[col, row] = A[row, col]
+                A[col, row] = adjoint(A[row, col])
             end
         end
     end

test/runtests.jl

@@ -10,7 +10,7 @@ include("test_exp_and_gram.jl")
 include("test_adaptive_exp_and_gram.jl")
 include("test_autodiff.jl")
 
-numeric_types = (Float64,)
+numeric_types = (Float64, ComplexF64)
 
 @testset "FiniteHorizonGramians.jl" begin
 

test/test_adaptive_exp_and_gram.jl

@@ -17,6 +17,7 @@ function test_adaptive_exp_and_gram(T)
         Φ, G = exp_and_gram(A, B, t, method)
         @test isapprox(err(Φ, Φgt), zero(T), atol = tol)
         @test isapprox(err(G, Ggt), zero(T), atol = tol)
+        @test ishermitian(G)
     end
     Φgt, Ggt = mf_exp_and_gram(A, B)
     Φ, G = exp_and_gram(A, B, method)
@@ -40,10 +41,10 @@ function test_adaptive_exp_and_gram(T)
 
     # test covering the case when initial cholesky factor is non-square
     q = 13
-    m = 1 
-    n = m * (q + 2) 
-    A = - tril(ones(n, n))
-    B = ones(n, m)
+    m = 1
+    n = m * (q + 2)
+    A = - tril(ones(T, n, n))
+    B = ones(T, n, m)
     @test_nowarn exp_and_gram(A, B, method)
 
 end

test/test_exp_and_gram.jl

@@ -4,8 +4,6 @@ function test_exp_and_gram(T)
     A = randn(T, n, n)
     B = randn(T, n, m)
 
-    tol = 1e-11
-
     ts = [0.5, 1.0, 1.5]
 
     qs = [3, 5, 7, 9, 13]
@@ -23,6 +21,7 @@ function test_exp_and_gram(T)
                 Φ, G = exp_and_gram(A, B, t, method)
                 @test isapprox(err(Φ, Φgt), zero(T), atol = tols[i])
                 @test isapprox(err(G, Ggt), zero(T), atol = tols[i])
+                @test ishermitian(G)
             end
             Φgt, Ggt = mf_exp_and_gram(A, B)
             Φ, G = exp_and_gram(A, B, method)
@@ -49,8 +48,8 @@ function test_exp_and_gram(T)
             # test covering the case when initial cholesky factor is non-square
             m = 1
             n = m * (q + 2)
-            A = - tril(ones(n, n))
-            B = ones(n, m)
+            A = - tril(ones(T, n, n))
+            B = ones(T, n, m)
             @test_nowarn exp_and_gram_chol(A, B, method)
         end