From 18df4e127e0f54baef0ba9ecb0c20efd8e098395 Mon Sep 17 00:00:00 2001
From: Marco Congedo <47817916+Marco-Congedo@users.noreply.github.com>
Date: Tue, 4 Feb 2020 18:03:55 +0100
Subject: [PATCH 01/15] Add files via upload

This unit implements a Gaussian elimination recursion to compute the Cholesky factorization and its inverse in one pass. For small matrices this is faster then using `cholesky` and inversing the Cholesky factor with LinearAlgebra.jl. The main function creates two Cholesky factorizations, inheriting all the methods from LinearAlgebra.
---
 src/choleskyinv.jl | 152 +++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 152 insertions(+)
 create mode 100644 src/choleskyinv.jl

diff --git a/src/choleskyinv.jl b/src/choleskyinv.jl
new file mode 100644
index 0000000..f9d5261
--- /dev/null
+++ b/src/choleskyinv.jl
@@ -0,0 +1,152 @@
+using LinearAlgebra
+
+import Base: show, iterate
+
+"""
+Cholesky factorization and inverse Cholesky factorization, accessible
+in fields `.c` ans `.ci`, respectively.
+"""
+struct Choleskyinv{T<:Factorization}
+    c::T
+	ci::T
+end
+
+# iteration for destructuring into components
+Base.iterate(Choleskyinv) = (C.c, C.ci)
+
+
+"""
+    choleskyinv(P::AbstractMatrix{T};
+		kind::Symbol = :LLt, tol::Real = √eps(T)) where T<:Union{Real, Complex}
+
+ Compute the Cholesky factorization of a dense positive definite
+ matrix P and return a `Choleskyinv` object, holding in field `.c`
+ the Cholesky factorization and in field `ci` the inverse of the Cholesky
+ factorization.
+
+ The two factorizations are obtained in one pass and this is faster
+ then calling Julia's [chelosky](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.cholesky)
+ function and inverting the lower factor for small matrices.
+
+ Input matrix `P` may be of type `Matrix` or `Hermitian`. Since only the
+ lower triangle is used, `P` may also be a `LowerTriangular` view of a
+ positive definite matrix.
+ If `P` is real, it can also be of the `Symmetric` type.
+
+ The algorithm is a *multiplicative Gaussian elimination*.
+
+ **Notes:**
+ The iverse Cholesky factor ``L^{-1}``, obtained as `.ci.L`,
+ is an inverse square root (whitening matrix) of `P`, since
+ ``L^{-1}PL^{-H}=I``. It therefore yields the inversion of ``P`` as
+ ``P^{-1}=L^{-H}L^{-1}``. This is the fastest whitening matrix to be computed,
+ however it yields poor numerical precision, especially for large matrices.
+
+ The following relations holds:
+ - ``L=PL^{-H}``,
+ - ``L^{H}=L^{-1}P``,
+ - ``L^{-H}=P^{-1}L``
+ - ``L^{-1}=L^{H}P^{-1}``.
+
+ We also have ``L^{H}L=L^{-1}P^{2}L^{-H}=UPU^H``, with ``U`` orthogonal
+ (see below) and ``L^{-1}L^{-H}=L^{H}P^{-2}L=UP^{-1}U^H``.
+ ``LL^{H}`` and ``L^{H}L`` are unitarily similar, that is,
+
+ ``ULL^{H}U^H=L^{H}L``,
+
+ where ``U=L^{-1}P^{1/2}``, with ``P^{1/2}=H`` the *principal* (unique symmetric)
+ square root of ``P``. This is seen writing
+ ``PP^{-1}=HHL^{-H}L^{-1}``; multiplying both sides on the left by ``L^{-1}``
+ and on the right by ``L`` we obtain
+ ``L^{-1}PP^{-1}L=L^{-1}HHL^{-H}=I=(L^{-1}H)(L^{-1}H)^H`` and since
+ ``L^{-1}H`` is square it must be unitary.
+
+ From these expressions we have
+ - ``H=LU=U^HL^H,
+ - ``L=HU^H``,
+ - ``H^{-1}=U^HL^{-1}``
+ - ``L^{-1}=UH^{-1}``.
+
+ ``U`` is the *polar factor* of ``L^{H}``, *i.e.*, ``L^{H}=UH``,
+ since ``LL^{H}=HU^HUH^H=H^2=P``.
+ From ``L^{H}L=UCU^H`` we have ``L^{H}LU=UC=ULL^{H}`` and from
+ ``U=L^{-1}H`` we have ``L=HU^H``.
+
+ ## Examples
+	using PosDefManifold, Test
+	n = 40
+	etol = 1e-9
+	Y=randP(n)
+	Yi=inv(Y)
+	a=cholesky(Y)
+
+	C=choleskyinv!(copy(Matrix(Y)))
+	@test(norm(C.c.L*C.c.U-Y)/√n < etol)
+	@test(norm(C.ci.U*C.ci.L-Yi)/√n < etol)
+
+	# repeat the test for complex matrices
+	Y=randP(ComplexF64, n)
+	Yi=inv(Y)
+	C=choleskyinv!(copy(Matrix(Y)))
+	@test(norm(C.c.L*C.c.U-Y)/√n < etol)
+	@test(norm(C.ci.U*C.ci.L-Yi)/√n < etol)
+
+	# Benchmark
+
+	# computing the Cholesky factor and its inverse using LinearAlgebra
+	function linearAlgebraway(P)
+		C=cholesky(P)
+		Q=inv(C.L)
+	end
+
+	Y=randP(n)
+	@benchmark(choleskyinv(Y))
+	@benchmark(linearAlgebraway(Y))
+
+
+"""
+choleskyinv(P::AbstractMatrix{T};
+	   		check::Bool = true,
+	   		tol::Real = √eps(T)) where T<:Union{Real, Complex} =
+    choleskyinv!(P isa Hermitian || P isa Symmetric ? copy(Matrix(P)) : copy(P);
+			check=check, tol=tol)
+
+"""
+    choleskyinv!(P::AbstractMatrix{T};
+		kind::Symbol = :LLt, tol::Real = √eps(T)) where T<:RealOrComplex
+ The same thing as [`choleskyinv`](@ref), but destroys the input matrix.
+ This function does nt require copying the input matrix,
+ thus it is slightly faster.
+"""
+function choleskyinv!(	P::AbstractMatrix{T};
+			  	 		check::Bool = true,
+						tol::Real = √eps(real(T))) where T<:Union{Real, Complex}
+	P isa Matrix || P isa LowerTriangular || throw(ArgumentError("function choleskyinv!: input matrix must be of the Matrix or LowerTriangular type Call `choleskyinv` instead"))
+	#require_one_based_indexing(P)
+	n = LinearAlgebra.checksquare(P)
+	L₁ 	= LowerTriangular(Matrix{T}(I, n, n))
+	U₁⁻¹= UpperTriangular(Matrix{T}(I, n, n))
+
+	@inbounds begin
+		for j=1:n-1
+			check && abs2(P[j, j])<tol && throw(LinearAlgebra.PosDefException(1))
+			for i=j+1:n
+				θ = conj(P[i, j] / -P[j, j])
+				for k=i:n P[k, i] += θ * P[k, j] end # update A and write D
+				L₁[i, j] = conj(-θ) # write Cholesky factor
+				for k=1:j-1 U₁⁻¹[k, i] += θ * U₁⁻¹[k, j] end # write inv Cholesky factor
+				U₁⁻¹[j, i] = θ # write inv Cholesky factor
+			end
+		end
+	end
+
+	D=sqrt.(Diagonal(P))
+	return Choleskyinv(Cholesky(L₁*D, :L, 0), Cholesky(LowerTriangular(Matrix((U₁⁻¹*inv(D))')), :L, 0))
+end
+
+
+function show(io::IO, ::MIME{Symbol("text/plain")}, C::Choleskyinv)
+    println(io, "Cholesky (.c) and inverse Cholesky (.ci) factorizations")
+	#show(io, C.c)
+	#show(io, C.ci)
+end

From a8050fc9db3fb9ea620a690cf0ffbe3256bfbd14 Mon Sep 17 00:00:00 2001
From: Marco Congedo <47817916+Marco-Congedo@users.noreply.github.com>
Date: Tue, 4 Feb 2020 18:04:36 +0100
Subject: [PATCH 02/15] Update choleskyinv.jl

---
 src/choleskyinv.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/choleskyinv.jl b/src/choleskyinv.jl
index f9d5261..3732f84 100644
--- a/src/choleskyinv.jl
+++ b/src/choleskyinv.jl
@@ -122,7 +122,7 @@ function choleskyinv!(	P::AbstractMatrix{T};
 			  	 		check::Bool = true,
 						tol::Real = √eps(real(T))) where T<:Union{Real, Complex}
 	P isa Matrix || P isa LowerTriangular || throw(ArgumentError("function choleskyinv!: input matrix must be of the Matrix or LowerTriangular type Call `choleskyinv` instead"))
-	#require_one_based_indexing(P)
+	require_one_based_indexing(P)
 	n = LinearAlgebra.checksquare(P)
 	L₁ 	= LowerTriangular(Matrix{T}(I, n, n))
 	U₁⁻¹= UpperTriangular(Matrix{T}(I, n, n))

From 3b2b9ec838e0cc26871f37fda6d79d460abdbf54 Mon Sep 17 00:00:00 2001
From: Marco Congedo <47817916+Marco-Congedo@users.noreply.github.com>
Date: Tue, 4 Feb 2020 18:06:03 +0100
Subject: [PATCH 03/15] added include for unit `choleskyinv.jl`

---
 src/MatrixFactorizations.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/MatrixFactorizations.jl b/src/MatrixFactorizations.jl
index b57fd7c..9397023 100644
--- a/src/MatrixFactorizations.jl
+++ b/src/MatrixFactorizations.jl
@@ -35,5 +35,6 @@ export ql, ql!, qrunblocked, qrunblocked!, QL
 
 include("qr.jl")
 include("ql.jl")
+include("choleskyinv.jl")
 
 end #module

From b38095a98a2aa08c6cdde85d83118af41c34e4fe Mon Sep 17 00:00:00 2001
From: Marco Congedo <47817916+Marco-Congedo@users.noreply.github.com>
Date: Tue, 4 Feb 2020 18:09:39 +0100
Subject: [PATCH 04/15] updated docs

---
 src/choleskyinv.jl | 11 +++++------
 1 file changed, 5 insertions(+), 6 deletions(-)

diff --git a/src/choleskyinv.jl b/src/choleskyinv.jl
index 3732f84..9ddc3c1 100644
--- a/src/choleskyinv.jl
+++ b/src/choleskyinv.jl
@@ -17,7 +17,7 @@ Base.iterate(Choleskyinv) = (C.c, C.ci)
 
 """
     choleskyinv(P::AbstractMatrix{T};
-		kind::Symbol = :LLt, tol::Real = √eps(T)) where T<:Union{Real, Complex}
+		check::Bool=true, tol::Real = √eps(T)) where T<:Union{Real, Complex}
 
  Compute the Cholesky factorization of a dense positive definite
  matrix P and return a `Choleskyinv` object, holding in field `.c`
@@ -36,7 +36,7 @@ Base.iterate(Choleskyinv) = (C.c, C.ci)
  The algorithm is a *multiplicative Gaussian elimination*.
 
  **Notes:**
- The iverse Cholesky factor ``L^{-1}``, obtained as `.ci.L`,
+ The inverse Cholesky factor ``L^{-1}``, obtained as `.ci.L`,
  is an inverse square root (whitening matrix) of `P`, since
  ``L^{-1}PL^{-H}=I``. It therefore yields the inversion of ``P`` as
  ``P^{-1}=L^{-H}L^{-1}``. This is the fastest whitening matrix to be computed,
@@ -92,16 +92,17 @@ Base.iterate(Choleskyinv) = (C.c, C.ci)
 	@test(norm(C.ci.U*C.ci.L-Yi)/√n < etol)
 
 	# Benchmark
+	using BenchmarkTools
 
 	# computing the Cholesky factor and its inverse using LinearAlgebra
-	function linearAlgebraway(P)
+	function linearAlgebraWay(P)
 		C=cholesky(P)
 		Q=inv(C.L)
 	end
 
 	Y=randP(n)
 	@benchmark(choleskyinv(Y))
-	@benchmark(linearAlgebraway(Y))
+	@benchmark(linearAlgebraWay(Y))
 
 
 """
@@ -115,8 +116,6 @@ choleskyinv(P::AbstractMatrix{T};
     choleskyinv!(P::AbstractMatrix{T};
 		kind::Symbol = :LLt, tol::Real = √eps(T)) where T<:RealOrComplex
  The same thing as [`choleskyinv`](@ref), but destroys the input matrix.
- This function does nt require copying the input matrix,
- thus it is slightly faster.
 """
 function choleskyinv!(	P::AbstractMatrix{T};
 			  	 		check::Bool = true,

From e88b2708820a24ed1dbcb5525ff6f4c6510cf782 Mon Sep 17 00:00:00 2001
From: Marco Congedo <47817916+Marco-Congedo@users.noreply.github.com>
Date: Wed, 5 Feb 2020 10:50:36 +0100
Subject: [PATCH 05/15] added export of choelskyinv and choleskyinv!

---
 src/MatrixFactorizations.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/MatrixFactorizations.jl b/src/MatrixFactorizations.jl
index 9397023..b7558e6 100644
--- a/src/MatrixFactorizations.jl
+++ b/src/MatrixFactorizations.jl
@@ -30,7 +30,7 @@ else
     import Base: require_one_based_indexing    
 end                            
 
-export ql, ql!, qrunblocked, qrunblocked!, QL                        
+export ql, ql!, qrunblocked, qrunblocked!, QL, choleskyinv!, choleskyinv                        
 
 
 include("qr.jl")

From 412ba0611b8c39ed265948327cdcd929f822d14b Mon Sep 17 00:00:00 2001
From: Marco Congedo <47817916+Marco-Congedo@users.noreply.github.com>
Date: Thu, 6 Feb 2020 11:51:32 +0100
Subject: [PATCH 06/15] Update choleskyinv.jl

---
 src/choleskyinv.jl | 78 ++++++++++++++++++++++------------------------
 1 file changed, 38 insertions(+), 40 deletions(-)

diff --git a/src/choleskyinv.jl b/src/choleskyinv.jl
index 9ddc3c1..91bbbea 100644
--- a/src/choleskyinv.jl
+++ b/src/choleskyinv.jl
@@ -1,32 +1,32 @@
-using LinearAlgebra
+# choleskyinv.jl
+# Cholesky factor and its inverse (complex-conjugate) transposed in one pass.
+# For small matrices this is faster than computing the Cholesky factor
+# using the `cholesky` function in LinearAlgebra and then invert it.
 
-import Base: show, iterate
 
 """
 Cholesky factorization and inverse Cholesky factorization, accessible
-in fields `.c` ans `.ci`, respectively.
+via fields `.c` ans `.ci`, respectively.
 """
-struct Choleskyinv{T<:Factorization}
-    c::T
-	ci::T
+struct CholeskyInv{T, F<:Factorization{T}} <: Factorization{T}
+     c::F
+     ci::F
 end
 
-# iteration for destructuring into components
-Base.iterate(Choleskyinv) = (C.c, C.ci)
-
 
 """
-    choleskyinv(P::AbstractMatrix{T};
-		check::Bool=true, tol::Real = √eps(T)) where T<:Union{Real, Complex}
+choleskyinv(P::Union{Hermitian, Symmetric, Matrix, LowerTriangular};
+		check::Bool = true,
+   		tol::Real = √eps(real(eltype(P))))
 
  Compute the Cholesky factorization of a dense positive definite
  matrix P and return a `Choleskyinv` object, holding in field `.c`
  the Cholesky factorization and in field `ci` the inverse of the Cholesky
  factorization.
 
- The two factorizations are obtained in one pass and this is faster
+ The two factorizations are obtained in one pass and for small matrices this is faster
  then calling Julia's [chelosky](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.cholesky)
- function and inverting the lower factor for small matrices.
+ function and inverting the lower factor.
 
  Input matrix `P` may be of type `Matrix` or `Hermitian`. Since only the
  lower triangle is used, `P` may also be a `LowerTriangular` view of a
@@ -104,47 +104,45 @@ Base.iterate(Choleskyinv) = (C.c, C.ci)
 	@benchmark(choleskyinv(Y))
 	@benchmark(linearAlgebraWay(Y))
 
-
 """
-choleskyinv(P::AbstractMatrix{T};
-	   		check::Bool = true,
-	   		tol::Real = √eps(T)) where T<:Union{Real, Complex} =
-    choleskyinv!(P isa Hermitian || P isa Symmetric ? copy(Matrix(P)) : copy(P);
-			check=check, tol=tol)
+choleskyinv(P::Union{Hermitian, Symmetric, Matrix, LowerTriangular};
+	   	check::Bool = true,
+	   	tol::Real = √eps(real(eltype(P)))) =
+    choleskyinv!(copy(Matrix(P)); check=check, tol=tol)
+
+
 
 """
     choleskyinv!(P::AbstractMatrix{T};
-		kind::Symbol = :LLt, tol::Real = √eps(T)) where T<:RealOrComplex
+		 kind::Symbol = :LLt, tol::Real = √eps(T)) where T<:RealOrComplex
  The same thing as [`choleskyinv`](@ref), but destroys the input matrix.
 """
-function choleskyinv!(	P::AbstractMatrix{T};
-			  	 		check::Bool = true,
-						tol::Real = √eps(real(T))) where T<:Union{Real, Complex}
-	P isa Matrix || P isa LowerTriangular || throw(ArgumentError("function choleskyinv!: input matrix must be of the Matrix or LowerTriangular type Call `choleskyinv` instead"))
-	require_one_based_indexing(P)
+function choleskyinv!(	P::Matrix{T};
+			check::Bool = true,
+			tol::Real = √eps(real(T))) where T<:Union{Real, Complex}
+	LinearAlgebra.require_one_based_indexing(P)
 	n = LinearAlgebra.checksquare(P)
-	L₁ 	= LowerTriangular(Matrix{T}(I, n, n))
-	U₁⁻¹= UpperTriangular(Matrix{T}(I, n, n))
-
-	@inbounds begin
-		for j=1:n-1
-			check && abs2(P[j, j])<tol && throw(LinearAlgebra.PosDefException(1))
-			for i=j+1:n
-				θ = conj(P[i, j] / -P[j, j])
-				for k=i:n P[k, i] += θ * P[k, j] end # update A and write D
-				L₁[i, j] = conj(-θ) # write Cholesky factor
-				for k=1:j-1 U₁⁻¹[k, i] += θ * U₁⁻¹[k, j] end # write inv Cholesky factor
-				U₁⁻¹[j, i] = θ # write inv Cholesky factor
-			end
+	L₁ = LowerTriangular(Matrix{T}(I, n, n))
+	U₁⁻¹ = UpperTriangular(Matrix{T}(I, n, n))
+
+	@inbounds 
+	for j=1:n-1
+		check && abs2(P[j, j])<tol && throw(LinearAlgebra.PosDefException(1))
+		for i=j+1:n
+			θ = conj(P[i, j] / -P[j, j])
+			for k=i:n P[k, i] += θ * P[k, j] end # update A and write D
+			L₁[i, j] = conj(-θ) # write Cholesky factor
+			for k=1:j-1 U₁⁻¹[k, i] += θ * U₁⁻¹[k, j] end # write inv Cholesky factor
+			U₁⁻¹[j, i] = θ # write inv Cholesky factor
 		end
 	end
 
 	D=sqrt.(Diagonal(P))
-	return Choleskyinv(Cholesky(L₁*D, :L, 0), Cholesky(LowerTriangular(Matrix((U₁⁻¹*inv(D))')), :L, 0))
+	return CholeskyInv(Cholesky(L₁*D, :L, 0), Cholesky(LowerTriangular(Matrix((U₁⁻¹*inv(D))')), :L, 0))
 end
 
 
-function show(io::IO, ::MIME{Symbol("text/plain")}, C::Choleskyinv)
+function show(io::IO, ::MIME{Symbol("text/plain")}, C::CholeskyInv)
     println(io, "Cholesky (.c) and inverse Cholesky (.ci) factorizations")
 	#show(io, C.c)
 	#show(io, C.ci)

From bac102a10ff1635382b04f8b6e1d47b72e61db65 Mon Sep 17 00:00:00 2001
From: Marco Congedo <47817916+Marco-Congedo@users.noreply.github.com>
Date: Thu, 6 Feb 2020 17:31:07 +0100
Subject: [PATCH 07/15] Update choleskyinv.jl

---
 src/choleskyinv.jl | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/choleskyinv.jl b/src/choleskyinv.jl
index 91bbbea..45801e1 100644
--- a/src/choleskyinv.jl
+++ b/src/choleskyinv.jl
@@ -13,6 +13,8 @@ struct CholeskyInv{T, F<:Factorization{T}} <: Factorization{T}
      ci::F
 end
 
+# destructuring
+Base.iterate(C::CholeskyInv) = interate((C.c,C.ci))
 
 """
 choleskyinv(P::Union{Hermitian, Symmetric, Matrix, LowerTriangular};

From 84f3ae895a25ae3acd1383911adf9ede5662f27a Mon Sep 17 00:00:00 2001
From: Marco Congedo <47817916+Marco-Congedo@users.noreply.github.com>
Date: Thu, 6 Feb 2020 17:51:45 +0100
Subject: [PATCH 08/15] Update runtests.jl

Added tests for `choleskyinv!` and `choleskyinv`
---
 test/runtests.jl | 29 ++++++++++++++++++++++++++++-
 1 file changed, 28 insertions(+), 1 deletion(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index cb1654e..68902de 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -378,4 +378,31 @@ end
         @test rmul!(copy(c), Q) ≈ c*Matrix(Q)
         @test rmul!(copy(c), Q') ≈ c*Matrix(Q')
     end
-end
\ No newline at end of file
+end
+
+# choleskyinv.jl
+@testset "choleskyinv" begin
+    etol = 1e-9
+    n = 20
+	# real matrices
+    A = randn(n, n)
+    P = A*A'
+    Pi = inv(P)
+    C = choleskyinv!(copy(Matrix(P)))
+	@test(norm(C.c.L*C.c.U-P)/√n < etol)
+	@test(norm(C.ci.U*C.ci.L-Pi)/√n < etol)
+    C = choleskyinv(P)
+	@test(norm(C.c.L*C.c.U-P)/√n < etol)
+	@test(norm(C.ci.U*C.ci.L-Pi)/√n < etol)
+    
+	# repeat the test for complex matrices
+    A=randn(ComplexF64, n, n)
+    P=A*A'
+    Pi=inv(P)
+	C=choleskyinv!(copy(Matrix(P)))
+	@test(norm(C.c.L*C.c.U-P)/√n < etol)
+	@test(norm(C.ci.U*C.ci.L-Pi)/√n < etol)
+    C = choleskyinv(P)
+	@test(norm(C.c.L*C.c.U-P)/√n < etol)
+	@test(norm(C.ci.U*C.ci.L-Pi)/√n < etol)
+end

From 96f408693d6ca5e67bba661e57ad12d6a8f4ea00 Mon Sep 17 00:00:00 2001
From: Marco Congedo <47817916+Marco-Congedo@users.noreply.github.com>
Date: Fri, 7 Feb 2020 11:53:23 +0100
Subject: [PATCH 09/15] Update choleskyinv.jl

---
 src/choleskyinv.jl | 51 +++++++++++++++++++++-------------------------
 1 file changed, 23 insertions(+), 28 deletions(-)

diff --git a/src/choleskyinv.jl b/src/choleskyinv.jl
index 45801e1..5140161 100644
--- a/src/choleskyinv.jl
+++ b/src/choleskyinv.jl
@@ -6,29 +6,25 @@
 
 """
 Cholesky factorization and inverse Cholesky factorization, accessible
-via fields `.c` ans `.ci`, respectively.
+in fields `.c` ans `.ci`, respectively.
 """
 struct CholeskyInv{T, F<:Factorization{T}} <: Factorization{T}
      c::F
      ci::F
 end
 
-# destructuring
-Base.iterate(C::CholeskyInv) = interate((C.c,C.ci))
-
 """
-choleskyinv(P::Union{Hermitian, Symmetric, Matrix, LowerTriangular};
-		check::Bool = true,
-   		tol::Real = √eps(real(eltype(P))))
+    choleskyinv(P::AbstractMatrix{T};
+		check::Bool=true, tol::Real = √eps(T)) where T<:Union{Real, Complex}
 
  Compute the Cholesky factorization of a dense positive definite
  matrix P and return a `Choleskyinv` object, holding in field `.c`
  the Cholesky factorization and in field `ci` the inverse of the Cholesky
  factorization.
 
- The two factorizations are obtained in one pass and for small matrices this is faster
+ The two factorizations are obtained in one pass and this is faster
  then calling Julia's [chelosky](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.cholesky)
- function and inverting the lower factor.
+ function and inverting the lower factor for small matrices.
 
  Input matrix `P` may be of type `Matrix` or `Hermitian`. Since only the
  lower triangle is used, `P` may also be a `LowerTriangular` view of a
@@ -108,34 +104,33 @@ choleskyinv(P::Union{Hermitian, Symmetric, Matrix, LowerTriangular};
 
 """
 choleskyinv(P::Union{Hermitian, Symmetric, Matrix, LowerTriangular};
-	   	check::Bool = true,
-	   	tol::Real = √eps(real(eltype(P)))) =
+	   		check::Bool = true,
+	   		tol::Real = √eps(real(eltype(P)))) =
     choleskyinv!(copy(Matrix(P)); check=check, tol=tol)
 
-
-
 """
     choleskyinv!(P::AbstractMatrix{T};
-		 kind::Symbol = :LLt, tol::Real = √eps(T)) where T<:RealOrComplex
+		kind::Symbol = :LLt, tol::Real = √eps(T)) where T<:RealOrComplex
  The same thing as [`choleskyinv`](@ref), but destroys the input matrix.
 """
 function choleskyinv!(	P::Matrix{T};
-			check::Bool = true,
-			tol::Real = √eps(real(T))) where T<:Union{Real, Complex}
+			  	 		check::Bool = true,
+						tol::Real = √eps(real(T))) where T<:Union{Real, Complex}
 	LinearAlgebra.require_one_based_indexing(P)
 	n = LinearAlgebra.checksquare(P)
-	L₁ = LowerTriangular(Matrix{T}(I, n, n))
-	U₁⁻¹ = UpperTriangular(Matrix{T}(I, n, n))
-
-	@inbounds 
-	for j=1:n-1
-		check && abs2(P[j, j])<tol && throw(LinearAlgebra.PosDefException(1))
-		for i=j+1:n
-			θ = conj(P[i, j] / -P[j, j])
-			for k=i:n P[k, i] += θ * P[k, j] end # update A and write D
-			L₁[i, j] = conj(-θ) # write Cholesky factor
-			for k=1:j-1 U₁⁻¹[k, i] += θ * U₁⁻¹[k, j] end # write inv Cholesky factor
-			U₁⁻¹[j, i] = θ # write inv Cholesky factor
+	L₁ 	= LowerTriangular(Matrix{T}(I, n, n))
+	U₁⁻¹= UpperTriangular(Matrix{T}(I, n, n))
+
+	@inbounds begin
+		for j=1:n-1
+			check && abs2(P[j, j])<tol && throw(LinearAlgebra.PosDefException(1))
+			for i=j+1:n
+				θ = conj(P[i, j] / -P[j, j])
+				for k=i:n P[k, i] += θ * P[k, j] end # update A and write D
+				L₁[i, j] = conj(-θ) # write Cholesky factor
+				for k=1:j-1 U₁⁻¹[k, i] += θ * U₁⁻¹[k, j] end # write inv Cholesky factor
+				U₁⁻¹[j, i] = θ # write inv Cholesky factor
+			end
 		end
 	end
 

From bc39db757ee85a6b557d54c988ab9c05b09790f2 Mon Sep 17 00:00:00 2001
From: Marco Congedo <47817916+Marco-Congedo@users.noreply.github.com>
Date: Fri, 7 Feb 2020 13:22:29 +0100
Subject: [PATCH 10/15] Update MatrixFactorizations.jl

---
 src/MatrixFactorizations.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/MatrixFactorizations.jl b/src/MatrixFactorizations.jl
index b7558e6..0c9702c 100644
--- a/src/MatrixFactorizations.jl
+++ b/src/MatrixFactorizations.jl
@@ -8,7 +8,7 @@ import LinearAlgebra.LAPACK: liblapack, chkuplo, chktrans
 import LinearAlgebra: cholesky, cholesky!, norm, diag, eigvals!, eigvals, eigen!, eigen,
             qr, axpy!, ldiv!, mul!, lu, lu!, ldlt, ldlt!, AbstractTriangular,
             chkstride1, kron, lmul!, rmul!, factorize, StructuredMatrixStyle, logabsdet,
-            AbstractQ, _zeros, _cut_B, _ret_size
+            AbstractQ, _zeros, _cut_B, _ret_size, require_one_based_indexing, checksquare
 
 import Base: getindex, setindex!, *, +, -, ==, <, <=, >,
                 >=, /, ^, \, transpose, showerror, reindex, checkbounds, @propagate_inbounds

From 4a0451e9f8e7b1d2b53b5d9fd1ece2b8a7de0bdf Mon Sep 17 00:00:00 2001
From: Marco Congedo <47817916+Marco-Congedo@users.noreply.github.com>
Date: Fri, 7 Feb 2020 13:46:33 +0100
Subject: [PATCH 11/15] Update runtests.jl

---
 test/runtests.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index 68902de..fef2eb2 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,5 +1,5 @@
 using MatrixFactorizations, LinearAlgebra, Random, Test
-using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, rmul!, lmul!
+using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, rmul!, lmul!, require_one_based_indexing, checksquare
 
 n = 10
 

From bfadd51407c811a242433f95462b2714c1156f6b Mon Sep 17 00:00:00 2001
From: Marco Congedo <47817916+Marco-Congedo@users.noreply.github.com>
Date: Fri, 7 Feb 2020 17:41:07 +0100
Subject: [PATCH 12/15] Improved `choleskyinv` algorithm

Now it does not require any extra memory:
- the ``L_1`` factor overwrites the strictly lower triangle of the input matrix
- its complex-conjugate transpose inverse overwrites the stricly upper triangle of the input matrix
- ``D`` overwrites the diagonal part of the input matrix
---
 src/choleskyinv.jl | 31 ++++++++++++++-----------------
 1 file changed, 14 insertions(+), 17 deletions(-)

diff --git a/src/choleskyinv.jl b/src/choleskyinv.jl
index 5140161..89fb4ae 100644
--- a/src/choleskyinv.jl
+++ b/src/choleskyinv.jl
@@ -113,29 +113,26 @@ choleskyinv(P::Union{Hermitian, Symmetric, Matrix, LowerTriangular};
 		kind::Symbol = :LLt, tol::Real = √eps(T)) where T<:RealOrComplex
  The same thing as [`choleskyinv`](@ref), but destroys the input matrix.
 """
-function choleskyinv!(	P::Matrix{T};
-			  	 		check::Bool = true,
-						tol::Real = √eps(real(T))) where T<:Union{Real, Complex}
+function choleskyinv!(P::Matrix{T};
+			check::Bool = true,
+			tol::Real = √eps(real(T))) where T<:Union{Real, Complex}
 	LinearAlgebra.require_one_based_indexing(P)
 	n = LinearAlgebra.checksquare(P)
-	L₁ 	= LowerTriangular(Matrix{T}(I, n, n))
-	U₁⁻¹= UpperTriangular(Matrix{T}(I, n, n))
-
-	@inbounds begin
-		for j=1:n-1
-			check && abs2(P[j, j])<tol && throw(LinearAlgebra.PosDefException(1))
-			for i=j+1:n
-				θ = conj(P[i, j] / -P[j, j])
-				for k=i:n P[k, i] += θ * P[k, j] end # update A and write D
-				L₁[i, j] = conj(-θ) # write Cholesky factor
-				for k=1:j-1 U₁⁻¹[k, i] += θ * U₁⁻¹[k, j] end # write inv Cholesky factor
-				U₁⁻¹[j, i] = θ # write inv Cholesky factor
-			end
+	
+	@inbounds for j=1:n-1
+		check && abs2(P[j, j])<tol && throw(LinearAlgebra.PosDefException(1))
+		for i=j+1:n
+			θ = conj(P[i, j] / -P[j, j])
+			for k=i:n P[k, i] += θ * P[k, j] end # update P and write D
+			P[i, j] = conj(-θ) # write Cholesky factor
+			for k=1:j-1 P[k, i] += θ * P[k, j] end # write inv Cholesky factor
+			P[j, i] = θ # write inv Cholesky factor
 		end
 	end
 
+	LT, ULT, UUT = LowerTriangular, UnitLowerTriangular, UnitUpperTriangular
 	D=sqrt.(Diagonal(P))
-	return CholeskyInv(Cholesky(L₁*D, :L, 0), Cholesky(LowerTriangular(Matrix((U₁⁻¹*inv(D))')), :L, 0))
+	return CholeskyInv(Cholesky(LT(ULT(P)*D), :L, 0), Cholesky(LT(Matrix((UUT(P)*inv(D))')) , :L, 0 ))
 end
 
 

From 75b67d3c43674880dbc1291e69df6fdd4fbf9606 Mon Sep 17 00:00:00 2001
From: Marco Congedo <47817916+Marco-Congedo@users.noreply.github.com>
Date: Fri, 7 Feb 2020 17:42:00 +0100
Subject: [PATCH 13/15] Creted v 0.3.0

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 465e94b..c26f9f6 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "MatrixFactorizations"
 uuid = "a3b82374-2e81-5b9e-98ce-41277c0e4c87"
-version = "0.2.1"
+version = "0.3.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

From 390143509f77262ba985b831a16f0f69b43c9a60 Mon Sep 17 00:00:00 2001
From: Marco Congedo <47817916+Marco-Congedo@users.noreply.github.com>
Date: Fri, 7 Feb 2020 18:00:19 +0100
Subject: [PATCH 14/15] removed Travis builds for Julia v1.0 and v1.1

---
 .travis.yml | 2 --
 1 file changed, 2 deletions(-)

diff --git a/.travis.yml b/.travis.yml
index cc5829a..da02d42 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -5,8 +5,6 @@ os:
   - osx
   - windows
 julia:
-  - 1.0
-  - 1.1
   - 1.2
   - 1.3
   - nightly

From 1f87a971664aa8203189683394bc0631a70baf4b Mon Sep 17 00:00:00 2001
From: Marco Congedo <47817916+Marco-Congedo@users.noreply.github.com>
Date: Sat, 8 Feb 2020 12:03:59 +0100
Subject: [PATCH 15/15] Not needed anymore since only versions >=1.2 are
 supported

---
 src/MatrixFactorizations.jl | 7 +------
 1 file changed, 1 insertion(+), 6 deletions(-)

diff --git a/src/MatrixFactorizations.jl b/src/MatrixFactorizations.jl
index 0c9702c..577086d 100644
--- a/src/MatrixFactorizations.jl
+++ b/src/MatrixFactorizations.jl
@@ -23,12 +23,7 @@ import Base: convert, size, view, unsafe_indices,
 
 import ArrayLayouts: reflector!, reflectorApply!
 
-if VERSION < v"1.2-"
-    import Base: has_offset_axes
-    require_one_based_indexing(A...) = !has_offset_axes(A...) || throw(ArgumentError("offset arrays are not supported but got an array with index other than 1"))
-else
-    import Base: require_one_based_indexing    
-end                            
+                        
 
 export ql, ql!, qrunblocked, qrunblocked!, QL, choleskyinv!, choleskyinv