1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
|
// 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/floats"
"gonum.org/v1/gonum/lapack"
)
type Dgesc2er interface {
Dgesc2(n int, a []float64, lda int, rhs []float64, ipiv, jpiv []int) (scale float64)
Dgetc2er
}
func Dgesc2Test(t *testing.T, impl Dgesc2er) {
rnd := rand.New(rand.NewSource(1))
for _, n := range []int{0, 1, 2, 3, 4, 5, 10, 20, 50} {
for _, lda := range []int{n, n + 3} {
testDgesc2(t, impl, rnd, n, lda, false)
testDgesc2(t, impl, rnd, n, lda, true)
}
}
}
func testDgesc2(t *testing.T, impl Dgesc2er, rnd *rand.Rand, n, lda int, big bool) {
const tol = 1e-14
name := fmt.Sprintf("n=%v,lda=%v,big=%v", n, lda, big)
// Generate random general matrix.
a := randomGeneral(n, n, max(1, lda), rnd)
// Generate a random right hand side vector.
b := randomGeneral(n, 1, 1, rnd)
if big {
for i := 0; i < n; i++ {
b.Data[i] *= bignum
}
}
// Compute the LU factorization of A with full pivoting.
lu := cloneGeneral(a)
ipiv := make([]int, n)
jpiv := make([]int, n)
impl.Dgetc2(n, lu.Data, lu.Stride, ipiv, jpiv)
// Make copies of const input to Dgesc2.
luCopy := cloneGeneral(lu)
ipivCopy := make([]int, len(ipiv))
copy(ipivCopy, ipiv)
jpivCopy := make([]int, len(jpiv))
copy(jpivCopy, jpiv)
// Call Dgesc2 to solve A*x = scale*b.
x := cloneGeneral(b)
scale := impl.Dgesc2(n, lu.Data, lu.Stride, x.Data, ipiv, jpiv)
if n == 0 {
return
}
// Check that const input to Dgesc2 hasn't been modified.
if !floats.Same(lu.Data, luCopy.Data) {
t.Errorf("%v: unexpected modification in lu", name)
}
if !intsEqual(ipiv, ipivCopy) {
t.Errorf("%v: unexpected modification in ipiv", name)
}
if !intsEqual(jpiv, jpivCopy) {
t.Errorf("%v: unexpected modification in jpiv", name)
}
if scale <= 0 || 1 < scale {
t.Errorf("%v: scale %v out of bounds (0,1]", name, scale)
}
if !big && scale != 1 {
t.Errorf("%v: unexpected scaling, scale=%v", name, scale)
}
// Compute the difference rhs := A*x - scale*b.
diff := b
for i := 0; i < n; i++ {
diff.Data[i] *= scale
}
blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, a, x, -1, diff)
// Compute the residual |A*x - scale*b| / |x|.
xnorm := dlange(lapack.MaxColumnSum, n, 1, x.Data, 1)
resid := dlange(lapack.MaxColumnSum, n, 1, diff.Data, 1) / xnorm
if resid > tol || math.IsNaN(resid) {
t.Errorf("%v: unexpected result; resid=%v, want<=%v", name, resid, tol)
}
}
|
