blob: fa4fce0e630afb35de436b2d7da81490fff00a25 [file] [log] [blame]
 // Copyright ©2015 The Gonum Authors. All rights reserved. // Use of this source code is governed by a BSD-style // license that can be found in the LICENSE file. package testlapack import ( "math" "testing" "golang.org/x/exp/rand" "gonum.org/v1/gonum/blas" "gonum.org/v1/gonum/blas/blas64" "gonum.org/v1/gonum/floats" "gonum.org/v1/gonum/lapack" ) type Dlasrer interface { Dlasr(side blas.Side, pivot lapack.Pivot, direct lapack.Direct, m, n int, c, s, a []float64, lda int) } func DlasrTest(t *testing.T, impl Dlasrer) { rnd := rand.New(rand.NewSource(1)) for _, side := range []blas.Side{blas.Left, blas.Right} { for _, pivot := range []lapack.Pivot{lapack.Variable, lapack.Top, lapack.Bottom} { for _, direct := range []lapack.Direct{lapack.Forward, lapack.Backward} { for _, test := range []struct { m, n, lda int }{ {5, 5, 0}, {5, 10, 0}, {10, 5, 0}, {5, 5, 20}, {5, 10, 20}, {10, 5, 20}, } { m := test.m n := test.n lda := test.lda if lda == 0 { lda = n } // Allocate n×n matrix A and fill it with random numbers. a := make([]float64, m*lda) for i := range a { a[i] = rnd.Float64() } // Allocate slices for implicitly // represented rotation matrices. var s, c []float64 if side == blas.Left { s = make([]float64, m-1) c = make([]float64, m-1) } else { s = make([]float64, n-1) c = make([]float64, n-1) } for k := range s { // Generate a random number in [0,2*pi). theta := rnd.Float64() * 2 * math.Pi s[k] = math.Sin(theta) c[k] = math.Cos(theta) } aCopy := make([]float64, len(a)) copy(a, aCopy) // Apply plane a sequence of plane // rotation in s and c to the matrix A. impl.Dlasr(side, pivot, direct, m, n, c, s, a, lda) // Compute a reference solution by multiplying A // by explicitly formed rotation matrix P. pSize := m if side == blas.Right { pSize = n } // Allocate matrix P. p := blas64.General{ Rows: pSize, Cols: pSize, Stride: pSize, Data: make([]float64, pSize*pSize), } // Allocate matrix P_k. pk := blas64.General{ Rows: pSize, Cols: pSize, Stride: pSize, Data: make([]float64, pSize*pSize), } ptmp := blas64.General{ Rows: pSize, Cols: pSize, Stride: pSize, Data: make([]float64, pSize*pSize), } // Initialize P to the identity matrix. for i := 0; i < pSize; i++ { p.Data[i*p.Stride+i] = 1 ptmp.Data[i*p.Stride+i] = 1 } // Iterate over the sequence of plane rotations. for k := range s { // Set P_k to the identity matrix. for i := range p.Data { pk.Data[i] = 0 } for i := 0; i < pSize; i++ { pk.Data[i*p.Stride+i] = 1 } // Set the corresponding elements of P_k. switch pivot { case lapack.Variable: pk.Data[k*p.Stride+k] = c[k] pk.Data[k*p.Stride+k+1] = s[k] pk.Data[(k+1)*p.Stride+k] = -s[k] pk.Data[(k+1)*p.Stride+k+1] = c[k] case lapack.Top: pk.Data[0] = c[k] pk.Data[k+1] = s[k] pk.Data[(k+1)*p.Stride] = -s[k] pk.Data[(k+1)*p.Stride+k+1] = c[k] case lapack.Bottom: pk.Data[(pSize-1-k)*p.Stride+pSize-k-1] = c[k] pk.Data[(pSize-1-k)*p.Stride+pSize-1] = s[k] pk.Data[(pSize-1)*p.Stride+pSize-1-k] = -s[k] pk.Data[(pSize-1)*p.Stride+pSize-1] = c[k] } // Compute P <- P_k * P or P <- P * P_k. if direct == lapack.Forward { blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, pk, ptmp, 0, p) } else { blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, ptmp, pk, 0, p) } copy(ptmp.Data, p.Data) } aMat := blas64.General{ Rows: m, Cols: n, Stride: lda, Data: make([]float64, m*lda), } copy(a, aCopy) newA := blas64.General{ Rows: m, Cols: n, Stride: lda, Data: make([]float64, m*lda), } // Compute P * A or A * P. if side == blas.Left { blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, p, aMat, 0, newA) } else { blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, aMat, p, 0, newA) } // Compare the result from Dlasr with the reference solution. if !floats.EqualApprox(newA.Data, a, 1e-12) { t.Errorf("A update mismatch") } } } } } }