lapack/{gonum,testlapack}: clean up Dorgr2 and its test

lapack/gonum/dorgr2.go

@@ -9,39 +9,34 @@ import (
 	"gonum.org/v1/gonum/blas/blas64"
 )
 
-// Dorgr2 generates an m×n real matrix Q with orthonormal rows,
+// Dorgr2 generates an m×n real matrix Q with orthonormal rows, which is defined
-// which is defined as the last m rows of a product of k elementary
+// as the last m rows of a product of k elementary reflectors of order n
-// reflectors of order n
 //  Q = H_0 * H_1 * ... * H_{k-1}
 // as returned by Dgerqf.
 //
-// Each entry of tau contains the scalar factor of the elementary reflector H_i
+// On entry, the (m-k+i)-th row of A must contain the vector which defines the
-// and the length of tau must be at least k. The length of work must be at least m.
+// elementary reflector H_i, for i = 0,1,...,k, as returned by Dgerqf. On
-// a is a matrix of dimensions (n,lda). On entry the [m-k+i-1]-th row must contain (counting from zero)
+// return, A will contain the m×n matrix Q.
-// the vector which defines the elementary reflector H_i, for i = 0,1,2,...,k-1, as
-// returned by Dgerqf in the last k rows of its array argument A.
-// On exit, the m×n matrix Q.
-//  n >= m >= k >= 0
 //
