mat/tridiag_test.go

// Copyright ©2021 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 (
	"fmt"
	"reflect"
	"testing"

	"golang.org/x/exp/rand"

	"gonum.org/v1/gonum/lapack/lapack64"
)

func TestNewTridiag(t *testing.T) {
	for i, test := range []struct {
		n         int
		dl, d, du []float64
		panics    bool
		want      *Tridiag
		dense     *Dense
	}{
		{
			n:      1,
			dl:     nil,
			d:      []float64{1.2},
			du:     nil,
			panics: false,
			want: &Tridiag{
				mat: lapack64.Tridiagonal{
					N:  1,
					DL: nil,
					D:  []float64{1.2},
					DU: nil,
				},
			},
			dense: NewDense(1, 1, []float64{1.2}),
		},
		{
			n:      1,
			dl:     []float64{},
			d:      []float64{1.2},
			du:     []float64{},
			panics: false,
			want: &Tridiag{
				mat: lapack64.Tridiagonal{
					N:  1,
					DL: []float64{},
					D:  []float64{1.2},
					DU: []float64{},
				},
			},
			dense: NewDense(1, 1, []float64{1.2}),
		},
		{
			n:      4,
			dl:     []float64{1.2, 2.3, 3.4},
			d:      []float64{4.5, 5.6, 6.7, 7.8},
			du:     []float64{8.9, 9.0, 0.1},
			panics: false,
			want: &Tridiag{
				mat: lapack64.Tridiagonal{
					N:  4,
					DL: []float64{1.2, 2.3, 3.4},
					D:  []float64{4.5, 5.6, 6.7, 7.8},
					DU: []float64{8.9, 9.0, 0.1},
				},
			},
			dense: NewDense(4, 4, []float64{
				4.5, 8.9, 0, 0,
				1.2, 5.6, 9.0, 0,
				0, 2.3, 6.7, 0.1,
				0, 0, 3.4, 7.8,
			}),
		},
		{
			n:      4,
			dl:     nil,
			d:      nil,
			du:     nil,
			panics: false,
			want: &Tridiag{
				mat: lapack64.Tridiagonal{
					N:  4,
					DL: []float64{0, 0, 0},
					D:  []float64{0, 0, 0, 0},
					DU: []float64{0, 0, 0},
				},
			},
			dense: NewDense(4, 4, nil),
		},
		{
			n:      -1,
			panics: true,
		},
		{
			n:      0,
			panics: true,
		},
		{
			n:      1,
			dl:     []float64{1.2},
			d:      nil,
			du:     nil,
			panics: true,
		},
		{
			n:      1,
			dl:     nil,
			d:      []float64{1.2, 2.3},
			du:     nil,
			panics: true,
		},
		{
			n:      1,
			dl:     []float64{},
			d:      nil,
			du:     []float64{},
			panics: true,
		},
		{
			n:      4,
			dl:     []float64{1.2},
			d:      nil,
			du:     nil,
			panics: true,
		},
		{
			n:      4,
			dl:     []float64{1.2, 2.3, 3.4},
			d:      []float64{4.5, 5.6, 6.7, 7.8, 1.2},
			du:     []float64{8.9, 9.0, 0.1},
			panics: true,
		},
	} {
		var a *Tridiag
		panicked, msg := panics(func() {
			a = NewTridiag(test.n, test.dl, test.d, test.du)
		})
		if panicked {
			if !test.panics {
				t.Errorf("Case %d: unexpected panic: %s", i, msg)
			}
			continue
		}
		if test.panics {
			t.Errorf("Case %d: expected panic", i)
			continue
		}

		r, c := a.Dims()
		if r != test.n {
			t.Errorf("Case %d: unexpected number of rows: got=%d want=%d", i, r, test.n)
		}
		if c != test.n {
			t.Errorf("Case %d: unexpected number of columns: got=%d want=%d", i, c, test.n)
		}

		kl, ku := a.Bandwidth()
		if kl != 1 || ku != 1 {
			t.Errorf("Case %d: unexpected bandwidth: got=%d,%d want=1,1", i, kl, ku)
		}

		if !reflect.DeepEqual(a, test.want) {
			t.Errorf("Case %d: unexpected value via reflect: got=%v, want=%v", i, a, test.want)
		}
		if !Equal(a, test.want) {
			t.Errorf("Case %d: unexpected value via mat.Equal: got=%v, want=%v", i, a, test.want)
		}
		if !Equal(a, test.dense) {
			t.Errorf("Case %d: unexpected value via mat.Equal(Tridiag,Dense):\ngot:\n% v\nwant:\n% v", i, Formatted(a), Formatted(test.dense))
		}
	}
}

