vendor/gonum.org/v1/gonum/lapack/testlapack/dgecon.go

   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         "gonum.org/v1/gonum/floats"
  12         "gonum.org/v1/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 }