mirror of
https://github.com/gonum/gonum.git
synced 2025-10-05 15:16:59 +08:00
testlapack: add test for Dtrtrs
This commit is contained in:
Vladimir Chalupecky
committed by
Vladimír Chalupecký
Vladimír Chalupecký
parent
971fc50f31
commit
aacdd93956
@@ -573,6 +573,11 @@ func TestDtrtri(t *testing.T) {
|
||||
testlapack.DtrtriTest(t, impl)
|
||||
}
|
||||
|
||||
func TestDtrtrs(t *testing.T) {
|
||||
t.Parallel()
|
||||
testlapack.DtrtrsTest(t, impl)
|
||||
}
|
||||
|
||||
func TestIladlc(t *testing.T) {
|
||||
t.Parallel()
|
||||
testlapack.IladlcTest(t, impl)
|
||||
|
||||
138
lapack/testlapack/dtrtrs.go
Normal file
138
lapack/testlapack/dtrtrs.go
Normal file
@@ -0,0 +1,138 @@
|
||||
// Copyright ©2020 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/floats"
|
||||
"gonum.org/v1/gonum/lapack"
|
||||
)
|
||||
|
||||
type Dtrtrser interface {
|
||||
Dtrtrs(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n, nrhs int, a []float64, lda int, b []float64, ldb int) bool
|
||||
}
|
||||
|
||||
func DtrtrsTest(t *testing.T, impl Dtrtrser) {
|
||||
rnd := rand.New(rand.NewSource(1))
|
||||
|
||||
for _, trans := range []blas.Transpose{blas.NoTrans, blas.Trans, blas.ConjTrans} {
|
||||
name := transToString(trans)
|
||||
t.Run(name, func(t *testing.T) {
|
||||
for _, uplo := range []blas.Uplo{blas.Upper, blas.Lower} {
|
||||
for _, diag := range []blas.Diag{blas.Unit, blas.NonUnit} {
|
||||
for _, n := range []int{0, 1, 2, 3, 4, 5, 10} {
|
||||
for _, nrhs := range []int{0, 1, 2, 3, 4, 5, 10, 15} {
|
||||
for _, lda := range []int{max(1, n), n + 3} {
|
||||
for _, ldb := range []int{max(1, nrhs), nrhs + 3} {
|
||||
if diag == blas.Unit {
|
||||
dtrtrsTest(t, impl, rnd, uplo, trans, diag, n, nrhs, lda, ldb, false)
|
||||
} else {
|
||||
dtrtrsTest(t, impl, rnd, uplo, trans, diag, n, nrhs, lda, ldb, true)
|
||||
dtrtrsTest(t, impl, rnd, uplo, trans, diag, n, nrhs, lda, ldb, false)
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
})
|
||||
}
|
||||
}
|
||||
|
||||
func dtrtrsTest(t *testing.T, impl Dtrtrser, rnd *rand.Rand, uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n, nrhs int, lda, ldb int, singular bool) {
|
||||
if singular && diag == blas.Unit {
|
||||
panic("blas.Unit triangular matrix cannot be singular")
|
||||
}
|
||||
|
||||
const tol = 1e-14
|
||||
|
||||
if n == 0 {
|
||||
singular = false
|
||||
}
|
||||
name := fmt.Sprintf("uplo=%v,diag=%v,n=%v,nrhs=%v,lda=%v,ldb=%v,sing=%v", string(uplo), string(diag), n, nrhs, lda, ldb, singular)
|
||||
|
||||
// Generate a random triangular matrix A. One of its triangles won't be
|
||||
// referenced.
|
||||
a := make([]float64, n*lda)
|
||||
for i := range a {
|
||||
a[i] = rnd.NormFloat64()
|
||||
}
|
||||
if singular {
|
||||
i := rnd.Intn(n)
|
||||
a[i*lda+i] = 0
|
||||
}
|
||||
aCopy := make([]float64, len(a))
|
||||
copy(aCopy, a)
|
||||
|
||||
// Generate a random solution matrix X.
|
||||
x := make([]float64, n*ldb)
|
||||
for i := range x {
|
||||
x[i] = rnd.NormFloat64()
|
||||
}
|
||||
|
||||
// Generate the right-hand side as A * X or Aᵀ * X.
|
||||
b := make([]float64, len(x))
|
||||
copy(b, x)
|
||||
bi := blas64.Implementation()
|
||||
bi.Dtrmm(blas.Left, uplo, trans, diag, n, nrhs, 1, a, lda, b, ldb)
|
||||
|
||||
got := make([]float64, len(b))
|
||||
copy(got, b)
|
||||
|
||||
ok := impl.Dtrtrs(uplo, trans, diag, n, nrhs, a, lda, got, ldb)
|
||||
|
||||
if !floats.Equal(a, aCopy) {
|
||||
t.Errorf("%v: unexpected modification of A", name)
|
||||
}
|
||||
|
||||
if ok == singular {
|
||||
t.Errorf("%v: misdetected singular matrix, ok=%v", name, ok)
|
||||
}
|
||||
|
||||
if !ok {
|
||||
if !floats.Equal(got, b) {
|
||||
t.Errorf("%v: unexpected modification of B when singular", name)
|
||||
}
|
||||
return
|
||||
}
|
||||
|
||||
if n == 0 || nrhs == 0 {
|
||||
return
|
||||
}
|
||||
|
||||
work := make([]float64, n)
|
||||
|
||||
// Compute the 1-norm of A or Aᵀ.
|
||||
var aNorm float64
|
||||
if trans == blas.NoTrans {
|
||||
aNorm = dlantr(lapack.MaxColumnSum, uplo, diag, n, n, a, lda, work)
|
||||
} else {
|
||||
aNorm = dlantr(lapack.MaxRowSum, uplo, diag, n, n, a, lda, work)
|
||||
}
|
||||
|
||||
// Compute the maximum over the number of right-hand sides of
|
||||
// |op(A)*x-b| / (|op(A)| * |x|)
|
||||
var resid float64
|
||||
for j := 0; j < nrhs; j++ {
|
||||
bi.Dcopy(n, got[j:], ldb, work, 1)
|
||||
bi.Dtrmv(uplo, trans, diag, n, a, lda, work, 1)
|
||||
bi.Daxpy(n, -1, b[j:], ldb, work, 1)
|
||||
rjNorm := bi.Dasum(n, work, 1)
|
||||
xNorm := bi.Dasum(n, got[j:], ldb)
|
||||
resid = math.Max(resid, rjNorm/aNorm/xNorm)
|
||||
}
|
||||
t.Log(name, resid)
|
||||
if resid > tol {
|
||||
t.Errorf("%v: unexpected result; resid=%v,want<=%v", name, resid, tol)
|
||||
}
|
||||
}
|
||||
@@ -110,6 +110,19 @@ func sideToString(side blas.Side) string {
|
||||
}
|
||||
}
|
||||
|
||||
func transToString(trans blas.Transpose) string {
|
||||
switch trans {
|
||||
case blas.NoTrans:
|
||||
return "NoTrans"
|
||||
case blas.Trans:
|
||||
return "Trans"
|
||||
case blas.ConjTrans:
|
||||
return "ConjTrans"
|
||||
default:
|
||||
panic("invalid trans")
|
||||
}
|
||||
}
|
||||
|
||||
// nanSlice allocates a new slice of length n filled with NaN.
|
||||
func nanSlice(n int) []float64 {
|
||||
s := make([]float64, n)
|
||||
@@ -1565,3 +1578,191 @@ func dlansb(norm lapack.MatrixNorm, uplo blas.Uplo, n, kd int, ab []float64, lda
|
||||
}
|
||||
return value
|
||||
}
|
||||
|
||||
func dlantr(norm lapack.MatrixNorm, uplo blas.Uplo, diag blas.Diag, m, n int, a []float64, lda int, work []float64) float64 {
|
||||
// Quick return if possible.
|
||||
minmn := min(m, n)
|
||||
if minmn == 0 {
|
||||
return 0
|
||||
}
|
||||
switch norm {
|
||||
case lapack.MaxAbs:
|
||||
if diag == blas.Unit {
|
||||
value := 1.0
|
||||
if uplo == blas.Upper {
|
||||
for i := 0; i < m; i++ {
|
||||
for j := i + 1; j < n; j++ {
|
||||
tmp := math.Abs(a[i*lda+j])
|
||||
if math.IsNaN(tmp) {
|
||||
return tmp
|
||||
}
|
||||
if tmp > value {
|
||||
value = tmp
|
||||
}
|
||||
}
|
||||
}
|
||||
return value
|
||||
}
|
||||
for i := 1; i < m; i++ {
|
||||
for j := 0; j < min(i, n); j++ {
|
||||
tmp := math.Abs(a[i*lda+j])
|
||||
if math.IsNaN(tmp) {
|
||||
return tmp
|
||||
}
|
||||
if tmp > value {
|
||||
value = tmp
|
||||
}
|
||||
}
|
||||
}
|
||||
return value
|
||||
}
|
||||
var value float64
|
||||
if uplo == blas.Upper {
|
||||
for i := 0; i < m; i++ {
|
||||
for j := i; j < n; j++ {
|
||||
tmp := math.Abs(a[i*lda+j])
|
||||
if math.IsNaN(tmp) {
|
||||
return tmp
|
||||
}
|
||||
if tmp > value {
|
||||
value = tmp
|
||||
}
|
||||
}
|
||||
}
|
||||
return value
|
||||
}
|
||||
for i := 0; i < m; i++ {
|
||||
for j := 0; j <= min(i, n-1); j++ {
|
||||
tmp := math.Abs(a[i*lda+j])
|
||||
if math.IsNaN(tmp) {
|
||||
return tmp
|
||||
}
|
||||
if tmp > value {
|
||||
value = tmp
|
||||
}
|
||||
}
|
||||
}
|
||||
return value
|
||||
case lapack.MaxColumnSum:
|
||||
if diag == blas.Unit {
|
||||
for i := 0; i < minmn; i++ {
|
||||
work[i] = 1
|
||||
}
|
||||
for i := minmn; i < n; i++ {
|
||||
work[i] = 0
|
||||
}
|
||||
if uplo == blas.Upper {
|
||||
for i := 0; i < m; i++ {
|
||||
for j := i + 1; j < n; j++ {
|
||||
work[j] += math.Abs(a[i*lda+j])
|
||||
}
|
||||
}
|
||||
} else {
|
||||
for i := 1; i < m; i++ {
|
||||
for j := 0; j < min(i, n); j++ {
|
||||
work[j] += math.Abs(a[i*lda+j])
|
||||
}
|
||||
}
|
||||
}
|
||||
} else {
|
||||
for i := 0; i < n; i++ {
|
||||
work[i] = 0
|
||||
}
|
||||
if uplo == blas.Upper {
|
||||
for i := 0; i < m; i++ {
|
||||
for j := i; j < n; j++ {
|
||||
work[j] += math.Abs(a[i*lda+j])
|
||||
}
|
||||
}
|
||||
} else {
|
||||
for i := 0; i < m; i++ {
|
||||
for j := 0; j <= min(i, n-1); j++ {
|
||||
work[j] += math.Abs(a[i*lda+j])
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
var max float64
|
||||
for _, v := range work[:n] {
|
||||
if math.IsNaN(v) {
|
||||
return math.NaN()
|
||||
}
|
||||
if v > max {
|
||||
max = v
|
||||
}
|
||||
}
|
||||
return max
|
||||
case lapack.MaxRowSum:
|
||||
var maxsum float64
|
||||
if diag == blas.Unit {
|
||||
if uplo == blas.Upper {
|
||||
for i := 0; i < m; i++ {
|
||||
var sum float64
|
||||
if i < minmn {
|
||||
sum = 1
|
||||
}
|
||||
for j := i + 1; j < n; j++ {
|
||||
sum += math.Abs(a[i*lda+j])
|
||||
}
|
||||
if math.IsNaN(sum) {
|
||||
return math.NaN()
|
||||
}
|
||||
if sum > maxsum {
|
||||
maxsum = sum
|
||||
}
|
||||
}
|
||||
return maxsum
|
||||
} else {
|
||||
for i := 0; i < m; i++ {
|
||||
var sum float64
|
||||
if i < minmn {
|
||||
sum = 1
|
||||
}
|
||||
for j := 0; j < min(i, n); j++ {
|
||||
sum += math.Abs(a[i*lda+j])
|
||||
}
|
||||
if math.IsNaN(sum) {
|
||||
return math.NaN()
|
||||
}
|
||||
if sum > maxsum {
|
||||
maxsum = sum
|
||||
}
|
||||
}
|
||||
return maxsum
|
||||
}
|
||||
} else {
|
||||
if uplo == blas.Upper {
|
||||
for i := 0; i < m; i++ {
|
||||
var sum float64
|
||||
for j := i; j < n; j++ {
|
||||
sum += math.Abs(a[i*lda+j])
|
||||
}
|
||||
if math.IsNaN(sum) {
|
||||
return sum
|
||||
}
|
||||
if sum > maxsum {
|
||||
maxsum = sum
|
||||
}
|
||||
}
|
||||
return maxsum
|
||||
} else {
|
||||
for i := 0; i < m; i++ {
|
||||
var sum float64
|
||||
for j := 0; j <= min(i, n-1); j++ {
|
||||
sum += math.Abs(a[i*lda+j])
|
||||
}
|
||||
if math.IsNaN(sum) {
|
||||
return sum
|
||||
}
|
||||
if sum > maxsum {
|
||||
maxsum = sum
|
||||
}
|
||||
}
|
||||
return maxsum
|
||||
}
|
||||
}
|
||||
default:
|
||||
// lapack.Frobenius
|
||||
panic("not implemented")
|
||||
}
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user