Adddorgbr

cgo/lapack.go

@@ -17,6 +17,7 @@ const (
 	badD            = "lapack: d has insufficient length"
 	badDecompUpdate = "lapack: bad decomp update"
 	badDiag         = "lapack: bad diag"
+	badDims         = "lapack: bad input dimensions"
 	badDirect       = "lapack: bad direct"
 	badE            = "lapack: e has insufficient length"
 	badIpiv         = "lapack: insufficient permutation length"
@@ -542,6 +543,47 @@ func (impl Implementation) Dgetrs(trans blas.Transpose, n, nrhs int, a []float64
 	clapack.Dgetrs(trans, n, nrhs, a, lda, ipiv32, b, ldb)
 }
 
+// Dorgbr generates one of the matrices Q or P^T computed by Dgebrd
+// computed from the decomposition Dgebrd. See Dgebd2 for the description of
+// Q and P^T.
+//
+// If vect == lapack.ApplyQ, then a is assumed to have been an m×k matrix and
+// Q is of order m. If m >= k, then Dorgbr returns the first n columns of Q
+// where m >= n >= k. If m < k, then Dorgbr returns Q as an m×m matrix.
+//
+// If vect == lapack.ApplyP, then A is assumed to have been a k×n matrix, and
+// P^T is of order n. If k < n, then Dorgbr returns the first m rows of P^T,
+// where n >= m >= k. If k >= n, then Dorgbr returns P^T as an n×n matrix.
+func (impl Implementation) Dorgbr(vect lapack.DecompUpdate, m, n, k int, a []float64, lda int, tau, work []float64, lwork int) {
+	mn := min(m, n)
+	wantq := vect == lapack.ApplyQ
+	if wantq {
+		if m < n || n < min(m, k) || m < min(m, k) {
+			panic(badDims)
+		}
+	} else {
+		if n < m || m < min(n, k) || n < min(n, k) {
+			panic(badDims)
+		}
+	}
+	if wantq {
+		checkMatrix(m, k, a, lda)
+	} else {
+		checkMatrix(k, n, a, lda)
+	}
+	if lwork == -1 {
+		work[0] = float64(mn)
+		return
+	}
+	if len(work) < lwork {
+		panic(badWork)
+	}
+	if lwork < mn {
+		panic(badWork)
+	}
+	clapack.Dorgbr(byte(vect), m, n, k, a, lda, tau)
+}
+
 // Dorglq generates an m×n matrix Q with orthonormal rows defined by the
 // product of elementary reflectors as computed by Dgelqf.
 //  Q = H(0) * H(2) * ... * H(k-1)

cgo/lapack_test.go

@@ -127,6 +127,10 @@ func TestDormbr(t *testing.T) {
 }
 */
 
+func TestDorgbr(t *testing.T) {
+	testlapack.DorgbrTest(t, blockedTranslate{impl})
+}
+
 func TestDormqr(t *testing.T) {
 	testlapack.Dorm2rTest(t, blockedTranslate{impl})
 }

native/dorgbr.go

@@ -0,0 +1,114 @@
+// 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/lapack"
+
+// Dorgbr generates one of the matrices Q or P^T computed by Dgebrd
+// computed from the decomposition Dgebrd. See Dgebd2 for the description of
+// Q and P^T.
+//
+// If vect == lapack.ApplyQ, then a is assumed to have been an m×k matrix and
+// Q is of order m. If m >= k, then Dorgbr returns the first n columns of Q
+// where m >= n >= k. If m < k, then Dorgbr returns Q as an m×m matrix.
+//
+// If vect == lapack.ApplyP, then A is assumed to have been a k×n matrix, and
+// P^T is of order n. If k < n, then Dorgbr returns the first m rows of P^T,
+// where n >= m >= k. If k >= n, then Dorgbr returns P^T as an n×n matrix.
+func (impl Implementation) Dorgbr(vect lapack.DecompUpdate, m, n, k int, a []float64, lda int, tau, work []float64, lwork int) {
+	mn := min(m, n)
+	wantq := vect == lapack.ApplyQ
+	if wantq {
+		if m < n || n < min(m, k) || m < min(m, k) {
+			panic(badDims)
+		}
+	} else {
+		if n < m || m < min(n, k) || n < min(n, k) {
+			panic(badDims)
+		}
+	}
+	if wantq {
+		checkMatrix(m, k, a, lda)
+	} else {
+		checkMatrix(k, n, a, lda)
+	}
+	work[0] = 1
+	if wantq {
+		if m >= k {
+			impl.Dorgqr(m, n, k, a, lda, tau, work, -1)
+		} else if m > 1 {
+			impl.Dorgqr(m-1, m-1, m-1, a[lda+1:], lda, tau, work, -1)
+		}
+	} else {
+		if k < n {
+			impl.Dorglq(m, n, k, a, lda, tau, work, -1)
+		} else if n > 1 {
+			impl.Dorglq(n-1, n-1, n-1, a[lda+1:], lda, tau, work, -1)
+		}
+	}
+	lworkopt := int(work[0])
+	lworkopt = max(lworkopt, mn)
+	if lwork == -1 {
+		work[0] = float64(lworkopt)
+		return
+	}
+	if len(work) < lwork {
+		panic(badWork)
+	}
+	if lwork < mn {
+		panic(badWork)
+	}
+	if m == 0 || n == 0 {
+		work[0] = 1
+		return
+	}
+	if wantq {
+		// Form Q, determined by a call to Dgebrd to reduce an m×k matrix.
+		if m >= k {
+			impl.Dorgqr(m, n, k, a, lda, tau, work, lwork)
+		} else {
+			// Shift the vectors which define the elementary reflectors one
+			// column to the right, and set the first row and column of Q to
+			// those of the unit matrix.
+			for j := m - 2; j >= 0; j-- {
+				a[j] = 0
+				for i := j; i < m; i++ {
+					a[i*lda+j] = a[i*lda+j-1]
+				}
+			}
+			a[0] = 1
+			for i := 1; i < m; i++ {
+				a[i*lda] = 0
+			}
+			if m > 1 {
+				// Form Q[1:m-1, 1:m-1]
+				impl.Dorgqr(m-1, m-1, m-1, a[lda+1:], lda, tau, work, lwork)
+			}
+		}
+	} else {
+		// Form P^T, determined by a call to Dgebrd to reduce a k×n matrix.
+		if k < n {
+			impl.Dorglq(m, n, k, a, lda, tau, work, lwork)
+		} else {
+			// Shift the vectors which define the elementary reflectors one
+			// row downward, and set the first row and column of P^T to
+			// those of the unit matrix.
+			a[0] = 1
+			for i := 1; i < n; i++ {
+				a[i*lda] = 0
+			}
+			for j := 1; j < n; j++ {
+				for i := j - 1; i >= 1; i-- {
+					a[i*lda+j] = a[(i-1)*lda+j]
+				}
+				a[j] = 0
+			}
+			if n > 1 {
+				impl.Dorglq(n-1, n-1, n-1, a[lda+1:], lda, tau, work, lwork)
+			}
+		}
+	}
+	work[0] = float64(lworkopt)
+}

native/dorgl2.go

@@ -35,7 +35,7 @@ func (impl Implementation) Dorgl2(m, n, k int, a []float64, lda int, tau, work [
 		return
 	}
 	bi := blas64.Implementation()
-	if k < m-1 {
+	if k < m {
 		for i := k; i < m; i++ {
 			for j := 0; j < n; j++ {
 				a[i*lda+j] = 0

native/general.go

@@ -23,6 +23,7 @@ const (
 	badD            = "lapack: d has insufficient length"
 	badDecompUpdate = "lapack: bad decomp update"
 	badDiag         = "lapack: bad diag"
+	badDims         = "lapack: bad input dimensions"
 	badDirect       = "lapack: bad direct"
 	badE            = "lapack: e has insufficient length"
 	badIpiv         = "lapack: insufficient permutation length"

native/lapack_test.go

@@ -140,6 +140,10 @@ func TestDorg2r(t *testing.T) {
 	testlapack.Dorg2rTest(t, impl)
 }
 
+func TestDorgbr(t *testing.T) {
+	testlapack.DorgbrTest(t, impl)
+}
+
 func TestDorgl2(t *testing.T) {
 	testlapack.Dorgl2Test(t, impl)
 }

testlapack/dorgbr.go

@@ -0,0 +1,145 @@
+// 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/blas/blas64"
+	"github.com/gonum/floats"
+	"github.com/gonum/lapack"
+)
+
+type Dorgbrer interface {
+	Dorgbr(vect lapack.DecompUpdate, m, n, k int, a []float64, lda int, tau, work []float64, lwork int)
+	Dgebrder
+}
+
+func DorgbrTest(t *testing.T, impl Dorgbrer) {
+	for _, vect := range []lapack.DecompUpdate{lapack.ApplyQ, lapack.ApplyP} {
+		for _, test := range []struct {
+			m, n, k, lda int
+		}{
+			{5, 5, 5, 0},
+			{5, 5, 3, 0},
+			{5, 3, 5, 0},
+			{3, 5, 5, 0},
+			{3, 4, 5, 0},
+			{3, 5, 4, 0},
+			{4, 3, 5, 0},
+			{4, 5, 3, 0},
+			{5, 3, 4, 0},
+			{5, 4, 3, 0},
+
+			{5, 5, 5, 10},
+			{5, 5, 3, 10},
+			{5, 3, 5, 10},
+			{3, 5, 5, 10},
+			{3, 4, 5, 10},
+			{3, 5, 4, 10},
+			{4, 3, 5, 10},
+			{4, 5, 3, 10},
+			{5, 3, 4, 10},
+			{5, 4, 3, 10},
+		} {
+			m := test.m
+			n := test.n
+			k := test.k
+			lda := test.lda
+			// Filter out bad tests
+			if vect == lapack.ApplyQ {
+				if m < n || n < min(m, k) || m < min(m, k) {
+					continue
+				}
+			} else {
+				if n < m || m < min(n, k) || n < min(n, k) {
+					continue
+				}
+			}
+			// Sizes for Dorgbr.
+			var ma, na int
+			if vect == lapack.ApplyQ {
+				ma = m
+				na = k
+			} else {
+				ma = k
+				na = n
+			}
+			// a eventually needs to store either P or Q, so it must be
+			// sufficiently big.
+			var a []float64
+			if vect == lapack.ApplyQ {
+				lda = max(m, lda)
+				a = make([]float64, m*lda)
+			} else {
+				lda = max(n, lda)
+				a = make([]float64, n*lda)
+			}
+			for i := range a {
+				a[i] = rand.NormFloat64()
+			}
+
+			nTau := min(ma, na)
+			tauP := make([]float64, nTau)
+			tauQ := make([]float64, nTau)
+			d := make([]float64, nTau)
+			e := make([]float64, nTau)
+			lwork := -1
+			work := make([]float64, 1)
+			impl.Dgebrd(ma, na, a, lda, d, e, tauQ, tauP, work, lwork)
+			work = make([]float64, int(work[0]))
+			lwork = len(work)
+			impl.Dgebrd(ma, na, a, lda, d, e, tauQ, tauP, work, lwork)
+
+			aCopy := make([]float64, len(a))
+			copy(aCopy, a)
+
+			var tau []float64
+			if vect == lapack.ApplyQ {
+				tau = tauQ
+			} else {
+				tau = tauP
+			}
+
+			impl.Dorgbr(vect, m, n, k, a, lda, tau, work, -1)
+			work = make([]float64, int(work[0]))
+			lwork = len(work)
+			impl.Dorgbr(vect, m, n, k, a, lda, tau, work, lwork)
+
+			var ans blas64.General
+			var nRows, nCols int
+			equal := true
+			if vect == lapack.ApplyQ {
+				nRows = m
+				nCols = m
+				if m >= k {
+					nCols = n
+				}
+				ans = constructQPBidiagonal(vect, ma, na, min(m, k), aCopy, lda, tau)
+			} else {
+				nRows = n
+				if k < n {
+					nRows = m
+				}
+				nCols = n
+				ansTmp := constructQPBidiagonal(vect, ma, na, min(k, n), aCopy, lda, tau)
+				// Dorgbr actually computes P^T
+				ans = transposeGeneral(ansTmp)
+			}
+			for i := 0; i < nRows; i++ {
+				for j := 0; j < nCols; j++ {
+					if !floats.EqualWithinAbsOrRel(a[i*lda+j], ans.Data[i*ans.Stride+j], 1e-8, 1e-8) {
+						equal = false
+					}
+				}
+			}
+			if !equal {
+				applyQ := vect == lapack.ApplyQ
+				t.Errorf("Extracted matrix mismatch. applyQ: %v, m = %v, n = %v, k = %v", applyQ, m, n, k)
+			}
+		}
+	}
+}

testlapack/dorml2.go

@@ -19,6 +19,9 @@ type Dorml2er interface {
 }
 
 func Dorml2Test(t *testing.T, impl Dorml2er) {
+	// TODO(btracey): This test is not complete, because it
+	// doesn't test individual values of m, n, and k, instead only testing
+	// a specific subset of possible k values.
 	for _, side := range []blas.Side{blas.Left, blas.Right} {
 		for _, trans := range []blas.Transpose{blas.NoTrans, blas.Trans} {
 			for _, test := range []struct {

testlapack/general.go

@@ -37,6 +37,23 @@ func nanSlice(n int) []float64 {
 	return s
 }
 
+// transposeGeneral returns a new general matrix that is the transpose of the
+// input. Nothing is done with data outside the {rows, cols} limit of the general.
+func transposeGeneral(a blas64.General) blas64.General {
+	ans := blas64.General{
+		Rows:   a.Cols,
+		Cols:   a.Rows,
+		Stride: a.Rows,
+		Data:   make([]float64, a.Cols*a.Rows),
+	}
+	for i := 0; i < a.Rows; i++ {
+		for j := 0; j < a.Cols; j++ {
+			ans.Data[j*ans.Stride+i] = a.Data[i*a.Stride+j]
+		}
+	}
+	return ans
+}
+
 // extractVMat collects the single reflectors from a into a matrix.
 func extractVMat(m, n int, a []float64, lda int, direct lapack.Direct, store lapack.StoreV) blas64.General {
 	k := min(m, n)