func TestTridiagAtSet(t *testing.T) {
	t.Parallel()
	for _, n := range []int{1, 2, 3, 4, 7, 10} {
		tri, ref := newTestTridiag(n)

		name := fmt.Sprintf("Case n=%v", n)

		// Check At explicitly with all valid indices.
		for i := 0; i < n; i++ {
			for j := 0; j < n; j++ {
				if tri.At(i, j) != ref.At(i, j) {
					t.Errorf("%v: unexpected value for At(%d,%d): got %v, want %v",
						name, i, j, tri.At(i, j), ref.At(i, j))
				}
			}
		}
		// Check At via a call to Equal.
		if !Equal(tri, ref) {
			t.Errorf("%v: unexpected value:\ngot: % v\nwant:% v",
				name, Formatted(tri, Prefix("     ")), Formatted(ref, Prefix("     ")))
		}

		// Check At out of bounds.
		for _, i := range []int{-1, n, n + 1} {
			for j := 0; j < n; j++ {
				panicked, message := panics(func() { tri.At(i, j) })
				if !panicked || message != ErrRowAccess.Error() {
					t.Errorf("%v: expected panic for invalid row access at (%d,%d)", name, i, j)
				}
			}
		}
		for _, j := range []int{-1, n, n + 1} {
			for i := 0; i < n; i++ {
				panicked, message := panics(func() { tri.At(i, j) })
				if !panicked || message != ErrColAccess.Error() {
					t.Errorf("%v: expected panic for invalid column access at (%d,%d)", name, i, j)
				}
			}
		}

		// Check SetBand out of bounds.
		for _, i := range []int{-1, n, n + 1} {
			for j := 0; j < n; j++ {
				panicked, message := panics(func() { tri.SetBand(i, j, 1.2) })
				if !panicked || message != ErrRowAccess.Error() {
					t.Errorf("%v: expected panic for invalid row access at (%d,%d)", name, i, j)
				}
			}
		}
		for _, j := range []int{-1, n, n + 1} {
			for i := 0; i < n; i++ {
				panicked, message := panics(func() { tri.SetBand(i, j, 1.2) })
				if !panicked || message != ErrColAccess.Error() {
					t.Errorf("%v: expected panic for invalid column access at (%d,%d)", name, i, j)
				}
			}
		}
		for i := 0; i < n; i++ {
			for j := 0; j <= i-2; j++ {
				panicked, message := panics(func() { tri.SetBand(i, j, 1.2) })
				if !panicked || message != ErrBandSet.Error() {
					t.Errorf("%v: expected panic for invalid access at (%d,%d)", name, i, j)
				}
			}
			for j := i + 2; j < n; j++ {
				panicked, message := panics(func() { tri.SetBand(i, j, 1.2) })
				if !panicked || message != ErrBandSet.Error() {
					t.Errorf("%v: expected panic for invalid access at (%d,%d)", name, i, j)
				}
			}
		}

		// Check SetBand within bandwidth.
		for i := 0; i < n; i++ {
			for j := max(0, i-1); j <= min(i+1, n-1); j++ {
				want := float64(i*n + j + 100)
				tri.SetBand(i, j, want)
				if got := tri.At(i, j); got != want {
					t.Errorf("%v: unexpected value at (%d,%d) after SetBand: got %v, want %v", name, i, j, got, want)
				}
			}
		}
	}
}

func newTestTridiag(n int) (*Tridiag, *Dense) {
	var dl, d, du []float64
	d = make([]float64, n)
	if n > 1 {
		dl = make([]float64, n-1)
		du = make([]float64, n-1)
	}
	for i := range d {
		d[i] = float64(i*n + i + 1)
	}
	for j := range dl {
		i := j + 1
		dl[j] = float64(i*n + j + 1)
	}
	for i := range du {
		j := i + 1
		du[i] = float64(i*n + j + 1)
	}
	dense := make([]float64, n*n)
	for i := 0; i < n; i++ {
		for j := max(0, i-1); j <= min(i+1, n-1); j++ {
			dense[i*n+j] = float64(i*n + j + 1)
		}
	}
	return NewTridiag(n, dl, d, du), NewDense(n, n, dense)
}

func TestTridiagReset(t *testing.T) {
	t.Parallel()
	for _, n := range []int{1, 2, 3, 4, 7, 10} {
		a, _ := newTestTridiag(n)
		if a.IsEmpty() {
			t.Errorf("Case n=%d: matrix is empty", n)
		}
		a.Reset()
		if !a.IsEmpty() {
			t.Errorf("Case n=%d: matrix is not empty after Reset", n)
		}
	}
}

func TestTridiagDiagView(t *testing.T) {
	t.Parallel()
	for _, n := range []int{1, 2, 3, 4, 7, 10} {
		a, _ := newTestTridiag(n)
		testDiagView(t, n, a)
	}
}

func TestTridiagZero(t *testing.T) {
	t.Parallel()
	for _, n := range []int{1, 2, 3, 4, 7, 10} {
		a, _ := newTestTridiag(n)
		a.Zero()
		for i := 0; i < n; i++ {
			for j := 0; j < n; j++ {
				if a.At(i, j) != 0 {
					t.Errorf("Case n=%d: unexpected non-zero at (%d,%d): got %f", n, i, j, a.At(i, j))
				}
			}
		}
	}
}

