From 9d4132a30c4b2d156dee905d04606b6400bbbc2c Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Fri, 4 Jan 2019 14:33:43 +0100
Subject: [PATCH] mat: use lapack64.Potri in Cholesky.InverseTo
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Also,
 - make tricky shadowing conversion from a SymDense to a TriDense more
 explicit in Cholesky.ToSym,
 - add a benchmark for InverseTo:

name                        old time/op  new time/op  delta
CholeskyInverseTo/n=10-4    10.8µs ± 2%   3.2µs ± 1%  -70.16%  (p=0.008 n=5+5)
CholeskyInverseTo/n=100-4   1.07ms ± 2%  0.51ms ± 2%  -52.06%  (p=0.008 n=5+5)
CholeskyInverseTo/n=1000-4   713ms ± 1%   315ms ± 1%  -55.83%  (p=0.008 n=5+5)
---
 mat/cholesky.go      | 51 +++++++++++++++++++++++++++++++-------------
 mat/cholesky_test.go | 28 ++++++++++++++++++++++++
 2 files changed, 64 insertions(+), 15 deletions(-)

diff --git a/mat/cholesky.go b/mat/cholesky.go
index 90e24c2a..03e70b0b 100644
--- a/mat/cholesky.go
+++ b/mat/cholesky.go
@@ -291,10 +291,19 @@ func (c *Cholesky) ToSym(dst *SymDense) *SymDense {
 	} else {
 		dst.reuseAs(n)
 	}
-	// Copy c.chol into a TriDense U with dst's backing slice.
-	// Operations on u are reflected in dst.
-	u := *c.chol
-	u.mat.Data = dst.mat.Data
+	// Create a TriDense representing the Cholesky factor U with dst's
+	// backing slice.
+	// Operations on u are reflected in s.
+	u := &TriDense{
+		mat: blas64.Triangular{
+			Uplo:   blas.Upper,
+			Diag:   blas.NonUnit,
+			N:      n,
+			Data:   dst.mat.Data,
+			Stride: dst.mat.Stride,
+		},
+		cap: n,
+	}
 	u.Copy(c.chol)
 	// Compute the product U^T*U using the algorithm from LAPACK/TESTING/LIN/dpot01.f
 	a := u.mat.Data
@@ -318,18 +327,30 @@ func (c *Cholesky) InverseTo(s *SymDense) error {
 	if !c.valid() {
 		panic(badCholesky)
 	}
-	// TODO(btracey): Replace this code with a direct call to Dpotri when it
-	// is available.
 	s.reuseAs(c.chol.mat.N)
-	// If:
-	//  chol(A) = U^T * U
-	// Then:
-	//  chol(A^-1) = S * S^T
-	// where S = U^-1
-	var t TriDense
-	err := t.InverseTri(c.chol)
-	s.SymOuterK(1, &t)
-	return err
+	// Create a TriDense representing the Cholesky factor U with the backing
+	// slice from s.
+	// Operations on u are reflected in s.
+	u := &TriDense{
+		mat: blas64.Triangular{
+			Uplo:   blas.Upper,
+			Diag:   blas.NonUnit,
+			N:      s.mat.N,
+			Data:   s.mat.Data,
+			Stride: s.mat.Stride,
+		},
+		cap: s.mat.N,
+	}
+	u.Copy(c.chol)
+
+	_, ok := lapack64.Potri(u.mat)
+	if !ok {
+		return Condition(math.Inf(1))
+	}
+	if c.cond > ConditionTolerance {
+		return Condition(c.cond)
+	}
+	return nil
 }
 
 // Scale multiplies the original matrix A by a positive constant using
diff --git a/mat/cholesky_test.go b/mat/cholesky_test.go
index 61623268..c8674543 100644
--- a/mat/cholesky_test.go
+++ b/mat/cholesky_test.go
@@ -616,3 +616,31 @@ func BenchmarkCholeskyToSym(b *testing.B) {
 		})
 	}
 }
+
+func BenchmarkCholeskyInverseTo(b *testing.B) {
+	for _, n := range []int{10, 100, 1000} {
+		b.Run("n="+strconv.Itoa(n), func(b *testing.B) {
+			rnd := rand.New(rand.NewSource(1))
+
+			data := make([]float64, n*n)
+			for i := range data {
+				data[i] = rnd.NormFloat64()
+			}
+			var a SymDense
+			a.SymOuterK(1, NewDense(n, n, data))
+
+			var chol Cholesky
+			ok := chol.Factorize(&a)
+			if !ok {
+				panic("not positive definite")
+			}
+
+			dst := NewSymDense(n, nil)
+
+			b.ResetTimer()
+			for i := 0; i < b.N; i++ {
+				chol.InverseTo(dst)
+			}
+		})
+	}
+}