+++ /dev/null
-// Copyright ©2013 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
-
-import (
- "testing"
-
- "golang.org/x/exp/rand"
-)
-
-func TestLUD(t *testing.T) {
- for _, n := range []int{1, 5, 10, 11, 50} {
- a := NewDense(n, n, nil)
- for i := 0; i < n; i++ {
- for j := 0; j < n; j++ {
- a.Set(i, j, rand.NormFloat64())
- }
- }
- var want Dense
- want.Clone(a)
-
- var lu LU
- lu.Factorize(a)
-
- l := lu.LTo(nil)
- u := lu.UTo(nil)
- var p Dense
- pivot := lu.Pivot(nil)
- p.Permutation(n, pivot)
- var got Dense
- got.Product(&p, l, u)
- if !EqualApprox(&got, &want, 1e-12) {
- t.Errorf("PLU does not equal original matrix.\nWant: %v\n Got: %v", want, got)
- }
- }
-}
-
-func TestLURankOne(t *testing.T) {
- for _, pivoting := range []bool{true} {
- for _, n := range []int{3, 10, 50} {
- // Construct a random LU factorization
- lu := &LU{}
- lu.lu = NewDense(n, n, nil)
- for i := 0; i < n; i++ {
- for j := 0; j < n; j++ {
- lu.lu.Set(i, j, rand.Float64())
- }
- }
- lu.pivot = make([]int, n)
- for i := range lu.pivot {
- lu.pivot[i] = i
- }
- if pivoting {
- // For each row, randomly swap with itself or a row after (like is done)
- // in the actual LU factorization.
- for i := range lu.pivot {
- idx := i + rand.Intn(n-i)
- lu.pivot[i], lu.pivot[idx] = lu.pivot[idx], lu.pivot[i]
- }
- }
- // Apply a rank one update. Ensure the update magnitude is larger than
- // the equal tolerance.
- alpha := rand.Float64() + 1
- x := NewVecDense(n, nil)
- y := NewVecDense(n, nil)
- for i := 0; i < n; i++ {
- x.setVec(i, rand.Float64()+1)
- y.setVec(i, rand.Float64()+1)
- }
- a := luReconstruct(lu)
- a.RankOne(a, alpha, x, y)
-
- var luNew LU
- luNew.RankOne(lu, alpha, x, y)
- lu.RankOne(lu, alpha, x, y)
-
- aR1New := luReconstruct(&luNew)
- aR1 := luReconstruct(lu)
-
- if !Equal(aR1, aR1New) {
- t.Error("Different answer when new receiver")
- }
- if !EqualApprox(aR1, a, 1e-10) {
- t.Errorf("Rank one mismatch, pivot %v.\nWant: %v\nGot:%v\n", pivoting, a, aR1)
- }
- }
- }
-}
-
-// luReconstruct reconstructs the original A matrix from an LU decomposition.
-func luReconstruct(lu *LU) *Dense {
- var L, U TriDense
- lu.LTo(&L)
- lu.UTo(&U)
- var P Dense
- pivot := lu.Pivot(nil)
- P.Permutation(len(pivot), pivot)
-
- var a Dense
- a.Mul(&L, &U)
- a.Mul(&P, &a)
- return &a
-}
-
-func TestSolveLU(t *testing.T) {
- for _, test := range []struct {
- n, bc int
- }{
- {5, 5},
- {5, 10},
- {10, 5},
- } {
- n := test.n
- bc := test.bc
- a := NewDense(n, n, nil)
- for i := 0; i < n; i++ {
- for j := 0; j < n; j++ {
- a.Set(i, j, rand.NormFloat64())
- }
- }
- b := NewDense(n, bc, nil)
- for i := 0; i < n; i++ {
- for j := 0; j < bc; j++ {
- b.Set(i, j, rand.NormFloat64())
- }
- }
- var lu LU
- lu.Factorize(a)
- var x Dense
- if err := lu.Solve(&x, false, b); err != nil {
- continue
- }
- var got Dense
- got.Mul(a, &x)
- if !EqualApprox(&got, b, 1e-12) {
- t.Errorf("Solve mismatch for non-singular matrix. n = %v, bc = %v.\nWant: %v\nGot: %v", n, bc, b, got)
- }
- }
- // TODO(btracey): Add testOneInput test when such a function exists.
-}
-
-func TestSolveLUCond(t *testing.T) {
- for _, test := range []*Dense{
- NewDense(2, 2, []float64{1, 0, 0, 1e-20}),
- } {
- m, _ := test.Dims()
- var lu LU
- lu.Factorize(test)
- b := NewDense(m, 2, nil)
- var x Dense
- if err := lu.Solve(&x, false, b); err == nil {
- t.Error("No error for near-singular matrix in matrix solve.")
- }
-
- bvec := NewVecDense(m, nil)
- var xvec VecDense
- if err := lu.SolveVec(&xvec, false, bvec); err == nil {
- t.Error("No error for near-singular matrix in matrix solve.")
- }
- }
-}
-
-func TestSolveLUVec(t *testing.T) {
- for _, n := range []int{5, 10} {
- a := NewDense(n, n, nil)
- for i := 0; i < n; i++ {
- for j := 0; j < n; j++ {
- a.Set(i, j, rand.NormFloat64())
- }
- }
- b := NewVecDense(n, nil)
- for i := 0; i < n; i++ {
- b.SetVec(i, rand.NormFloat64())
- }
- var lu LU
- lu.Factorize(a)
- var x VecDense
- if err := lu.SolveVec(&x, false, b); err != nil {
- continue
- }
- var got VecDense
- got.MulVec(a, &x)
- if !EqualApprox(&got, b, 1e-12) {
- t.Errorf("Solve mismatch n = %v.\nWant: %v\nGot: %v", n, b, got)
- }
- }
- // TODO(btracey): Add testOneInput test when such a function exists.
-}
OSDN Git Service
