2245

Committed 23 Jul 2018 - 10:21 coverage increased (+3.5%) to 75.543%

Build # 2245

Build Type

Pull #556

travis-ci

Committed by

web-flow

Commit Message

lapack/testlapack: add test for randomOrthogonal

... and remove the orthogonality assertion in randomOrthogonal.

Pull Request Pull Request #556: lapack/testlapack: clean up and use isOrthogonal helper

Run Details

42344 of 56053 relevant lines covered (75.54%)

13.9 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

98.44

/lapack/gonum/dlatrs.go

// Copyright ©2015 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 gonum

import (
        "math"

        "gonum.org/v1/gonum/blas"
        "gonum.org/v1/gonum/blas/blas64"
)

// Dlatrs solves a triangular system of equations scaled to prevent overflow. It
// solves
//  A * x = scale * b if trans == blas.NoTrans
//  A^T * x = scale * b if trans == blas.Trans
// where the scale s is set for numeric stability.
//
// A is an n×n triangular matrix. On entry, the slice x contains the values of
// of b, and on exit it contains the solution vector x.
//
// If normin == true, cnorm is an input and cnorm[j] contains the norm of the off-diagonal
// part of the j^th column of A. If trans == blas.NoTrans, cnorm[j] must be greater
// than or equal to the infinity norm, and greater than or equal to the one-norm
// otherwise. If normin == false, then cnorm is treated as an output, and is set
// to contain the 1-norm of the off-diagonal part of the j^th column of A.
//
// Dlatrs is an internal routine. It is exported for testing purposes.
func (impl Implementation) Dlatrs(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, normin bool, n int, a []float64, lda int, x []float64, cnorm []float64) (scale float64) {
        if uplo != blas.Upper && uplo != blas.Lower {
                panic(badUplo)
        }
        if trans != blas.Trans && trans != blas.NoTrans {
                panic(badTrans)
        }
        if diag != blas.Unit && diag != blas.NonUnit {
                panic(badDiag)
        }
        upper := uplo == blas.Upper
        noTrans := trans == blas.NoTrans
        nonUnit := diag == blas.NonUnit

        if n < 0 {
                panic(nLT0)
        }
        checkMatrix(n, n, a, lda)
        checkVector(n, x, 1)
        checkVector(n, cnorm, 1)

        if n == 0 {
                return 0
        }
        smlnum := dlamchS / dlamchP
        bignum := 1 / smlnum
        scale = 1
        bi := blas64.Implementation()
        if !normin {
                if upper {
                        cnorm[0] = 0
                        for j := 1; j < n; j++ {
                                cnorm[j] = bi.Dasum(j, a[j:], lda)
                        }
                } else {
                        for j := 0; j < n-1; j++ {
                                cnorm[j] = bi.Dasum(n-j-1, a[(j+1)*lda+j:], lda)
                        }
                        cnorm[n-1] = 0
                }
        }
        // Scale the column norms by tscal if the maximum element in cnorm is greater than bignum.
        imax := bi.Idamax(n, cnorm, 1)
        tmax := cnorm[imax]
        var tscal float64
        if tmax <= bignum {
                tscal = 1
        } else {
                tscal = 1 / (smlnum * tmax)
                bi.Dscal(n, tscal, cnorm, 1)
        }

        // Compute a bound on the computed solution vector to see if bi.Dtrsv can be used.
        j := bi.Idamax(n, x, 1)
        xmax := math.Abs(x[j])
        xbnd := xmax
        var grow float64
        var jfirst, jlast, jinc int
        if noTrans {
                if upper {
                        jfirst = n - 1
                        jlast = -1
                        jinc = -1
                } else {
                        jfirst = 0
                        jlast = n
                        jinc = 1
                }
                // Compute the growth in A * x = b.
                if tscal != 1 {
                        grow = 0
                        goto Solve
                }
                if nonUnit {
                        grow = 1 / math.Max(xbnd, smlnum)
                        xbnd = grow
                        for j := jfirst; j != jlast; j += jinc {
                                if grow <= smlnum {
                                        goto Solve
                                }
                                tjj := math.Abs(a[j*lda+j])
                                xbnd = math.Min(xbnd, math.Min(1, tjj)*grow)
                                if tjj+cnorm[j] >= smlnum {
                                        grow *= tjj / (tjj + cnorm[j])
                                } else {
                                        grow = 0
                                }
                        }
                        grow = xbnd
                } else {
                        grow = math.Min(1, 1/math.Max(xbnd, smlnum))
                        for j := jfirst; j != jlast; j += jinc {
                                if grow <= smlnum {
                                        goto Solve
                                }
                                grow *= 1 / (1 + cnorm[j])
                        }
                }
        } else {
                if upper {
                        jfirst = 0
                        jlast = n
                        jinc = 1
                } else {
                        jfirst = n - 1
                        jlast = -1
                        jinc = -1
                }
                if tscal != 1 {
                        grow = 0
                        goto Solve
                }
                if nonUnit {
                        grow = 1 / (math.Max(xbnd, smlnum))
                        xbnd = grow
                        for j := jfirst; j != jlast; j += jinc {
                                if grow <= smlnum {
                                        goto Solve
                                }
                                xj := 1 + cnorm[j]
                                grow = math.Min(grow, xbnd/xj)
                                tjj := math.Abs(a[j*lda+j])
                                if xj > tjj {
                                        xbnd *= tjj / xj
                                }
                        }
                        grow = math.Min(grow, xbnd)
                } else {
                        grow = math.Min(1, 1/math.Max(xbnd, smlnum))
                        for j := jfirst; j != jlast; j += jinc {
                                if grow <= smlnum {
                                        goto Solve
                                }
                                xj := 1 + cnorm[j]
                                grow /= xj
                        }
                }
        }

Solve:
        if grow*tscal > smlnum {
                // Use the Level 2 BLAS solve if the reciprocal of the bound on
                // elements of X is not too small.
                bi.Dtrsv(uplo, trans, diag, n, a, lda, x, 1)
                if tscal != 1 {
                        bi.Dscal(n, 1/tscal, cnorm, 1)
                }
                return scale
        }

        // Use a Level 1 BLAS solve, scaling intermediate results.
        if xmax > bignum {
                scale = bignum / xmax
                bi.Dscal(n, scale, x, 1)
                xmax = bignum
        }
        if noTrans {
                for j := jfirst; j != jlast; j += jinc {
                        xj := math.Abs(x[j])
                        var tjj, tjjs float64
                        if nonUnit {
                                tjjs = a[j*lda+j] * tscal
                        } else {
                                tjjs = tscal
                                if tscal == 1 {
                                        goto Skip1
                                }
                        }
                        tjj = math.Abs(tjjs)
                        if tjj > smlnum {
                                if tjj < 1 {
                                        if xj > tjj*bignum {
                                                rec := 1 / xj
                                                bi.Dscal(n, rec, x, 1)
                                                scale *= rec
                                                xmax *= rec
                                        }
                                }
                                x[j] /= tjjs
                                xj = math.Abs(x[j])
                        } else if tjj > 0 {
                                if xj > tjj*bignum {
                                        rec := (tjj * bignum) / xj
                                        if cnorm[j] > 1 {
                                                rec /= cnorm[j]
                                        }
                                        bi.Dscal(n, rec, x, 1)
                                        scale *= rec
                                        xmax *= rec
                                }
                                x[j] /= tjjs
                                xj = math.Abs(x[j])
                        } else {
                                for i := 0; i < n; i++ {
                                        x[i] = 0
                                }
                                x[j] = 1
                                xj = 1
                                scale = 0
                                xmax = 0
                        }
                Skip1:
                        if xj > 1 {
                                rec := 1 / xj
                                if cnorm[j] > (bignum-xmax)*rec {
                                        rec *= 0.5
                                        bi.Dscal(n, rec, x, 1)
                                        scale *= rec
                                }
                        } else if xj*cnorm[j] > bignum-xmax {
                                bi.Dscal(n, 0.5, x, 1)
                                scale *= 0.5
                        }
                        if upper {
                                if j > 0 {
                                        bi.Daxpy(j, -x[j]*tscal, a[j:], lda, x, 1)
                                        i := bi.Idamax(j, x, 1)
                                        xmax = math.Abs(x[i])
                                }
                        } else {
                                if j < n-1 {
                                        bi.Daxpy(n-j-1, -x[j]*tscal, a[(j+1)*lda+j:], lda, x[j+1:], 1)
                                        i := j + bi.Idamax(n-j-1, x[j+1:], 1)
                                        xmax = math.Abs(x[i])
                                }
                        }
                }
        } else {
                for j := jfirst; j != jlast; j += jinc {
                        xj := math.Abs(x[j])
                        uscal := tscal
                        rec := 1 / math.Max(xmax, 1)
                        var tjjs float64
                        if cnorm[j] > (bignum-xj)*rec {
                                rec *= 0.5
                                if nonUnit {
                                        tjjs = a[j*lda+j] * tscal
                                } else {
                                        tjjs = tscal
                                }
                                tjj := math.Abs(tjjs)
                                if tjj > 1 {
                                        rec = math.Min(1, rec*tjj)
                                        uscal /= tjjs
                                }
                                if rec < 1 {
                                        bi.Dscal(n, rec, x, 1)
                                        scale *= rec
                                        xmax *= rec
                                }
                        }
                        var sumj float64
                        if uscal == 1 {
                                if upper {
                                        sumj = bi.Ddot(j, a[j:], lda, x, 1)
                                } else if j < n-1 {
                                        sumj = bi.Ddot(n-j-1, a[(j+1)*lda+j:], lda, x[j+1:], 1)
                                }
                        } else {
                                if upper {
                                        for i := 0; i < j; i++ {
                                                sumj += (a[i*lda+j] * uscal) * x[i]
                                        }
                                } else if j < n {
                                        for i := j + 1; i < n; i++ {
                                                sumj += (a[i*lda+j] * uscal) * x[i]
                                        }
                                }
                        }
                        if uscal == tscal {
                                x[j] -= sumj
                                xj := math.Abs(x[j])
                                var tjjs float64
                                if nonUnit {
                                        tjjs = a[j*lda+j] * tscal
                                } else {
                                        tjjs = tscal
                                        if tscal == 1 {
                                                goto Skip2
                                        }
                                }
                                tjj := math.Abs(tjjs)
                                if tjj > smlnum {
                                        if tjj < 1 {
                                                if xj > tjj*bignum {
                                                        rec = 1 / xj
                                                        bi.Dscal(n, rec, x, 1)
                                                        scale *= rec
                                                        xmax *= rec
                                                }
                                        }
                                        x[j] /= tjjs
                                } else if tjj > 0 {
                                        if xj > tjj*bignum {
                                                rec = (tjj * bignum) / xj
                                                bi.Dscal(n, rec, x, 1)
                                                scale *= rec
                                                xmax *= rec
                                        }
                                        x[j] /= tjjs
                                } else {
                                        for i := 0; i < n; i++ {
                                                x[i] = 0
                                        }
                                        x[j] = 1
                                        scale = 0
                                        xmax = 0
                                }
                        } else {
                                x[j] = x[j]/tjjs - sumj
                        }
                Skip2:
                        xmax = math.Max(xmax, math.Abs(x[j]))
                }
        }
        scale /= tscal
        if tscal != 1 {
                bi.Dscal(n, 1/tscal, cnorm, 1)
        }
        return scale
}

