gonum.org/v1/gonum@v0.15.1-0.20240517103525-f853624cb1bb/lapack/testlapack/dgecon.go (about)

     1  // Copyright ©2015 The Gonum Authors. All rights reserved.
     2  // Use of this source code is governed by a BSD-style
     3  // license that can be found in the LICENSE file.
     4  
     5  package testlapack
     6  
     7  import (
     8  	"fmt"
     9  	"math"
    10  	"testing"
    11  
    12  	"golang.org/x/exp/rand"
    13  
    14  	"gonum.org/v1/gonum/floats"
    15  	"gonum.org/v1/gonum/lapack"
    16  )
    17  
    18  type Dgeconer interface {
    19  	Dgecon(norm lapack.MatrixNorm, n int, a []float64, lda int, anorm float64, work []float64, iwork []int) float64
    20  
    21  	Dgetrier
    22  	Dlanger
    23  }
    24  
    25  func DgeconTest(t *testing.T, impl Dgeconer) {
    26  	rnd := rand.New(rand.NewSource(1))
    27  	for _, n := range []int{0, 1, 2, 3, 4, 5, 10, 50} {
    28  		for _, lda := range []int{max(1, n), n + 3} {
    29  			dgeconTest(t, impl, rnd, n, lda)
    30  		}
    31  	}
    32  }
    33  
    34  func dgeconTest(t *testing.T, impl Dgeconer, rnd *rand.Rand, n, lda int) {
    35  	const ratioThresh = 10
    36  
    37  	// Generate a random square matrix A with elements uniformly in [-1,1).
    38  	a := make([]float64, max(0, (n-1)*lda+n))
    39  	for i := range a {
    40  		a[i] = 2*rnd.Float64() - 1
    41  	}
    42  
    43  	// Allocate work slices.
    44  	iwork := make([]int, n)
    45  	work := make([]float64, max(1, 4*n))
    46  
    47  	// Compute the LU factorization of A.
    48  	aFac := make([]float64, len(a))
    49  	copy(aFac, a)
    50  	ipiv := make([]int, n)
    51  	ok := impl.Dgetrf(n, n, aFac, lda, ipiv)
    52  	if !ok {
    53  		t.Fatalf("n=%v,lda=%v: bad matrix, Dgetrf failed", n, lda)
    54  	}
    55  	aFacCopy := make([]float64, len(aFac))
    56  	copy(aFacCopy, aFac)
    57  
    58  	// Compute the inverse A^{-1} from the LU factorization.
    59  	aInv := make([]float64, len(aFac))
    60  	copy(aInv, aFac)
    61  	ok = impl.Dgetri(n, aInv, lda, ipiv, work, len(work))
    62  	if !ok {
    63  		t.Fatalf("n=%v,lda=%v: bad matrix, Dgetri failed", n, lda)
    64  	}
    65  
    66  	for _, norm := range []lapack.MatrixNorm{lapack.MaxColumnSum, lapack.MaxRowSum} {
    67  		name := fmt.Sprintf("norm=%v,n=%v,lda=%v", string(norm), n, lda)
    68  
    69  		// Compute the norm of A and A^{-1}.
    70  		aNorm := impl.Dlange(norm, n, n, a, lda, work)
    71  		aInvNorm := impl.Dlange(norm, n, n, aInv, lda, work)
    72  
    73  		// Compute a good estimate of the condition number
    74  		//  rcondWant := 1/(norm(A) * norm(inv(A)))
    75  		rcondWant := 1.0
    76  		if aNorm > 0 && aInvNorm > 0 {
    77  			rcondWant = 1 / aNorm / aInvNorm
    78  		}
    79  
    80  		// Compute an estimate of rcond using the LU factorization and Dgecon.
    81  		rcondGot := impl.Dgecon(norm, n, aFac, lda, aNorm, work, iwork)
    82  		if !floats.Equal(aFac, aFacCopy) {
    83  			t.Errorf("%v: unexpected modification of aFac", name)
    84  		}
    85  
    86  		ratio := rCondTestRatio(rcondGot, rcondWant)
    87  		if ratio >= ratioThresh {
    88  			t.Errorf("%v: unexpected value of rcond; got=%v, want=%v (ratio=%v)",
    89  				name, rcondGot, rcondWant, ratio)
    90  		}
    91  
    92  		// Check for corner-case values of anorm.
    93  		for _, anorm := range []float64{0, math.Inf(1), math.NaN()} {
    94  			rcondGot = impl.Dgecon(norm, n, aFac, lda, anorm, work, iwork)
    95  			if n == 0 {
    96  				if rcondGot != 1 {
    97  					t.Errorf("%v: unexpected rcond when anorm=%v: got=%v, want=1", name, anorm, rcondGot)
    98  				}
    99  				continue
   100  			}
   101  			if math.IsNaN(anorm) {
   102  				if !math.IsNaN(rcondGot) {
   103  					t.Errorf("%v: NaN not propagated when anorm=NaN: got=%v", name, rcondGot)
   104  				}
   105  				continue
   106  			}
   107  			if rcondGot != 0 {
   108  				t.Errorf("%v: unexpected rcond when anorm=%v: got=%v, want=0", name, anorm, rcondGot)
   109  			}
   110  		}
   111  	}
   112  }