gonum/lapack/testlapack/dgebal.go

// Copyright ©2016 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 testlapack

import (
	"fmt"
	"testing"

	"golang.org/x/exp/rand"

	"gonum.org/v1/gonum/blas"
	"gonum.org/v1/gonum/blas/blas64"
	"gonum.org/v1/gonum/lapack"
)

type Dgebaler interface {
	Dgebal(job lapack.BalanceJob, n int, a []float64, lda int, scale []float64) (int, int)
}

func DgebalTest(t *testing.T, impl Dgebaler) {
	rnd := rand.New(rand.NewSource(1))

	for _, job := range []lapack.BalanceJob{lapack.BalanceNone, lapack.Permute, lapack.Scale, lapack.PermuteScale} {
		for _, n := range []int{0, 1, 2, 3, 4, 5, 6, 10, 18, 31, 53, 100} {
			for _, extra := range []int{0, 11} {
				for cas := 0; cas < 100; cas++ {
					a := unbalancedSparseGeneral(n, n, n+extra, 2*n, rnd)
					testDgebal(t, impl, job, a)
				}
			}
		}
	}
}

func testDgebal(t *testing.T, impl Dgebaler, job lapack.BalanceJob, a blas64.General) {
	const tol = 1e-14

	n := a.Rows
	extra := a.Stride - n

	var scale []float64
	if n > 0 {
		scale = nanSlice(n)
	}

	want := cloneGeneral(a)

	ilo, ihi := impl.Dgebal(job, n, a.Data, a.Stride, scale)

	prefix := fmt.Sprintf("Case job=%c, n=%v, extra=%v", job, n, extra)

	if !generalOutsideAllNaN(a) {
		t.Errorf("%v: out-of-range write to A\n%v", prefix, a.Data)
	}

	if n == 0 {
		if ilo != 0 {
			t.Errorf("%v: unexpected ilo when n=0. Want 0, got %v", prefix, ilo)
		}
		if ihi != -1 {
			t.Errorf("%v: unexpected ihi when n=0. Want -1, got %v", prefix, ihi)
		}
		return
	}

	if job == lapack.BalanceNone {
		if ilo != 0 {
			t.Errorf("%v: unexpected ilo when job=BalanceNone. Want 0, got %v", prefix, ilo)
		}
		if ihi != n-1 {
			t.Errorf("%v: unexpected ihi when job=BalanceNone. Want %v, got %v", prefix, n-1, ihi)
		}
		k := -1
		for i := range scale {
			if scale[i] != 1 {
				k = i
				break
			}
		}
		if k != -1 {
			t.Errorf("%v: unexpected scale[%v] when job=BalanceNone. Want 1, got %v", prefix, k, scale[k])
		}
		if !equalApproxGeneral(a, want, 0) {
			t.Errorf("%v: unexpected modification of A when job=BalanceNone", prefix)
		}
		return
	}

	if ilo < 0 || ihi < ilo || n <= ihi {
		t.Errorf("%v: invalid ordering of ilo=%v and ihi=%v", prefix, ilo, ihi)
	}

	if ilo >= 2 && !isUpperTriangular(blas64.General{Rows: ilo - 1, Cols: ilo - 1, Data: a.Data, Stride: a.Stride}) {
		t.Errorf("%v: T1 is not upper triangular", prefix)
	}
	m := n - ihi - 1 // Order of T2.
	k := ihi + 1
	if m >= 2 && !isUpperTriangular(blas64.General{Rows: m, Cols: m, Data: a.Data[k*a.Stride+k:], Stride: a.Stride}) {
		t.Errorf("%v: T2 is not upper triangular", prefix)
	}

	if job == lapack.Permute || job == lapack.PermuteScale {
		// Check that all rows in [ilo:ihi+1] have at least one nonzero
		// off-diagonal element.
		zeroRow := -1
		for i := ilo; i <= ihi; i++ {
			onlyZeros := true
			for j := ilo; j <= ihi; j++ {
				if i != j && a.Data[i*a.Stride+j] != 0 {
					onlyZeros = false
					break
				}
			}
			if onlyZeros {
				zeroRow = i
				break
			}
		}
		if zeroRow != -1 && ilo != ihi {
			t.Errorf("%v: row %v has only zero off-diagonal elements, ilo=%v, ihi=%v", prefix, zeroRow, ilo, ihi)
		}
		// Check that all columns in [ilo:ihi+1] have at least one nonzero
		// off-diagonal element.
		zeroCol := -1
		for j := ilo; j <= ihi; j++ {
			onlyZeros := true
			for i := ilo; i <= ihi; i++ {
				if i != j && a.Data[i*a.Stride+j] != 0 {
					onlyZeros = false
					break
				}
			}
			if onlyZeros {
				zeroCol = j
				break
			}
		}
		if zeroCol != -1 && ilo != ihi {
			t.Errorf("%v: column %v has only zero off-diagonal elements, ilo=%v, ihi=%v", prefix, zeroCol, ilo, ihi)
		}

		// Create the permutation matrix P.
		p := eye(n, n)
		for j := n - 1; j > ihi; j-- {
			blas64.Swap(blas64.Vector{N: n, Data: p.Data[j:], Inc: p.Stride},
				blas64.Vector{N: n, Data: p.Data[int(scale[j]):], Inc: p.Stride})
		}
		for j := 0; j < ilo; j++ {
			blas64.Swap(blas64.Vector{N: n, Data: p.Data[j:], Inc: p.Stride},
				blas64.Vector{N: n, Data: p.Data[int(scale[j]):], Inc: p.Stride})
		}
		// Compute Pᵀ*A*P and store into want.
		ap := zeros(n, n, n)
		blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, want, p, 0, ap)
		blas64.Gemm(blas.Trans, blas.NoTrans, 1, p, ap, 0, want)
	}
	if job == lapack.Scale || job == lapack.PermuteScale {
		// Modify want by D and D^{-1}.
		d := eye(n, n)
		dinv := eye(n, n)
		for i := ilo; i <= ihi; i++ {
			d.Data[i*d.Stride+i] = scale[i]
			dinv.Data[i*dinv.Stride+i] = 1 / scale[i]
		}
		ad := zeros(n, n, n)
		blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, want, d, 0, ad)
		blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, dinv, ad, 0, want)
	}
	if !equalApproxGeneral(want, a, tol) {
		t.Errorf("%v: unexpected value of A, ilo=%v, ihi=%v", prefix, ilo, ihi)
	}
}