lapack/gonum: fix inf-norm computation for unit lower triangular matrices
MODIFY
lapack/gonum/dlantr.go
lapack/gonum/dlantr.go
@@ -171,7 +171,7 @@ func (impl Implementation) Dlantr(norm lapack.MatrixNorm, uplo blas.Uplo, diag b
}
return maxsum
} else {
- for i := 1; i < m; i++ {
+ for i := 0; i < m; i++ {
var sum float64
if i < minmn {
sum = 1

MODIFY
lapack/testlapack/dlantr.go
lapack/testlapack/dlantr.go
package testlapack
import (
+ "fmt"
"math"
"testing"
"golang.org/x/exp/rand"
"gonum.org/v1/gonum/blas"
+ "gonum.org/v1/gonum/floats"
"gonum.org/v1/gonum/lapack"
)
@@ -21,70 +23,115 @@ type Dlantrer interface {
func DlantrTest(t *testing.T, impl Dlantrer) {
rnd := rand.New(rand.NewSource(1))
- for _, norm := range []lapack.MatrixNorm{lapack.MaxAbs, lapack.MaxColumnSum, lapack.MaxRowSum, lapack.Frobenius} {
- for _, diag := range []blas.Diag{blas.NonUnit, blas.Unit} {
+ for _, m := range []int{0, 1, 2, 3, 4, 5, 10} {
+ for _, n := range []int{0, 1, 2, 3, 4, 5, 10} {
for _, uplo := range []blas.Uplo{blas.Lower, blas.Upper} {
- for _, test := range []struct {
- m, n, lda int
- }{
- {3, 3, 0},
- {3, 5, 0},
- {10, 5, 0},
-
- {5, 5, 11},
- {5, 10, 11},
- {10, 5, 11},
- } {
- // Do a couple of random trials since the values change.
- for trial := 0; trial < 100; trial++ {
- m := test.m
- n := test.n
- lda := test.lda
- if lda == 0 {
- lda = n
- }
- a := make([]float64, m*lda)
- if trial == 0 {
- for i := range a {
- a[i] = float64(i)
+ for _, diag := range []blas.Diag{blas.NonUnit, blas.Unit} {
+ for _, lda := range []int{max(1, n), n + 3} {
+ dlantrTest(t, impl, rnd, uplo, diag, m, n, lda)
}
- } else {
+ }
+ }
+ }
+ }
+}
+
+func dlantrTest(t *testing.T, impl Dlantrer, rnd *rand.Rand, uplo blas.Uplo, diag blas.Diag, m, n, lda int) {
+ const tol = 1e-14
+
+ // Generate a random triangular matrix. If the matrix has unit diagonal,
+ // don't set the diagonal elements to 1.
+ a := make([]float64, max(0, (m-1)*lda+n))
for i := range a {
a[i] = rnd.NormFloat64()
}
+ rowsum := make([]float64, m)
+ colsum := make([]float64, n)
+ var frobWant, maxabsWant float64
+ if diag == blas.Unit {
+ // Account for the unit diagonal.
+ for i := 0; i < min(m, n); i++ {
+ rowsum[i] = 1
+ colsum[i] = 1
}
- aDense := make([]float64, len(a))
- if uplo == blas.Lower {
- for i := 0; i < m; i++ {
- for j := 0; j <= min(i, n-1); j++ {
- aDense[i*lda+j] = a[i*lda+j]
+ frobWant = float64(min(m, n))
+ if min(m, n) > 0 {
+ maxabsWant = 1
}
}
- } else {
- for i := 0; i < m; i++ {
- for j := i; j < n; j++ {
- aDense[i*lda+j] = a[i*lda+j]
+ if uplo == blas.Upper {
+ for i := 0; i < min(m, n); i++ {
+ start := i
+ if diag == blas.Unit {
+ start = i + 1
}
+ for j := start; j < n; j++ {
+ aij := 2*rnd.Float64() - 1
+ a[i*lda+j] = aij
+ rowsum[i] += math.Abs(aij)
+ colsum[j] += math.Abs(aij)
+ maxabsWant = math.Max(maxabsWant, math.Abs(aij))
+ frobWant += aij * aij
}
}
+ } else {
+ for i := 0; i < m; i++ {
+ end := i
if diag == blas.Unit {
- for i := 0; i < min(m, n); i++ {
- aDense[i*lda+i] = 1
+ end = i - 1
+ }
+ for j := 0; j <= min(end, n-1); j++ {
+ aij := 2*rnd.Float64() - 1
+ a[i*lda+j] = aij
+ rowsum[i] += math.Abs(aij)
+ colsum[j] += math.Abs(aij)
+ maxabsWant = math.Max(maxabsWant, math.Abs(aij))
+ frobWant += aij * aij
}
}
- work := make([]float64, n+6)
- for i := range work {
- work[i] = rnd.Float64()
}
- got := impl.Dlantr(norm, uplo, diag, m, n, a, lda, work)
- want := impl.Dlange(norm, m, n, aDense, lda, work)
- if math.Abs(got-want) > 1e-13 {
- t.Errorf("Norm mismatch. norm = %c, unitdiag = %v, upper = %v, m = %v, n = %v, lda = %v, Want %v, got %v.",
- norm, diag == blas.Unit, uplo == blas.Upper, m, n, lda, got, want)
+ frobWant = math.Sqrt(frobWant)
+ var maxcolsumWant, maxrowsumWant float64
+ if n > 0 {
+ maxcolsumWant = floats.Max(colsum)
}
+ if m > 0 {
+ maxrowsumWant = floats.Max(rowsum)
}
+
+ aCopy := make([]float64, len(a))
+ copy(aCopy, a)
+
+ work := make([]float64, n)
+ for _, norm := range []lapack.MatrixNorm{lapack.MaxAbs, lapack.MaxColumnSum, lapack.MaxRowSum, lapack.Frobenius} {
+ name := fmt.Sprintf("norm=%v,uplo=%v,diag=%v,m=%v,n=%v,lda=%v", string(norm), string(uplo), string(diag), m, n, lda)
+
+ normGot := impl.Dlantr(norm, uplo, diag, m, n, a, lda, work)
+
+ if !floats.Equal(a, aCopy) {
+ t.Fatalf("%v: unexpected modification of a", name)
+ }
+
+ if norm == lapack.MaxAbs {
+ // MaxAbs norm involves no floating-point computation,
+ // so we expect exact equality here.
+ if normGot != maxabsWant {
+ t.Errorf("%v: unexpected result; got %v, want %v", name, normGot, maxabsWant)
}
+ continue
}
+
+ var normWant float64
+ switch norm {
+ case lapack.MaxColumnSum:
+ normWant = maxcolsumWant
+ case lapack.MaxRowSum:
+ normWant = maxrowsumWant
+ case lapack.Frobenius:
+ normWant = frobWant
+ }
+ if math.Abs(normGot-normWant) >= tol {
+ t.Errorf("%v: unexpected result; got %v, want %v", name, normGot, normWant)
}
}
}

lapack/gonum: fix inf-norm computation for unit lower triangular matrices