Merge fa212df into e56ddb0

mat/cholesky.go

@@ -158,13 +158,13 @@ func (c *Cholesky) LogDet() float64 {
 	return det
 }
 
-// SolveCholesky finds the matrix m that solves A * m = b where A is represented
-// by the Cholesky decomposition, placing the result in the receiver.
-func (m *Dense) SolveCholesky(chol *Cholesky, b Matrix) error {
-	if !chol.valid() {
+// Solve finds the matrix m that solves A * m = b where A is represented
+// by the Cholesky decomposition, placing the result in m.
+func (c *Cholesky) Solve(m *Dense, b Matrix) error {
+	if !c.valid() {
 		panic(badCholesky)
 	}
-	n := chol.chol.mat.N
+	n := c.chol.mat.N
 	bm, bn := b.Dims()
 	if n != bm {
 		panic(ErrShape)
@@ -174,10 +174,10 @@ func (m *Dense) SolveCholesky(chol *Cholesky, b Matrix) error {
 	if b != m {
 		m.Copy(b)
 	}
-	blas64.Trsm(blas.Left, blas.Trans, 1, chol.chol.mat, m.mat)
-	blas64.Trsm(blas.Left, blas.NoTrans, 1, chol.chol.mat, m.mat)
-	if chol.cond > ConditionTolerance {
-		return Condition(chol.cond)
+	blas64.Trsm(blas.Left, blas.Trans, 1, c.chol.mat, m.mat)
+	blas64.Trsm(blas.Left, blas.NoTrans, 1, c.chol.mat, m.mat)
+	if c.cond > ConditionTolerance {
+		return Condition(c.cond)
 	}
 	return nil
 }
@@ -205,13 +205,13 @@ func (m *Dense) solveTwoChol(a, b *Cholesky) error {
 	return nil
 }
 
-// SolveCholeskyVec finds the vector v that solves A * v = b where A is represented
-// by the Cholesky decomposition, placing the result in the receiver.
-func (v *Vector) SolveCholeskyVec(chol *Cholesky, b *Vector) error {
-	if !chol.valid() {
+// SolveVec finds the vector v that solves A * v = b where A is represented
+// by the Cholesky decomposition, placing the result in v.
+func (c *Cholesky) SolveVec(v, b *Vector) error {
+	if !c.valid() {
 		panic(badCholesky)
 	}
-	n := chol.chol.mat.N
+	n := c.chol.mat.N
 	vn := b.Len()
 	if vn != n {
 		panic(ErrShape)
@@ -223,69 +223,87 @@ func (v *Vector) SolveCholeskyVec(chol *Cholesky, b *Vector) error {
 	if v != b {
 		v.CopyVec(b)
 	}
-	blas64.Trsv(blas.Trans, chol.chol.mat, v.mat)
-	blas64.Trsv(blas.NoTrans, chol.chol.mat, v.mat)
-	if chol.cond > ConditionTolerance {
-		return Condition(chol.cond)
+	blas64.Trsv(blas.Trans, c.chol.mat, v.mat)
+	blas64.Trsv(blas.NoTrans, c.chol.mat, v.mat)
+	if c.cond > ConditionTolerance {
+		return Condition(c.cond)
 	}
 	return nil
 
 }
 
-// UFromCholesky extracts the n×n upper triangular matrix U from a Cholesky
-// decomposition
+// UTo extracts the n×n upper triangular matrix U from a Cholesky
+// decomposition into dst and returns the result. If dst is nil a new
+// TriDense is allocated.
 //  A = U^T * U.
-func (t *TriDense) UFromCholesky(chol *Cholesky) {
-	if !chol.valid() {
+func (c *Cholesky) UTo(dst *TriDense) *TriDense {
+	if !c.valid() {
 		panic(badCholesky)
 	}
-	n := chol.chol.mat.N
-	t.reuseAs(n, Upper)
-	t.Copy(chol.chol)
+	n := c.chol.mat.N
+	if dst == nil {
+		dst = NewTriDense(n, Upper, make([]float64, n*n))
+	} else {
+		dst.reuseAs(n, Upper)
+	}
+	dst.Copy(c.chol)
+	return dst
 }
 
-// LFromCholesky extracts the n×n lower triangular matrix L from a Cholesky
-// decomposition
+// LTo extracts the n×n lower triangular matrix L from a Cholesky
+// decomposition into dst and returns the result. If dst is nil a new
+// TriDense is allocated.
 //  A = L * L^T.
-func (t *TriDense) LFromCholesky(chol *Cholesky) {
-	if !chol.valid() {
+func (c *Cholesky) LTo(dst *TriDense) *TriDense {
+	if !c.valid() {
 		panic(badCholesky)
 	}
-	n := chol.chol.mat.N
-	t.reuseAs(n, Lower)
-	t.Copy(chol.chol.TTri())
+	n := c.chol.mat.N
+	if dst == nil {
+		dst = NewTriDense(n, Lower, make([]float64, n*n))
+	} else {
+		dst.reuseAs(n, Lower)
+	}
+	dst.Copy(c.chol.TTri())
+	return dst
 }
 
-// FromCholesky reconstructs the original positive definite matrix given its
-// Cholesky decomposition.
-func (s *SymDense) FromCholesky(chol *Cholesky) {
-	if !chol.valid() {
+// To reconstructs the original positive definite matrix given its
+// Cholesky decomposition into dst and returns the result. If dst is nil
+// a new SymDense is allocated.
+func (c *Cholesky) To(dst *SymDense) *SymDense {
+	if !c.valid() {
 		panic(badCholesky)
 	}
-	n := chol.chol.mat.N
-	s.reuseAs(n)
-	s.SymOuterK(1, chol.chol.T())
+	n := c.chol.mat.N
+	if dst == nil {
+		dst = NewSymDense(n, make([]float64, n*n))
+	} else {
+		dst.reuseAs(n)
+	}
+	dst.SymOuterK(1, c.chol.T())
+	return dst
 }
 
-// InverseCholesky computes the inverse of the matrix represented by its Cholesky
-// factorization and stores the result into the receiver. If the factorized
+// InverseTo computes the inverse of the matrix represented by its Cholesky
+// factorization and stores the result into s. If the factorized
 // matrix is ill-conditioned, a Condition error will be returned.
 // Note that matrix inversion is numerically unstable, and should generally be
 // avoided where possible, for example by using the Solve routines.
-func (s *SymDense) InverseCholesky(chol *Cholesky) error {
-	if !chol.valid() {
+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(chol.chol.mat.N)
+	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(chol.chol)
+	err := t.InverseTri(c.chol)
 	s.SymOuterK(1, &t)
 	return err
 }

mat/cholesky_example_test.go

@@ -35,17 +35,16 @@ func ExampleCholesky() {
 	// Use the factorization to solve the system of equations a * x = b.
 	b := mat.NewVector(4, []float64{1, 2, 3, 4})
 	var x mat.Vector
-	if err := x.SolveCholeskyVec(&chol, b); err != nil {
+	if err := chol.SolveVec(&x, b); err != nil {
 		fmt.Println("Matrix is near singular: ", err)
 	}
 	fmt.Println("Solve a * x = b")
 	fmt.Printf("x = %0.4v\n", mat.Formatted(&x, mat.Prefix("    ")))
 
 	// Extract the factorization and check that it equals the original matrix.
-	var t mat.TriDense
-	t.LFromCholesky(&chol)
+	t := chol.LTo(nil)
 	var test mat.Dense
-	test.Mul(&t, t.T())
+	test.Mul(t, t.T())
 	fmt.Println()
 	fmt.Printf("L * L^T = %0.4v\n", mat.Formatted(&a, mat.Prefix("          ")))
 
@@ -92,13 +91,12 @@ func ExampleCholesky_SymRankOne() {
 	// Rank-1 update the matrix a.
 	a.SymRankOne(a, 1, x)
 
-	var au mat.SymDense
-	au.FromCholesky(&chol)
+	au := chol.To(nil)
 
 	// Print the matrix that was updated directly.
 	fmt.Printf("\nA' =        %0.4v\n", mat.Formatted(a, mat.Prefix("            ")))
 	// Print the matrix recovered from the factorization.
-	fmt.Printf("\nU'^T * U' = %0.4v\n", mat.Formatted(&au, mat.Prefix("            ")))
+	fmt.Printf("\nU'^T * U' = %0.4v\n", mat.Formatted(au, mat.Prefix("            ")))
 
 	// Output:
 	// A = ⎡ 1   1   1   1⎤

mat/cholesky_test.go

@@ -53,18 +53,16 @@ func TestCholesky(t *testing.T) {
 			if math.Abs(test.cond-chol.cond) > 1e-13 {
 				t.Errorf("Condition number mismatch: Want %v, got %v", test.cond, chol.cond)
 			}
-			var U TriDense
-			U.UFromCholesky(chol)
+			U := chol.UTo(nil)
 			aCopy := DenseCopyOf(test.a)
 			var a Dense
-			a.Mul(U.TTri(), &U)
+			a.Mul(U.TTri(), U)
 			if !EqualApprox(&a, aCopy, 1e-14) {
 				t.Error("unexpected Cholesky factor product")
 			}
 
-			var L TriDense
-			L.LFromCholesky(chol)
-			a.Mul(&L, L.TTri())
+			L := chol.LTo(nil)
+			a.Mul(L, L.TTri())
 			if !EqualApprox(&a, aCopy, 1e-14) {
 				t.Error("unexpected Cholesky factor product")
 			}
@@ -103,7 +101,7 @@ func TestCholeskySolve(t *testing.T) {
 		}
 
 		var x Dense
-		x.SolveCholesky(&chol, test.b)
+		chol.Solve(&x, test.b)
 		if !EqualApprox(&x, test.ans, 1e-12) {
 			t.Error("incorrect Cholesky solve solution")
 		}
@@ -208,7 +206,7 @@ func TestCholeskySolveVec(t *testing.T) {
 		}
 
 		var x Vector
-		x.SolveCholeskyVec(&chol, test.b)
+		chol.SolveVec(&x, test.b)
 		if !EqualApprox(&x, test.ans, 1e-12) {
 			t.Error("incorrect Cholesky solve solution")
 		}
@@ -221,7 +219,7 @@ func TestCholeskySolveVec(t *testing.T) {
 	}
 }
 
-func TestFromCholesky(t *testing.T) {
+func TestCholeskyTo(t *testing.T) {
 	for _, test := range []*SymDense{
 		NewSymDense(3, []float64{
 			53, 59, 37,
@@ -234,11 +232,10 @@ func TestFromCholesky(t *testing.T) {
 		if !ok {
 			t.Fatal("unexpected Cholesky factorization failure: not positive definite")
 		}
-		var s SymDense
-		s.FromCholesky(&chol)
+		s := chol.To(nil)
 
-		if !EqualApprox(&s, test, 1e-12) {
-			t.Errorf("Cholesky reconstruction not equal to original matrix.\nWant:\n% v\nGot:\n% v\n", Formatted(test), Formatted(&s))
+		if !EqualApprox(s, test, 1e-12) {
+			t.Errorf("Cholesky reconstruction not equal to original matrix.\nWant:\n% v\nGot:\n% v\n", Formatted(test), Formatted(s))
 		}
 	}
 }
@@ -279,7 +276,7 @@ func TestCloneCholesky(t *testing.T) {
 	}
 }
 
-func TestInverseCholesky(t *testing.T) {
+func TestCholeskyInverseTo(t *testing.T) {
 	for _, n := range []int{1, 3, 5, 9} {
 		data := make([]float64, n*n)
 		for i := range data {
@@ -295,7 +292,7 @@ func TestInverseCholesky(t *testing.T) {
 		}
 
 		var sInv SymDense
-		sInv.InverseCholesky(&chol)
+		chol.InverseTo(&sInv)
 
 		var ans Dense
 		ans.Mul(&sInv, &s)
@@ -341,7 +338,7 @@ func TestCholeskySymRankOne(t *testing.T) {
 			a.SymRankOne(&a, alpha, x)
 
 			var achol SymDense
-			achol.FromCholesky(&chol)
+			chol.To(&achol)
 			if !EqualApprox(&achol, &a, 1e-13) {
 				t.Errorf("n=%v, alpha=%v: mismatch between updated matrix and from Cholesky:\nupdated:\n%v\nfrom Cholesky:\n%v",
 					n, alpha, Formatted(&a), Formatted(&achol))
@@ -416,7 +413,7 @@ func TestCholeskySymRankOne(t *testing.T) {
 		a.SymRankOne(a, test.alpha, x)
 
 		var achol SymDense
-		achol.FromCholesky(&chol)
+		chol.To(&achol)
 		if !EqualApprox(&achol, a, 1e-13) {
 			t.Errorf("Case %v: mismatch between updated matrix and from Cholesky:\nupdated:\n%v\nfrom Cholesky:\n%v",
 				i, Formatted(a), Formatted(&achol))

mat/lq.go

@@ -121,16 +121,16 @@ func (lq *LQ) QTo(dst *Dense) *Dense {
 	return dst
 }
 
-// SolveLQ finds a minimum-norm solution to a system of linear equations defined
+// Solve finds a minimum-norm solution to a system of linear equations defined
 // by the matrices A and b, where A is an m×n matrix represented in its LQ factorized
 // form. If A is singular or near-singular a Condition error is returned. Please
 // see the documentation for Condition for more information.
 //
 // The minimization problem solved depends on the input parameters.
 //  If trans == false, find the minimum norm solution of A * X = b.
 //  If trans == true, find X such that ||A*X - b||_2 is minimized.
-// The solution matrix, X, is stored in place into the receiver.
-func (m *Dense) SolveLQ(lq *LQ, trans bool, b Matrix) error {
+// The solution matrix, X, is stored in place into m.
+func (lq *LQ) Solve(m *Dense, trans bool, b Matrix) error {
 	r, c := lq.lq.Dims()
 	br, bc := b.Dims()
 
@@ -188,9 +188,9 @@ func (m *Dense) SolveLQ(lq *LQ, trans bool, b Matrix) error {
 	return nil
 }
 
-// SolveLQVec finds a minimum-norm solution to a system of linear equations.
-// Please see Dense.SolveLQ for the full documentation.
-func (v *Vector) SolveLQVec(lq *LQ, trans bool, b *Vector) error {
+// SolveVec finds a minimum-norm solution to a system of linear equations.
+// Please see LQ.Solve for the full documentation.
+func (lq *LQ) SolveVec(v *Vector, trans bool, b *Vector) error {
 	if v != b {
 		v.checkOverlap(b.mat)
 	}
@@ -202,5 +202,5 @@ func (v *Vector) SolveLQVec(lq *LQ, trans bool, b *Vector) error {
 	} else {
 		v.reuseAs(c)
 	}
-	return v.asDense().SolveLQ(lq, trans, b.asDense())
+	return lq.Solve(v.asDense(), trans, b.asDense())
 }

mat/lq_test.go

@@ -77,7 +77,7 @@ func TestSolveLQ(t *testing.T) {
 			var x Dense
 			lq := &LQ{}
 			lq.Factorize(a)
-			x.SolveLQ(lq, trans, b)
+			lq.Solve(&x, trans, b)
 
 			// Test that the normal equations hold.
 			// A^T * A * x = A^T * b if !trans
@@ -130,7 +130,7 @@ func TestSolveLQVec(t *testing.T) {
 			var x Vector
 			lq := &LQ{}
 			lq.Factorize(a)
-			x.SolveLQVec(lq, trans, b)
+			lq.SolveVec(&x, trans, b)
 
 			// Test that the normal equations hold.
 			// A^T * A * x = A^T * b if !trans
@@ -166,13 +166,13 @@ func TestSolveLQCond(t *testing.T) {
 		lq.Factorize(test)
 		b := NewDense(m, 2, nil)
 		var x Dense
-		if err := x.SolveLQ(&lq, false, b); err == nil {
+		if err := lq.Solve(&x, false, b); err == nil {
 			t.Error("No error for near-singular matrix in matrix solve.")
 		}
 
 		bvec := NewVector(m, nil)
 		var xvec Vector
-		if err := xvec.SolveLQVec(&lq, false, bvec); err == nil {
+		if err := lq.SolveVec(&xvec, false, bvec); err == nil {
 			t.Error("No error for near-singular matrix in matrix solve.")
 		}
 	}