github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/testlapack/dormr2.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  	"math/rand"
     9  	"testing"
    10  
    11  	"github.com/gonum/blas"
    12  	"github.com/gonum/blas/blas64"
    13  	"github.com/gonum/floats"
    14  )
    15  
    16  type Dormr2er interface {
    17  	Dgerqf(m, n int, a []float64, lda int, tau, work []float64, lwork int)
    18  	Dormr2(side blas.Side, trans blas.Transpose, m, n, k int, a []float64, lda int, tau, c []float64, ldc int, work []float64)
    19  }
    20  
    21  func Dormr2Test(t *testing.T, impl Dormr2er) {
    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 _, test := range []struct {
    26  				common, adim, cdim, lda, ldc int
    27  			}{
    28  				{3, 4, 5, 0, 0},
    29  				{3, 5, 4, 0, 0},
    30  				{4, 3, 5, 0, 0},
    31  				{4, 5, 3, 0, 0},
    32  				{5, 3, 4, 0, 0},
    33  				{5, 4, 3, 0, 0},
    34  				{3, 4, 5, 6, 20},
    35  				{3, 5, 4, 6, 20},
    36  				{4, 3, 5, 6, 20},
    37  				{4, 5, 3, 6, 20},
    38  				{5, 3, 4, 6, 20},
    39  				{5, 4, 3, 6, 20},
    40  				{3, 4, 5, 20, 6},
    41  				{3, 5, 4, 20, 6},
    42  				{4, 3, 5, 20, 6},
    43  				{4, 5, 3, 20, 6},
    44  				{5, 3, 4, 20, 6},
    45  				{5, 4, 3, 20, 6},
    46  			} {
    47  				ma := test.adim
    48  				na := test.common
    49  				var mc, nc int
    50  				if side == blas.Left {
    51  					mc = test.common
    52  					nc = test.cdim
    53  				} else {
    54  					mc = test.cdim
    55  					nc = test.common
    56  				}
    57  
    58  				// Generate a random matrix
    59  				lda := test.lda
    60  				if lda == 0 {
    61  					lda = na
    62  				}
    63  				a := make([]float64, ma*lda)
    64  				for i := range a {
    65  					a[i] = rnd.Float64()
    66  				}
    67  				ldc := test.ldc
    68  				if ldc == 0 {
    69  					ldc = nc
    70  				}
    71  				// Compute random C matrix
    72  				c := make([]float64, mc*ldc)
    73  				for i := range c {
    74  					c[i] = rnd.Float64()
    75  				}
    76  
    77  				// Compute RQ
    78  				k := min(ma, na)
    79  				tau := make([]float64, k)
    80  				work := make([]float64, 1)
    81  				impl.Dgerqf(ma, na, a, lda, tau, work, -1)
    82  				work = make([]float64, int(work[0]))
    83  				impl.Dgerqf(ma, na, a, lda, tau, work, len(work))
    84  
    85  				// Build Q from result
    86  				q := constructQ("RQ", ma, na, a, lda, tau)
    87  
    88  				cMat := blas64.General{
    89  					Rows:   mc,
    90  					Cols:   nc,
    91  					Stride: ldc,
    92  					Data:   make([]float64, len(c)),
    93  				}
    94  				copy(cMat.Data, c)
    95  				cMatCopy := blas64.General{
    96  					Rows:   cMat.Rows,
    97  					Cols:   cMat.Cols,
    98  					Stride: cMat.Stride,
    99  					Data:   make([]float64, len(cMat.Data)),
   100  				}
   101  				copy(cMatCopy.Data, cMat.Data)
   102  				switch {
   103  				default:
   104  					panic("bad test")
   105  				case side == blas.Left && trans == blas.NoTrans:
   106  					blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, q, cMatCopy, 0, cMat)
   107  				case side == blas.Left && trans == blas.Trans:
   108  					blas64.Gemm(blas.Trans, blas.NoTrans, 1, q, cMatCopy, 0, cMat)
   109  				case side == blas.Right && trans == blas.NoTrans:
   110  					blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, cMatCopy, q, 0, cMat)
   111  				case side == blas.Right && trans == blas.Trans:
   112  					blas64.Gemm(blas.NoTrans, blas.Trans, 1, cMatCopy, q, 0, cMat)
   113  				}
   114  				// Do Dorm2r ard compare
   115  				if side == blas.Left {
   116  					work = make([]float64, nc)
   117  				} else {
   118  					work = make([]float64, mc)
   119  				}
   120  				aCopy := make([]float64, len(a))
   121  				copy(aCopy, a)
   122  				tauCopy := make([]float64, len(tau))
   123  				copy(tauCopy, tau)
   124  				impl.Dormr2(side, trans, mc, nc, k, a[(ma-k)*lda:], lda, tau, c, ldc, work)
   125  				if !floats.Equal(a, aCopy) {
   126  					t.Errorf("a changed in call")
   127  				}
   128  				if !floats.Equal(tau, tauCopy) {
   129  					t.Errorf("tau changed in call")
   130  				}
   131  				if !floats.EqualApprox(cMat.Data, c, 1e-14) {
   132  					t.Errorf("Multiplication mismatch.\n Want %v \n got %v.", cMat.Data, c)
   133  				}
   134  			}
   135  		}
   136  	}
   137  }