1	// Copyright ©2015 The Gonum Authors. All rights reserved.
2	// Use of this source code is governed by a BSD-style
3	// license that can be found in the LICENSE file.
4
5	package gonum
6
7	import (
8	"math"
9
10	"gonum.org/v1/gonum/blas"
11	"gonum.org/v1/gonum/blas/blas64"
12	)
13
14	// Dlatrs solves a triangular system of equations scaled to prevent overflow. It
15	// solves
16	// A * x = scale * b if trans == blas.NoTrans
17	// A^T * x = scale * b if trans == blas.Trans
18	// where the scale s is set for numeric stability.
19	//
20	// A is an n×n triangular matrix. On entry, the slice x contains the values of
21	// of b, and on exit it contains the solution vector x.
22	//
23	// If normin == true, cnorm is an input and cnorm[j] contains the norm of the off-diagonal
24	// part of the j^th column of A. If trans == blas.NoTrans, cnorm[j] must be greater
25	// than or equal to the infinity norm, and greater than or equal to the one-norm
26	// otherwise. If normin == false, then cnorm is treated as an output, and is set
27	// to contain the 1-norm of the off-diagonal part of the j^th column of A.
28	//
29	// Dlatrs is an internal routine. It is exported for testing purposes.
30	func (impl Implementation) Dlatrs(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, normin bool, n int, a []float64, lda int, x []float64, cnorm []float64) (scale float64) {	16×
31	if uplo != blas.Upper && uplo != blas.Lower {	16×
32	panic(badUplo)	!
33	}
34	if trans != blas.Trans && trans != blas.NoTrans {	16×
35	panic(badTrans)	!
36	}
37	if diag != blas.Unit && diag != blas.NonUnit {	16×
38	panic(badDiag)	!
39	}
40	upper := uplo == blas.Upper	16×
41	noTrans := trans == blas.NoTrans	16×
42	nonUnit := diag == blas.NonUnit	16×
43		16×
44	if n < 0 {	16×
45	panic(nLT0)	!
46	}
47	checkMatrix(n, n, a, lda)	16×
48	checkVector(n, x, 1)	16×
49	checkVector(n, cnorm, 1)	16×
50		16×
51	if n == 0 {	32×
52	return 0	16×
53	}	16×
54	smlnum := dlamchS / dlamchP	16×
55	bignum := 1 / smlnum	16×
56	scale = 1	16×
57	bi := blas64.Implementation()	16×
58	if !normin {	32×
59	if upper {	32×
60	cnorm[0] = 0	16×
61	for j := 1; j < n; j++ {	32×
62	cnorm[j] = bi.Dasum(j, a[j:], lda)	16×
63	}	16×
64	} else {	16×
65	for j := 0; j < n-1; j++ {	32×
66	cnorm[j] = bi.Dasum(n-j-1, a[(j+1)*lda+j:], lda)	16×
67	}	16×
68	cnorm[n-1] = 0	16×
69	}
70	}
71	// Scale the column norms by tscal if the maximum element in cnorm is greater than bignum.
72	imax := bi.Idamax(n, cnorm, 1)	16×
73	tmax := cnorm[imax]	16×
74	var tscal float64	16×
75	if tmax <= bignum {	32×
76	tscal = 1	16×
77	} else {	32×
78	tscal = 1 / (smlnum * tmax)	16×
79	bi.Dscal(n, tscal, cnorm, 1)	16×
80	}	16×
81
82	// Compute a bound on the computed solution vector to see if bi.Dtrsv can be used.
83	j := bi.Idamax(n, x, 1)	16×
84	xmax := math.Abs(x[j])	16×
85	xbnd := xmax	16×
86	var grow float64	16×
87	var jfirst, jlast, jinc int	16×
88	if noTrans {	32×
89	if upper {	32×
90	jfirst = n - 1	16×
91	jlast = -1	16×
92	jinc = -1	16×
93	} else {	32×
94	jfirst = 0	16×
95	jlast = n	16×
96	jinc = 1	16×
97	}	16×
98	// Compute the growth in A * x = b.
99	if tscal != 1 {	32×
100	grow = 0	16×
101	goto Solve	16×
102	}
103	if nonUnit {	32×
104	grow = 1 / math.Max(xbnd, smlnum)	16×
105	xbnd = grow	16×
106	for j := jfirst; j != jlast; j += jinc {	32×
107	if grow <= smlnum {	32×
108	goto Solve	16×
109	}
110	tjj := math.Abs(a[j*lda+j])	16×
111	xbnd = math.Min(xbnd, math.Min(1, tjj)*grow)	16×
112	if tjj+cnorm[j] >= smlnum {	32×
113	grow *= tjj / (tjj + cnorm[j])	16×
114	} else {	32×
115	grow = 0	16×
116	}	16×
117	}
118	grow = xbnd	16×
119	} else {	16×
120	grow = math.Min(1, 1/math.Max(xbnd, smlnum))	16×
121	for j := jfirst; j != jlast; j += jinc {	32×
122	if grow <= smlnum {	32×
123	goto Solve	16×
124	}
125	grow *= 1 / (1 + cnorm[j])	16×
126	}
127	}
128	} else {	16×
129	if upper {	32×
130	jfirst = 0	16×
131	jlast = n	16×
132	jinc = 1	16×
133	} else {	32×
134	jfirst = n - 1	16×
135	jlast = -1	16×
136	jinc = -1	16×
137	}	16×
138	if tscal != 1 {	32×
139	grow = 0	16×
140	goto Solve	16×
141	}
142	if nonUnit {	32×
143	grow = 1 / (math.Max(xbnd, smlnum))	16×
144	xbnd = grow	16×
145	for j := jfirst; j != jlast; j += jinc {	32×
146	if grow <= smlnum {	32×
147	goto Solve	16×
148	}
149	xj := 1 + cnorm[j]	16×
150	grow = math.Min(grow, xbnd/xj)	16×
151	tjj := math.Abs(a[j*lda+j])	16×
152	if xj > tjj {	32×
153	xbnd *= tjj / xj	16×
154	}	16×
155	}
156	grow = math.Min(grow, xbnd)	16×
157	} else {	16×
158	grow = math.Min(1, 1/math.Max(xbnd, smlnum))	16×
159	for j := jfirst; j != jlast; j += jinc {	32×
160	if grow <= smlnum {	32×
161	goto Solve	16×
162	}
163	xj := 1 + cnorm[j]	16×
164	grow /= xj	16×
165	}
166	}
167	}
168
169	Solve:
170	if grow*tscal > smlnum {	24×
171	// Use the Level 2 BLAS solve if the reciprocal of the bound on
172	// elements of X is not too small.
173	bi.Dtrsv(uplo, trans, diag, n, a, lda, x, 1)
174	if tscal != 1 {
UNCOV 175	bi.Dscal(n, 1/tscal, cnorm, 1)
UNCOV 176	}
177	return scale
178	}
179
180	// Use a Level 1 BLAS solve, scaling intermediate results.
181	if xmax > bignum {	32×
182	scale = bignum / xmax	16×
183	bi.Dscal(n, scale, x, 1)	16×
184	xmax = bignum	16×
185	}	16×
186	if noTrans {	32×
187	for j := jfirst; j != jlast; j += jinc {	32×
188	xj := math.Abs(x[j])	16×
189	var tjj, tjjs float64	16×
190	if nonUnit {	32×
191	tjjs = a[jlda+j] tscal	16×
192	} else {	32×
193	tjjs = tscal	16×
194	if tscal == 1 {	32×
195	goto Skip1	16×
196	}
197	}
198	tjj = math.Abs(tjjs)	16×
199	if tjj > smlnum {	32×
200	if tjj < 1 {	32×
201	if xj > tjj*bignum {	32×
202	rec := 1 / xj	16×
203	bi.Dscal(n, rec, x, 1)	16×
204	scale *= rec	16×
205	xmax *= rec	16×
206	}	16×
207	}
208	x[j] /= tjjs	16×
209	xj = math.Abs(x[j])	16×
210	} else if tjj > 0 {	32×
211	if xj > tjj*bignum {	32×
212	rec := (tjj * bignum) / xj	16×
213	if cnorm[j] > 1 {	32×
214	rec /= cnorm[j]	16×
215	}	16×
216	bi.Dscal(n, rec, x, 1)	16×
217	scale *= rec	16×
218	xmax *= rec	16×
219	}
220	x[j] /= tjjs	16×
221	xj = math.Abs(x[j])	16×
222	} else {	16×
223	for i := 0; i < n; i++ {	32×
224	x[i] = 0	16×
225	}	16×
226	x[j] = 1	16×
227	xj = 1	16×
228	scale = 0	16×
229	xmax = 0	16×
230	}
231	Skip1:
232	if xj > 1 {	32×
233	rec := 1 / xj	16×
234	if cnorm[j] > (bignum-xmax)*rec {	32×
235	rec *= 0.5	16×
236	bi.Dscal(n, rec, x, 1)	16×
237	scale *= rec	16×
238	}	16×
239	} else if xj*cnorm[j] > bignum-xmax {	32×
240	bi.Dscal(n, 0.5, x, 1)	16×
241	scale *= 0.5	16×
242	}	16×
243	if upper {	24×
244	if j > 0 {
245	bi.Daxpy(j, -x[j]*tscal, a[j:], lda, x, 1)
246	i := bi.Idamax(j, x, 1)
247	xmax = math.Abs(x[i])
248	}
249	} else {
250	if j < n-1 {
251	bi.Daxpy(n-j-1, -x[j]tscal, a[(j+1)lda+j:], lda, x[j+1:], 1)
252	i := j + bi.Idamax(n-j-1, x[j+1:], 1)
253	xmax = math.Abs(x[i])
254	}
255	}
256	}
257	} else {	16×
258	for j := jfirst; j != jlast; j += jinc {	32×
259	xj := math.Abs(x[j])	16×
260	uscal := tscal	16×
261	rec := 1 / math.Max(xmax, 1)	16×
262	var tjjs float64	16×
263	if cnorm[j] > (bignum-xj)*rec {	32×
264	rec *= 0.5	16×
265	if nonUnit {	32×
266	tjjs = a[jlda+j] tscal	16×
267	} else {	32×
268	tjjs = tscal	16×
269	}	16×
270	tjj := math.Abs(tjjs)	16×
271	if tjj > 1 {	32×
272	rec = math.Min(1, rec*tjj)	16×
273	uscal /= tjjs	16×
274	}	16×
275	if rec < 1 {	32×
276	bi.Dscal(n, rec, x, 1)	16×
277	scale *= rec	16×
278	xmax *= rec	16×
279	}	16×
280	}
281	var sumj float64	16×
282	if uscal == 1 {	32×
283	if upper {	32×
284	sumj = bi.Ddot(j, a[j:], lda, x, 1)	16×
285	} else if j < n-1 {	48×
286	sumj = bi.Ddot(n-j-1, a[(j+1)*lda+j:], lda, x[j+1:], 1)	16×
287	}	16×
288	} else {	16×
289	if upper {	32×
290	for i := 0; i < j; i++ {	32×
291	sumj += (a[ilda+j] uscal) * x[i]	16×
292	}	16×
293	} else if j < n {	32×
294	for i := j + 1; i < n; i++ {	32×
295	sumj += (a[ilda+j] uscal) * x[i]	16×
296	}	16×
297	}
298	}
299	if uscal == tscal {	32×
300	x[j] -= sumj	16×
301	xj := math.Abs(x[j])	16×
302	var tjjs float64	16×
303	if nonUnit {	32×
304	tjjs = a[jlda+j] tscal	16×
305	} else {	32×
306	tjjs = tscal	16×
307	if tscal == 1 {	32×
308	goto Skip2	16×
309	}
310	}
311	tjj := math.Abs(tjjs)	16×
312	if tjj > smlnum {	32×
313	if tjj < 1 {	32×
314	if xj > tjj*bignum {	32×
315	rec = 1 / xj	16×
316	bi.Dscal(n, rec, x, 1)	16×
317	scale *= rec	16×
318	xmax *= rec	16×
319	}	16×
320	}
321	x[j] /= tjjs	16×
322	} else if tjj > 0 {	32×
323	if xj > tjj*bignum {	32×
324	rec = (tjj * bignum) / xj	16×
325	bi.Dscal(n, rec, x, 1)	16×
326	scale *= rec	16×
327	xmax *= rec	16×
328	}	16×
329	x[j] /= tjjs	16×
330	} else {	16×
331	for i := 0; i < n; i++ {	32×
332	x[i] = 0	16×
333	}	16×
334	x[j] = 1	16×
335	scale = 0	16×
336	xmax = 0	16×
337	}
338	} else {	16×
339	x[j] = x[j]/tjjs - sumj	16×
340	}	16×
341	Skip2:
342	xmax = math.Max(xmax, math.Abs(x[j]))	16×
343	}
344	}
345	scale /= tscal	16×
346	if tscal != 1 {	32×
347	bi.Dscal(n, 1/tscal, cnorm, 1)	16×
348	}	16×
349	return scale	16×
350	}

gonum / gonum / 2245

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous