File: dlascl.go

package info (click to toggle)
golang-gonum-v1-gonum 0.15.1-1
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid, trixie
  • size: 18,792 kB
  • sloc: asm: 6,252; fortran: 5,271; sh: 377; ruby: 211; makefile: 98
file content (122 lines) | stat: -rw-r--r-- 2,949 bytes parent folder | download
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
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
// Copyright ©2016 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/lapack"
)

type Dlascler interface {
	Dlascl(kind lapack.MatrixType, kl, ku int, cfrom, cto float64, m, n int, a []float64, lda int)
}

func DlasclTest(t *testing.T, impl Dlascler) {
	const tol = 1e-15

	rnd := rand.New(rand.NewSource(1))
	for ti, test := range []struct {
		m, n int
	}{
		{0, 0},
		{1, 1},
		{1, 10},
		{10, 1},
		{2, 2},
		{2, 11},
		{11, 2},
		{3, 3},
		{3, 11},
		{11, 3},
		{11, 11},
		{11, 100},
		{100, 11},
	} {
		m := test.m
		n := test.n
		for _, extra := range []int{0, 11} {
			for _, kind := range []lapack.MatrixType{lapack.General, lapack.UpperTri, lapack.LowerTri} {
				a := randomGeneral(m, n, n+extra, rnd)
				aCopy := cloneGeneral(a)
				cfrom := rnd.NormFloat64()
				cto := rnd.NormFloat64()
				scale := cto / cfrom

				impl.Dlascl(kind, -1, -1, cfrom, cto, m, n, a.Data, a.Stride)

				prefix := fmt.Sprintf("Case #%v: kind=%v,m=%v,n=%v,extra=%v", ti, kind, m, n, extra)
				if !generalOutsideAllNaN(a) {
					t.Errorf("%v: out-of-range write to A", prefix)
				}
				switch kind {
				case lapack.UpperTri:
					var mod bool
				loopLower:
					for i := 0; i < m; i++ {
						for j := 0; j < min(i, n); j++ {
							if a.Data[i*a.Stride+j] != aCopy.Data[i*aCopy.Stride+j] {
								mod = true
								break loopLower
							}
						}
					}
					if mod {
						t.Errorf("%v: unexpected modification in lower triangle of A", prefix)
					}
				case lapack.LowerTri:
					var mod bool
				loopUpper:
					for i := 0; i < m; i++ {
						for j := i + 1; j < n; j++ {
							if a.Data[i*a.Stride+j] != aCopy.Data[i*aCopy.Stride+j] {
								mod = true
								break loopUpper
							}
						}
					}
					if mod {
						t.Errorf("%v: unexpected modification in upper triangle of A", prefix)
					}
				}

				var resid float64
				switch kind {
				case lapack.General:
					for i := 0; i < m; i++ {
						for j := 0; j < n; j++ {
							want := scale * aCopy.Data[i*aCopy.Stride+j]
							got := a.Data[i*a.Stride+j]
							resid = math.Max(resid, math.Abs(want-got))
						}
					}
				case lapack.UpperTri:
					for i := 0; i < m; i++ {
						for j := i; j < n; j++ {
							want := scale * aCopy.Data[i*aCopy.Stride+j]
							got := a.Data[i*a.Stride+j]
							resid = math.Max(resid, math.Abs(want-got))
						}
					}
				case lapack.LowerTri:
					for i := 0; i < m; i++ {
						for j := 0; j <= min(i, n-1); j++ {
							want := scale * aCopy.Data[i*aCopy.Stride+j]
							got := a.Data[i*a.Stride+j]
							resid = math.Max(resid, math.Abs(want-got))
						}
					}
				}
				if resid > tol*float64(max(m, n)) {
					t.Errorf("%v: unexpected result; residual=%v, want<=%v", prefix, resid, tol*float64(max(m, n)))
				}
			}
		}
	}
}