dtrcon.go

// 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 testlapack

import (
	"fmt"
	"testing"

	"golang.org/x/exp/rand"

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

type Dtrconer interface {
	Dtrcon(norm lapack.MatrixNorm, uplo blas.Uplo, diag blas.Diag, n int, a []float64, lda int, work []float64, iwork []int) float64

	Dtrtri(uplo blas.Uplo, diag blas.Diag, n int, a []float64, lda int) bool
	Dlantr(norm lapack.MatrixNorm, uplo blas.Uplo, diag blas.Diag, m, n int, a []float64, lda int, work []float64) float64
}

func DtrconTest(t *testing.T, impl Dtrconer) {
	rnd := rand.New(rand.NewSource(1))
	for _, n := range []int{0, 1, 2, 3, 4, 5, 10, 50} {
		for _, uplo := range []blas.Uplo{blas.Lower, blas.Upper} {
			for _, diag := range []blas.Diag{blas.NonUnit, blas.Unit} {
				for _, lda := range []int{max(1, n), n + 3} {
					for _, mattype := range []int{0, 1, 2} {
						dtrconTest(t, impl, rnd, uplo, diag, n, lda, mattype)
					}
				}
			}
		}
	}
}

func dtrconTest(t *testing.T, impl Dtrconer, rnd *rand.Rand, uplo blas.Uplo, diag blas.Diag, n, lda, mattype int) {
	const ratioThresh = 10

	a := make([]float64, max(0, (n-1)*lda+n))
	for i := range a {
		a[i] = rnd.Float64()
	}
	switch mattype {
	default:
		panic("bad mattype")
	case 0:
		// Matrix filled with consecutive integer values.
		// For lapack.MaxRowSum norm (infinity-norm) these matrices
		// sometimes lead to a slightly inaccurate estimate of the condition
		// number.
		c := 2.0
		for i := 0; i < n; i++ {
			for j := 0; j < n; j++ {
				a[i*lda+j] = c
				c += 1
			}
		}
	case 1:
		// Identity matrix.
		if uplo == blas.Upper {
			for i := 0; i < n; i++ {
				for j := i + 1; j < n; j++ {
					a[i*lda+j] = 0
				}
			}
		} else {
			for i := 0; i < n; i++ {
				for j := 0; j < i; j++ {
					a[i*lda+j] = 0
				}
			}
		}
		if diag == blas.NonUnit {
			for i := 0; i < n; i++ {
				a[i*lda+i] = 1
			}
		}
	case 2:
		// Matrix filled with random values uniformly in [-1,1).
		// These matrices often lead to a slightly inaccurate estimate
		// of the condition number.
		for i := 0; i < n; i++ {
			for j := 0; j < n; j++ {
				a[i*lda+j] = 2*rnd.Float64() - 1
			}
		}
	}
	aCopy := make([]float64, len(a))
	copy(aCopy, a)

	// Compute the inverse A^{-1}.
	aInv := make([]float64, len(a))
	copy(aInv, a)
	ok := impl.Dtrtri(uplo, diag, n, aInv, lda)
	if !ok {
		t.Fatalf("uplo=%v,diag=%v,n=%v,lda=%v,mattype=%v: bad matrix, Dtrtri failed", string(uplo), string(diag), n, lda, mattype)
	}

	work := make([]float64, 3*n)
	iwork := make([]int, n)
	for _, norm := range []lapack.MatrixNorm{lapack.MaxColumnSum, lapack.MaxRowSum} {
		name := fmt.Sprintf("norm=%v,uplo=%v,diag=%v,n=%v,lda=%v,mattype=%v", string(norm), string(uplo), string(diag), n, lda, mattype)

		// Compute the norm of A and A^{-1}.
		aNorm := impl.Dlantr(norm, uplo, diag, n, n, a, lda, work)
		aInvNorm := impl.Dlantr(norm, uplo, diag, n, n, aInv, lda, work)

		// Compute a good estimate of the condition number
		//  rcondWant := 1/(norm(A) * norm(inv(A)))
		rcondWant := 1.0
		if aNorm > 0 && aInvNorm > 0 {
			rcondWant = 1 / aNorm / aInvNorm
		}

		// Compute an estimate of rcond using Dtrcon.
		rcondGot := impl.Dtrcon(norm, uplo, diag, n, a, lda, work, iwork)
		if !floats.Equal(a, aCopy) {
			t.Errorf("%v: unexpected modification of a", name)
		}

		ratio := rCondTestRatio(rcondGot, rcondWant)
		if ratio >= ratioThresh {
			t.Errorf("%v: unexpected value of rcond; got=%v, want=%v (ratio=%v)",
				name, rcondGot, rcondWant, ratio)
		}
	}
}