vendor/gonum.org/v1/gonum/lapack/testlapack/dgelqf.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/floats"
  13 )
  14
  15 type Dgelqfer interface {
  16         Dgelq2er
  17         Dgelqf(m, n int, a []float64, lda int, tau, work []float64, lwork int)
  18 }
  19
  20 func DgelqfTest(t *testing.T, impl Dgelqfer) {
  21         rnd := rand.New(rand.NewSource(1))
  22         for c, test := range []struct {
  23                 m, n, lda int
  24         }{
  25                 {10, 5, 0},
  26                 {5, 10, 0},
  27                 {10, 10, 0},
  28                 {300, 5, 0},
  29                 {3, 500, 0},
  30                 {200, 200, 0},
  31                 {300, 200, 0},
  32                 {204, 300, 0},
  33                 {1, 3000, 0},
  34                 {3000, 1, 0},
  35                 {10, 5, 30},
  36                 {5, 10, 30},
  37                 {10, 10, 30},
  38                 {300, 5, 500},
  39                 {3, 500, 600},
  40                 {200, 200, 300},
  41                 {300, 200, 300},
  42                 {204, 300, 400},
  43                 {1, 3000, 4000},
  44                 {3000, 1, 4000},
  45         } {
  46                 m := test.m
  47                 n := test.n
  48                 lda := test.lda
  49                 if lda == 0 {
  50                         lda = n
  51                 }
  52                 a := make([]float64, m*lda)
  53                 for i := 0; i < m; i++ {
  54                         for j := 0; j < n; j++ {
  55                                 a[i*lda+j] = rnd.Float64()
  56                         }
  57                 }
  58                 tau := make([]float64, n)
  59                 for i := 0; i < n; i++ {
  60                         tau[i] = rnd.Float64()
  61                 }
  62                 aCopy := make([]float64, len(a))
  63                 copy(aCopy, a)
  64                 ans := make([]float64, len(a))
  65                 copy(ans, a)
  66                 work := make([]float64, m)
  67                 for i := range work {
  68                         work[i] = rnd.Float64()
  69                 }
  70                 // Compute unblocked QR.
  71                 impl.Dgelq2(m, n, ans, lda, tau, work)
  72                 // Compute blocked QR with small work.
  73                 impl.Dgelqf(m, n, a, lda, tau, work, len(work))
  74                 if !floats.EqualApprox(ans, a, 1e-12) {
  75                         t.Errorf("Case %v, mismatch small work.", c)
  76                 }
  77                 // Try the full length of work.
  78                 impl.Dgelqf(m, n, a, lda, tau, work, -1)
  79                 lwork := int(work[0])
  80                 work = make([]float64, lwork)
  81                 copy(a, aCopy)
  82                 impl.Dgelqf(m, n, a, lda, tau, work, lwork)
  83                 if !floats.EqualApprox(ans, a, 1e-12) {
  84                         t.Errorf("Case %v, mismatch large work.", c)
  85                 }
  86
  87                 // Try a slightly smaller version of work to test blocking code.
  88                 if len(work) <= m {
  89                         continue
  90                 }
  91                 work = work[1:]
  92                 lwork--
  93                 copy(a, aCopy)
  94                 impl.Dgelqf(m, n, a, lda, tau, work, lwork)
  95                 if !floats.EqualApprox(ans, a, 1e-12) {
  96                         t.Errorf("Case %v, mismatch large work.", c)
  97                 }
  98         }
  99 }