From ff3e32092eee436eec078fb361dd6ba3b08bc6e9 Mon Sep 17 00:00:00 2001 From: Vladimir Chalupecky Date: Thu, 21 Sep 2023 06:20:21 +0200 Subject: [PATCH] lapack/lapack64: add Geqp3 and clean up docs --- lapack/gonum/dgeqp3.go | 41 +++++++++++++++++++++---------------- lapack/lapack.go | 1 + lapack/lapack64/lapack64.go | 41 +++++++++++++++++++++++++++++++++++++ 3 files changed, 65 insertions(+), 18 deletions(-) diff --git a/lapack/gonum/dgeqp3.go b/lapack/gonum/dgeqp3.go index f96f03be..9b720a78 100644 --- a/lapack/gonum/dgeqp3.go +++ b/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 -// m×n matrix A: A*P = Q*R using Level 3 BLAS. +// Dgeqp3 computes a QR factorization with column pivoting of the m×n matrix A: // -// 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; -// v[i:m] is stored on exit in A[i:m, i], and tau in tau[i]. +// 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. 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, Dgeqp3 -// 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 Dgeqp3 will panic. +// 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, Dgeqp3 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 Dgeqp3 will panic. // -// tau holds the scalar factors of the elementary reflectors. -// It must have length min(m, n), otherwise Dgeqp3 will panic. +// tau holds the scalar factors of the elementary reflectors. It must have +// 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 -// 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. +// 3*n+1, otherwise Dgeqp3 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 Dgeqp3, only the optimal value of lwork // will be stored in work[0]. diff --git a/lapack/lapack.go b/lapack/lapack.go index 18790462..a08324d9 100644 --- a/lapack/lapack.go +++ b/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) diff --git a/lapack/lapack64/lapack64.go b/lapack/lapack64/lapack64.go index 0cf94bfc..19a40cfc 100644 --- a/lapack/lapack64/lapack64.go +++ b/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