diff --git a/cgo/lapack.go b/cgo/lapack.go index 2e051ed5..2b11de89 100644 --- a/cgo/lapack.go +++ b/cgo/lapack.go @@ -28,6 +28,9 @@ const ( badUplo = "lapack: illegal triangle" badWork = "lapack: insufficient working memory" badWorkStride = "lapack: insufficient working array stride" + kGTM = "lapack: k > m" + kGTN = "lapack: k > n" + mLTN = "lapack: m < n" negDimension = "lapack: negative matrix dimension" nLT0 = "lapack: n < 0" shortWork = "lapack: working array shorter than declared" @@ -462,6 +465,37 @@ func (impl Implementation) Dgetrs(trans blas.Transpose, n, nrhs int, a []float64 clapack.Dgetrs(trans, n, nrhs, a, lda, ipiv32, b, ldb) } +func (impl Implementation) Dorglq(m, n, k int, a []float64, lda int, tau, work []float64) { + checkMatrix(m, n, a, lda) + if len(tau) < k { + panic(badTau) + } + if k > m { + panic(kGTM) + } + if k > m { + panic(kGTM) + } + clapack.Dorglq(m, n, k, a, lda, tau) +} + +func (impl Implementation) Dorgqr(m, n, k int, a []float64, lda int, tau, work []float64) { + checkMatrix(m, n, a, lda) + if len(tau) < k { + panic(badTau) + } + if len(work) < n { + panic(badWork) + } + if k > n { + panic(kGTN) + } + if n > m { + panic(mLTN) + } + clapack.Dorgqr(m, n, k, a, lda, tau) +} + // Dormlq multiplies the matrix C by the othogonal matrix Q defined by the // slices a and tau. A and tau are as returned from Dgelqf. // C = Q * C if side == blas.Left and trans == blas.NoTrans diff --git a/cgo/lapack_test.go b/cgo/lapack_test.go index 6ab342ba..5934c7eb 100644 --- a/cgo/lapack_test.go +++ b/cgo/lapack_test.go @@ -84,6 +84,22 @@ func (d blockedTranslate) Dorml2(side blas.Side, trans blas.Transpose, m, n, k i impl.Dormlq(side, trans, m, n, k, a, lda, tau, c, ldc, work, len(work)) } +func (d blockedTranslate) Dorg2r(m, n, k int, a []float64, lda int, tau, work []float64) { + impl.Dorgqr(m, n, k, a, lda, tau, work) +} + +func (d blockedTranslate) Dorgl2(m, n, k int, a []float64, lda int, tau, work []float64) { + impl.Dorglq(m, n, k, a, lda, tau, work) +} + +func TestDorglq(t *testing.T) { + testlapack.Dorgl2Test(t, blockedTranslate{impl}) +} + +func TestDorgqr(t *testing.T) { + testlapack.Dorg2rTest(t, blockedTranslate{impl}) +} + func TestDormqr(t *testing.T) { testlapack.Dorm2rTest(t, blockedTranslate{impl}) } diff --git a/native/dgebd2.go b/native/dgebd2.go new file mode 100644 index 00000000..4f84ee5b --- /dev/null +++ b/native/dgebd2.go @@ -0,0 +1,26 @@ +// 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 native + +// Dgebd2 reduces an m×n matrix A to upper or lower bidiagonal form by an orthogonal +// transformation. +// Q^T * A * P = B +// if m >= n, B is upper diagonal, otherwise B is lower bidiagonal. +// d is the diagonal, len = min(m,n) +// e is the off-diagonal len = min(m,n)-1 +func (impl Implementation) Dgebd2(m, n int, a []float64, lda int, d, e, tauQ, tauP, work []float64) { + checkMatrix(m, n, a, lda) + if len(d) < min(m, n) { + panic("lapack: insufficient d") + } + if len(e) < min(m, n)-1 { + panic("lapack: insufficient e") + } + if m > n { + for i := 0; i < n; i++ { + impl.Dlarfg(m-i, a[i*lda+i], a[min(i+1, m-1)*lda+i:], 1) + } + } +} diff --git a/native/dorg2r.go b/native/dorg2r.go new file mode 100644 index 00000000..7342ddc9 --- /dev/null +++ b/native/dorg2r.go @@ -0,0 +1,61 @@ +// 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 native + +import ( + "github.com/gonum/blas" + "github.com/gonum/blas/blas64" +) + +// Dorg2r generates an m×n matrix Q with orthonormal columns defined by the +// product of elementary reflectors as computed by Dgeqrf. +// Q = H(0) * H(2) * ... * H(k-1) +// The length of tau must be equal to k, and the length of work must be at least n. +// It also must be that 0 <= k <= n. m >= n >= 0. Dorg2r will panic if these +// conditions are not met. +func (impl Implementation) Dorg2r(m, n, k int, a []float64, lda int, tau []float64, work []float64) { + checkMatrix(m, n, a, lda) + if len(tau) < k { + panic(badTau) + } + if len(work) < n { + panic(badWork) + } + if k > n { + panic(kGTN) + } + if n > m { + panic(mLTN) + } + if n == 0 { + return + } + bi := blas64.Implementation() + // Initialize columns k+1:n to columns of the unit matrix. + for l := 0; l < m; l++ { + for j := k; j < n; j++ { + a[l*lda+j] = 0 + } + } + for j := k; j < n; j++ { + a[j*lda+j] = 1 + } + for i := k - 1; i >= 0; i-- { + for i := range work { + work[i] = 0 + } + if i < n-1 { + a[i*lda+i] = 1 + impl.Dlarf(blas.Left, m-i, n-i-1, a[i*lda+i:], lda, tau[i], a[i*lda+i+1:], lda, work) + } + if i < m-1 { + bi.Dscal(m-i-1, -tau[i], a[(i+1)*lda+i:], lda) + } + a[i*lda+i] = 1 - tau[i] + for l := 0; l < i; l++ { + a[l*lda+i] = 0 + } + } +} diff --git a/native/dorgl2.go b/native/dorgl2.go new file mode 100644 index 00000000..e64f4033 --- /dev/null +++ b/native/dorgl2.go @@ -0,0 +1,56 @@ +// 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 native + +import ( + "github.com/gonum/blas" + "github.com/gonum/blas/blas64" +) + +// Dorgl2 generates an m×n matrix Q with orthonormal rows defined by the +// first m rows product of elementary reflectors as computed by Dgelqf. +// Q = H(0) * H(2) * ... * H(k-1) +// The length of tau must be equal to k, and the length of work must be at least n. +// It also must be that 0 <= k <= m. n >= m >= 0. Dorgl2 will panic if these +// conditions are not met. +func (impl Implementation) Dorgl2(m, n, k int, a []float64, lda int, tau, work []float64) { + checkMatrix(m, n, a, lda) + if len(tau) < k { + panic(badTau) + } + if k > m { + panic(kGTM) + } + if k > m { + panic(kGTM) + } + if m == 0 { + return + } + bi := blas64.Implementation() + if k < m-1 { + for i := k; i < m; i++ { + for j := 0; j < n; j++ { + a[i*lda+j] = 0 + } + } + for j := k; j < m; j++ { + a[j*lda+j] = 1 + } + } + for i := k - 1; i >= 0; i-- { + if i < n-1 { + if i < m-1 { + a[i*lda+i] = 1 + impl.Dlarf(blas.Right, m-i-1, n-i, a[i*lda+i:], 1, tau[i], a[(i+1)*lda+i:], lda, work) + } + bi.Dscal(n-i-1, -tau[i], a[i*lda+i+1:], 1) + } + a[i*lda+i] = 1 - tau[i] + for l := 0; l < i; l++ { + a[i*lda+l] = 0 + } + } +} diff --git a/native/general.go b/native/general.go index eefb6f3d..d2bb83a3 100644 --- a/native/general.go +++ b/native/general.go @@ -34,6 +34,9 @@ const ( badUplo = "lapack: illegal triangle" badWork = "lapack: insufficient working memory" badWorkStride = "lapack: insufficient working array stride" + kGTM = "lapack: k > m" + kGTN = "lapack: k > n" + mLTN = "lapack: m < n" negDimension = "lapack: negative matrix dimension" nLT0 = "lapack: n < 0" shortWork = "lapack: working array shorter than declared" diff --git a/native/lapack_test.go b/native/lapack_test.go index a6d54459..929179e7 100644 --- a/native/lapack_test.go +++ b/native/lapack_test.go @@ -120,6 +120,14 @@ func TestDorm2r(t *testing.T) { testlapack.Dorm2rTest(t, impl) } +func TestDorg2r(t *testing.T) { + testlapack.Dorg2rTest(t, impl) +} + +func TestDorgl2(t *testing.T) { + testlapack.Dorgl2Test(t, impl) +} + func TestDpocon(t *testing.T) { testlapack.DpoconTest(t, impl) } diff --git a/testlapack/dorg2r.go b/testlapack/dorg2r.go new file mode 100644 index 00000000..a7a951f8 --- /dev/null +++ b/testlapack/dorg2r.go @@ -0,0 +1,70 @@ +// 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" + "math/rand" + "testing" + + "github.com/gonum/floats" +) + +type Dorg2rer interface { + Dgeqrfer + Dorg2r(m, n, k int, a []float64, lda int, tau []float64, work []float64) +} + +func Dorg2rTest(t *testing.T, impl Dorg2rer) { + for _, test := range []struct { + m, n, lda int + }{ + {3, 3, 0}, + {4, 3, 0}, + + {5, 5, 20}, + {10, 5, 20}, + } { + m := test.m + n := test.n + lda := test.lda + if lda == 0 { + lda = test.n + } + a := make([]float64, m*lda) + for i := range a { + a[i] = rand.NormFloat64() + } + k := min(m, n) + tau := make([]float64, k) + work := make([]float64, 1) + impl.Dgeqrf(m, n, a, lda, tau, work, -1) + work = make([]float64, int(work[0])) + impl.Dgeqrf(m, n, a, lda, tau, work, len(work)) + + q := constructQ("QR", m, n, a, lda, tau) + + impl.Dorg2r(m, n, k, a, lda, tau, work) + + // Check that the first n columns match. + same := true + for i := 0; i < m; i++ { + for j := 0; j < n; j++ { + if !floats.EqualWithinAbsOrRel(q.Data[i*q.Stride+j], a[i*lda+j], 1e-12, 1e-12) { + same = false + break + } + } + } + if !same { + fmt.Println() + fmt.Println("a =") + printRowise(a, m, n, lda, false) + fmt.Println("q =") + printRowise(q.Data, q.Rows, q.Cols, q.Stride, false) + t.Errorf("Q mismatch") + } + } +} diff --git a/testlapack/dorgl2.go b/testlapack/dorgl2.go new file mode 100644 index 00000000..be07b676 --- /dev/null +++ b/testlapack/dorgl2.go @@ -0,0 +1,64 @@ +// 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 ( + "math/rand" + "testing" + + "github.com/gonum/floats" +) + +type Dorgl2er interface { + Dgelqfer + Dorgl2(m, n, k int, a []float64, lda int, tau []float64, work []float64) +} + +func Dorgl2Test(t *testing.T, impl Dorgl2er) { + for _, test := range []struct { + m, n, lda int + }{ + {3, 3, 0}, + {3, 4, 0}, + + {5, 5, 20}, + {5, 10, 20}, + } { + m := test.m + n := test.n + lda := test.lda + if lda == 0 { + lda = test.n + } + a := make([]float64, m*lda) + for i := range a { + a[i] = rand.NormFloat64() + } + k := min(m, n) + tau := make([]float64, k) + work := make([]float64, 1) + impl.Dgelqf(m, n, a, lda, tau, work, -1) + work = make([]float64, int(work[0])) + impl.Dgelqf(m, n, a, lda, tau, work, len(work)) + + q := constructQ("LQ", m, n, a, lda, tau) + + impl.Dorgl2(m, n, k, a, lda, tau, work) + + // Check that the first m rows match. + same := true + for i := 0; i < m; i++ { + for j := 0; j < n; j++ { + if !floats.EqualWithinAbsOrRel(q.Data[i*q.Stride+j], a[i*lda+j], 1e-12, 1e-12) { + same = false + break + } + } + } + if !same { + t.Errorf("Q mismatch") + } + } +}