cholesky renamed

src/cholesky.jl

@@ -1,26 +1,26 @@
 
 struct SkewCholesky{T,R<:UpperTriangular{<:T},J<:JMatrix{<:T},P<:AbstractVector{<:Integer}}
-    Rm::R #Uppertriangular matrix
-    Jm::J # Block diagonal skew-symmetric matrix of type JMatrix
-    Pv::P #Permutation vector
+    R::R #Uppertriangular matrix
+    J::J # Block diagonal skew-symmetric matrix of type JMatrix
+    p::P #Permutation vector
 
-    function SkewCholesky{T,R,J,P}(Rm,Jm,Pv) where {T,R,J,P}
+    function SkewCholesky{T,R,J,P}(Rm,Jm,pv) where {T,R,J,P}
         LA.require_one_based_indexing(Rm)
-        new{T,R,J,P}(Rm,Jm,Pv)
+        new{T,R,J,P}(Rm,Jm,pv)
     end
 end
 
 """
-    SkewCholesky(Rm,Pv)
+    SkewCholesky(R,p)
 
 Construct a `SkewCholesky` structure from the `UpperTriangular`
-matrix `Rm` and the permutation vector `Pv`. A matrix `Jm` of type `JMatrix`
+matrix `R` and the permutation vector `p`. A matrix `J` of type `JMatrix`
 is build calling this function.
-The `SkewCholesky` structure has three arguments: `Rm`,`Jm` and `Pv`.
+The `SkewCholesky` structure has three arguments: `R`,`J` and `p`.
 """
-function SkewCholesky(Rm::UpperTriangular{<:T},Pv::AbstractVector{<:Integer}) where {T<:Real}
-    n = size(Rm, 1)
-    return SkewCholesky{T,typeof(Rm),JMatrix{T,+1},typeof(Pv)}(Rm, JMatrix{T,+1}(n), Pv)
+function SkewCholesky(R::UpperTriangular{<:T},p::AbstractVector{<:Integer}) where {T<:Real}
+    n = size(R, 1)
+    return SkewCholesky{T,typeof(R),JMatrix{T,+1},typeof(p)}(R, JMatrix{T,+1}(n), p)
 
 end
 
@@ -93,9 +93,9 @@ 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 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]`
+`R`,`J`,`p`. `R` is `UpperTriangular`, `J` is a `JMatrix`,
+`p` is an array of integers. Let `S` be the returned structure, then the factorization
+is such that `S.R'*S.J*S.R = A[S.p,S.p]`
 
 This factorization (and the underlying algorithm) is described in from P. Benner et al,
 "[Cholesky-like factorizations of skew-symmetric matrices](https://etna.ricam.oeaw.ac.at/vol.11.2000/pp85-93.dir/pp85-93.pdf)"(2000).

src/tridiag.jl

