github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/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  	"testing"
    10  
    11  	"golang.org/x/exp/rand"
    12  
    13  	"github.com/jingcheng-WU/gonum/floats"
    14  	"github.com/jingcheng-WU/gonum/lapack"
    15  )
    16  
    17  type Dgeconer interface {
    18  	Dgecon(norm lapack.MatrixNorm, n int, a []float64, lda int, anorm float64, work []float64, iwork []int) float64
    19  
    20  	Dgetrier
    21  	Dlanger
    22  }
    23  
    24  func DgeconTest(t *testing.T, impl Dgeconer) {
    25  	rnd := rand.New(rand.NewSource(1))
    26  	for _, n := range []int{0, 1, 2, 3, 4, 5, 10, 50} {
    27  		for _, lda := range []int{max(1, n), n + 3} {
    28  			dgeconTest(t, impl, rnd, n, lda)
    29  		}
    30  	}
    31  }
    32  
    33  func dgeconTest(t *testing.T, impl Dgeconer, rnd *rand.Rand, n, lda int) {
    34  	const ratioThresh = 10
    35  
    36  	// Generate a random square matrix A with elements uniformly in [-1,1).
    37  	a := make([]float64, max(0, (n-1)*lda+n))
    38  	for i := range a {
    39  		a[i] = 2*rnd.Float64() - 1
    40  	}
    41  
    42  	// Allocate work slices.
    43  	iwork := make([]int, n)
    44  	work := make([]float64, max(1, 4*n))
    45  
    46  	// Compute the LU factorization of A.
    47  	aFac := make([]float64, len(a))
    48  	copy(aFac, a)
    49  	ipiv := make([]int, n)
    50  	ok := impl.Dgetrf(n, n, aFac, lda, ipiv)
    51  	if !ok {
    52  		t.Fatalf("n=%v,lda=%v: bad matrix, Dgetrf failed", n, lda)
    53  	}
    54  	aFacCopy := make([]float64, len(aFac))
    55  	copy(aFacCopy, aFac)
    56  
    57  	// Compute the inverse A^{-1} from the LU factorization.
    58  	aInv := make([]float64, len(aFac))
    59  	copy(aInv, aFac)
    60  	ok = impl.Dgetri(n, aInv, lda, ipiv, work, len(work))
    61  	if !ok {
    62  		t.Fatalf("n=%v,lda=%v: bad matrix, Dgetri failed", n, lda)
    63  	}
    64  
    65  	for _, norm := range []lapack.MatrixNorm{lapack.MaxColumnSum, lapack.MaxRowSum} {
    66  		name := fmt.Sprintf("norm=%v,n=%v,lda=%v", string(norm), n, lda)
    67  
    68  		// Compute the norm of A and A^{-1}.
    69  		aNorm := impl.Dlange(norm, n, n, a, lda, work)
    70  		aInvNorm := impl.Dlange(norm, n, n, aInv, lda, work)
    71  
    72  		// Compute a good estimate of the condition number
    73  		//  rcondWant := 1/(norm(A) * norm(inv(A)))
    74  		rcondWant := 1.0
    75  		if aNorm > 0 && aInvNorm > 0 {
    76  			rcondWant = 1 / aNorm / aInvNorm
    77  		}
    78  
    79  		// Compute an estimate of rcond using the LU factorization and Dgecon.
    80  		rcondGot := impl.Dgecon(norm, n, aFac, lda, aNorm, work, iwork)
    81  		if !floats.Equal(aFac, aFacCopy) {
    82  			t.Errorf("%v: unexpected modification of aFac", name)
    83  		}
    84  
    85  		ratio := rCondTestRatio(rcondGot, rcondWant)
    86  		if ratio >= ratioThresh {
    87  			t.Errorf("%v: unexpected value of rcond; got=%v, want=%v (ratio=%v)",
    88  				name, rcondGot, rcondWant, ratio)
    89  		}
    90  	}
    91  }