Updated for Julia 0.7

Author: andreasnoack

README.md

@@ -28,13 +28,13 @@ julia> using GenericLinearAlgebra
 
 julia> A = randn(1000,1000); A = A'A;
 
-julia> cholfact(A);
+julia> cholesky(A);
 
-julia> @time cholfact(A);
+julia> @time cholesky(A);
   0.013036 seconds (16 allocations: 7.630 MB)
 
-julia> GenericLinearAlgebra.CholeskyModule.cholRecursive!(copy(A), Val{:L});
+julia> GenericLinearAlgebra.cholRecursive!(copy(A), Val{:L});
 
-julia> @time GenericLinearAlgebra.CholeskyModule.cholRecursive!(copy(A), Val{:L});
+julia> @time GenericLinearAlgebra.cholRecursive!(copy(A), Val{:L});
   0.012098 seconds (7.00 k allocations: 7.934 MB)
 ```

src/GenericLinearAlgebra.jl

@@ -2,6 +2,8 @@ __precompile__()
 
 module GenericLinearAlgebra
 
+import LinearAlgebra: mul!, ldiv!
+
 include("juliaBLAS.jl")
 include("lapack.jl")
 include("cholesky.jl")

src/cholesky.jl

@@ -1,56 +1,51 @@
-module CholeskyModule
-
-    using ..JuliaBLAS
-    using Base.LinAlg: A_rdiv_Bc!
-
-    function cholUnblocked!(A::AbstractMatrix{T}, ::Type{Val{:L}}) where T<:Number
-        n = LinAlg.checksquare(A)
-        A[1,1] = sqrt(A[1,1])
-        if n > 1
-            a21 = view(A, 2:n, 1)
-            scale!(a21, inv(real(A[1,1])))
-
-            A22 = view(A, 2:n, 2:n)
-            rankUpdate!(-one(real(T)), a21, Hermitian(A22, :L))
-            cholUnblocked!(A22, Val{:L})
-        end
-        A
+using LinearAlgebra: rdiv!
+
+function cholUnblocked!(A::AbstractMatrix{T}, ::Type{Val{:L}}) where T<:Number
+    n = LinearAlgebra.checksquare(A)
+    A[1,1] = sqrt(A[1,1])
+    if n > 1
+        a21 = view(A, 2:n, 1)
+        rmul!(a21, inv(real(A[1,1])))
+
+        A22 = view(A, 2:n, 2:n)
+        rankUpdate!(-one(real(T)), a21, Hermitian(A22, :L))
+        cholUnblocked!(A22, Val{:L})
     end
-
-    function cholBlocked!(A::AbstractMatrix{T}, ::Type{Val{:L}}, blocksize::Integer) where T<:Number
-        n = LinAlg.checksquare(A)
-        mnb = min(n, blocksize)
-        A11 = view(A, 1:mnb, 1:mnb)
-        cholUnblocked!(A11, Val{:L})
-        if n > blocksize
-            A21 = view(A, (blocksize + 1):n, 1:blocksize)
-            A_rdiv_Bc!(A21, LowerTriangular(A11))
-
-            A22 = view(A, (blocksize + 1):n, (blocksize + 1):n)
-            rankUpdate!(-one(real(T)), A21, Hermitian(A22, :L))
-            cholBlocked!(A22, Val{:L}, blocksize)
-        end
-        A
+    A
+end
+
+function cholBlocked!(A::AbstractMatrix{T}, ::Type{Val{:L}}, blocksize::Integer) where T<:Number
+    n = LinearAlgebra.checksquare(A)
+    mnb = min(n, blocksize)
+    A11 = view(A, 1:mnb, 1:mnb)
+    cholUnblocked!(A11, Val{:L})
+    if n > blocksize
+        A21 = view(A, (blocksize + 1):n, 1:blocksize)
+        rdiv!(A21, LowerTriangular(A11)')
+
+        A22 = view(A, (blocksize + 1):n, (blocksize + 1):n)
+        rankUpdate!(-one(real(T)), A21, Hermitian(A22, :L))
+        cholBlocked!(A22, Val{:L}, blocksize)
     end
+    A
+end
 
-    function cholRecursive!(A::StridedMatrix{T}, ::Type{Val{:L}}, cutoff = 1) where T
-        n = LinAlg.checksquare(A)
-        if n == 1
-            A[1,1] = sqrt(A[1,1])
-        elseif n < cutoff
-            cholUnblocked!(A, Val{:L})
-        else
-            n2 = div(n, 2)
-            A11 = view(A, 1:n2, 1:n2)
-            cholRecursive!(A11, Val{:L})
-            A21 = view(A, n2 + 1:n, 1:n2)
-            A_rdiv_Bc!(A21, LowerTriangular(A11))
-
-            A22 = view(A, n2 + 1:n, n2 + 1:n)
-            rankUpdate!(-real(one(T)), A21, Hermitian(A22, :L))
-            cholRecursive!(A22, Val{:L}, cutoff)
-        end
-        return LowerTriangular(A)
+function cholRecursive!(A::StridedMatrix{T}, ::Type{Val{:L}}, cutoff = 1) where T
+    n = LinearAlgebra.checksquare(A)
+    if n == 1
+        A[1,1] = sqrt(A[1,1])
+    elseif n < cutoff
+        cholUnblocked!(A, Val{:L})
+    else
+        n2 = div(n, 2)
+        A11 = view(A, 1:n2, 1:n2)
+        cholRecursive!(A11, Val{:L})
+        A21 = view(A, n2 + 1:n, 1:n2)
+        rdiv!(A21, LowerTriangular(A11)')
+
+        A22 = view(A, n2 + 1:n, n2 + 1:n)
+        rankUpdate!(-real(one(T)), A21, Hermitian(A22, :L))
+        cholRecursive!(A22, Val{:L}, cutoff)
     end
-
-end
+    return LowerTriangular(A)
+end