Add dlansy and test

cgo/lapack.go

@@ -88,6 +88,25 @@ func (impl Implementation) Dlange(norm lapack.MatrixNorm, m, n int, a []float64,
 	return clapack.Dlange(byte(norm), m, n, a, lda)
 }
 
+// Dlansy computes the specified norm of an n×n symmetric matrix. If
+// norm == lapack.MaxColumnSum or norm == lapackMaxRowSum work must have length
+// at least n, otherwise work is unused.
+func (impl Implementation) Dlansy(norm lapack.MatrixNorm, uplo blas.Uplo, n int, a []float64, lda int, work []float64) float64 {
+	checkMatrix(n, n, a, lda)
+	switch norm {
+	case lapack.MaxRowSum, lapack.MaxColumnSum, lapack.NormFrob, lapack.MaxAbs:
+	default:
+		panic(badNorm)
+	}
+	if (norm == lapack.MaxColumnSum || norm == lapack.MaxRowSum) && len(work) < n {
+		panic(badWork)
+	}
+	if uplo != blas.Upper && uplo != blas.Lower {
+		panic(badUplo)
+	}
+	return clapack.Dlansy(byte(norm), uplo, n, a, lda)
+}
+
 // Dlantr computes the specified norm of an m×n trapezoidal matrix A. If
 // norm == lapack.MaxColumnSum work must have length at least n, otherwise work
 // is unused.

native/dlansy.go

@@ -0,0 +1,121 @@
+package native
+
+import (
+	"math"
+
+	"github.com/gonum/blas"
+	"github.com/gonum/lapack"
+)
+
+// Dlansy computes the specified norm of an n×n symmetric matrix. If
+// norm == lapack.MaxColumnSum or norm == lapackMaxRowSum work must have length
+// at least n, otherwise work is unused.
+func (impl Implementation) Dlansy(norm lapack.MatrixNorm, uplo blas.Uplo, n int, a []float64, lda int, work []float64) float64 {
+	checkMatrix(n, n, a, lda)
+	switch norm {
+	case lapack.MaxRowSum, lapack.MaxColumnSum, lapack.NormFrob, lapack.MaxAbs:
+	default:
+		panic(badNorm)
+	}
+	if (norm == lapack.MaxColumnSum || norm == lapack.MaxRowSum) && len(work) < n {
+		panic(badWork)
+	}
+	if uplo != blas.Upper && uplo != blas.Lower {
+		panic(badUplo)
+	}
+
+	if n == 0 {
+		return 0
+	}
+	switch norm {
+	default:
+		panic("unreachable")
+	case lapack.MaxAbs:
+		if uplo == blas.Upper {
+			var max float64
+			for i := 0; i < n; i++ {
+				for j := i; j < n; j++ {
+					v := math.Abs(a[i*lda+j])
+					if math.IsNaN(v) {
+						return math.NaN()
+					}
+					if v > max {
+						max = v
+					}
+				}
+			}
+			return max
+		}
+		var max float64
+		for i := 0; i < n; i++ {
+			for j := 0; j <= i; j++ {
+				v := math.Abs(a[i*lda+j])
+				if math.IsNaN(v) {
+					return math.NaN()
+				}
+				if v > max {
+					max = v
+				}
+			}
+		}
+		return max
+	case lapack.MaxRowSum, lapack.MaxColumnSum:
+		// A symmetric matrix has the same 1-norm and ∞-norm.
+		for i := 0; i < n; i++ {
+			work[i] = 0
+		}
+		if uplo == blas.Upper {
+			for i := 0; i < n; i++ {
+				work[i] += math.Abs(a[i*lda+i])
+				for j := i + 1; j < n; j++ {
+					v := math.Abs(a[i*lda+j])
+					work[i] += v
+					work[j] += v
+				}
+			}
+		} else {
+			for i := 0; i < n; i++ {
+				for j := 0; j < i; j++ {
+					v := math.Abs(a[i*lda+j])
+					work[i] += v
+					work[j] += v
+				}
+				work[i] += math.Abs(a[i*lda+i])
+			}
+		}
+		var max float64
+		for i := 0; i < n; i++ {
+			v := work[i]
+			if math.IsNaN(v) {
+				return math.NaN()
+			}
+			if v > max {
+				max = v
+			}
+		}
+		return max
+	case lapack.NormFrob:
+		if uplo == blas.Upper {
+			var sum float64
+			for i := 0; i < n; i++ {
+				v := a[i*lda+i]
+				sum += v * v
+				for j := i + 1; j < n; j++ {
+					v := a[i*lda+j]
+					sum += 2 * v * v
+				}
+			}
+			return math.Sqrt(sum)
+		}
+		var sum float64
+		for i := 0; i < n; i++ {
+			for j := 0; j < i; j++ {
+				v := a[i*lda+j]
+				sum += 2 * v * v
+			}
+			v := a[i*lda+i]
+			sum += v * v
+		}
+		return math.Sqrt(sum)
+	}
+}

native/lapack_test.go

@@ -48,6 +48,10 @@ func TestDlange(t *testing.T) {
 	testlapack.DlangeTest(t, impl)
 }
 
+func TestDlansy(t *testing.T) {
+	testlapack.DlansyTest(t, impl)
+}
+
 func TestDlantr(t *testing.T) {
 	testlapack.DlantrTest(t, impl)
 }

testlapack/dlansy.go

@@ -0,0 +1,75 @@
+package testlapack
+
+import (
+	"math"
+	"math/rand"
+	"testing"
+
+	"github.com/gonum/blas"
+	"github.com/gonum/lapack"
+)
+
+type Dlansyer interface {
+	Dlanger
+	Dlansy(norm lapack.MatrixNorm, uplo blas.Uplo, n int, a []float64, lda int, work []float64) float64
+}
+
+func DlansyTest(t *testing.T, impl Dlansyer) {
+	for _, norm := range []lapack.MatrixNorm{lapack.MaxAbs, lapack.MaxColumnSum, lapack.MaxRowSum, lapack.NormFrob} {
+		for _, uplo := range []blas.Uplo{blas.Lower, blas.Upper} {
+			for _, test := range []struct {
+				n, lda int
+			}{
+				{1, 0},
+				{3, 0},
+
+				{1, 10},
+				{3, 10},
+			} {
+				for trial := 0; trial < 100; trial++ {
+					n := test.n
+					lda := test.lda
+					if lda == 0 {
+						lda = n
+					}
+					a := make([]float64, lda*n)
+					if trial == 0 {
+						for i := range a {
+							a[i] = float64(i)
+						}
+					} else {
+						for i := range a {
+							a[i] = rand.NormFloat64()
+						}
+					}
+
+					aDense := make([]float64, n*n)
+					if uplo == blas.Upper {
+						for i := 0; i < n; i++ {
+							for j := i; j < n; j++ {
+								v := a[i*lda+j]
+								aDense[i*n+j] = v
+								aDense[j*n+i] = v
+							}
+						}
+					} else {
+						for i := 0; i < n; i++ {
+							for j := 0; j <= i; j++ {
+								v := a[i*lda+j]
+								aDense[i*n+j] = v
+								aDense[j*n+i] = v
+							}
+						}
+					}
+					work := make([]float64, n)
+					got := impl.Dlansy(norm, uplo, n, a, lda, work)
+					want := impl.Dlange(norm, n, n, aDense, n, work)
+					if math.Abs(want-got) > 1e-14 {
+						t.Errorf("Norm mismatch. norm = %c, upper = %v, n = %v, lda = %v, want %v, got %v.",
+							norm, uplo == blas.Upper, n, lda, got, want)
+					}
+				}
+			}
+		}
+	}
+}