Add Dorg2r and test

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

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})
 }

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)
+		}
+	}
+}

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
+		}
+	}
+}

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
+		}
+	}
+}

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"

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)
 }

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")
+		}
+	}
+}

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")
+		}
+	}
+}