github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/lapack/testlapack/dorml2.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/blas/blas64"
    14  	"github.com/jingcheng-WU/gonum/floats"
    15  )
    16  
    17  type Dorml2er interface {
    18  	Dgelqfer
    19  	Dorml2(side blas.Side, trans blas.Transpose, m, n, k int, a []float64, lda int, tau, c []float64, ldc int, work []float64)
    20  }
    21  
    22  func Dorml2Test(t *testing.T, impl Dorml2er) {
    23  	rnd := rand.New(rand.NewSource(1))
    24  	// TODO(btracey): This test is not complete, because it
    25  	// doesn't test individual values of m, n, and k, instead only testing
    26  	// a specific subset of possible k values.
    27  	for _, side := range []blas.Side{blas.Left, blas.Right} {
    28  		for _, trans := range []blas.Transpose{blas.NoTrans, blas.Trans} {
    29  			for _, test := range []struct {
    30  				common, adim, cdim, lda, ldc int
    31  			}{
    32  				{3, 4, 5, 0, 0},
    33  				{3, 5, 4, 0, 0},
    34  				{4, 3, 5, 0, 0},
    35  				{4, 5, 3, 0, 0},
    36  				{5, 3, 4, 0, 0},
    37  				{5, 4, 3, 0, 0},
    38  
    39  				{3, 4, 5, 6, 20},
    40  				{3, 5, 4, 6, 20},
    41  				{4, 3, 5, 6, 20},
    42  				{4, 5, 3, 6, 20},
    43  				{5, 3, 4, 6, 20},
    44  				{5, 4, 3, 6, 20},
    45  				{3, 4, 5, 20, 6},
    46  				{3, 5, 4, 20, 6},
    47  				{4, 3, 5, 20, 6},
    48  				{4, 5, 3, 20, 6},
    49  				{5, 3, 4, 20, 6},
    50  				{5, 4, 3, 20, 6},
    51  			} {
    52  				var ma, na, mc, nc int
    53  				if side == blas.Left {
    54  					ma = test.adim
    55  					na = test.common
    56  					mc = test.common
    57  					nc = test.cdim
    58  				} else {
    59  					ma = test.adim
    60  					na = test.common
    61  					mc = test.cdim
    62  					nc = test.common
    63  				}
    64  				// Generate a random matrix
    65  				lda := test.lda
    66  				if lda == 0 {
    67  					lda = na
    68  				}
    69  				a := make([]float64, ma*lda)
    70  				for i := range a {
    71  					a[i] = rnd.Float64()
    72  				}
    73  				ldc := test.ldc
    74  				if ldc == 0 {
    75  					ldc = nc
    76  				}
    77  				// Compute random C matrix
    78  				c := make([]float64, mc*ldc)
    79  				for i := range c {
    80  					c[i] = rnd.Float64()
    81  				}
    82  
    83  				// Compute LQ
    84  				k := min(ma, na)
    85  				tau := make([]float64, k)
    86  				work := make([]float64, 1)
    87  				impl.Dgelqf(ma, na, a, lda, tau, work, -1)
    88  				work = make([]float64, int(work[0]))
    89  				impl.Dgelqf(ma, na, a, lda, tau, work, len(work))
    90  
    91  				// Build Q from result
    92  				q := constructQ("LQ", ma, na, a, lda, tau)
    93  
    94  				cMat := blas64.General{
    95  					Rows:   mc,
    96  					Cols:   nc,
    97  					Stride: ldc,
    98  					Data:   make([]float64, len(c)),
    99  				}
   100  				copy(cMat.Data, c)
   101  				cMatCopy := blas64.General{
   102  					Rows:   cMat.Rows,
   103  					Cols:   cMat.Cols,
   104  					Stride: cMat.Stride,
   105  					Data:   make([]float64, len(cMat.Data)),
   106  				}
   107  				copy(cMatCopy.Data, cMat.Data)
   108  				switch {
   109  				default:
   110  					panic("bad test")
   111  				case side == blas.Left && trans == blas.NoTrans:
   112  					blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, q, cMatCopy, 0, cMat)
   113  				case side == blas.Left && trans == blas.Trans:
   114  					blas64.Gemm(blas.Trans, blas.NoTrans, 1, q, cMatCopy, 0, cMat)
   115  				case side == blas.Right && trans == blas.NoTrans:
   116  					blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, cMatCopy, q, 0, cMat)
   117  				case side == blas.Right && trans == blas.Trans:
   118  					blas64.Gemm(blas.NoTrans, blas.Trans, 1, cMatCopy, q, 0, cMat)
   119  				}
   120  				// Do Dorm2r ard compare
   121  				if side == blas.Left {
   122  					work = make([]float64, nc)
   123  				} else {
   124  					work = make([]float64, mc)
   125  				}
   126  				aCopy := make([]float64, len(a))
   127  				copy(aCopy, a)
   128  				tauCopy := make([]float64, len(tau))
   129  				copy(tauCopy, tau)
   130  				impl.Dorml2(side, trans, mc, nc, k, a, lda, tau, c, ldc, work)
   131  				if !floats.Equal(a, aCopy) {
   132  					t.Errorf("a changed in call")
   133  				}
   134  				if !floats.Equal(tau, tauCopy) {
   135  					t.Errorf("tau changed in call")
   136  				}
   137  				if !floats.EqualApprox(cMat.Data, c, 1e-14) {
   138  					isLeft := side == blas.Left
   139  					isTrans := trans == blas.Trans
   140  					t.Errorf("Multiplication mismatch. IsLeft = %v. IsTrans = %v", isLeft, isTrans)
   141  				}
   142  			}
   143  		}
   144  	}
   145  }