@@ -198,12 +198,12 @@ Base.conj(M::SkewHermTridiagonal{<:Complex}) = SkewHermTridiagonal(conj.(M.ev),(
 Base.copy(M::SkewHermTridiagonal{<:Real}) = SkewHermTridiagonal(copy(M.ev))
 Base.copy(M::SkewHermTridiagonal{<:Complex}) = SkewHermTridiagonal(copy(M.ev), (M.dvim !==nothing ? copy(M.dvim) : nothing))
 
-function Base.imag(M::SkewHermTridiagonal)
+function Base.imag(M::SkewHermTridiagonal{T}) where T
     if M.dvim !== nothing
         LA.SymTridiagonal(M.dvim, imag.(M.ev))
     else
         n=size(M,1)
-        LA.SymTridiagonal(zeros(eltype(imag(M.ev[1])), n), imag.(M.ev))
+        LA.SymTridiagonal(zeros(real(T), n), imag.(M.ev))
     end
 end
 Base.real(M::SkewHermTridiagonal) = SkewHermTridiagonal(real.(M.ev))
@@ -370,15 +370,13 @@ end
     if n != size(C, 2)
         throw(DimensionMismatch("second dimension of B, $n, doesn't match second dimension of C, $(size(C,2))"))
     end
-
     if m == 0
         return C
     elseif iszero(_add.alpha)
         return LA._rmul_or_fill!(C, _add.beta)
     end
     α = S.dvim
     β = S.ev
-
     if α === nothing
         @inbounds begin
             for j = 1:n
@@ -391,7 +389,7 @@ end
                     β₋, β₀ = β₀, β[i]
                     LA._modify!(_add, β₋*x₋ - adjoint(β₀) * x₊, C, (i, j))
                 end
-                LA._modify!(_add, β[m-1] * x₀ , C, (m, j))
+                LA._modify!(_add, β₀  * x₀ , C, (m, j))
             end
         end
     else
@@ -406,7 +404,7 @@ end
                     β₋, β₀ = β₀, β[i]
                     LA._modify!(_add, β₋*x₋ +α[i]*x₀*1im -adjoint(β₀)*x₊, C, (i, j))
                 end
-                LA._modify!(_add, β[m-1]*x₀+α[m]*x₊*1im , C, (m, j))
+                LA._modify!(_add, β₀*x₀+α[m]*x₊*1im , C, (m, j))
             end
         end
     end
@@ -426,8 +424,8 @@ function LA.dot(x::AbstractVector, S::SkewHermTridiagonal, y::AbstractVector)
     dv = S.dvim
     ev = S.ev
     x₀ = x[1]
-    x₊ = x[2]
-    sub = ev[1]
+    x₊ = (nx > 1 ? x[2] : zero(x[1]))
+    sub = (nx > 1 ? ev[1] : zero(x[1]))
     if dv !== nothing
         r = dot( adjoint(sub)*x₊+complex(zero(dv[1]),-dv[1])*x₀, y[1])
         @inbounds for j in 2:nx-1

test/runtests.jl

@@ -22,7 +22,7 @@ Random.seed!(314159) # use same pseudorandom stream for every test
     @test eigvals(A, 1:3) ≈ [iλ₁,iλ₂,-iλ₂]*im
     @test svdvals(A) ≈ [iλ₁,iλ₁,iλ₂,iλ₂]
     C = SLA.skewchol(A)
-    @test transpose(C.Rm)*C.Jm*C.Rm≈A[C.Pv,C.Pv]
+    @test transpose(C.R)*C.J*C.R≈A[C.p,C.p]
 end
 
 @testset "SkewLinearAlgebra.jl" begin
@@ -81,7 +81,6 @@ end
             @test A[n, n-1] === T(3)
             @test A[n-1, n] === T(-3)
         end
-
         x = rand(T, n)
         y = zeros(T, n)
         mul!(y, A, x, T(2), T(0))
@@ -236,7 +235,7 @@ end
 end
 
 @testset "tridiag.jl" begin
-    for T in (Int32,Float32,Float64,ComplexF32), n in [ 2, 10, 11]
+    for T in (Int32,Float32,Float64,ComplexF32), n in [1, 2, 10, 11]
         if T<:Integer
             C = SLA.skewhermitian(rand(convert(Array{T},-20:20), n, n) * T(2))
         else
@@ -286,7 +285,9 @@ end
         @test Matrix(A)/2 == Matrix(A / 2)
         @test Matrix(A + A) == Matrix(A * 2)
         @test Matrix(A- 2 * A) == Matrix(-A)
-        @test dot(x, A, y) ≈ dot(x, Matrix(A), y)
+        if n>1
+            @test dot(x, A, y) ≈ dot(x, Matrix(A), y)
+        end
         if T<:Complex
             z = rand(T)
             @test A * z ≈ Tridiagonal(A) * z
@@ -306,9 +307,6 @@ end
         @test A * x ≈ B * x
         @test yb * A ≈ yb * B
         @test B * A ≈ A * B ≈ B * B
-        if n>1
-            @test A[1,2] == B[1,2]
-        end
         @test size(A,1) == n
 
         EA = eigen(A)
@@ -327,17 +325,17 @@ end
         for f in (real, imag)
             @test Matrix(f(A)) == f(B)
         end
-
-        A[2,1] = 2
-        @test A[2,1] === T(2) === -A[1,2]'
+        if n > 1
+            A[2,1] = 2
+            @test A[2,1] === T(2) === -A[1,2]'
+        end
     end
-
     B = SLA.SkewHermTridiagonal([3,4,5])
     @test B == [0 -3 0 0; 3 0 -4 0; 0 4 0 -5; 0 0 5 0]
     @test repr("text/plain", B) == "4×4 SkewLinearAlgebra.SkewHermTridiagonal{$Int, Vector{$Int}, Nothing}:\n ⋅  -3   ⋅   ⋅\n 3   ⋅  -4   ⋅\n ⋅   4   ⋅  -5\n ⋅   ⋅   5   ⋅"
     C = SLA.SkewHermTridiagonal(complex.([3,4,5]), [6,7,8,9])
     @test C == [6im -3 0 0; 3 7im -4 0; 0 4 8im -5; 0 0 5 9im]
-    @test repr("text/plain", C) == "4×4 SkewLinearAlgebra.SkewHermTridiagonal{Complex{$Int}, Vector{Complex{$Int}}, Vector{$Int}}:\n 0+6im  -3+0im     ⋅       ⋅  \n 3+0im   0+7im  -4+0im     ⋅  \n   ⋅     4+0im   0+8im  -5+0im\n   ⋅       ⋅     5+0im   0+9im"
+    @test repr("text/plain", C) == "4×4 SkewLinearAlgebra.SkewHermTridiagonal{Complex{$Int}, Vector{Complex{$Int}}, Vector{$Int}}:\n 0+6im  -3+0im     ⋅       ⋅  \n 3+0im   0+7im  -4+0im     ⋅  \n   ⋅     4+0im   0+8im  -5+0im\n   ⋅       ⋅     5+0im   0+9im"    
 end
 
 @testset "pfaffian.jl" begin
@@ -364,10 +362,10 @@ end
             A = SLA.skewhermitian(randn(T, n, n))
         end
         C = SLA.skewchol(A)
-        @test transpose(C.Rm) * C.Jm *C.Rm ≈ A.data[C.Pv, C.Pv]
+        @test transpose(C.R) * C.J *C.R ≈ A.data[C.p, C.p]
         B = Matrix(A)
         C = SLA.skewchol(B)
-        @test transpose(C.Rm)* C.Jm *C.Rm ≈ B[C.Pv, C.Pv]
+        @test transpose(C.R)* C.J *C.R ≈ B[C.p, C.p]
     end
 end