X-Git-Url: http://git.osdn.net/view?a=blobdiff_plain;ds=sidebyside;f=vendor%2Fgonum.org%2Fv1%2Fgonum%2Flapack%2Ftestlapack%2Fdormhr.go;fp=vendor%2Fgonum.org%2Fv1%2Fgonum%2Flapack%2Ftestlapack%2Fdormhr.go;h=af13d21b903dfb686c7dcdc8482e9457bd7d51d8;hb=db158dcf09436b003defd333f1a665e7e051d820;hp=0000000000000000000000000000000000000000;hpb=d09b7a78d44dc259725902b8141cdba0d716b121;p=bytom%2Fvapor.git diff --git a/vendor/gonum.org/v1/gonum/lapack/testlapack/dormhr.go b/vendor/gonum.org/v1/gonum/lapack/testlapack/dormhr.go new file mode 100644 index 00000000..af13d21b --- /dev/null +++ b/vendor/gonum.org/v1/gonum/lapack/testlapack/dormhr.go @@ -0,0 +1,133 @@ +// 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" + "testing" + + "golang.org/x/exp/rand" + + "gonum.org/v1/gonum/blas" + "gonum.org/v1/gonum/blas/blas64" +) + +type Dormhrer interface { + Dormhr(side blas.Side, trans blas.Transpose, m, n, ilo, ihi int, a []float64, lda int, tau, c []float64, ldc int, work []float64, lwork int) + + Dgehrder +} + +func DormhrTest(t *testing.T, impl Dormhrer) { + 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 _, m := range []int{1, 2, 3, 4, 5, 8, 9, 10, 23} { + for _, n := range []int{1, 2, 3, 4, 5, 8, 9, 10, 23} { + for _, extra := range []int{0, 1, 13} { + for cas := 0; cas < 10; cas++ { + nq := m + if side == blas.Right { + nq = n + } + ilo := rnd.Intn(nq) + ihi := rnd.Intn(nq) + if ilo > ihi { + ilo, ihi = ihi, ilo + } + testDormhr(t, impl, side, trans, m, n, ilo, ihi, extra, true, rnd) + testDormhr(t, impl, side, trans, m, n, ilo, ihi, extra, false, rnd) + } + } + } + } + } + } + for _, side := range []blas.Side{blas.Left, blas.Right} { + for _, trans := range []blas.Transpose{blas.NoTrans, blas.Trans} { + testDormhr(t, impl, side, trans, 0, 0, 0, -1, 0, true, rnd) + testDormhr(t, impl, side, trans, 0, 0, 0, -1, 0, false, rnd) + } + } +} + +func testDormhr(t *testing.T, impl Dormhrer, side blas.Side, trans blas.Transpose, m, n, ilo, ihi, extra int, optwork bool, rnd *rand.Rand) { + const tol = 1e-14 + + var nq, nw int + switch side { + case blas.Left: + nq = m + nw = n + case blas.Right: + nq = n + nw = m + } + + // Compute the elementary reflectors and tau. + a := randomGeneral(nq, nq, nq+extra, rnd) + var tau []float64 + if nq > 1 { + tau = nanSlice(nq - 1) + } + work := nanSlice(max(1, nq)) // Minimum work for Dgehrd. + impl.Dgehrd(nq, ilo, ihi, a.Data, a.Stride, tau, work, len(work)) + + // Construct Q from the elementary reflectors in a and from tau. + q := eye(nq, nq) + qCopy := eye(nq, nq) + for j := ilo; j < ihi; j++ { + h := eye(nq, nq) + v := blas64.Vector{ + Inc: 1, + Data: make([]float64, nq), + } + 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) + } + + c := randomGeneral(m, n, n+extra, rnd) + + // Compute the product of Q and C explicitly. + qc := randomGeneral(m, n, n+extra, rnd) + if side == blas.Left { + blas64.Gemm(trans, blas.NoTrans, 1, q, c, 0, qc) + } else { + blas64.Gemm(blas.NoTrans, trans, 1, c, q, 0, qc) + } + + // Compute the product of Q and C using Dormhr. + if optwork { + work = nanSlice(1) + impl.Dormhr(side, trans, m, n, ilo, ihi, nil, a.Stride, nil, nil, c.Stride, work, -1) + work = nanSlice(int(work[0])) + } else { + work = nanSlice(max(1, nw)) + } + impl.Dormhr(side, trans, m, n, ilo, ihi, a.Data, a.Stride, tau, c.Data, c.Stride, work, len(work)) + + // Compare the two answers. + prefix := fmt.Sprintf("Case side=%v, trans=%v, m=%v, n=%v, ilo=%v, ihi=%v, extra=%v, optwork=%v", + side, trans, m, n, ilo, ihi, extra, optwork) + if !generalOutsideAllNaN(c) { + t.Errorf("%v: out-of-range write to C\n%v", prefix, c.Data) + } + for i := 0; i < m; i++ { + for j := 0; j < n; j++ { + cij := c.Data[i*c.Stride+j] + qcij := qc.Data[i*qc.Stride+j] + if math.Abs(cij-qcij) > tol { + t.Errorf("%v: unexpected value of the QC product at [%v,%v]: want %v, got %v", prefix, i, j, qcij, cij) + } + } + } +}