github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/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  	"log"
     9  	"testing"
    10  
    11  	"github.com/gonum/floats"
    12  	"github.com/gonum/lapack"
    13  )
    14  
    15  type Dgeconer interface {
    16  	Dlanger
    17  	Dgetrfer
    18  	Dgecon(norm lapack.MatrixNorm, n int, a []float64, lda int, anorm float64, work []float64, iwork []int) float64
    19  }
    20  
    21  func DgeconTest(t *testing.T, impl Dgeconer) {
    22  	for _, test := range []struct {
    23  		m       int
    24  		n       int
    25  		a       []float64
    26  		condOne float64
    27  		condInf float64
    28  	}{
    29  		{
    30  			a: []float64{
    31  				8, 1, 6,
    32  				3, 5, 7,
    33  				4, 9, 2,
    34  			},
    35  			m:       3,
    36  			n:       3,
    37  			condOne: 3.0 / 16,
    38  			condInf: 3.0 / 16,
    39  		},
    40  		{
    41  			a: []float64{
    42  				2, 9, 3, 2,
    43  				10, 9, 9, 3,
    44  				1, 1, 5, 2,
    45  				8, 4, 10, 2,
    46  			},
    47  			m:       4,
    48  			n:       4,
    49  			condOne: 0.024740155174938,
    50  			condInf: 0.012034465570035,
    51  		},
    52  		// Dgecon does not match Dpocon for this case. https://github.com/xianyi/OpenBLAS/issues/664.
    53  		{
    54  			a: []float64{
    55  				2.9995576045549965, -2.0898894566158663, 3.965560740124006,
    56  				-2.0898894566158663, 1.9634729526261008, -2.8681002706874104,
    57  				3.965560740124006, -2.8681002706874104, 5.502416670471008,
    58  			},
    59  			m:       3,
    60  			n:       3,
    61  			condOne: 0.024054837369015203,
    62  			condInf: 0.024054837369015203,
    63  		},
    64  	} {
    65  		m := test.m
    66  		n := test.n
    67  		lda := n
    68  		a := make([]float64, len(test.a))
    69  		copy(a, test.a)
    70  		ipiv := make([]int, min(m, n))
    71  
    72  		// Find the norms of the original matrix.
    73  		work := make([]float64, 4*n)
    74  		oneNorm := impl.Dlange(lapack.MaxColumnSum, m, n, a, lda, work)
    75  		infNorm := impl.Dlange(lapack.MaxRowSum, m, n, a, lda, work)
    76  
    77  		// Compute LU factorization of a.
    78  		impl.Dgetrf(m, n, a, lda, ipiv)
    79  
    80  		// Compute the condition number
    81  		iwork := make([]int, n)
    82  		condOne := impl.Dgecon(lapack.MaxColumnSum, n, a, lda, oneNorm, work, iwork)
    83  		condInf := impl.Dgecon(lapack.MaxRowSum, n, a, lda, infNorm, work, iwork)
    84  
    85  		// Error if not the same order, otherwise log the difference.
    86  		if !floats.EqualWithinAbsOrRel(condOne, test.condOne, 1e0, 1e0) {
    87  			t.Errorf("One norm mismatch. Want %v, got %v.", test.condOne, condOne)
    88  		} else if !floats.EqualWithinAbsOrRel(condOne, test.condOne, 1e-14, 1e-14) {
    89  			log.Printf("Dgecon one norm mismatch. Want %v, got %v.", test.condOne, condOne)
    90  		}
    91  		if !floats.EqualWithinAbsOrRel(condInf, test.condInf, 1e0, 1e0) {
    92  			t.Errorf("One norm mismatch. Want %v, got %v.", test.condInf, condInf)
    93  		} else if !floats.EqualWithinAbsOrRel(condInf, test.condInf, 1e-14, 1e-14) {
    94  			log.Printf("Dgecon one norm mismatch. Want %v, got %v.", test.condInf, condInf)
    95  		}
    96  	}
    97  }