+++ /dev/null
-// 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 (
- "fmt"
- "testing"
-
- "golang.org/x/exp/rand"
-
- "gonum.org/v1/gonum/blas"
- "gonum.org/v1/gonum/floats"
-)
-
-type Dormqrer interface {
- Dorm2rer
- Dormqr(side blas.Side, trans blas.Transpose, m, n, k int, a []float64, lda int, tau, c []float64, ldc int, work []float64, lwork int)
-}
-
-func DormqrTest(t *testing.T, impl Dormqrer) {
- rnd := rand.New(rand.NewSource(1))
- for _, side := range []blas.Side{blas.Left, blas.Right} {
- for _, trans := range []blas.Transpose{blas.NoTrans, blas.Trans} {
- for _, test := range []struct {
- common, adim, cdim, lda, ldc int
- }{
- {6, 7, 8, 0, 0},
- {6, 8, 7, 0, 0},
- {7, 6, 8, 0, 0},
- {7, 8, 6, 0, 0},
- {8, 6, 7, 0, 0},
- {8, 7, 6, 0, 0},
- {100, 200, 300, 0, 0},
- {100, 300, 200, 0, 0},
- {200, 100, 300, 0, 0},
- {200, 300, 100, 0, 0},
- {300, 100, 200, 0, 0},
- {300, 200, 100, 0, 0},
- {100, 200, 300, 400, 500},
- {100, 300, 200, 400, 500},
- {200, 100, 300, 400, 500},
- {200, 300, 100, 400, 500},
- {300, 100, 200, 400, 500},
- {300, 200, 100, 400, 500},
- {100, 200, 300, 500, 400},
- {100, 300, 200, 500, 400},
- {200, 100, 300, 500, 400},
- {200, 300, 100, 500, 400},
- {300, 100, 200, 500, 400},
- {300, 200, 100, 500, 400},
- } {
- var ma, na, mc, nc int
- if side == blas.Left {
- ma = test.common
- na = test.adim
- mc = test.common
- nc = test.cdim
- } else {
- ma = test.common
- na = test.adim
- mc = test.cdim
- nc = test.common
- }
- // Generate a random matrix
- lda := test.lda
- if lda == 0 {
- lda = na
- }
- a := make([]float64, ma*lda)
- for i := range a {
- a[i] = rnd.Float64()
- }
- // Compute random C matrix
- ldc := test.ldc
- if ldc == 0 {
- ldc = nc
- }
- c := make([]float64, mc*ldc)
- for i := range c {
- c[i] = rnd.Float64()
- }
-
- // Compute QR
- k := min(ma, na)
- tau := make([]float64, k)
- work := make([]float64, 1)
- impl.Dgeqrf(ma, na, a, lda, tau, work, -1)
- work = make([]float64, int(work[0]))
- impl.Dgeqrf(ma, na, a, lda, tau, work, len(work))
-
- cCopy := make([]float64, len(c))
- copy(cCopy, c)
- ans := make([]float64, len(c))
- copy(ans, cCopy)
-
- if side == blas.Left {
- work = make([]float64, nc)
- } else {
- work = make([]float64, mc)
- }
- impl.Dorm2r(side, trans, mc, nc, k, a, lda, tau, ans, ldc, work)
-
- // Make sure Dorm2r and Dormqr match with small work
- for i := range work {
- work[i] = rnd.Float64()
- }
- copy(c, cCopy)
- impl.Dormqr(side, trans, mc, nc, k, a, lda, tau, c, ldc, work, len(work))
- if !floats.EqualApprox(c, ans, 1e-12) {
- t.Errorf("Dormqr and Dorm2r mismatch for small work")
- }
-
- // Try with the optimum amount of work
- copy(c, cCopy)
- impl.Dormqr(side, trans, mc, nc, k, nil, lda, nil, nil, ldc, work, -1)
- work = make([]float64, int(work[0]))
- for i := range work {
- work[i] = rnd.Float64()
- }
- impl.Dormqr(side, trans, mc, nc, k, a, lda, tau, c, ldc, work, len(work))
- if !floats.EqualApprox(c, ans, 1e-12) {
- t.Errorf("Dormqr and Dorm2r mismatch for full work")
- fmt.Println("ccopy")
- for i := 0; i < mc; i++ {
- fmt.Println(cCopy[i*ldc : (i+1)*ldc])
- }
- fmt.Println("ans =")
- for i := 0; i < mc; i++ {
- fmt.Println(ans[i*ldc : (i+1)*ldc])
- }
- fmt.Println("c =")
- for i := 0; i < mc; i++ {
- fmt.Println(c[i*ldc : (i+1)*ldc])
- }
- }
-
- // Try with amount of work that is less than
- // optimal but still long enough to use the
- // blocked code.
- copy(c, cCopy)
- if side == blas.Left {
- work = make([]float64, 3*nc)
- } else {
- work = make([]float64, 3*mc)
- }
- impl.Dormqr(side, trans, mc, nc, k, a, lda, tau, c, ldc, work, len(work))
- if !floats.EqualApprox(c, ans, 1e-12) {
- t.Errorf("Dormqr and Dorm2r mismatch for medium work")
- }
- }
- }
- }
-}
OSDN Git Service
