mat: rewrite Cholesky.ToSym without calling SymOuterK

name                    old time/op  new time/op  delta
CholeskyToSym/n=10-4    2.65µs ± 2%  0.78µs ± 0%  -70.75%  (p=0.029 n=4+4)
CholeskyToSym/n=100-4    329µs ± 1%   124µs ± 0%  -62.33%  (p=0.029 n=4+4)
CholeskyToSym/n=1000-4   241ms ± 2%   175ms ± 0%  -27.34%  (p=0.029 n=4+4)

mat/cholesky.go

@@ -287,11 +287,24 @@ func (c *Cholesky) ToSym(dst *SymDense) *SymDense {
 	}
 	n := c.chol.mat.N
 	if dst == nil {
-		dst = NewSymDense(n, make([]float64, n*n))
+		dst = NewSymDense(n, nil)
 	} else {
 		dst.reuseAs(n)
 	}
-	dst.SymOuterK(1, c.chol.T())
+	// Copy c.chol into a TriDense U with dst's backing slice.
+	u := *c.chol
+	u.mat.Data = dst.mat.Data
+	u.Copy(c.chol)
+	// Compute the product U^T*U using the algorithm from LAPACK/TESTING/LIN/dpot01.f
+	bi := blas64.Implementation()
+	a := dst.mat.Data
+	lda := dst.mat.Stride
+	for k := n - 1; k >= 0; k-- {
+		a[k*lda+k] = bi.Ddot(k+1, a[k:], lda, a[k:], lda)
+		if k > 0 {
+			bi.Dtrmv(blas.Upper, blas.Trans, blas.NonUnit, k, a, lda, a[k:], lda)
+		}
+	}
 	return dst
 }