lapack/lapack64: add Geqp3 and clean up docs

lapack/gonum/dgeqp3.go

@@ -9,35 +9,40 @@ import (
 	"gonum.org/v1/gonum/blas/blas64"
 )
 
-// Dgeqp3 computes a QR factorization with column pivoting of the
+// Dgeqp3 computes a QR factorization with column pivoting of the m×n matrix A:
-// m×n matrix A: A*P = Q*R using Level 3 BLAS.
 //
-// The matrix Q is represented as a product of elementary reflectors
+//	A*P = Q*R
 //
-//	Q = H_0 H_1 . . . H_{k-1}, where k = min(m,n).
+// where P is a permutation matrix, Q is an orthogonal matrix and R is a
+// min(m,n)×n upper trapezoidal matrix.
+//
+// On return, the upper triangle of A contains the matrix R. The elements below
+// the diagonal together with tau represent the matrix Q as a product of
+// elementary reflectors
+//
+//	Q = H_0 * H_1 * ... * H_{k-1}, where k = min(m,n).
 //
 // Each H_i has the form
 //
 //	H_i = I - tau * v * vᵀ
 //
-// where tau and v are real vectors with v[0:i-1] = 0 and v[i] = 1;
+// where tau is a scalar and v is a vector with v[0:i] = 0 and v[i] = 1;
-// v[i:m] is stored on exit in A[i:m, i], and tau in tau[i].
+// v[i+1:m] is stored on exit in A[i+1:m,i], and tau in tau[i].
 //
-// jpvt specifies a column pivot to be applied to A. If
+// jpvt specifies a column pivot to be applied to A. On entry, if jpvt[j] is at
-// jpvt[j] is at least zero, the jth column of A is permuted
+// least zero, the jth column of A is permuted to the front of A*P (a leading
-// to the front of A*P (a leading column), if jpvt[j] is -1
+// column), if jpvt[j] is -1 the jth column of A is a free column. If jpvt[j] <
-// the jth column of A is a free column. If jpvt[j] < -1, Dgeqp3
+// -1, Dgeqp3 will panic. On return, jpvt holds the permutation that was
-// will panic. On return, jpvt holds the permutation that was
+// applied; the jth column of A*P was the jpvt[j] column of A. jpvt must have
-// applied; the jth column of A*P was the jpvt[j] column of A.
+// length n or Dgeqp3 will panic.
-// jpvt must have length n or Dgeqp3 will panic.
 //
-// tau holds the scalar factors of the elementary reflectors.
+// tau holds the scalar factors of the elementary reflectors. It must have
-// It must have length min(m, n), otherwise Dgeqp3 will panic.
+// length min(m,n), otherwise Dgeqp3 will panic.
 //
 // work must have length at least max(1,lwork), and lwork must be at least
-// 3*n+1, otherwise Dgeqp3 will panic. For optimal performance lwork must
+// 3*n+1, otherwise Dgeqp3 will panic. For optimal performance lwork must be at
-// be at least 2*n+(n+1)*nb, where nb is the optimal blocksize. On return,
+// least 2*n+(n+1)*nb, where nb is the optimal blocksize. On return, work[0]
-// work[0] will contain the optimal value of lwork.
+// will contain the optimal value of lwork.
 //
 // If lwork == -1, instead of performing Dgeqp3, only the optimal value of lwork
 // will be stored in work[0].

lapack/lapack.go

@@ -15,6 +15,7 @@ type Float64 interface {
 	Dgeev(jobvl LeftEVJob, jobvr RightEVJob, n int, a []float64, lda int, wr, wi []float64, vl []float64, ldvl int, vr []float64, ldvr int, work []float64, lwork int) (first int)
 	Dgels(trans blas.Transpose, m, n, nrhs int, a []float64, lda int, b []float64, ldb int, work []float64, lwork int) bool
 	Dgelqf(m, n int, a []float64, lda int, tau, work []float64, lwork int)
+	Dgeqp3(m, n int, a []float64, lda int, jpvt []int, tau, work []float64, lwork int)
 	Dgeqrf(m, n int, a []float64, lda int, tau, work []float64, lwork int)
 	Dgesvd(jobU, jobVT SVDJob, m, n int, a []float64, lda int, s, u []float64, ldu int, vt []float64, ldvt int, work []float64, lwork int) (ok bool)
 	Dgetrf(m, n int, a []float64, lda int, ipiv []int) (ok bool)

lapack/lapack64/lapack64.go

@@ -222,6 +222,47 @@ func Gels(trans blas.Transpose, a blas64.General, b blas64.General, work []float
 	return lapack64.Dgels(trans, a.Rows, a.Cols, b.Cols, a.Data, max(1, a.Stride), b.Data, max(1, b.Stride), work, lwork)
 }
 
+// Geqp3 computes a QR factorization with column pivoting of the m×n matrix A:
+//
+//	A*P = Q*R
+//
+// where P is a permutation matrix, Q is an orthogonal matrix and R is a
+// min(m,n)×n upper trapezoidal matrix.
+//
+// On return, the upper triangle of A contains the matrix R. The elements below
+// the diagonal together with tau represent the matrix Q as a product of
+// elementary reflectors
+//
+//	Q = H_0 * H_1 * ... * H_{k-1}, where k = min(m,n).
+//
+// Each H_i has the form
+//
+//	H_i = I - tau * v * vᵀ
+//
+// where tau is a scalar and v is a vector with v[0:i] = 0 and v[i] = 1;
+// v[i+1:m] is stored on exit in A[i+1:m,i], and tau in tau[i].
+//
+// jpvt specifies a column pivot to be applied to A. On entry, if jpvt[j] is at
+// least zero, the jth column of A is permuted to the front of A*P (a leading
+// column), if jpvt[j] is -1 the jth column of A is a free column. If jpvt[j] <
+// -1, Geqp3 will panic. On return, jpvt holds the permutation that was applied;
+// the jth column of A*P was the jpvt[j] column of A. jpvt must have length n or
+// Geqp3 will panic.
+//
+// tau holds the scalar factors of the elementary reflectors. It must have
+// length min(m,n), otherwise Geqp3 will panic.
+//
+// work must have length at least max(1,lwork), and lwork must be at least
+// 3*n+1, otherwise Geqp3 will panic. For optimal performance lwork must be at
+// least 2*n+(n+1)*nb, where nb is the optimal blocksize. On return, work[0]
+// will contain the optimal value of lwork.
+//
+// If lwork == -1, instead of performing Geqp3, only the optimal value of lwork
+// will be stored in work[0].
+func Geqp3(a blas64.General, jpvt []int, tau, work []float64, lwork int) {
+	lapack64.Dgeqp3(a.Rows, a.Cols, a.Data, max(1, a.Stride), jpvt, tau, work, lwork)
+}
+
 // Geqrf computes the QR factorization of the m×n matrix A using a blocked
 // algorithm. A is modified to contain the information to construct Q and R.
 // The upper triangle of a contains the matrix R. The lower triangular elements