mat: use lapack64.Potri in Cholesky.InverseTo

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

@@ -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

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)
+			}
+		})
+	}
+}