Correct algorithm

cgo/lapack.go

@@ -30,6 +30,7 @@ const (
 	badWorkStride = "lapack: insufficient working array stride"
 	kGTM          = "lapack: k > m"
 	kGTN          = "lapack: k > n"
+	kLT0          = "lapack: k < 0"
 	mLTN          = "lapack: m < n"
 	negDimension  = "lapack: negative matrix dimension"
 	nLT0          = "lapack: n < 0"
@@ -479,20 +480,46 @@ func (impl Implementation) Dorglq(m, n, k int, a []float64, lda int, tau, work [
 	clapack.Dorglq(m, n, k, a, lda, tau)
 }
 
+// Dorgqr 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)
+// Dorgqr is the blocked version of dorg2r that makes greater use of level-3 BLAS
+// routines.
+//
+// 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 and 0 <= n <= m.
+//
+// Work is temporary storage, and lwork specifies the usable memory length.
+// The C interface does not support providing temporary storage. To provide compatibility
+// with native, lwork == -1 will not run Dorgqr but will instead write the minimum
+// work necessary to work[0]. If len(work) < lwork, Dgeqrf will panic, and at minimum
+// lwork >= n.
+//
+// Dorgqr will panic if the conditions on input values are not met.
 func (impl Implementation) Dorgqr(m, n, k int, a []float64, lda int, tau, work []float64, lwork int) {
-	checkMatrix(m, n, a, lda)
-	if len(tau) < k {
-		panic(badTau)
+	if lwork == -1 {
+		work[0] = float64(n)
+		return
 	}
-	if len(work) < n {
-		panic(badWork)
+	checkMatrix(m, n, a, lda)
+	if k < 0 {
+		panic(kLT0)
 	}
 	if k > n {
 		panic(kGTN)
 	}
 	if n > m {
 		panic(mLTN)
 	}
+	if len(tau) < k {
+		panic(badTau)
+	}
+	if len(work) < lwork {
+		panic(shortWork)
+	}
+	if lwork < n {
+		panic(badWork)
+	}
 	clapack.Dorgqr(m, n, k, a, lda, tau)
 }
 

native/dorgqr.go

@@ -20,13 +20,15 @@ import (
 // Dorgqr will panic if the conditions on input values are not met.
 func (impl Implementation) Dorgqr(m, n, k int, a []float64, lda int, tau, work []float64, lwork int) {
 	nb := impl.Ilaenv(1, "DORGQR", " ", m, n, k, -1)
-	nb = 2
 	// work is treated as an n×nb matrix
 	if lwork == -1 {
-		work[0] = float64(n * nb)
+		work[0] = float64(max(1, n) * nb)
 		return
 	}
 	checkMatrix(m, n, a, lda)
+	if k < 0 {
+		panic(kLT0)
+	}
 	if k > n {
 		panic(kGTN)
 	}
@@ -39,6 +41,9 @@ func (impl Implementation) Dorgqr(m, n, k int, a []float64, lda int, tau, work [
 	if len(work) < lwork {
 		panic(shortWork)
 	}
+	if lwork < n {
+		panic(badWork)
+	}
 	if n == 0 {
 		return
 	}
@@ -48,7 +53,6 @@ func (impl Implementation) Dorgqr(m, n, k int, a []float64, lda int, tau, work [
 	var ldwork int
 	if nb > 1 && nb < k {
 		nx = max(0, impl.Ilaenv(3, "DORGQR", " ", m, n, k, -1))
-		nx = 2
 		if nx < k {
 			ldwork = nb
 			iws = n * ldwork
@@ -83,8 +87,18 @@ func (impl Implementation) Dorgqr(m, n, k int, a []float64, lda int, tau, work [
 	for i := ki; i >= 0; i -= nb {
 		ib := min(nb, k-i)
 		if i+ib < n {
-			impl.Dlarft(lapack.Forward, lapack.ColumnWise, m-i, ib, a[i*lda+i:], lda, tau[i:], work, ldwork)
-			impl.Dlarfb(blas.Left, blas.NoTrans, lapack.Forward, lapack.ColumnWise, m-i, n-i-ib, ib, a[i*lda+i:], lda, work, ldwork, a[i*lda+i+ib:], lda, work[ib:], ldwork)
+			impl.Dlarft(lapack.Forward, lapack.ColumnWise,
+				m-i, ib,
+				a[i*lda+i:], lda,
+				tau[i:],
+				work, ldwork)
+
+			impl.Dlarfb(blas.Left, blas.NoTrans, lapack.Forward, lapack.ColumnWise,
+				m-i, n-i-ib, ib,
+				a[i*lda+i:], lda,
+				work, ldwork,
+				a[i*lda+i+ib:], lda,
+				work[ib*ldwork:], ldwork)
 		}
 		impl.Dorg2r(m-i, ib, ib, a[i*lda+i:], lda, tau[i:], work)
 		// Set rows 0:i-1 of current block to zero

native/general.go

@@ -36,6 +36,7 @@ const (
 	badWorkStride = "lapack: insufficient working array stride"
 	kGTM          = "lapack: k > m"
 	kGTN          = "lapack: k > n"
+	kLT0          = "lapack: k < 0"
 	mLTN          = "lapack: m < n"
 	negDimension  = "lapack: negative matrix dimension"
 	nLT0          = "lapack: n < 0"

testlapack/dorgqr.go

@@ -1,7 +1,6 @@
 package testlapack
 
 import (
-	"fmt"
 	"math"
 	"math/rand"
 	"testing"
@@ -17,20 +16,28 @@ type Dorgqrer interface {
 func DorgqrTest(t *testing.T, impl Dorgqrer) {
 	// TODO(btracey): Base tests off of nb and nx.
 	for _, test := range []struct{ m, n, k, lda int }{
-		/*
-			{10, 10, 10, 0},
-			{10, 10, 10, 20},
-			{30, 10, 10, 0},
-			{30, 20, 10, 0},
-			{100, 100, 100, 0},
-		*/
+		{10, 10, 10, 0},
+		{10, 10, 10, 20},
+		{30, 10, 10, 0},
+		{30, 20, 10, 0},
 
-		//{400, 400, 200, 500},
-		{4, 4, 4, 10},
-		//{8, 8, 3, 0},
-		//{4, 4, 1, 0},
-		//{11, 11, 5, 0},
-		//{2, 2, 2, 0},
+		{100, 100, 100, 0},
+		{100, 100, 50, 0},
+		{130, 100, 100, 0},
+		{130, 100, 50, 0},
+		{100, 100, 100, 150},
+		{100, 100, 50, 150},
+		{130, 100, 100, 150},
+		{130, 100, 50, 150},
+
+		{200, 200, 200, 0},
+		{200, 200, 150, 0},
+		{230, 200, 200, 0},
+		{230, 200, 150, 0},
+		{200, 200, 200, 250},
+		{200, 200, 150, 250},
+		{230, 200, 200, 250},
+		{230, 200, 150, 250},
 	} {
 		m := test.m
 		n := test.n
@@ -66,10 +73,6 @@ func DorgqrTest(t *testing.T, impl Dorgqrer) {
 		impl.Dorgqr(m, n, k, a, lda, tau, work, len(work))
 		if !floats.EqualApprox(a, aUnblocked, 1e-10) {
 			t.Errorf("Q Mismatch. m = %d, n = %d, k = %d, lda = %d", m, n, k, lda)
-			fmt.Println("orgqr")
-			printRowise(a, m, n, lda, false)
-			fmt.Println("org2r")
-			printRowise(aUnblocked, m, n, lda, false)
 		}
 	}
 }