1 // Copyright ©2013 The Gonum Authors. All rights reserved.
2 // Use of this source code is governed by a BSD-style
3 // license that can be found in the LICENSE file.
10 "golang.org/x/exp/rand"
13 func TestLUD(t *testing.T) {
14 for _, n := range []int{1, 5, 10, 11, 50} {
15 a := NewDense(n, n, nil)
16 for i := 0; i < n; i++ {
17 for j := 0; j < n; j++ {
18 a.Set(i, j, rand.NormFloat64())
30 pivot := lu.Pivot(nil)
31 p.Permutation(n, pivot)
34 if !EqualApprox(&got, &want, 1e-12) {
35 t.Errorf("PLU does not equal original matrix.\nWant: %v\n Got: %v", want, got)
40 func TestLURankOne(t *testing.T) {
41 for _, pivoting := range []bool{true} {
42 for _, n := range []int{3, 10, 50} {
43 // Construct a random LU factorization
45 lu.lu = NewDense(n, n, nil)
46 for i := 0; i < n; i++ {
47 for j := 0; j < n; j++ {
48 lu.lu.Set(i, j, rand.Float64())
51 lu.pivot = make([]int, n)
52 for i := range lu.pivot {
56 // For each row, randomly swap with itself or a row after (like is done)
57 // in the actual LU factorization.
58 for i := range lu.pivot {
59 idx := i + rand.Intn(n-i)
60 lu.pivot[i], lu.pivot[idx] = lu.pivot[idx], lu.pivot[i]
63 // Apply a rank one update. Ensure the update magnitude is larger than
64 // the equal tolerance.
65 alpha := rand.Float64() + 1
66 x := NewVecDense(n, nil)
67 y := NewVecDense(n, nil)
68 for i := 0; i < n; i++ {
69 x.setVec(i, rand.Float64()+1)
70 y.setVec(i, rand.Float64()+1)
72 a := luReconstruct(lu)
73 a.RankOne(a, alpha, x, y)
76 luNew.RankOne(lu, alpha, x, y)
77 lu.RankOne(lu, alpha, x, y)
79 aR1New := luReconstruct(&luNew)
80 aR1 := luReconstruct(lu)
82 if !Equal(aR1, aR1New) {
83 t.Error("Different answer when new receiver")
85 if !EqualApprox(aR1, a, 1e-10) {
86 t.Errorf("Rank one mismatch, pivot %v.\nWant: %v\nGot:%v\n", pivoting, a, aR1)
92 // luReconstruct reconstructs the original A matrix from an LU decomposition.
93 func luReconstruct(lu *LU) *Dense {
98 pivot := lu.Pivot(nil)
99 P.Permutation(len(pivot), pivot)
107 func TestSolveLU(t *testing.T) {
108 for _, test := range []struct {
117 a := NewDense(n, n, nil)
118 for i := 0; i < n; i++ {
119 for j := 0; j < n; j++ {
120 a.Set(i, j, rand.NormFloat64())
123 b := NewDense(n, bc, nil)
124 for i := 0; i < n; i++ {
125 for j := 0; j < bc; j++ {
126 b.Set(i, j, rand.NormFloat64())
132 if err := lu.Solve(&x, false, b); err != nil {
137 if !EqualApprox(&got, b, 1e-12) {
138 t.Errorf("Solve mismatch for non-singular matrix. n = %v, bc = %v.\nWant: %v\nGot: %v", n, bc, b, got)
141 // TODO(btracey): Add testOneInput test when such a function exists.
144 func TestSolveLUCond(t *testing.T) {
145 for _, test := range []*Dense{
146 NewDense(2, 2, []float64{1, 0, 0, 1e-20}),
151 b := NewDense(m, 2, nil)
153 if err := lu.Solve(&x, false, b); err == nil {
154 t.Error("No error for near-singular matrix in matrix solve.")
157 bvec := NewVecDense(m, nil)
159 if err := lu.SolveVec(&xvec, false, bvec); err == nil {
160 t.Error("No error for near-singular matrix in matrix solve.")
165 func TestSolveLUVec(t *testing.T) {
166 for _, n := range []int{5, 10} {
167 a := NewDense(n, n, nil)
168 for i := 0; i < n; i++ {
169 for j := 0; j < n; j++ {
170 a.Set(i, j, rand.NormFloat64())
173 b := NewVecDense(n, nil)
174 for i := 0; i < n; i++ {
175 b.SetVec(i, rand.NormFloat64())
180 if err := lu.SolveVec(&x, false, b); err != nil {
185 if !EqualApprox(&got, b, 1e-12) {
186 t.Errorf("Solve mismatch n = %v.\nWant: %v\nGot: %v", n, b, got)
189 // TODO(btracey): Add testOneInput test when such a function exists.