lapack/testlapack/dgecon.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"
	"math"
	"testing"

	"golang.org/x/exp/rand"

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

type Dgeconer interface {
	Dgecon(norm lapack.MatrixNorm, n int, a []float64, lda int, anorm float64, work []float64, iwork []int) float64

	Dgetrier
	Dlanger
}

func DgeconTest(t *testing.T, impl Dgeconer) {
	rnd := rand.New(rand.NewSource(1))
	for _, n := range []int{0, 1, 2, 3, 4, 5, 10, 50} {
		for _, lda := range []int{max(1, n), n + 3} {
			dgeconTest(t, impl, rnd, n, lda)
		}
	}
}

func dgeconTest(t *testing.T, impl Dgeconer, rnd *rand.Rand, n, lda int) {
	const ratioThresh = 10

	// Generate a random square matrix A with elements uniformly in [-1,1).
	a := make([]float64, max(0, (n-1)*lda+n))
	for i := range a {
		a[i] = 2*rnd.Float64() - 1
	}

	// Allocate work slices.
	iwork := make([]int, n)
	work := make([]float64, max(1, 4*n))

	// Compute the LU factorization of A.
	aFac := make([]float64, len(a))
	copy(aFac, a)
	ipiv := make([]int, n)
	ok := impl.Dgetrf(n, n, aFac, lda, ipiv)
	if !ok {
		t.Fatalf("n=%v,lda=%v: bad matrix, Dgetrf failed", n, lda)
	}
	aFacCopy := make([]float64, len(aFac))
	copy(aFacCopy, aFac)

	// Compute the inverse A^{-1} from the LU factorization.
	aInv := make([]float64, len(aFac))
	copy(aInv, aFac)
	ok = impl.Dgetri(n, aInv, lda, ipiv, work, len(work))
	if !ok {
		t.Fatalf("n=%v,lda=%v: bad matrix, Dgetri failed", n, lda)
	}

	for _, norm := range []lapack.MatrixNorm{lapack.MaxColumnSum, lapack.MaxRowSum} {
		name := fmt.Sprintf("norm=%v,n=%v,lda=%v", string(norm), n, lda)

		// Compute the norm of A and A^{-1}.
		aNorm := impl.Dlange(norm, n, n, a, lda, work)
		aInvNorm := impl.Dlange(norm, 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 the LU factorization and Dgecon.
		rcondGot := impl.Dgecon(norm, n, aFac, lda, aNorm, work, iwork)
		if !floats.Equal(aFac, aFacCopy) {
			t.Errorf("%v: unexpected modification of aFac", 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)
		}

		// Check for corner-case values of anorm.
		for _, anorm := range []float64{0, math.Inf(1), math.NaN()} {
			rcondGot = impl.Dgecon(norm, n, aFac, lda, anorm, work, iwork)
			if n == 0 {
				if rcondGot != 1 {
					t.Errorf("%v: unexpected rcond when anorm=%v: got=%v, want=1", name, anorm, rcondGot)
				}
				continue
			}
			if math.IsNaN(anorm) {
				if !math.IsNaN(rcondGot) {
					t.Errorf("%v: NaN not propagated when anorm=NaN: got=%v", name, rcondGot)
				}
				continue
			}
			if rcondGot != 0 {
				t.Errorf("%v: unexpected rcond when anorm=%v: got=%v, want=0", name, anorm, rcondGot)
			}
		}
	}
}