From fb46a9be5ed23643ec24ae7dca1e19f9cfb996ca Mon Sep 17 00:00:00 2001 From: Vladimir Chalupecky Date: Mon, 30 May 2016 17:28:12 +0900 Subject: [PATCH] testlapack: add test for Dgehrd --- native/lapack_test.go | 4 + testlapack/dgehrd.go | 230 ++++++++++++++++++++++++++++++++++++++++++ testlapack/general.go | 6 ++ 3 files changed, 240 insertions(+) create mode 100644 testlapack/dgehrd.go diff --git a/native/lapack_test.go b/native/lapack_test.go index 8b6640ad..d9506cf8 100644 --- a/native/lapack_test.go +++ b/native/lapack_test.go @@ -32,6 +32,10 @@ func TestDgehd2(t *testing.T) { testlapack.Dgehd2Test(t, impl) } +func TestDgehrd(t *testing.T) { + testlapack.DgehrdTest(t, impl) +} + func TestDgelqf(t *testing.T) { testlapack.DgelqfTest(t, impl) } diff --git a/testlapack/dgehrd.go b/testlapack/dgehrd.go new file mode 100644 index 00000000..e0bbeae3 --- /dev/null +++ b/testlapack/dgehrd.go @@ -0,0 +1,230 @@ +// Copyright ©2016 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" + "math" + "math/rand" + "testing" + + "github.com/gonum/blas" + "github.com/gonum/blas/blas64" +) + +type Dgehrder interface { + Dgehrd(n, ilo, ihi int, a []float64, lda int, tau, work []float64, lwork int) +} + +func DgehrdTest(t *testing.T, impl Dgehrder) { + rnd := rand.New(rand.NewSource(1)) + + // Randomized tests for small matrix sizes that will most likely + // use the unblocked algorithm. + for _, n := range []int{1, 2, 3, 4, 5, 10, 34} { + for _, extra := range []int{0, 1, 13} { + for _, optwork := range []bool{true, false} { + for cas := 0; cas < 100; cas++ { + ilo := rnd.Intn(n) + ihi := rnd.Intn(n) + if ilo > ihi { + ilo, ihi = ihi, ilo + } + testDgehrd(t, impl, n, ilo, ihi, extra, optwork, rnd) + } + } + } + } + + // These are selected tests for larger matrix sizes to test the blocked + // algorithm. Use sizes around several powers of two because that is + // where the blocked path will most likely start to be taken. For + // example, at present the blocked algorithm is used for sizes larger + // than 129. + for _, test := range []struct { + n, ilo, ihi int + }{ + {0, 0, -1}, + + {68, 0, 63}, + {68, 0, 64}, + {68, 0, 65}, + {68, 0, 66}, + {68, 0, 67}, + + {132, 2, 129}, + {132, 1, 129}, // Size = 129, unblocked. + {132, 0, 129}, // Size = 130, blocked. + {132, 1, 130}, + {132, 0, 130}, + {132, 1, 131}, + {132, 0, 131}, + + {260, 2, 257}, + {260, 1, 257}, + {260, 0, 257}, + {260, 0, 258}, + {260, 0, 259}, + } { + for _, extra := range []int{0, 1, 13} { + for _, optwork := range []bool{true, false} { + testDgehrd(t, impl, test.n, test.ilo, test.ihi, extra, optwork, rnd) + } + } + } +} + +func testDgehrd(t *testing.T, impl Dgehrder, n, ilo, ihi, extra int, optwork bool, rnd *rand.Rand) { + a := randomGeneral(n, n, n+extra, rnd) + aCopy := a + aCopy.Data = make([]float64, len(a.Data)) + copy(aCopy.Data, a.Data) + + var tau []float64 + if n > 1 { + tau = nanSlice(n - 1) + } + + var work []float64 + if optwork { + work = nanSlice(1) + impl.Dgehrd(n, ilo, ihi, a.Data, a.Stride, tau, work, -1) + work = nanSlice(int(work[0])) + } else { + work = nanSlice(max(1, n)) + } + + impl.Dgehrd(n, ilo, ihi, a.Data, a.Stride, tau, work, len(work)) + + if n == 0 { + // Just make sure there is no panic. + return + } + + prefix := fmt.Sprintf("Case n=%v, ilo=%v, ihi=%v, extra=%v", n, ilo, ihi, extra) + + // Check any invalid modifications of a. + if !generalOutsideAllNaN(a) { + t.Errorf("%v: out-of-range write to A\n%v", prefix, a.Data) + } + for i := ilo; i <= ihi; i++ { + for j := 0; j < min(ilo, i); j++ { + if a.Data[i*a.Stride+j] != aCopy.Data[i*aCopy.Stride+j] { + t.Errorf("%v: unexpected modification of A[%v,%v]", prefix, i, j) + } + } + } + for i := ihi + 1; i < n; i++ { + for j := 0; j < i; j++ { + if a.Data[i*a.Stride+j] != aCopy.Data[i*aCopy.Stride+j] { + t.Errorf("%v: unexpected modification of A[%v,%v]", prefix, i, j) + } + } + } + for i := 0; i <= ilo; i++ { + for j := i; j < ilo+1; j++ { + if a.Data[i*a.Stride+j] != aCopy.Data[i*aCopy.Stride+j] { + t.Errorf("%v: unexpected modification at A[%v,%v]", prefix, i, j) + } + } + for j := ihi + 1; j < n; j++ { + if a.Data[i*a.Stride+j] != aCopy.Data[i*aCopy.Stride+j] { + t.Errorf("%v: unexpected modification at A[%v,%v]", prefix, i, j) + } + } + } + for i := ihi + 1; i < n; i++ { + for j := i; j < n; j++ { + if a.Data[i*a.Stride+j] != aCopy.Data[i*aCopy.Stride+j] { + t.Errorf("%v: unexpected modification at A[%v,%v]", prefix, i, j) + } + } + } + + // Check that tau has been assigned properly. + for i, v := range tau { + if math.IsNaN(v) { + t.Errorf("%v: unexpected NaN at tau[%v]", prefix, i) + } + } + + // Extract Q and check that it is orthogonal. + q := blas64.General{ + Rows: n, + Cols: n, + Stride: n, + Data: make([]float64, n*n), + } + for i := 0; i < q.Rows; i++ { + q.Data[i*q.Stride+i] = 1 + } + qCopy := q + qCopy.Data = make([]float64, len(q.Data)) + for j := ilo; j < ihi; j++ { + h := blas64.General{ + Rows: n, + Cols: n, + Stride: n, + Data: make([]float64, n*n), + } + for i := 0; i < h.Rows; i++ { + h.Data[i*h.Stride+i] = 1 + } + v := blas64.Vector{ + Inc: 1, + Data: make([]float64, n), + } + v.Data[j+1] = 1 + for i := j + 2; i < ihi+1; i++ { + v.Data[i] = a.Data[i*a.Stride+j] + } + blas64.Ger(-tau[j], v, v, h) + copy(qCopy.Data, q.Data) + blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, qCopy, h, 0, q) + } + if !isOrthonormal(q) { + t.Errorf("%v: Q is not orthogonal\nQ=%v", prefix, q) + } + + // Construct Q^T * AOrig * Q and check that it is upper Hessenberg. + aq := blas64.General{ + Rows: n, + Cols: n, + Stride: n, + Data: make([]float64, n*n), + } + blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, aCopy, q, 0, aq) + qaq := blas64.General{ + Rows: n, + Cols: n, + Stride: n, + Data: make([]float64, n*n), + } + blas64.Gemm(blas.Trans, blas.NoTrans, 1, q, aq, 0, qaq) + for i := 0; i <= ilo; i++ { + for j := ilo + 1; j <= ihi; j++ { + qaqij := qaq.Data[i*qaq.Stride+j] + diff := qaqij - a.Data[i*a.Stride+j] + if math.Abs(diff) > 1e-13 { + t.Errorf("%v: Q^T*AOrig*Q and A are not equal, diff at [%v,%v]=%v", prefix, i, j, diff) + } + } + } + for i := ilo + 1; i <= ihi; i++ { + for j := ilo; j < n; j++ { + qaqij := qaq.Data[i*qaq.Stride+j] + if j < i-1 { + if math.Abs(qaqij) > 1e-13 { + t.Errorf("%v: Q^T*AOrig*Q is not upper Hessenberg, [%v,%v]=%v", prefix, i, j, qaqij) + } + continue + } + diff := qaqij - a.Data[i*a.Stride+j] + if math.Abs(diff) > 1e-13 { + t.Errorf("%v: Q^T*AOrig*Q and A are not equal, diff at [%v,%v]=%v", prefix, i, j, diff) + } + } + } +} diff --git a/testlapack/general.go b/testlapack/general.go index e5639ee6..e4dc27ab 100644 --- a/testlapack/general.go +++ b/testlapack/general.go @@ -40,6 +40,12 @@ func nanSlice(n int) []float64 { // nanGeneral allocates a new r×c general matrix filled with NaN values. func nanGeneral(r, c, stride int) blas64.General { + if r < 0 || c < 0 { + panic("bad matrix size") + } + if r == 0 || c == 0 { + return blas64.General{} + } return blas64.General{ Rows: r, Cols: c,