func TestTridiagSolveTo(t *testing.T) {
	t.Parallel()

	const tol = 1e-13

	rnd := rand.New(rand.NewSource(1))
	random := func(n int) []float64 {
		d := make([]float64, n)
		for i := range d {
			d[i] = rnd.NormFloat64()
		}
		return d
	}

	for _, n := range []int{1, 2, 3, 4, 7, 10} {
		a := NewTridiag(n, random(n-1), random(n), random(n-1))
		var aDense Dense
		aDense.CloneFrom(a)
		for _, trans := range []bool{false, true} {
			for _, nrhs := range []int{1, 2, 5} {
				const (
					denseB = iota
					rawB
					basicB
				)
				for _, bType := range []int{denseB, rawB, basicB} {
					const (
						emptyDst = iota
						shapedDst
						bIsDst
					)
					for _, dstType := range []int{emptyDst, shapedDst, bIsDst} {
						if dstType == bIsDst && bType != denseB {
							continue
						}

						var b Matrix
						switch bType {
						case denseB:
							b = NewDense(n, nrhs, random(n*nrhs))
						case rawB:
							b = &rawMatrix{asBasicMatrix(NewDense(n, nrhs, random(n*nrhs)))}
						case basicB:
							b = asBasicMatrix(NewDense(n, nrhs, random(n*nrhs)))
						default:
							panic("bad bType")
						}

						var dst *Dense
						switch dstType {
						case emptyDst:
							dst = new(Dense)
						case shapedDst:
							dst = NewDense(n, nrhs, random(n*nrhs))
						case bIsDst:
							dst = b.(*Dense)
						default:
							panic("bad dstType")
						}

						name := fmt.Sprintf("n=%d,nrhs=%d,trans=%t,dstType=%d,bType=%d", n, nrhs, trans, dstType, bType)

						var want Dense
						var err error
						if !trans {
							err = want.Solve(&aDense, b)
						} else {
							err = want.Solve(aDense.T(), b)
						}
						if err != nil {
							t.Fatalf("%v: unexpected failure when computing reference solution: %v", name, err)
						}

						err = a.SolveTo(dst, trans, b)
						if err != nil {
							t.Fatalf("%v: unexpected failure from Tridiag.SolveTo: %v", name, err)
						}

						var diff Dense
						diff.Sub(dst, &want)
						if resid := Norm(&diff, 1); resid > tol*float64(n) {
							t.Errorf("%v: unexpected result; resid=%v, want<=%v", name, resid, tol*float64(n))
						}
					}
				}
			}
		}
	}
}

func TestTridiagSolveVecTo(t *testing.T) {
	t.Parallel()

	const tol = 1e-13

	rnd := rand.New(rand.NewSource(1))
	random := func(n int) []float64 {
		d := make([]float64, n)
		for i := range d {
			d[i] = rnd.NormFloat64()
		}
		return d
	}

	for _, n := range []int{1, 2, 3, 4, 7, 10} {
		a := NewTridiag(n, random(n-1), random(n), random(n-1))
		var aDense Dense
		aDense.CloneFrom(a)
		for _, trans := range []bool{false, true} {
			const (
				denseB = iota
				rawB
				basicB
			)
			for _, bType := range []int{denseB, rawB, basicB} {
				const (
					emptyDst = iota
					shapedDst
					bIsDst
				)
				for _, dstType := range []int{emptyDst, shapedDst, bIsDst} {
					if dstType == bIsDst && bType != denseB {
						continue
					}

					var b Vector
					switch bType {
					case denseB:
						b = NewVecDense(n, random(n))
					case rawB:
						b = &rawVector{asBasicVector(NewVecDense(n, random(n)))}
					case basicB:
						b = asBasicVector(NewVecDense(n, random(n)))
					default:
						panic("bad bType")
					}

					var dst *VecDense
					switch dstType {
					case emptyDst:
						dst = new(VecDense)
					case shapedDst:
						dst = NewVecDense(n, random(n))
					case bIsDst:
						dst = b.(*VecDense)
					default:
						panic("bad dstType")
					}

					name := fmt.Sprintf("n=%d,trans=%t,dstType=%d,bType=%d", n, trans, dstType, bType)

					var want VecDense
					var err error
					if !trans {
						err = want.SolveVec(&aDense, b)
					} else {
						err = want.SolveVec(aDense.T(), b)
					}
					if err != nil {
						t.Fatalf("%v: unexpected failure when computing reference solution: %v", name, err)
					}

					err = a.SolveVecTo(dst, trans, b)
					if err != nil {
						t.Fatalf("%v: unexpected failure from Tridiag.SolveTo: %v", name, err)
					}

					var diff Dense
					diff.Sub(dst, &want)
					if resid := Norm(&diff, 1); resid > tol*float64(n) {
						t.Errorf("%v: unexpected result; resid=%v, want<=%v", name, resid, tol*float64(n))
					}
				}
			}
		}
	}
}