Add the linear solve routines (Dgetrs, Dgels) to the lapack64 interface

cgo/lapack.go

@@ -36,6 +36,13 @@ func min(m, n int) int {
 	return n
 }
 
+func max(m, n int) int {
+	if m < n {
+		return n
+	}
+	return m
+}
+
 // checkMatrix verifies the parameters of a matrix input.
 // Copied from lapack/native. Keep in sync.
 func checkMatrix(m, n int, a []float64, lda int) {
@@ -193,6 +200,50 @@ func (impl Implementation) Dgeqrf(m, n int, a []float64, lda int, tau, work []fl
 	clapack.Dgeqrf(m, n, a, lda, tau)
 }
 
+// Dgels finds a minimum-norm solution based on the matrices A and B using the
+// QR or LQ factorization. Dgels returns false if the matrix
+// A is singular, and true if this solution was successfully found.
+//
+// The minimization problem solved depends on the input parameters.
+//
+//  1. If m >= n and trans == blas.NoTrans, Dgels finds X such that || A*X - B||_2
+//     is minimized.
+//  2. If m < n and trans == blas.NoTrans, Dgels finds the minimum norm solution of
+//     A * X = B.
+//  3. If m >= n and trans == blas.Trans, Dgels finds the minimum norm solution of
+//     A^T * X = B.
+//  4. If m < n and trans == blas.Trans, Dgels finds X such that || A*X - B||_2
+//     is minimized.
+// Note that the least-squares solutions (cases 1 and 3) perform the minimization
+// per column of B. This is not the same as finding the minimum-norm matrix.
+//
+// The matrix A is a general matrix of size m×n and is modified during this call.
+// The input matrix B is of size max(m,n)×nrhs, and serves two purposes. On entry,
+// the elements of b specify the input matrix B. B has size m×nrhs if
+// trans == blas.NoTrans, and n×nrhs if trans == blas.Trans. On exit, the
+// leading submatrix of b contains the solution vectors X. If trans == blas.NoTrans,
+// this submatrix is of size n×nrhs, and of size m×nrhs otherwise.
+//
+// The C interface does not support providing temporary storage. To provide compatibility
+// with native, lwork == -1 will not run Dgeqrf but will instead write the minimum
+// work necessary to work[0]. If len(work) < lwork, Dgeqrf will panic.
+func (impl Implementation) Dgels(trans blas.Transpose, m, n, nrhs int, a []float64, lda int, b []float64, ldb int, work []float64, lwork int) bool {
+	mn := min(m, n)
+	if lwork == -1 {
+		work[0] = float64(mn + max(mn, nrhs))
+		return true
+	}
+	checkMatrix(m, n, a, lda)
+	checkMatrix(mn, nrhs, b, ldb)
+	if len(work) < lwork {
+		panic(shortWork)
+	}
+	if lwork < mn+max(mn, nrhs) {
+		panic(badWork)
+	}
+	return clapack.Dgels(trans, m, n, nrhs, a, lda, b, ldb)
+}
+
 // Dgetf2 computes the LU decomposition of the m×n matrix A.
 // The LU decomposition is a factorization of a into
 //  A = P * L * U
@@ -223,7 +274,7 @@ func (Implementation) Dgetf2(m, n int, a []float64, lda int, ipiv []int) (ok boo
 }
 
 // Dgetrf computes the LU decomposition of the m×n matrix A.
-// The LU decomposition is a factorization of a into
+// The LU decomposition is a factorization of A into
 //  A = P * L * U
 // where P is a permutation matrix, L is a unit lower triangular matrix, and
 // U is a (usually) non-unit upper triangular matrix. On exit, L and U are stored