From 09b8684471f991270432bc8ec8e0285b1361e926 Mon Sep 17 00:00:00 2001
From: Simon Mataigne <61558943+smataigne@users.noreply.github.com>
Date: Mon, 15 Aug 2022 15:08:19 -0400
Subject: [PATCH] cholesky.jl improved and fixed (#52)

* 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
---
 src/SkewLinearAlgebra.jl |   6 +-
 src/cholesky.jl          | 143 +++++++++++++++++++--------------------
 test/runtests.jl         |  26 ++++---
 3 files changed, 89 insertions(+), 86 deletions(-)

diff --git a/src/SkewLinearAlgebra.jl b/src/SkewLinearAlgebra.jl
index f5ad466..18b063d 100644
--- a/src/SkewLinearAlgebra.jl
+++ b/src/SkewLinearAlgebra.jl
@@ -11,18 +11,22 @@ export
     #Types
     SkewHermitian,
     SkewHermTridiagonal,
+    SkewCholesky,
     #functions
     isskewhermitian,
     skewhermitian,
     skewhermitian!,
     pfaffian,
-    pfaffian!
+    pfaffian!,
+    skewchol,
+    skewchol!
 
 include("skewhermitian.jl")
 include("tridiag.jl")
 include("hessenberg.jl")
 include("eigen.jl")
 include("exp.jl")
+include("cholesky.jl")
 include("pfaffian.jl")
 
 end
diff --git a/src/cholesky.jl b/src/cholesky.jl
index 992b783..1823e5a 100644
--- a/src/cholesky.jl
+++ b/src/cholesky.jl
@@ -1,62 +1,68 @@
-#import .SkewLinearAlgebra as SLA
-@views function skewchol!(A::SLA.SkewHermitian)
-    B = A.data
-    tol = 1e-15
+
+struct SkewCholesky{T,R<:UpperTriangular{<:T},J<:SkewHermTridiagonal{<:T},P<:AbstractVector{<:Integer}} 
+    Rm::R #Uppertriangular matrix
+    Jm::J # Block diagonal skew-symmetric matrix
+    Pv::P #Permutation vector
+
+    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)
+    end
+end
+#SkewCholesky(Rm::UpperTriangular{<:T},Jm::SkewHermTridiagonal{<:T},Pv::AbstractVector{<:Integer}) where {T<:Real} = SkewCholesky(Rm,Jm,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)
+
+end
+
+function _skewchol!(A::SkewHermitian)
+    @views B = A.data
+    tol = 1e-15*norm(B)
     m = size(B,1)
-    J2 = [0 1;-1 0]
+    J2 = similar(B,2,2)
+    J2[1,1] = 0; J2[2,1] = -1; J2[1,2] = 1; J2[2,2] = 0
     ii = 0; jj = 0; kk = 0
     P = Array(1:m)
-    temp = similar(B,m)
     tempM = similar(B,2,m-2)
     for j = 1:m÷2
         j2 = 2*j
-        M = maximum(B[j2-1:m,j2-1:m])
-        for i1 = j2-1:m
-            for i2 = j2-1:m
-                if B[i1,i2] == M
-                    ii = i1
-                    jj = i2
-                end
-            end
-        end
+        
+        M = findmax(B[j2-1:m,j2-1:m])
+        ii = M[2][1] + j2 - 2
+        jj = M[2][2] + j2 - 2
         if abs(B[ii,jj])<tol
             rank = j2-2
-            return P
+            return P, rank
         end
         if jj == j2-1
-            kk=ii
+            kk = ii
         else
             kk = jj
         end
         if ii != j2-1
-
-            I = Array(1:m)
-            I[ii] = j2-1
-            I[j2-1] = ii
-            
-            temp2 = P[ii]
-            P[ii] = P[j2-1]
-            P[j2-1] = temp2
+            P[ii],P[j2-1] = P[j2-1],P[ii]
+            for t = 1:m 
+                B[t,ii], B[t,j2-1] = B[t,j2-1], B[t,ii]
+            end
+            for t = 1:m 
+                B[ii,t], B[j2-1,t] = B[j2-1,t], B[ii,t]
+            end
             
-            Base.permutecols!!(B,I)
-            temp .= B[ii,:]
-            B[ii,:] .= B[j2-1,:]
-            B[j2-1,:] .= temp
         end
         if kk != j2
-
-            I = Array(1:m)
-            I[kk] = j2
-            I[j2] = kk
-            
-            temp3 = P[kk]
-            P[kk] = P[j2]
-            P[j2] = temp3
-            
-            Base.permutecols!!(B,I)
-            temp .= B[kk,:]
-            B[kk,:] .= B[j2,:]
-            B[j2,:] .= temp
+            P[kk],P[j2] = P[j2],P[kk]
+            for t = 1:m 
+                B[t,kk], B[t,j2] = B[t,j2], B[t,kk]
+            end
+            for t = 1:m 
+                B[kk,t], B[j2,t] = B[j2,t], B[kk,t]
+            end
         end
         
         l = m-j2
