vendor/gonum.org/v1/gonum/lapack/testlapack/dgetf2.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         "testing"
   9
  10         "golang.org/x/exp/rand"
  11
  12         "gonum.org/v1/gonum/blas"
  13         "gonum.org/v1/gonum/blas/blas64"
  14         "gonum.org/v1/gonum/floats"
  15 )
  16
  17 type Dgetf2er interface {
  18         Dgetf2(m, n int, a []float64, lda int, ipiv []int) bool
  19 }
  20
  21 func Dgetf2Test(t *testing.T, impl Dgetf2er) {
  22         rnd := rand.New(rand.NewSource(1))
  23         for _, test := range []struct {
  24                 m, n, lda int
  25         }{
  26                 {10, 10, 0},
  27                 {10, 5, 0},
  28                 {10, 5, 0},
  29
  30                 {10, 10, 20},
  31                 {5, 10, 20},
  32                 {10, 5, 20},
  33         } {
  34                 m := test.m
  35                 n := test.n
  36                 lda := test.lda
  37                 if lda == 0 {
  38                         lda = n
  39                 }
  40                 a := make([]float64, m*lda)
  41                 for i := range a {
  42                         a[i] = rnd.Float64()
  43                 }
  44                 aCopy := make([]float64, len(a))
  45                 copy(aCopy, a)
  46
  47                 mn := min(m, n)
  48                 ipiv := make([]int, mn)
  49                 for i := range ipiv {
  50                         ipiv[i] = rnd.Int()
  51                 }
  52                 ok := impl.Dgetf2(m, n, a, lda, ipiv)
  53                 checkPLU(t, ok, m, n, lda, ipiv, a, aCopy, 1e-14, true)
  54         }
  55
  56         // Test with singular matrices (random matrices are almost surely non-singular).
  57         for _, test := range []struct {
  58                 m, n, lda int
  59                 a         []float64
  60         }{
  61                 {
  62                         m:   2,
  63                         n:   2,
  64                         lda: 2,
  65                         a: []float64{
  66                                 1, 0,
  67                                 0, 0,
  68                         },
  69                 },
  70                 {
  71                         m:   2,
  72                         n:   2,
  73                         lda: 2,
  74                         a: []float64{
  75                                 1, 5,
  76                                 2, 10,
  77                         },
  78                 },
  79                 {
  80                         m:   3,
  81                         n:   3,
  82                         lda: 3,
  83                         // row 3 = row1 + 2 * row2
  84                         a: []float64{
  85                                 1, 5, 7,
  86                                 2, 10, -3,
  87                                 5, 25, 1,
  88                         },
  89                 },
  90                 {
  91                         m:   3,
  92                         n:   4,
  93                         lda: 4,
  94                         // row 3 = row1 + 2 * row2
  95                         a: []float64{
  96                                 1, 5, 7, 9,
  97                                 2, 10, -3, 11,
  98                                 5, 25, 1, 31,
  99                         },
 100                 },
 101         } {
 102                 if impl.Dgetf2(test.m, test.n, test.a, test.lda, make([]int, min(test.m, test.n))) {
 103                         t.Log("Returned ok with singular matrix.")
 104                 }
 105         }
 106 }
 107
 108 // checkPLU checks that the PLU factorization contained in factorize matches
 109 // the original matrix contained in original.
 110 func checkPLU(t *testing.T, ok bool, m, n, lda int, ipiv []int, factorized, original []float64, tol float64, print bool) {
 111         var hasZeroDiagonal bool
 112         for i := 0; i < min(m, n); i++ {
 113                 if factorized[i*lda+i] == 0 {
 114                         hasZeroDiagonal = true
 115                         break
 116                 }
 117         }
 118         if hasZeroDiagonal && ok {
 119                 t.Error("Has a zero diagonal but returned ok")
 120         }
 121         if !hasZeroDiagonal && !ok {
 122                 t.Error("Non-zero diagonal but returned !ok")
 123         }
 124
 125         // Check that the LU decomposition is correct.
 126         mn := min(m, n)
 127         l := make([]float64, m*mn)
 128         ldl := mn
 129         u := make([]float64, mn*n)
 130         ldu := n
 131         for i := 0; i < m; i++ {
 132                 for j := 0; j < n; j++ {
 133                         v := factorized[i*lda+j]
 134                         switch {
 135                         case i == j:
 136                                 l[i*ldl+i] = 1
 137                                 u[i*ldu+i] = v
 138                         case i > j:
 139                                 l[i*ldl+j] = v
 140                         case i < j:
 141                                 u[i*ldu+j] = v
 142                         }
 143                 }
 144         }
 145
 146         LU := blas64.General{
 147                 Rows:   m,
 148                 Cols:   n,
 149                 Stride: n,
 150                 Data:   make([]float64, m*n),
 151         }
 152         U := blas64.General{
 153                 Rows:   mn,
 154                 Cols:   n,
 155                 Stride: ldu,
 156                 Data:   u,
 157         }
 158         L := blas64.General{
 159                 Rows:   m,
 160                 Cols:   mn,
 161                 Stride: ldl,
 162                 Data:   l,
 163         }
 164         blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, L, U, 0, LU)
 165
 166         p := make([]float64, m*m)
 167         ldp := m
 168         for i := 0; i < m; i++ {
 169                 p[i*ldp+i] = 1
 170         }
 171         for i := len(ipiv) - 1; i >= 0; i-- {
 172                 v := ipiv[i]
 173                 blas64.Swap(m, blas64.Vector{Inc: 1, Data: p[i*ldp:]}, blas64.Vector{Inc: 1, Data: p[v*ldp:]})
 174         }
 175         P := blas64.General{
 176                 Rows:   m,
 177                 Cols:   m,
 178                 Stride: m,
 179                 Data:   p,
 180         }
 181         aComp := blas64.General{
 182                 Rows:   m,
 183                 Cols:   n,
 184                 Stride: lda,
 185                 Data:   make([]float64, m*lda),
 186         }
 187         copy(aComp.Data, factorized)
 188         blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, P, LU, 0, aComp)
 189         if !floats.EqualApprox(aComp.Data, original, tol) {
 190                 if print {
 191                         t.Errorf("PLU multiplication does not match original matrix.\nWant: %v\nGot: %v", original, aComp.Data)
 192                         return
 193                 }
 194                 t.Error("PLU multiplication does not match original matrix.")
 195         }
 196 }