github.com/gopherd/gonum@v0.0.4/lapack/testlapack/dgetri.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  	"testing"
     9  
    10  	"math/rand"
    11  
    12  	"github.com/gopherd/gonum/blas"
    13  	"github.com/gopherd/gonum/blas/blas64"
    14  )
    15  
    16  type Dgetrier interface {
    17  	Dgetrfer
    18  	Dgetri(n int, a []float64, lda int, ipiv []int, work []float64, lwork int) bool
    19  }
    20  
    21  func DgetriTest(t *testing.T, impl Dgetrier) {
    22  	const tol = 1e-13
    23  	rnd := rand.New(rand.NewSource(1))
    24  	bi := blas64.Implementation()
    25  	for _, test := range []struct {
    26  		n, lda int
    27  	}{
    28  		{5, 0},
    29  		{5, 8},
    30  		{45, 0},
    31  		{45, 50},
    32  		{63, 70},
    33  		{64, 70},
    34  		{65, 0},
    35  		{65, 70},
    36  		{66, 70},
    37  		{150, 0},
    38  		{150, 250},
    39  	} {
    40  		n := test.n
    41  		lda := test.lda
    42  		if lda == 0 {
    43  			lda = n
    44  		}
    45  		// Generate a random well conditioned matrix
    46  		perm := rnd.Perm(n)
    47  		a := make([]float64, n*lda)
    48  		for i := 0; i < n; i++ {
    49  			a[i*lda+perm[i]] = 1
    50  		}
    51  		for i := range a {
    52  			a[i] += 0.01 * rnd.Float64()
    53  		}
    54  		aCopy := make([]float64, len(a))
    55  		copy(aCopy, a)
    56  		ipiv := make([]int, n)
    57  		// Compute LU decomposition.
    58  		impl.Dgetrf(n, n, a, lda, ipiv)
    59  		// Test with various workspace sizes.
    60  		for _, wl := range []worklen{minimumWork, mediumWork, optimumWork} {
    61  			ainv := make([]float64, len(a))
    62  			copy(ainv, a)
    63  
    64  			var lwork int
    65  			switch wl {
    66  			case minimumWork:
    67  				lwork = max(1, n)
    68  			case mediumWork:
    69  				work := make([]float64, 1)
    70  				impl.Dgetri(n, ainv, lda, ipiv, work, -1)
    71  				lwork = max(int(work[0])-2*n, n)
    72  			case optimumWork:
    73  				work := make([]float64, 1)
    74  				impl.Dgetri(n, ainv, lda, ipiv, work, -1)
    75  				lwork = int(work[0])
    76  			}
    77  			work := make([]float64, lwork)
    78  
    79  			// Compute inverse.
    80  			ok := impl.Dgetri(n, ainv, lda, ipiv, work, lwork)
    81  			if !ok {
    82  				t.Errorf("Unexpected singular matrix.")
    83  			}
    84  
    85  			// Check that A(inv) * A = I.
    86  			ans := make([]float64, len(ainv))
    87  			bi.Dgemm(blas.NoTrans, blas.NoTrans, n, n, n, 1, aCopy, lda, ainv, lda, 0, ans, lda)
    88  			// The tolerance is so high because computing matrix inverses is very unstable.
    89  			dist := distFromIdentity(n, ans, lda)
    90  			if dist > tol {
    91  				t.Errorf("|Inv(A) * A - I|_inf = %v is too large. n = %v, lda = %v", dist, n, lda)
    92  			}
    93  		}
    94  	}
    95  }