@@ -64,42 +70,29 @@
         B[j2-1,j2-1] = r
         B[j2,j2] = r
         B[j2-1,j2] = 0
-        mul!(tempM[:,1:l],J2,B[j2-1:j2,j2+1:m])
+        @views mul!(tempM[:,1:l], J2, B[j2-1:j2,j2+1:m])
         B[j2-1:j2,j2+1:m] .= tempM[:,1:l]
         B[j2-1:j2,j2+1:m] .*= (-1/r)
-        mul!(tempM[:,1:l],J2,B[j2-1:j2,j2+1:m])
-        mul!(B[j2+1:m,j2+1:m],transpose(B[j2-1:j2,j2+1:m]),tempM[:,1:l],-1,1)
+        @views mul!(tempM[:,1:l], J2, B[j2-1:j2,j2+1:m])
+        @views mul!(B[j2+1:m,j2+1:m], transpose(B[j2-1:j2,j2+1:m]), tempM[:,1:l],-1,1)
+        
     end
-    r=2*(m÷2)
-    return P
+    rank=2*(m÷2)
+    return P,rank
+end
+copyeigtype(A::AbstractMatrix) = copyto!(similar(A, LA.eigtype(eltype(A))), A)
+@views function skewchol!(A::SkewHermitian)
+    P = _skewchol!(A)[1]
+    return SkewCholesky(UpperTriangular(A.data),P)
 end
 
-@views function skewsolve(A::SkewHermitian,b::AbstractVector)
-    
-    n = size(A,1)
-    P = skewchol!(A)
-    y1 = similar(A.data,n)
-    y2 = similar(A.data,n)
-    R = UpperTriangular(A.data)
-    vec = zeros(n-1)
-    for i = 1:2:n-1
-        vec[i] = 1
-    end
-    Jt = Tridiagonal(vec,zeros(n),-vec)
-    Base.permute!(b,P)
-    y1 .= transpose(R)\b
-    mul!(y2,Jt,y1)
-    y1 .= R\y2
-    Base.permute!(y1,P)
-    return y1
+skewchol(A::SkewHermitian) = skewchol!(copyeigtype(A))
+skewchol!(A::AbstractMatrix) = @views  skewchol!(SkewHermitian(A))
+
+function skewchol(A::AbstractMatrix)
+    isskewhermitian(A) || throw(ArgumentError("Pfaffian requires a skew-Hermitian matrix"))
+    return skewchol!(SkewHermitian(copyeigtype(A)))
 end
-#=
-A=SLA.skewhermitian(A)
-P=skewchol!(A)
-R=triu(A)
-dipslay(transpose(R)*)
-=#
-"""
-Q^TR^TJRQx=b
-=>RQx=J^T R^(-T) Q^Tb
-"""
+
+
+
diff --git a/test/runtests.jl b/test/runtests.jl
index b44766b..04fd70f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -265,16 +265,6 @@ end
     end
 end
 
-
-@testset "pfaffian.jl" begin
-    for n in [2,3,4,5,6,8,10,20,40]
-        A=SLA.skewhermitian(rand(-10:10,n,n)*2)
-        Abig = BigInt.(A.data)
-        @test SLA.pfaffian(A) ≈ SLA.pfaffian(Abig)  == SLA.pfaffian(SLA.SkewHermitian(Abig))
-        #@test SLA.pfaffian(Abig)^2 ≈ det(Abig) #ideally compare with == if recent Julia version
-    end
-
-end
 @testset "pfaffian.jl" begin
     for n in [2,3,4,5,6,8,10,20,40]
         A=SLA.skewhermitian(rand(-10:10,n,n)*2)
@@ -289,6 +279,22 @@ end
     @test SLA.pfaffian(big.([0 14 7 -10 0 10 0 -11; -14 0 -10 7 13 -9 -12 -13; -7 10 0 -4 6 -17 -1 18; 10 -7 4 0 -2 -4 0 11; 0 -13 -6 2 0 -8 -18 17; -10 9 17 4 8 0 -8 12; 0 12 1 0 18 8 0 0; 11 13 -18 -11 -17 -12 0 0])) == -119000
 end
 
+@testset "cholesky.jl" begin
+    for T in (Int32, Int64, Float32, Float64), n in [2,4,20,100]
+        if T<:Integer
+            A = SLA.skewhermitian(rand(convert(Array{T},-10:10),n,n)*2)
+        else
+            A = SLA.skewhermitian(randn(T,n,n))
+        end
+        C = SLA.skewchol(A)
+        @test transpose(C.Rm)*Matrix(C.Jm)*C.Rm ≈A.data[C.Pv,C.Pv]
+        B = Matrix(A)
+        C = SLA.skewchol(B)
+        @test transpose(C.Rm)*Matrix(C.Jm)*C.Rm ≈B[C.Pv,C.Pv]
+    end
+end
+
+
 
 
 #=