-// Dorgr2 will panic if the conditions on input values are not met.
+// The i-th element of tau must contain the scalar factor of the elementary
+// reflector H_i, as returned by Dgerqf.
+//
+// It must hold that
+//  n >= m >= k >= 0,
+// the length of tau must be k and the length of work must be m, otherwise
+// Dorgr2 will panic.
 //
 // Dorgr2 is an internal routine. It is exported for testing purposes.
 func (impl Implementation) Dorgr2(m, n, k int, a []float64, lda int, tau, work []float64) {
 	switch {
-	case m < 0:
-		panic(mLT0)
-	case n < 0:
-		panic(nLT0)
-	case m > n:
-		panic(mGTN)
 	case k < 0:
 		panic(kLT0)
-	case k > m:
+	case m < k:
 		panic(kGTM)
+	case n < m:
+		panic(mGTN)
 	case lda < max(1, n):
 		panic(badLdA)
-	case len(work) < m:
-		panic(shortWork)
 	}
 
 	// Quick return if possible.
@@ -54,17 +49,18 @@ func (impl Implementation) Dorgr2(m, n, k int, a []float64, lda int, tau, work [
 		panic(badLenTau)
 	case len(a) < (m-1)*lda+n:
 		panic(shortA)
+	case len(work) < m:
+		panic(shortWork)
 	}
 
-	if k < m {
 	// Initialise rows 0:m-k to rows of the unit matrix.
 	for l := 0; l < m-k; l++ {
-			for j := 0; j < n; j++ {
+		row := a[l*lda : l*lda+n]
-				a[l*lda+j] = 0
+		for j := range row {
+			row[j] = 0
 		}
 		a[l*lda+n-m+l] = 1
 	}
-	}
 	bi := blas64.Implementation()
 	for i := 0; i < k; i++ {
 		ii := m - k + i

lapack/testlapack/dorgr2.go

@@ -6,12 +6,14 @@ package testlapack
 
 import (
 	"fmt"
+	"math"
 	"testing"
 
 	"golang.org/x/exp/rand"
 
 	"gonum.org/v1/gonum/blas"
 	"gonum.org/v1/gonum/blas/blas64"
+	"gonum.org/v1/gonum/floats"
 	"gonum.org/v1/gonum/lapack"
 )
 
@@ -26,7 +28,7 @@ func Dorgr2Test(t *testing.T, impl Dorgr2er) {
 	for _, k := range []int{0, 1, 2, 5} {
 		for _, m := range []int{k, k + 1, k + 2, k + 4} {
 			for _, n := range []int{m, m + 1, m + 2, m + 4, m + 7} {
-				for _, lda := range []int{n, n + 5} {
+				for _, lda := range []int{max(1, n), n + 5} {
 					dorgr2Test(t, impl, rnd, m, n, k, lda)
 				}
 			}
@@ -35,60 +37,72 @@ func Dorgr2Test(t *testing.T, impl Dorgr2er) {
 }
 
 func dorgr2Test(t *testing.T, impl Dorgr2er, rnd *rand.Rand, m, n, k, lda int) {
-	const tol = 1e-12
+	const tol = 1e-14
+
 	name := fmt.Sprintf("m=%v,n=%v,k=%v,lda=%v", m, n, k, lda)
 
-	if lda == 0 {
+	// Generate a random m×n matrix A.
-		lda = n
-	}
-
 	a := randomGeneral(m, n, lda, rnd)
-	aCopy := cloneGeneral(a)
+
 	// Compute the RQ decomposition of A.
+	rq := cloneGeneral(a)
 	tau := make([]float64, m)
 	work := make([]float64, 1)
-	impl.Dgerqf(m, n, a.Data, a.Stride, tau, work, -1)
+	impl.Dgerqf(m, n, rq.Data, rq.Stride, tau, work, -1)
 	work = make([]float64, int(work[0]))
-	impl.Dgerqf(m, n, a.Data, a.Stride, tau, work, len(work))
+	impl.Dgerqf(m, n, rq.Data, rq.Stride, tau, work, len(work))
+
+	tauCopy := make([]float64, len(tau))
+	copy(tauCopy, tau)
 
-	// Generate the upper triangular matrix R from subarray A[m-k:m, n-m:n]
-	r := zeros(k, m, m)
-	for i := 0; i < k; i++ {
-		for j := 0; j < k; j++ {
-			ia := i + a.Rows - k
-			ja := j + a.Cols - k
-			jr := j + r.Cols - k
-			if i <= j {
-				r.Data[i*r.Stride+jr] = a.Data[ia*a.Stride+ja]
-			}
-		}
-	}
 	// Compute the matrix Q using Dorg2r.
-	impl.Dorgr2(m, n, k, a.Data, a.Stride, tau[m-k:m], work)
+	q := cloneGeneral(rq)
+	impl.Dorgr2(m, n, k, q.Data, q.Stride, tau[m-k:m], work)
+
 	if m == 0 {
 		return
 	}
-	q := a
+
-	// Test Q orthogonality.
+	// Check that tau hasn't been modified.
-	res := residualOrthogonal(q, true)
+	if !floats.Equal(tau, tauCopy) {
-	if res > tol {
+		t.Errorf("%v: unexpected modification in tau", name)
-		t.Errorf("%v: |I - Q * Qᵀ| residual too large (%g)", name, res)
 	}
+
+	// Check that Q has orthonormal rows.
+	res := residualOrthogonal(q, true)
+	if res > tol || math.IsNaN(res) {
+		t.Errorf("%v: residual |I - Q*Qᵀ| too large, got %v, want <= %v", name, res, tol)
+	}
+
 	if k == 0 {
 		return
 	}
 
-	// Reconstruct last rows of A.
+	// Extract the k×m upper triangular matrix R from RQ[m-k:m,n-k:n].
+	r := zeros(k, m, m)
+	for i := 0; i < k; i++ {
+		for j := 0; j < k; j++ {
+			ii := rq.Rows - k + i
+			jj := rq.Cols - k + j
+			jr := r.Cols - k + j
+			if i <= j {
+				r.Data[i*r.Stride+jr] = rq.Data[ii*rq.Stride+jj]
+			}
+		}
+	}
+
+	// Construct a view A[m-k:m,0:n] of the last k rows of A.
 	aRec := blas64.General{
 		Rows:   k,
-		Cols:   aCopy.Cols,
+		Cols:   n,
-		Data:   aCopy.Data[(a.Rows-k)*lda:],
+		Data:   a.Data[(m-k)*a.Stride:],
-		Stride: lda,
+		Stride: a.Stride,
 	}
-	// Test |A[m-k:m,0:n] - R[m-k:m,0:m] * Q| is small
+	// Compute A - R*Q.
-	blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, r, q, -1, aRec)
+	blas64.Gemm(blas.NoTrans, blas.NoTrans, -1, r, q, 1, aRec)
+	// Check that |A - R*Q| is small.
 	res = dlange(lapack.MaxColumnSum, aRec.Rows, aRec.Cols, aRec.Data, aRec.Stride)
-	if res > tol {
+	if res > tol || math.IsNaN(res) {
-		t.Errorf("%v: |A[m-k:m,0:n] - R[m-k:m,0:m] * Q| residual too large (%g)", name, res)
+		t.Errorf("%v: residual |A - R*Q| too large, got %v, want <= %v", name, res, tol)
 	}
 }