github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/lapack/testlapack/dormlq.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  	"golang.org/x/exp/rand"
    11  
    12  	"github.com/jingcheng-WU/gonum/blas"
    13  	"github.com/jingcheng-WU/gonum/floats"
    14  )
    15  
    16  type Dormlqer interface {
    17  	Dorml2er
    18  	Dormlq(side blas.Side, trans blas.Transpose, m, n, k int, a []float64, lda int, tau, c []float64, ldc int, work []float64, lwork int)
    19  }
    20  
    21  func DormlqTest(t *testing.T, impl Dormlqer) {
    22  	rnd := rand.New(rand.NewSource(1))
    23  	for _, side := range []blas.Side{blas.Left, blas.Right} {
    24  		for _, trans := range []blas.Transpose{blas.NoTrans, blas.Trans} {
    25  			for _, wl := range []worklen{minimumWork, mediumWork, optimumWork} {
    26  				for _, test := range []struct {
    27  					common, adim, cdim, lda, ldc int
    28  				}{
    29  					{0, 0, 0, 0, 0},
    30  					{6, 7, 8, 0, 0},
    31  					{6, 8, 7, 0, 0},
    32  					{7, 6, 8, 0, 0},
    33  					{7, 8, 6, 0, 0},
    34  					{8, 6, 7, 0, 0},
    35  					{8, 7, 6, 0, 0},
    36  					{100, 200, 300, 0, 0},
    37  					{100, 300, 200, 0, 0},
    38  					{200, 100, 300, 0, 0},
    39  					{200, 300, 100, 0, 0},
    40  					{300, 100, 200, 0, 0},
    41  					{300, 200, 100, 0, 0},
    42  					{100, 200, 300, 400, 500},
    43  					{100, 300, 200, 400, 500},
    44  					{200, 100, 300, 400, 500},
    45  					{200, 300, 100, 400, 500},
    46  					{300, 100, 200, 400, 500},
    47  					{300, 200, 100, 400, 500},
    48  					{100, 200, 300, 500, 400},
    49  					{100, 300, 200, 500, 400},
    50  					{200, 100, 300, 500, 400},
    51  					{200, 300, 100, 500, 400},
    52  					{300, 100, 200, 500, 400},
    53  					{300, 200, 100, 500, 400},
    54  				} {
    55  					var ma, na, mc, nc int
    56  					if side == blas.Left {
    57  						ma = test.adim
    58  						na = test.common
    59  						mc = test.common
    60  						nc = test.cdim
    61  					} else {
    62  						ma = test.adim
    63  						na = test.common
    64  						mc = test.cdim
    65  						nc = test.common
    66  					}
    67  					// Generate a random matrix
    68  					lda := test.lda
    69  					if lda == 0 {
    70  						lda = max(1, na)
    71  					}
    72  					a := make([]float64, ma*lda)
    73  					for i := range a {
    74  						a[i] = rnd.Float64()
    75  					}
    76  					// Compute random C matrix
    77  					ldc := test.ldc
    78  					if ldc == 0 {
    79  						ldc = nc
    80  					}
    81  					c := make([]float64, mc*ldc)
    82  					for i := range c {
    83  						c[i] = rnd.Float64()
    84  					}
    85  
    86  					// Compute LQ
    87  					k := min(ma, na)
    88  					tau := make([]float64, k)
    89  					work := make([]float64, 1)
    90  					impl.Dgelqf(ma, na, a, lda, tau, work, -1)
    91  					work = make([]float64, int(work[0]))
    92  					impl.Dgelqf(ma, na, a, lda, tau, work, len(work))
    93  
    94  					cCopy := make([]float64, len(c))
    95  					copy(cCopy, c)
    96  					ans := make([]float64, len(c))
    97  					copy(ans, cCopy)
    98  
    99  					var nw int
   100  					if side == blas.Left {
   101  						nw = nc
   102  					} else {
   103  						nw = mc
   104  					}
   105  					work = make([]float64, max(1, nw))
   106  					impl.Dorml2(side, trans, mc, nc, k, a, lda, tau, ans, ldc, work)
   107  
   108  					var lwork int
   109  					switch wl {
   110  					case minimumWork:
   111  						lwork = nw
   112  					case optimumWork:
   113  						impl.Dormlq(side, trans, mc, nc, k, a, lda, tau, c, ldc, work, -1)
   114  						lwork = int(work[0])
   115  					case mediumWork:
   116  						work := make([]float64, 1)
   117  						impl.Dormlq(side, trans, mc, nc, k, a, lda, tau, c, ldc, work, -1)
   118  						lwork = (int(work[0]) + nw) / 2
   119  					}
   120  					lwork = max(1, lwork)
   121  					work = make([]float64, lwork)
   122  
   123  					impl.Dormlq(side, trans, mc, nc, k, a, lda, tau, c, ldc, work, lwork)
   124  					if !floats.EqualApprox(c, ans, 1e-13) {
   125  						t.Errorf("Dormqr and Dorm2r results mismatch")
   126  					}
   127  				}
   128  			}
   129  		}
   130  	}
   131  }