Skip to content

Commit

Permalink
case n==1 removed from exp.jl and jmatrix added (#84)
Browse files Browse the repository at this point in the history
* hessenberg merged

* hessenberg merged

* pfaffian.jl

* pfaffian.jl

* non blas type

* non blas type

* non blas type

* eigen complex friendly

* eigen complex friendly

* cholesky.jl added

* complex trigo

* tridiag tests added, fixed copy,conj,getindex

* complex eigen etc tridiagonal, Int32 friendly

* complex eigen etc tridiagonal, Int32 friendly

* complex eigen etc tridiagonal, Int32 friendly

* Fix lapackerror 22and cleaned code

* Fix lapackerror 22and cleaned code

* size decreased

* hessenberg limit null case fixed

* lapack 22 friendly

* documentation

* trigo for tridiag

* Goodbye error22+little changes

* Goodbye error22+little changes second chance

* logpfaffian+rdiv,ldiv etc for vectors, ambiguity on *,/?... stiil to fix

* logpfaffian+rdiv,ldiv etc for vectors, ambiguity on *,/?... stiil to fix

* sincos(A)+ more tests

* Documenter

* Documenter step2

* Documenter step 3

* type stable exp file

* first jmatrix.jl file

* delete n==1 exp.jl and jmatrix.jl added
  • Loading branch information
smataigne authored Aug 22, 2022
1 parent 3a45df3 commit 38acc2a
Show file tree
Hide file tree
Showing 5 changed files with 129 additions and 54 deletions.
1 change: 1 addition & 0 deletions src/SkewLinearAlgebra.jl
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@ export

include("skewhermitian.jl")
include("tridiag.jl")
include("jmatrix.jl")
include("hessenberg.jl")
include("eigen.jl")
include("exp.jl")
Expand Down
14 changes: 5 additions & 9 deletions src/cholesky.jl
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@

struct SkewCholesky{T,R<:UpperTriangular{<:T},J<:SkewHermTridiagonal{<:T},P<:AbstractVector{<:Integer}}
struct SkewCholesky{T,R<:UpperTriangular{<:T},J<:JMatrix{<:T},P<:AbstractVector{<:Integer}}
Rm::R #Uppertriangular matrix
Jm::J # Block diagonal skew-symmetric matrix
Jm::J # Block diagonal skew-symmetric matrix of type JMatrix
Pv::P #Permutation vector

function SkewCholesky{T,R,J,P}(Rm,Jm,Pv) where {T,R,J,P}
Expand All @@ -14,17 +14,13 @@ end
SkewCholesky(Rm,Pv)
Construct a `SkewCholesky` structure from the `UpperTriangular`
matrix `Rm` and the permutation vector `Pv`. A matrix `Jm` of type `SkewHermTridiagonal`
matrix `Rm` and the permutation vector `Pv`. A matrix `Jm` of type `JMatrix`
is build calling this function.
The `SkewCholesky` structure has three arguments: `Rm`,`Jm` and `Pv`.
"""
function SkewCholesky(Rm::UpperTriangular{<:T},Pv::AbstractVector{<:Integer}) where {T<:Real}
n = size(Rm, 1)
vec = zeros(T, n - 1)
for i = 1 : 2 : n - 1
vec[i] = -1
end
return SkewCholesky{T,UpperTriangular{<:T},SkewHermTridiagonal{<:T},AbstractVector{<:Integer}}(Rm, SkewHermTridiagonal(vec), Pv)
return SkewCholesky{T,UpperTriangular{<:T},JMatrix{<:T},AbstractVector{<:Integer}}(Rm, JMatrix(T, n), Pv)

end

Expand Down Expand Up @@ -100,7 +96,7 @@ skewchol!(A::AbstractMatrix) = @views skewchol!(SkewHermitian(A))
Computes a Cholesky-like factorization of the real skew-symmetric matrix `A`.
The function returns a `SkewCholesky` structure composed of three fields:
`Rm`,`Jm`,`Pv`. `Rm` is `UpperTriangular`, `Jm` is `SkewHermTridiagonal`,
`Rm`,`Jm`,`Pv`. `Rm` is `UpperTriangular`, `Jm` is a `JMatrix`,
`Pv` is an array of integers. Let `S` be the returned structure, then the factorization
is such that `S.Rm'*S.Jm*S.Rm = A[S.Pv,S.Pv]`
Expand Down
44 changes: 0 additions & 44 deletions src/exp.jl
Original file line number Diff line number Diff line change
Expand Up @@ -4,10 +4,8 @@ function skewexp!(A::Union{SkewHermitian{T},SkewHermTridiagonal{T}}) where {T<:R
n = size(A, 1)

if typeof(A) <:SkewHermitian
n == 1 && return fill(T(1), 1, 1)
vals, Qr, Qim = skeweigen!(A)
else
n == 1 && return fill(T(1), 1, 1)
E = eigen!(A)
vals = E.values
Qr = real(E.vectors)
Expand Down Expand Up @@ -39,13 +37,6 @@ end

@views function skewexp!(A::Union{SkewHermitian{<:Complex},SkewHermTridiagonal{<:Complex}})
n = size(A, 1)
if n == 1
if typeof(A)<:SkewHermitian
return fill(exp(A.data[1,1]), 1, 1)
else
return fill(exp(complex(0, A.dvim[1])), 1, 1)
end
end
Eig = eigen!(A)
eig = exp.(Eig.values)
temp = similar(A, n, n)
Expand All @@ -59,7 +50,6 @@ Base.exp(A::Union{SkewHermitian,SkewHermTridiagonal}) = skewexp!(copyeigtype(A))

@views function skewcis!(A::Union{SkewHermitian{T},SkewHermTridiagonal{T}}) where {T<:Real}
n = size(A, 1)
n == 1 && Hermitian(fill(T(1), 1, 1))
Eig = eigen!(A)
Q = Eig.vectors
temp = similar(Q, n, n)
Expand All @@ -73,13 +63,6 @@ end

@views function skewcis!(A::Union{SkewHermitian{<:Complex},SkewHermTridiagonal{<:Complex}})
n = size(A,1)
if n == 1
if typeof(A)<:SkewHermitian
return Hermitian(fill(cis(A.data[1,1]), 1, 1))
else
return Hermitian(fill(cis(complex(0, A.dvim[1])), 1, 1))
end
end
Eig = eigen!(A)
eig = @. exp(-imag(Eig.values))
Cis = similar(A, n, n)
Expand All @@ -92,10 +75,8 @@ end
@views function skewcos!(A::Union{SkewHermitian{T},SkewHermTridiagonal{T}}) where {T<:Real}
n = size(A,1)
if typeof(A) <:SkewHermitian
n == 1 && return Symmetric(fill(T(1), 1, 1))
vals, Qr, Qim = skeweigen!(A)
else
n == 1 && return Symmetric(fill(T(1), 1, 1))
E = eigen!(A)
vals = E.values
Qr = real(E.vectors)
Expand All @@ -116,13 +97,6 @@ end

@views function skewcos!(A::Union{SkewHermitian{<:Complex},SkewHermTridiagonal{<:Complex}})
n = size(A,1)
if n == 1
if typeof(A)<:SkewHermitian
return Hermitian(fill(cos(A.data[1,1]), 1, 1))
else
return Hermitian(fill(cos(complex(0, A.dvim[1])), 1, 1))
end
end
Eig = eigen!(A)
eig1 = @. exp(-imag(Eig.values))
eig2 = @. exp(imag(Eig.values))
Expand All @@ -141,10 +115,8 @@ end
@views function skewsin!(A::Union{SkewHermitian{T},SkewHermTridiagonal{T}}) where {T<:Real}
n = size(A, 1)
if typeof(A) <:SkewHermitian
n == 1 && return fill(T(0), 1, 1)
vals, Qr, Qim = skeweigen!(A)
else
n == 1 && return fill(T(0), 1, 1)
E = eigen!(A)
vals = E.values
Qr = real(E.vectors)
Expand All @@ -165,13 +137,6 @@ end

@views function skewsin!(A::Union{SkewHermitian{<:Complex},SkewHermTridiagonal{<:Complex}})
n = size(A,1)
if n == 1
if typeof(A)<:SkewHermitian
return fill(sin(A.data[1,1]), 1, 1)
else
return fill(sin(complex(0, A.dvim[1])), 1, 1)
end
end
Eig = eigen!(A)
eig1 = @. exp(-imag(Eig.values))
eig2 = @. exp(imag(Eig.values))
Expand All @@ -195,10 +160,8 @@ Base.sin(A::Union{SkewHermitian,SkewHermTridiagonal}) = skewsin!(copyeigtype(A))
@views function skewsincos!(A::Union{SkewHermitian{T},SkewHermTridiagonal{T}}) where {T<:Real}
n = size(A,1)
if typeof(A) <:SkewHermitian
n == 1 && return fill(T(0), 1, 1), Symmetric(fill(T(1), 1, 1))
vals, Qr, Qim = skeweigen!(A)
else
n == 1 && return fill(T(0), 1, 1), Symmetric(fill(T(1), 1, 1))
E = eigen!(A)
vals = E.values
Qr = real(E.vectors)
Expand All @@ -224,13 +187,6 @@ Base.sin(A::Union{SkewHermitian,SkewHermTridiagonal}) = skewsin!(copyeigtype(A))
end
@views function skewsincos!(A::Union{SkewHermitian{<:Complex},SkewHermTridiagonal{<:Complex}})
n = size(A, 1)
if n == 1
if typeof(A)<:SkewHermitian
return fill(sin(A.data[1,1]), 1, 1), Hermitian(fill(cos(A.data[1,1]), 1, 1))
else
return fill(sin(A.data[1,1]), 1, 1), Hermitian(fill(cos(complex(0, A.dvim[1])), 1, 1))
end
end
Eig = eigen!(A)
eig1 = @. exp(-imag(Eig.values))
eig2 = @. exp(imag(Eig.values))
Expand Down
96 changes: 96 additions & 0 deletions src/jmatrix.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,96 @@

struct JMatrix{T, N<:Integer, SGN} <: AbstractMatrix{T}
n::N #size of the square matrix
sgn::SGN #+-1, allows to transpose,invert easily the matrix.
function JMatrix{T, N, SGN}(n, sgn) where {T, N<:Integer, SGN}
(sgn == T(1) || sgn == T(-1) ) || throw("sgn argument must be +-1")
new{T, N, SGN}(n,sgn)
end
end

JMatrix(T::Type, n::Integer) = JMatrix{T,typeof(n),Any}(n, T(1))
JMatrix(T::Type, n::Integer, sgn) = JMatrix{T,typeof(n), Any}(n, T(sgn))

Base.similar(J::JMatrix,::Type{T}) where {T} = JMatrix(T, J.n, J.sgn)
Base.similar(J::JMatrix, ::Type{T}, dims::Union{Dims{1},Dims{2}}) where {T} = zeros(T, dims...)
Base.copy(J::JMatrix{T}) where T = JMatrix(T, J.n, J.sgn)

Base.size(J::JMatrix) = (J.n, J.n)
function Base.size(J::JMatrix, n::Integer)
if n == 1 || n == 2
return J.n
else
return 1
end
end

function Base.Matrix(J::JMatrix{T}) where {T}
M = similar(J, T, J.n, J.n)
for i =1 : 2 : J.n-1
M[i+1,i] = -J.sgn
M[i,i+1] = J.sgn
end
return M
end

Base.Array(A::JMatrix) = Matrix(A)

function SkewHermTridiagonal(J::JMatrix{T}) where {T}
vec = zeros(T, J.n - 1)
for i = 1 : 2 : J.n - 1
vec[i] = -1
end
return SkewHermTridiagonal(vec)
end

Base.@propagate_inbounds function Base.getindex(J::JMatrix{T}, i::Integer, j::Integer) where T

@boundscheck checkbounds(J, i, j)
if i == j + 1 && i%2 == 0
return -J.sgn
elseif i + 1 == j && j%2 == 0
return J.sgn
else
return zero(T)
end
end
#=
function Base.:*(J::JMatrix,A::AbstractVecOrMat)
m, k = size(A, 1), size(1, 2)
if !(m == J.n)
throw(DimensionMismatch("A has first dimension $(size(S,1)), B has $(size(B,1)), C has $(size(C,1)) but all must match"))
end
B = similar(A,J.n, k)
for i = 1 : 2 : J.n-1
B[i,:] .= A[i+1,:]
B[i+1,:] .= -A[i,:]
end
return B
end
function Base.:*(A::AbstractVecOrMat, J::JMatrix)
m, k = size(A, 1), size(1, 2)
if !(k == J.n)
throw(DimensionMismatch("A has first dimension $(size(S,1)), B has $(size(B,1)), C has $(size(C,1)) but all must match"))
end
B = similar(A,m, J.n)
for i = 1 : 2 : J.n-1
B[:,i] .= A[:,i+1]
B[:, i+1] .= -A[:, i]
end
return B
end
\(J::JMatrix,A::AbstractVecOrMat) = - J * A
=#
Base.:-(J::JMatrix{T}) where T = JMatrix(T, J.n, -J.sgn)
LA.transpose(J::JMatrix) = -J
LA.adjoint(J::JMatrix) = -J
function LA.inv(J::JMatrix)
iszero(J.n %2) ||throw(SingularException)
return -J
end
LA.diag(J::JMatrix{T}) where T = zeros(T, J.n)
LA.tr(J::JMatrix{T}) where T = zero(T)



28 changes: 27 additions & 1 deletion test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ Random.seed!(314159) # use same pseudorandom stream for every test
@test eigvals(A, 0,15) [iλ₁,iλ₂]*im
@test eigvals(A, 1:3) [iλ₁,iλ₂,-iλ₂]*im
@test svdvals(A) [iλ₁,iλ₁,iλ₂,iλ₂]
C=SLA.skewchol(A)
C = SLA.skewchol(A)
@test transpose(C.Rm)*C.Jm*C.RmA[C.Pv,C.Pv]
end

Expand Down Expand Up @@ -363,3 +363,29 @@ end
@test transpose(C.Rm)* C.Jm *C.Rm B[C.Pv, C.Pv]
end
end

@testset "jmatrix.jl" begin
for T in (Int32, Int64, Float32, Float64), n in [2, 20, 153, 200]
A = rand(T,n,n)
J = SLA.JMatrix(T, n)
vec = zeros(T, n - 1)
for i = 1 : 2 : n - 1
vec[i] = -1
end
Jtest = SLA.SkewHermTridiagonal(vec)
@test size(J) == (n, n)
@test size(J, 1) == n
@test Matrix(J) == Matrix(Jtest)
@test A*Jtest A*J
Jtest2 = Matrix(J)
@test Matrix(-J) == -Jtest2
@test Matrix(transpose(J)) == -Jtest2
if iszero(n%2)
B = inv(J)
@test B == -Jtest2
@test B * A Matrix(J) \ A
end

end
end

0 comments on commit 38acc2a

Please sign in to comment.