+++ /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/blas/blas64"
- "gonum.org/v1/gonum/floats"
- "gonum.org/v1/gonum/lapack"
-)
-
-type Dgesvder interface {
- Dgesvd(jobU, jobVT lapack.SVDJob, m, n int, a []float64, lda int, s, u []float64, ldu int, vt []float64, ldvt int, work []float64, lwork int) (ok bool)
-}
-
-func DgesvdTest(t *testing.T, impl Dgesvder) {
- rnd := rand.New(rand.NewSource(1))
- // TODO(btracey): Add tests for all of the cases when the SVD implementation
- // is finished.
- // TODO(btracey): Add tests for m > mnthr and n > mnthr when other SVD
- // conditions are implemented. Right now mnthr is 5,000,000 which is too
- // large to create a square matrix of that size.
- for _, test := range []struct {
- m, n, lda, ldu, ldvt int
- }{
- {5, 5, 0, 0, 0},
- {5, 6, 0, 0, 0},
- {6, 5, 0, 0, 0},
- {5, 9, 0, 0, 0},
- {9, 5, 0, 0, 0},
-
- {5, 5, 10, 11, 12},
- {5, 6, 10, 11, 12},
- {6, 5, 10, 11, 12},
- {5, 5, 10, 11, 12},
- {5, 9, 10, 11, 12},
- {9, 5, 10, 11, 12},
-
- {300, 300, 0, 0, 0},
- {300, 400, 0, 0, 0},
- {400, 300, 0, 0, 0},
- {300, 600, 0, 0, 0},
- {600, 300, 0, 0, 0},
-
- {300, 300, 400, 450, 460},
- {300, 400, 500, 550, 560},
- {400, 300, 550, 550, 560},
- {300, 600, 700, 750, 760},
- {600, 300, 700, 750, 760},
- } {
- jobU := lapack.SVDAll
- jobVT := lapack.SVDAll
-
- m := test.m
- n := test.n
- lda := test.lda
- if lda == 0 {
- lda = n
- }
- ldu := test.ldu
- if ldu == 0 {
- ldu = m
- }
- ldvt := test.ldvt
- if ldvt == 0 {
- ldvt = n
- }
-
- a := make([]float64, m*lda)
- for i := range a {
- a[i] = rnd.NormFloat64()
- }
-
- u := make([]float64, m*ldu)
- for i := range u {
- u[i] = rnd.NormFloat64()
- }
-
- vt := make([]float64, n*ldvt)
- for i := range vt {
- vt[i] = rnd.NormFloat64()
- }
-
- uAllOrig := make([]float64, len(u))
- copy(uAllOrig, u)
- vtAllOrig := make([]float64, len(vt))
- copy(vtAllOrig, vt)
- aCopy := make([]float64, len(a))
- copy(aCopy, a)
-
- s := make([]float64, min(m, n))
-
- work := make([]float64, 1)
- impl.Dgesvd(jobU, jobVT, m, n, a, lda, s, u, ldu, vt, ldvt, work, -1)
-
- if !floats.Equal(a, aCopy) {
- t.Errorf("a changed during call to get work length")
- }
-
- work = make([]float64, int(work[0]))
- impl.Dgesvd(jobU, jobVT, m, n, a, lda, s, u, ldu, vt, ldvt, work, len(work))
-
- errStr := fmt.Sprintf("m = %v, n = %v, lda = %v, ldu = %v, ldv = %v", m, n, lda, ldu, ldvt)
- svdCheck(t, false, errStr, m, n, s, a, u, ldu, vt, ldvt, aCopy, lda)
- svdCheckPartial(t, impl, lapack.SVDAll, errStr, uAllOrig, vtAllOrig, aCopy, m, n, a, lda, s, u, ldu, vt, ldvt, work, false)
-
- // Test InPlace
- jobU = lapack.SVDInPlace
- jobVT = lapack.SVDInPlace
- copy(a, aCopy)
- copy(u, uAllOrig)
- copy(vt, vtAllOrig)
-
- impl.Dgesvd(jobU, jobVT, m, n, a, lda, s, u, ldu, vt, ldvt, work, len(work))
- svdCheck(t, true, errStr, m, n, s, a, u, ldu, vt, ldvt, aCopy, lda)
- svdCheckPartial(t, impl, lapack.SVDInPlace, errStr, uAllOrig, vtAllOrig, aCopy, m, n, a, lda, s, u, ldu, vt, ldvt, work, false)
- }
-}
-
-// svdCheckPartial checks that the singular values and vectors are computed when
-// not all of them are computed.
-func svdCheckPartial(t *testing.T, impl Dgesvder, job lapack.SVDJob, errStr string, uAllOrig, vtAllOrig, aCopy []float64, m, n int, a []float64, lda int, s, u []float64, ldu int, vt []float64, ldvt int, work []float64, shortWork bool) {
- rnd := rand.New(rand.NewSource(1))
- jobU := job
- jobVT := job
- // Compare the singular values when computed with {SVDNone, SVDNone.}
- sCopy := make([]float64, len(s))
- copy(sCopy, s)
- copy(a, aCopy)
- for i := range s {
- s[i] = rnd.Float64()
- }
- tmp1 := make([]float64, 1)
- tmp2 := make([]float64, 1)
- jobU = lapack.SVDNone
- jobVT = lapack.SVDNone
-
- impl.Dgesvd(jobU, jobVT, m, n, a, lda, s, tmp1, ldu, tmp2, ldvt, work, -1)
- work = make([]float64, int(work[0]))
- lwork := len(work)
- if shortWork {
- lwork--
- }
- ok := impl.Dgesvd(jobU, jobVT, m, n, a, lda, s, tmp1, ldu, tmp2, ldvt, work, lwork)
- if !ok {
- t.Errorf("Dgesvd did not complete successfully")
- }
- if !floats.EqualApprox(s, sCopy, 1e-10) {
- t.Errorf("Singular value mismatch when singular vectors not computed: %s", errStr)
- }
- // Check that the singular vectors are correctly computed when the other
- // is none.
- uAll := make([]float64, len(u))
- copy(uAll, u)
- vtAll := make([]float64, len(vt))
- copy(vtAll, vt)
-
- // Copy the original vectors so the data outside the matrix bounds is the same.
- copy(u, uAllOrig)
- copy(vt, vtAllOrig)
-
- jobU = job
- jobVT = lapack.SVDNone
- copy(a, aCopy)
- for i := range s {
- s[i] = rnd.Float64()
- }
- impl.Dgesvd(jobU, jobVT, m, n, a, lda, s, u, ldu, tmp2, ldvt, work, -1)
- work = make([]float64, int(work[0]))
- lwork = len(work)
- if shortWork {
- lwork--
- }
- impl.Dgesvd(jobU, jobVT, m, n, a, lda, s, u, ldu, tmp2, ldvt, work, len(work))
- if !floats.EqualApprox(uAll, u, 1e-10) {
- t.Errorf("U mismatch when VT is not computed: %s", errStr)
- }
- if !floats.EqualApprox(s, sCopy, 1e-10) {
- t.Errorf("Singular value mismatch when U computed VT not")
- }
- jobU = lapack.SVDNone
- jobVT = job
- copy(a, aCopy)
- for i := range s {
- s[i] = rnd.Float64()
- }
- impl.Dgesvd(jobU, jobVT, m, n, a, lda, s, tmp1, ldu, vt, ldvt, work, -1)
- work = make([]float64, int(work[0]))
- lwork = len(work)
- if shortWork {
- lwork--
- }
- impl.Dgesvd(jobU, jobVT, m, n, a, lda, s, tmp1, ldu, vt, ldvt, work, len(work))
- if !floats.EqualApprox(vtAll, vt, 1e-10) {
- t.Errorf("VT mismatch when U is not computed: %s", errStr)
- }
- if !floats.EqualApprox(s, sCopy, 1e-10) {
- t.Errorf("Singular value mismatch when VT computed U not")
- }
-}
-
-// svdCheck checks that the singular value decomposition correctly multiplies back
-// to the original matrix.
-func svdCheck(t *testing.T, thin bool, errStr string, m, n int, s, a, u []float64, ldu int, vt []float64, ldvt int, aCopy []float64, lda int) {
- sigma := blas64.General{
- Rows: m,
- Cols: n,
- Stride: n,
- Data: make([]float64, m*n),
- }
- for i := 0; i < min(m, n); i++ {
- sigma.Data[i*sigma.Stride+i] = s[i]
- }
-
- uMat := blas64.General{
- Rows: m,
- Cols: m,
- Stride: ldu,
- Data: u,
- }
- vTMat := blas64.General{
- Rows: n,
- Cols: n,
- Stride: ldvt,
- Data: vt,
- }
- if thin {
- sigma.Rows = min(m, n)
- sigma.Cols = min(m, n)
- uMat.Cols = min(m, n)
- vTMat.Rows = min(m, n)
- }
-
- tmp := blas64.General{
- Rows: m,
- Cols: n,
- Stride: n,
- Data: make([]float64, m*n),
- }
- ans := blas64.General{
- Rows: m,
- Cols: n,
- Stride: lda,
- Data: make([]float64, m*lda),
- }
- copy(ans.Data, a)
-
- blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, uMat, sigma, 0, tmp)
- blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, tmp, vTMat, 0, ans)
-
- if !floats.EqualApprox(ans.Data, aCopy, 1e-8) {
- t.Errorf("Decomposition mismatch. Trim = %v, %s", thin, errStr)
- }
-
- if !thin {
- // Check that U and V are orthogonal.
- for i := 0; i < uMat.Rows; i++ {
- for j := i + 1; j < uMat.Rows; j++ {
- dot := blas64.Dot(uMat.Cols,
- blas64.Vector{Inc: 1, Data: uMat.Data[i*uMat.Stride:]},
- blas64.Vector{Inc: 1, Data: uMat.Data[j*uMat.Stride:]},
- )
- if dot > 1e-8 {
- t.Errorf("U not orthogonal %s", errStr)
- }
- }
- }
- for i := 0; i < vTMat.Rows; i++ {
- for j := i + 1; j < vTMat.Rows; j++ {
- dot := blas64.Dot(vTMat.Cols,
- blas64.Vector{Inc: 1, Data: vTMat.Data[i*vTMat.Stride:]},
- blas64.Vector{Inc: 1, Data: vTMat.Data[j*vTMat.Stride:]},
- )
- if dot > 1e-8 {
- t.Errorf("V not orthogonal %s", errStr)
- }
- }
- }
- }
-}
OSDN Git Service
