--- /dev/null
+// Copyright ©2015 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package mat_test
+
+import (
+ "fmt"
+
+ "gonum.org/v1/gonum/mat"
+)
+
+func ExampleCholesky() {
+ // Construct a symmetric positive definite matrix.
+ tmp := mat.NewDense(4, 4, []float64{
+ 2, 6, 8, -4,
+ 1, 8, 7, -2,
+ 2, 2, 1, 7,
+ 8, -2, -2, 1,
+ })
+ var a mat.SymDense
+ a.SymOuterK(1, tmp)
+
+ fmt.Printf("a = %0.4v\n", mat.Formatted(&a, mat.Prefix(" ")))
+
+ // Compute the cholesky factorization.
+ var chol mat.Cholesky
+ if ok := chol.Factorize(&a); !ok {
+ fmt.Println("a matrix is not positive semi-definite.")
+ }
+
+ // Find the determinant.
+ fmt.Printf("\nThe determinant of a is %0.4g\n\n", chol.Det())
+
+ // Use the factorization to solve the system of equations a * x = b.
+ b := mat.NewVecDense(4, []float64{1, 2, 3, 4})
+ var x mat.VecDense
+ 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.
+ t := chol.LTo(nil)
+ var test mat.Dense
+ test.Mul(t, t.T())
+ fmt.Println()
+ fmt.Printf("L * L^T = %0.4v\n", mat.Formatted(&a, mat.Prefix(" ")))
+
+ // Output:
+ // a = ⎡120 114 -4 -16⎤
+ // ⎢114 118 11 -24⎥
+ // ⎢ -4 11 58 17⎥
+ // ⎣-16 -24 17 73⎦
+ //
+ // The determinant of a is 1.543e+06
+ //
+ // Solve a * x = b
+ // x = ⎡ -0.239⎤
+ // ⎢ 0.2732⎥
+ // ⎢-0.04681⎥
+ // ⎣ 0.1031⎦
+ //
+ // L * L^T = ⎡120 114 -4 -16⎤
+ // ⎢114 118 11 -24⎥
+ // ⎢ -4 11 58 17⎥
+ // ⎣-16 -24 17 73⎦
+}
+
+func ExampleCholesky_SymRankOne() {
+ a := mat.NewSymDense(4, []float64{
+ 1, 1, 1, 1,
+ 0, 2, 3, 4,
+ 0, 0, 6, 10,
+ 0, 0, 0, 20,
+ })
+ fmt.Printf("A = %0.4v\n", mat.Formatted(a, mat.Prefix(" ")))
+
+ // Compute the Cholesky factorization.
+ var chol mat.Cholesky
+ if ok := chol.Factorize(a); !ok {
+ fmt.Println("matrix a is not positive definite.")
+ }
+
+ x := mat.NewVecDense(4, []float64{0, 0, 0, 1})
+ fmt.Printf("\nx = %0.4v\n", mat.Formatted(x, mat.Prefix(" ")))
+
+ // Rank-1 update the factorization.
+ chol.SymRankOne(&chol, 1, x)
+ // Rank-1 update the matrix a.
+ a.SymRankOne(a, 1, x)
+
+ au := chol.ToSym(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(" ")))
+
+ // Output:
+ // A = ⎡ 1 1 1 1⎤
+ // ⎢ 1 2 3 4⎥
+ // ⎢ 1 3 6 10⎥
+ // ⎣ 1 4 10 20⎦
+ //
+ // x = ⎡0⎤
+ // ⎢0⎥
+ // ⎢0⎥
+ // ⎣1⎦
+ //
+ // A' = ⎡ 1 1 1 1⎤
+ // ⎢ 1 2 3 4⎥
+ // ⎢ 1 3 6 10⎥
+ // ⎣ 1 4 10 21⎦
+ //
+ // U'^T * U' = ⎡ 1 1 1 1⎤
+ // ⎢ 1 2 3 4⎥
+ // ⎢ 1 3 6 10⎥
+ // ⎣ 1 4 10 21⎦
+}