--- /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 (
+ "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 Dlarfter interface {
+ Dgeqr2er
+ Dlarft(direct lapack.Direct, store lapack.StoreV, n, k int, v []float64, ldv int, tau []float64, t []float64, ldt int)
+}
+
+func DlarftTest(t *testing.T, impl Dlarfter) {
+ rnd := rand.New(rand.NewSource(1))
+ for _, store := range []lapack.StoreV{lapack.ColumnWise, lapack.RowWise} {
+ for _, direct := range []lapack.Direct{lapack.Forward, lapack.Backward} {
+ for _, test := range []struct {
+ m, n, ldv, ldt int
+ }{
+ {6, 6, 0, 0},
+ {8, 6, 0, 0},
+ {6, 8, 0, 0},
+ {6, 6, 10, 15},
+ {8, 6, 10, 15},
+ {6, 8, 10, 15},
+ {6, 6, 15, 10},
+ {8, 6, 15, 10},
+ {6, 8, 15, 10},
+ } {
+ // Generate a matrix
+ m := test.m
+ n := test.n
+ lda := n
+ if lda == 0 {
+ lda = n
+ }
+
+ a := make([]float64, m*lda)
+ for i := 0; i < m; i++ {
+ for j := 0; j < lda; j++ {
+ a[i*lda+j] = rnd.Float64()
+ }
+ }
+ // Use dgeqr2 to find the v vectors
+ tau := make([]float64, n)
+ work := make([]float64, n)
+ impl.Dgeqr2(m, n, a, lda, tau, work)
+
+ // Construct H using these answers
+ vMatTmp := extractVMat(m, n, a, lda, lapack.Forward, lapack.ColumnWise)
+ vMat := constructVMat(vMatTmp, store, direct)
+ v := vMat.Data
+ ldv := vMat.Stride
+
+ h := constructH(tau, vMat, store, direct)
+
+ k := min(m, n)
+ ldt := test.ldt
+ if ldt == 0 {
+ ldt = k
+ }
+ // Find T from the actual function
+ tm := make([]float64, k*ldt)
+ for i := range tm {
+ tm[i] = 100 + rnd.Float64()
+ }
+ // The v data has been put into a.
+ impl.Dlarft(direct, store, m, k, v, ldv, tau, tm, ldt)
+
+ tData := make([]float64, len(tm))
+ copy(tData, tm)
+ if direct == lapack.Forward {
+ // Zero out the lower traingular portion.
+ for i := 0; i < k; i++ {
+ for j := 0; j < i; j++ {
+ tData[i*ldt+j] = 0
+ }
+ }
+ } else {
+ // Zero out the upper traingular portion.
+ for i := 0; i < k; i++ {
+ for j := i + 1; j < k; j++ {
+ tData[i*ldt+j] = 0
+ }
+ }
+ }
+
+ T := blas64.General{
+ Rows: k,
+ Cols: k,
+ Stride: ldt,
+ Data: tData,
+ }
+
+ vMatT := blas64.General{
+ Rows: vMat.Cols,
+ Cols: vMat.Rows,
+ Stride: vMat.Rows,
+ Data: make([]float64, vMat.Cols*vMat.Rows),
+ }
+ for i := 0; i < vMat.Rows; i++ {
+ for j := 0; j < vMat.Cols; j++ {
+ vMatT.Data[j*vMatT.Stride+i] = vMat.Data[i*vMat.Stride+j]
+ }
+ }
+ var comp blas64.General
+ if store == lapack.ColumnWise {
+ // H = I - V * T * V^T
+ tmp := blas64.General{
+ Rows: T.Rows,
+ Cols: vMatT.Cols,
+ Stride: vMatT.Cols,
+ Data: make([]float64, T.Rows*vMatT.Cols),
+ }
+ // T * V^T
+ blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, T, vMatT, 0, tmp)
+ comp = blas64.General{
+ Rows: vMat.Rows,
+ Cols: tmp.Cols,
+ Stride: tmp.Cols,
+ Data: make([]float64, vMat.Rows*tmp.Cols),
+ }
+ // V * (T * V^T)
+ blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, vMat, tmp, 0, comp)
+ } else {
+ // H = I - V^T * T * V
+ tmp := blas64.General{
+ Rows: T.Rows,
+ Cols: vMat.Cols,
+ Stride: vMat.Cols,
+ Data: make([]float64, T.Rows*vMat.Cols),
+ }
+ // T * V
+ blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, T, vMat, 0, tmp)
+ comp = blas64.General{
+ Rows: vMatT.Rows,
+ Cols: tmp.Cols,
+ Stride: tmp.Cols,
+ Data: make([]float64, vMatT.Rows*tmp.Cols),
+ }
+ // V^T * (T * V)
+ blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, vMatT, tmp, 0, comp)
+ }
+ // I - V^T * T * V
+ for i := 0; i < comp.Rows; i++ {
+ for j := 0; j < comp.Cols; j++ {
+ comp.Data[i*m+j] *= -1
+ if i == j {
+ comp.Data[i*m+j] += 1
+ }
+ }
+ }
+ if !floats.EqualApprox(comp.Data, h.Data, 1e-14) {
+ t.Errorf("T does not construct proper H. Store = %v, Direct = %v.\nWant %v\ngot %v.", string(store), string(direct), h.Data, comp.Data)
+ }
+ }
+ }
+ }
+}
OSDN Git Service
