diff --git a/lapack/testlapack/dgeqp3.go b/lapack/testlapack/dgeqp3.go index 6c211b17..4779a285 100644 --- a/lapack/testlapack/dgeqp3.go +++ b/lapack/testlapack/dgeqp3.go @@ -5,12 +5,14 @@ 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/lapack" ) type Dgeqp3er interface { @@ -20,120 +22,95 @@ type Dgeqp3er interface { func Dgeqp3Test(t *testing.T, impl Dgeqp3er) { rnd := rand.New(rand.NewSource(1)) - for c, test := range []struct { - m, n, lda int - }{ - {1, 1, 0}, - {2, 2, 0}, - {3, 2, 0}, - {2, 3, 0}, - {1, 12, 0}, - {2, 6, 0}, - {3, 4, 0}, - {4, 3, 0}, - {6, 2, 0}, - {12, 1, 0}, - {1, 1, 20}, - {2, 2, 20}, - {3, 2, 20}, - {2, 3, 20}, - {1, 12, 20}, - {2, 6, 20}, - {3, 4, 20}, - {4, 3, 20}, - {6, 2, 20}, - {12, 1, 20}, - {129, 256, 0}, - {256, 129, 0}, - {129, 256, 266}, - {256, 129, 266}, - } { - n := test.n - m := test.m - lda := test.lda - if lda == 0 { - lda = test.n - } - const ( - all = iota - some - none - ) - for _, free := range []int{all, some, none} { - // Allocate m×n matrix A and fill it with random numbers. - a := make([]float64, m*lda) - for i := range a { - a[i] = rnd.Float64() - } - // Store a copy of A for later comparison. - aCopy := make([]float64, len(a)) - copy(aCopy, a) - // Allocate a slice of column pivots. - jpvt := make([]int, n) - for j := range jpvt { - switch free { - case all: - // All columns are free. - jpvt[j] = -1 - case some: - // Some columns are free, some are leading columns. - jpvt[j] = rnd.Intn(2) - 1 // -1 or 0 - case none: - // All columns are leading. - jpvt[j] = 0 - default: - panic("bad freedom") - } - } - // Allocate a slice for scalar factors of elementary - // reflectors and fill it with random numbers. Dgeqp3 - // will overwrite them with valid data. - k := min(m, n) - tau := make([]float64, k) - for i := range tau { - tau[i] = rnd.Float64() - } - // Get optimal workspace size for Dgeqp3. - work := make([]float64, 1) - impl.Dgeqp3(m, n, a, lda, jpvt, tau, work, -1) - lwork := int(work[0]) - work = make([]float64, lwork) - for i := range work { - work[i] = rnd.Float64() - } - - // Compute a QR factorization of A with column pivoting. - impl.Dgeqp3(m, n, a, lda, jpvt, tau, work, lwork) - - // Compute Q based on the elementary reflectors stored in A. - q := constructQ("QR", m, n, a, lda, tau) - // Check that Q is orthogonal. - if !isOrthogonal(q) { - t.Errorf("Case %v, Q not orthogonal", c) - } - - // Copy the upper triangle of A into R. - r := blas64.General{ - Rows: m, - Cols: n, - Stride: n, - Data: make([]float64, m*n), - } - for i := 0; i < m; i++ { - for j := i; j < n; j++ { - r.Data[i*n+j] = a[i*lda+j] - } - } - // Compute Q * R. - got := nanGeneral(m, n, lda) - blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, q, r, 0, got) - // Compute A * P: rearrange the columns of A based on the permutation in jpvt. - want := blas64.General{Rows: m, Cols: n, Stride: lda, Data: aCopy} - impl.Dlapmt(true, want.Rows, want.Cols, want.Data, want.Stride, jpvt) - // Check that A * P = Q * R. - if !equalApproxGeneral(got, want, 1e-13) { - t.Errorf("Case %v, Q*R != A*P\nQ*R=%v\nA*P=%v", c, got, want) + for _, m := range []int{0, 1, 2, 3, 4, 5, 12, 23, 129} { + for _, n := range []int{0, 1, 2, 3, 4, 5, 12, 23, 129} { + for _, lda := range []int{max(1, n), n + 3} { + dgeqp3Test(t, impl, rnd, m, n, lda) } } } } + +func dgeqp3Test(t *testing.T, impl Dgeqp3er, rnd *rand.Rand, m, n, lda int) { + const ( + tol = 1e-14 + + all = iota + some + none + ) + for _, free := range []int{all, some, none} { + name := fmt.Sprintf("m=%d,n=%d,lda=%d,", m, n, lda) + + // Allocate m×n matrix A and fill it with random numbers. + a := randomGeneral(m, n, lda, rnd) + // Store a copy of A for later comparison. + aCopy := cloneGeneral(a) + // Allocate a slice of column pivots. + jpvt := make([]int, n) + for j := range jpvt { + switch free { + case all: + // All columns are free. + jpvt[j] = -1 + name += "free=all" + case some: + // Some columns are free, some are leading columns. + jpvt[j] = rnd.Intn(2) - 1 // -1 or 0 + name += "free=some" + case none: + // All columns are leading. + jpvt[j] = 0 + name += "free=none" + default: + panic("bad freedom") + } + } + // Allocate a slice for scalar factors of elementary + // reflectors and fill it with random numbers. Dgeqp3 + // will overwrite them with valid data. + k := min(m, n) + tau := make([]float64, k) + for i := range tau { + tau[i] = rnd.Float64() + } + // Get optimal workspace size for Dgeqp3. + work := make([]float64, 1) + impl.Dgeqp3(m, n, a.Data, a.Stride, jpvt, tau, work, -1) + lwork := int(work[0]) + work = make([]float64, lwork) + for i := range work { + work[i] = rnd.Float64() + } + + // Compute a QR factorization of A with column pivoting. + impl.Dgeqp3(m, n, a.Data, a.Stride, jpvt, tau, work, lwork) + + // Compute Q based on the elementary reflectors stored in A. + q := constructQ("QR", m, n, a.Data, a.Stride, tau) + // Check that Q is orthogonal. + if resid := residualOrthogonal(q, false); resid > tol*float64(max(m, n)) { + t.Errorf("Case %v: Q not orthogonal; resid=%v, want<=%v", name, resid, tol*float64(max(m, n))) + } + + // Copy the upper triangle of A into R. + r := zeros(m, n, lda) + for i := 0; i < m; i++ { + for j := i; j < n; j++ { + r.Data[i*r.Stride+j] = a.Data[i*a.Stride+j] + } + } + // Compute Q*R - A*P: + // 1. Rearrange the columns of A based on the permutation in jpvt. + qrap := cloneGeneral(aCopy) + impl.Dlapmt(true, qrap.Rows, qrap.Cols, qrap.Data, qrap.Stride, jpvt) + // Compute Q*R - A*P. + blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, q, r, -1, qrap) + // Check that |Q*R - A*P| is small. + resid := dlange(lapack.MaxColumnSum, qrap.Rows, qrap.Cols, qrap.Data, qrap.Stride) + if resid > tol*float64(max(m, n)) { + t.Errorf("Case %v: |Q*R - A*P|=%v, want<=%v", name, resid, tol*float64(max(m, n))) + + } + } +} diff --git a/lapack/testlapack/general.go b/lapack/testlapack/general.go index 18cc6f11..fd5baa26 100644 --- a/lapack/testlapack/general.go +++ b/lapack/testlapack/general.go @@ -582,7 +582,7 @@ func constructQK(kind string, m, n, k int, a []float64, lda int, tau []float64) q := blas64.General{ Rows: sz, Cols: sz, - Stride: sz, + Stride: max(1, sz), Data: make([]float64, sz*sz), } for i := 0; i < sz; i++ { @@ -598,7 +598,7 @@ func constructQK(kind string, m, n, k int, a []float64, lda int, tau []float64) h := blas64.General{ Rows: sz, Cols: sz, - Stride: sz, + Stride: max(1, sz), Data: make([]float64, sz*sz), } for j := 0; j < sz; j++ {