From 38acc2ae70f00cccfab87e114556c2bd8134d2db Mon Sep 17 00:00:00 2001 From: Simon Mataigne <61558943+smataigne@users.noreply.github.com> Date: Mon, 22 Aug 2022 11:45:15 -0400 Subject: [PATCH] case n==1 removed from exp.jl and jmatrix added (#84) * 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 --- src/SkewLinearAlgebra.jl | 1 + src/cholesky.jl | 14 +++--- src/exp.jl | 44 ------------------ src/jmatrix.jl | 96 ++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 28 +++++++++++- 5 files changed, 129 insertions(+), 54 deletions(-) create mode 100644 src/jmatrix.jl diff --git a/src/SkewLinearAlgebra.jl b/src/SkewLinearAlgebra.jl index 6437235..ba89ce5 100644 --- a/src/SkewLinearAlgebra.jl +++ b/src/SkewLinearAlgebra.jl @@ -26,6 +26,7 @@ export include("skewhermitian.jl") include("tridiag.jl") +include("jmatrix.jl") include("hessenberg.jl") include("eigen.jl") include("exp.jl") diff --git a/src/cholesky.jl b/src/cholesky.jl index 3fe2a6a..e4b3170 100644 --- a/src/cholesky.jl +++ b/src/cholesky.jl @@ -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} @@ -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 @@ -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]` diff --git a/src/exp.jl b/src/exp.jl index 0b52cd5..8fdbd90 100644 --- a/src/exp.jl +++ b/src/exp.jl @@ -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) @@ -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) @@ -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) @@ -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) @@ -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) @@ -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)) @@ -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) @@ -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)) @@ -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) @@ -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)) diff --git a/src/jmatrix.jl b/src/jmatrix.jl new file mode 100644 index 0000000..9e32563 --- /dev/null +++ b/src/jmatrix.jl @@ -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) + + + \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 6ba6f6e..020299f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -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.Rm≈A[C.Pv,C.Pv] end @@ -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 +