Repository: gonum/gonum
Pull request #1688 (2021-06-22 11:30:49)

lapack/gonum: add Dgesv
More for completeness than utility because Dgesv doesn't have a trans parameter to solve `A^T*X = B`.

Please take a look.

<!--
Checklist:

- API changes have been discussed
- code is goformated correctly (goimports)
- packages with generated code have had code generation run
- tests pass locally
- linked to relevant issues

Please make sure your commit message summary line and pull request
title match the Go convention; a one-line summary of the change,
prefixed by the primary affected package that should complete the
sentence, "This change modifies Gonum to _____."
-->


codecov-commenter: # [Codecov](https://codecov.io/gh/gonum/gonum/pull/1688?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum) Report
> Merging [#1688](https://codecov.io/gh/gonum/gonum/pull/1688?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum) (c4554c8) into [master](https://codecov.io/gh/gonum/gonum/commit/2600994ee4df66aa98cc884dd6eccd44bb60cba8?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum) (2600994) will **decrease** coverage by `0.02%`.
> The diff coverage is `36.36%`.

> :exclamation: Current head c4554c8 differs from pull request most recent head 3416c02. Consider uploading reports for the commit 3416c02 to get more accurate results
[![Impacted file tree graph](https://codecov.io/gh/gonum/gonum/pull/1688/graphs/tree.svg?width=650&height=150&src=pr&token=B7qHnqSPNL&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum)](https://codecov.io/gh/gonum/gonum/pull/1688?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum)

```diff
@@            Coverage Diff             @@
##           master    #1688      +/-   ##
==========================================
- Coverage   73.28%   73.25%   -0.03%     
==========================================
  Files         507      505       -2     
  Lines       60019    59982      -37     
==========================================
- Hits        43985    43940      -45     
- Misses      13553    13564      +11     
+ Partials     2481     2478       -3     
```


| [Impacted Files](https://codecov.io/gh/gonum/gonum/pull/1688?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum) | Coverage Δ | |
|---|---|---|
| [lapack/gonum/dgesv.go](https://codecov.io/gh/gonum/gonum/pull/1688/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-bGFwYWNrL2dvbnVtL2RnZXN2Lmdv) | `36.36% <36.36%> (ø)` | |
| [graph/internal/set/same\_appengine.go](https://codecov.io/gh/gonum/gonum/pull/1688/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-Z3JhcGgvaW50ZXJuYWwvc2V0L3NhbWVfYXBwZW5naW5lLmdv) | | |
| [mat/offset\_appengine.go](https://codecov.io/gh/gonum/gonum/pull/1688/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-bWF0L29mZnNldF9hcHBlbmdpbmUuZ28=) | | |
| [graph/iterator/nodes\_map\_safe.go](https://codecov.io/gh/gonum/gonum/pull/1688/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-Z3JhcGgvaXRlcmF0b3Ivbm9kZXNfbWFwX3NhZmUuZ28=) | | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/gonum/gonum/pull/1688?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/gonum/gonum/pull/1688?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum). Last update [2600994...3416c02](https://codecov.io/gh/gonum/gonum/pull/1688?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum).


Patch:
diff --git a/lapack/gonum/dgesv.go b/lapack/gonum/dgesv.go
new file mode 100644
index 00000000..53ba884d
--- /dev/null
+++ b/lapack/gonum/dgesv.go
@@ -0,0 +1,56 @@
+// Copyright ©2021 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package gonum
+
+import "gonum.org/v1/gonum/blas"
+
+// Dgesv computes the solution to a real system of linear equations
+//  A * X = B
+// where A is an n×n matrix and X and B are n×nrhs matrices.
+//
+// The LU decomposition with partial pivoting and row interchanges is used to
+// factor A as
+//  A = P * L * U
+// where P is a permutation matrix, L is unit lower triangular, and U is upper
+// triangular. On return, the factors L and U are stored in a; the unit diagonal
+// elements of L are not stored. The row pivot indices that define the
+// permutation matrix P are stored in ipiv.
+//
+// The factored form of A is then used to solve the system of equations A * X =
+// B. On entry, b contains the right hand side matrix B. On return, if ok is
+// true, b contains the solution matrix X.
+func (impl Implementation) Dgesv(n, nrhs int, a []float64, lda int, ipiv []int, b []float64, ldb int) (ok bool) {
+	switch {
+	case n < 0:
+		panic(nLT0)
+	case nrhs < 0:
+		panic(nrhsLT0)
+	case lda < max(1, n):
+		panic(badLdA)
+	case ldb < max(1, nrhs):
+		panic(badLdB)
+	}
+
+	// Quick return if possible.
+	if n == 0 || nrhs == 0 {
+		return true
+	}
+
+	switch {
+	case len(a) < (n-1)*lda+n:
+		panic(shortAB)
+	case len(ipiv) != n:
+		panic(badLenIpiv)
+	case len(b) < (n-1)*ldb+nrhs:
+		panic(shortB)
+	}
+
+	ok = impl.Dgetrf(n, n, a, lda, ipiv)
+	if ok {
+		impl.Dgetrs(blas.NoTrans, n, nrhs, a, lda, ipiv, b, ldb)
+	}
+
+	return ok
+}


Test patch:
diff --git a/lapack/gonum/lapack_test.go b/lapack/gonum/lapack_test.go
index a296d64f..97004ed4 100644
--- a/lapack/gonum/lapack_test.go
+++ b/lapack/gonum/lapack_test.go
@@ -117,6 +117,11 @@ func TestDgerqf(t *testing.T) {
 	testlapack.DgerqfTest(t, impl)
 }
 
+func TestDgesv(t *testing.T) {
+	t.Parallel()
+	testlapack.DgesvTest(t, impl)
+}
+
 func TestDgesvd(t *testing.T) {
 	t.Parallel()
 	const tol = 1e-13
diff --git a/lapack/testlapack/dgesv.go b/lapack/testlapack/dgesv.go
new file mode 100644
index 00000000..c41a37d1
--- /dev/null
+++ b/lapack/testlapack/dgesv.go
@@ -0,0 +1,164 @@
+// Copyright ©2021 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+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/lapack"
+)
+
+type Dgesver interface {
+	Dgesv(n, nrhs int, a []float64, lda int, ipiv []int, b []float64, ldb int) bool
+
+	Dgetri(n int, a []float64, lda int, ipiv []int, work []float64, lwork int) bool
+}
+
+func DgesvTest(t *testing.T, impl Dgesver) {
+	rnd := rand.New(rand.NewSource(1))
+	for _, n := range []int{0, 1, 2, 3, 4, 5, 10, 50, 100} {
+		for _, nrhs := range []int{0, 1, 2, 5} {
+			for _, lda := range []int{max(1, n), n + 5} {
+				for _, ldb := range []int{max(1, nrhs), nrhs + 5} {
+					dgesvTest(t, impl, rnd, n, nrhs, lda, ldb)
+				}
+			}
+		}
+	}
+}
+
+func dgesvTest(t *testing.T, impl Dgesver, rnd *rand.Rand, n, nrhs, lda, ldb int) {
+	const tol = 1e-15
+
+	name := fmt.Sprintf("n=%v,nrhs=%v,lda=%v,ldb=%v", n, nrhs, lda, ldb)
+
+	// Create a random system matrix A and the solution X.
+	a := randomGeneral(n, n, lda, rnd)
+	xWant := randomGeneral(n, nrhs, ldb, rnd)
+
+	// Compute the right hand side matrix B = A*X.
+	b := zeros(n, nrhs, ldb)
+	blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, a, xWant, 0, b)
+
+	// Allocate a slice for row pivots and fill it with invalid indices.
+	ipiv := make([]int, n)
+	for i := range ipiv {
+		ipiv[i] = -1
+	}
+
+	// Call Dgesv to solve A*X = B.
+	lu := cloneGeneral(a)
+	xGot := cloneGeneral(b)
+	ok := impl.Dgesv(n, nrhs, lu.Data, lu.Stride, ipiv, xGot.Data, xGot.Stride)
+
+	if !ok {
+		t.Errorf("%v: unexpected failure in Dgesv", name)
+		return
+	}
+
+	if n == 0 || nrhs == 0 {
+		return
+	}
+
+	// Check that all elements of ipiv have been set.
+	ipivSet := true
+	for _, ipv := range ipiv {
+		if ipv == -1 {
+			ipivSet = false
+			break
+		}
+	}
+	if !ipivSet {
+		t.Fatalf("%v: not all elements of ipiv set", name)
+		return
+	}
+
+	// Compute the reciprocal of the condition number of A from its LU
+	// decomposition before it's overwritten further below.
+	aInv := cloneGeneral(lu)
+	impl.Dgetri(n, aInv.Data, aInv.Stride, ipiv, make([]float64, n), n)
+	ainvnorm := dlange(lapack.MaxColumnSum, n, n, aInv.Data, aInv.Stride)
+	anorm := dlange(lapack.MaxColumnSum, n, n, a.Data, a.Stride)
+	rcond := 1 / anorm / ainvnorm
+
+	// Reconstruct matrix A from factors and compute residual.
+	//
+	// Extract L and U from lu.
+	l := zeros(n, n, n)
+	u := zeros(n, n, n)
+	for i := 0; i < n; i++ {
+		for j := 0; j < i; j++ {
+			l.Data[i*l.Stride+j] = lu.Data[i*lu.Stride+j]
+		}
+		l.Data[i*l.Stride+i] = 1
+		for j := i; j < n; j++ {
+			u.Data[i*u.Stride+j] = lu.Data[i*lu.Stride+j]
+		}
+	}
+	// Compute L*U.
+	blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, l, u, 0, lu)
+	// Apply P to L*U.
+	for i := n - 1; i >= 0; i-- {
+		ip := ipiv[i]
+		if ip == i {
+			continue
+		}
+		row1 := blas64.Vector{N: n, Data: lu.Data[i*lu.Stride:], Inc: 1}
+		row2 := blas64.Vector{N: n, Data: lu.Data[ip*lu.Stride:], Inc: 1}
+		blas64.Swap(row1, row2)
+	}
+	// Compute P*L*U - A.
+	for i := 0; i < n; i++ {
+		for j := 0; j < n; j++ {
+			lu.Data[i*lu.Stride+j] -= a.Data[i*a.Stride+j]
+		}
+	}
+	// Compute the residual |P*L*U - A|.
+	resid := dlange(lapack.MaxColumnSum, n, n, lu.Data, lu.Stride)
+	resid /= float64(n) * anorm
+	if resid > tol || math.IsNaN(resid) {
+		t.Errorf("%v: residual |P*L*U - A| is too large, got %v, want <= %v", name, resid, tol)
+	}
+
+	// Compute residual of the computed solution.
+	//
+	// Compute B - A*X.
+	blas64.Gemm(blas.NoTrans, blas.NoTrans, -1, a, xGot, 1, b)
+	// Compute the maximum over the number of right hand sides of |B - A*X| / (|A| * |X|).
+	resid = 0
+	for j := 0; j < nrhs; j++ {
+		bnorm := blas64.Asum(blas64.Vector{N: n, Data: b.Data[j:], Inc: b.Stride})
+		xnorm := blas64.Asum(blas64.Vector{N: n, Data: xGot.Data[j:], Inc: xGot.Stride})
+		resid = math.Max(resid, bnorm/anorm/xnorm)
+	}
+	if resid > tol || math.IsNaN(resid) {
+		t.Errorf("%v: residual |B - A*X| is too large, got %v, want <= %v", name, resid, tol)
+	}
+
+	// Compare the computed solution with the generated exact solution.
+	//
+	// Compute X - XWANT.
+	for i := 0; i < n; i++ {
+		for j := 0; j < nrhs; j++ {
+			xGot.Data[i*xGot.Stride+j] -= xWant.Data[i*xWant.Stride+j]
+		}
+	}
+	// Compute the maximum of |X - XWANT|/|XWANT| over all the vectors X and XWANT.
+	resid = 0
+	for j := 0; j < nrhs; j++ {
+		xnorm := dlange(lapack.MaxAbs, n, 1, xWant.Data[j:], xWant.Stride)
+		diff := dlange(lapack.MaxAbs, n, 1, xGot.Data[j:], xGot.Stride)
+		resid = math.Max(resid, diff/xnorm*rcond)
+	}
+	if resid > tol || math.IsNaN(resid) {
+		t.Errorf("%v: residual |X-XWANT| is too large, got %v, want <= %v", name, resid, tol)
+	}
+}