| // 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" |
| "testing" |
| |
| "golang.org/x/exp/rand" |
| |
| "gonum.org/v1/gonum/blas" |
| "gonum.org/v1/gonum/blas/blas64" |
| "gonum.org/v1/gonum/lapack" |
| ) |
| |
| type Dlaexcer interface { |
| Dlaexc(wantq bool, n int, t []float64, ldt int, q []float64, ldq int, j1, n1, n2 int, work []float64) bool |
| } |
| |
| func DlaexcTest(t *testing.T, impl Dlaexcer) { |
| rnd := rand.New(rand.NewSource(1)) |
| |
| for _, n := range []int{1, 2, 3, 4, 5, 6, 10, 18, 31, 53} { |
| for _, extra := range []int{0, 3} { |
| for cas := 0; cas < 100; cas++ { |
| testDlaexc(t, impl, rnd, n, extra) |
| } |
| } |
| } |
| } |
| |
| func testDlaexc(t *testing.T, impl Dlaexcer, rnd *rand.Rand, n, extra int) { |
| const tol = 1e-14 |
| |
| // Generate random T in Schur canonical form. |
| tmat, _, _ := randomSchurCanonical(n, n+extra, true, rnd) |
| tmatCopy := cloneGeneral(tmat) |
| |
| // Randomly pick the index of the first block. |
| j1 := rnd.Intn(n) |
| if j1 > 0 && tmat.Data[j1*tmat.Stride+j1-1] != 0 { |
| // Adjust j1 if it points to the second row of a 2x2 block. |
| j1-- |
| } |
| // Read sizes of the two blocks based on properties of T. |
| var n1, n2 int |
| switch j1 { |
| case n - 1: |
| n1, n2 = 1, 0 |
| case n - 2: |
| if tmat.Data[(j1+1)*tmat.Stride+j1] == 0 { |
| n1, n2 = 1, 1 |
| } else { |
| n1, n2 = 2, 0 |
| } |
| case n - 3: |
| if tmat.Data[(j1+1)*tmat.Stride+j1] == 0 { |
| n1, n2 = 1, 2 |
| } else { |
| n1, n2 = 2, 1 |
| } |
| default: |
| if tmat.Data[(j1+1)*tmat.Stride+j1] == 0 { |
| n1 = 1 |
| if tmat.Data[(j1+2)*tmat.Stride+j1+1] == 0 { |
| n2 = 1 |
| } else { |
| n2 = 2 |
| } |
| } else { |
| n1 = 2 |
| if tmat.Data[(j1+3)*tmat.Stride+j1+2] == 0 { |
| n2 = 1 |
| } else { |
| n2 = 2 |
| } |
| } |
| } |
| |
| name := fmt.Sprintf("Case n=%v,j1=%v,n1=%v,n2=%v,extra=%v", n, j1, n1, n2, extra) |
| |
| // 1. Test without accumulating Q. |
| |
| wantq := false |
| |
| work := nanSlice(n) |
| |
| ok := impl.Dlaexc(wantq, n, tmat.Data, tmat.Stride, nil, 1, j1, n1, n2, work) |
| |
| // 2. Test with accumulating Q. |
| |
| wantq = true |
| |
| tmat2 := cloneGeneral(tmatCopy) |
| q := eye(n, n+extra) |
| qCopy := cloneGeneral(q) |
| work = nanSlice(n) |
| |
| ok2 := impl.Dlaexc(wantq, n, tmat2.Data, tmat2.Stride, q.Data, q.Stride, j1, n1, n2, work) |
| |
| if !generalOutsideAllNaN(tmat) { |
| t.Errorf("%v: out-of-range write to T", name) |
| } |
| if !generalOutsideAllNaN(tmat2) { |
| t.Errorf("%v: out-of-range write to T2", name) |
| } |
| if !generalOutsideAllNaN(q) { |
| t.Errorf("%v: out-of-range write to Q", name) |
| } |
| |
| // Check that outputs from cases 1. and 2. are exactly equal, then check one of them. |
| if ok != ok2 { |
| t.Errorf("%v: ok != ok2", name) |
| } |
| if !equalGeneral(tmat, tmat2) { |
| t.Errorf("%v: T != T2", name) |
| } |
| |
| if !ok { |
| if n1 == 1 && n2 == 1 { |
| t.Errorf("%v: unexpected failure", name) |
| } else { |
| t.Logf("%v: Dlaexc returned false", name) |
| } |
| } |
| |
| if !ok || n1 == 0 || n2 == 0 || j1+n1 >= n { |
| // Check that T is not modified. |
| if !equalGeneral(tmat, tmatCopy) { |
| t.Errorf("%v: unexpected modification of T", name) |
| } |
| // Check that Q is not modified. |
| if !equalGeneral(q, qCopy) { |
| t.Errorf("%v: unexpected modification of Q", name) |
| } |
| return |
| } |
| |
| // Check that T is not modified outside of rows and columns [j1:j1+n1+n2]. |
| for i := 0; i < n; i++ { |
| if j1 <= i && i < j1+n1+n2 { |
| continue |
| } |
| for j := 0; j < n; j++ { |
| if j1 <= j && j < j1+n1+n2 { |
| continue |
| } |
| diff := tmat.Data[i*tmat.Stride+j] - tmatCopy.Data[i*tmatCopy.Stride+j] |
| if diff != 0 { |
| t.Errorf("%v: unexpected modification of T[%v,%v]", name, i, j) |
| } |
| } |
| } |
| |
| if !isSchurCanonicalGeneral(tmat) { |
| t.Errorf("%v: T is not in Schur canonical form", name) |
| } |
| |
| // Check that Q is orthogonal. |
| resid := residualOrthogonal(q, false) |
| if resid > tol { |
| t.Errorf("%v: Q is not orthogonal; resid=%v, want<=%v", name, resid, tol) |
| } |
| |
| // Check that Q is unchanged outside of columns [j1:j1+n1+n2]. |
| for i := 0; i < n; i++ { |
| for j := 0; j < n; j++ { |
| if j1 <= j && j < j1+n1+n2 { |
| continue |
| } |
| diff := q.Data[i*q.Stride+j] - qCopy.Data[i*qCopy.Stride+j] |
| if diff != 0 { |
| t.Errorf("%v: unexpected modification of Q[%v,%v]", name, i, j) |
| } |
| } |
| } |
| |
| // Check that Qᵀ * TOrig * Q == T |
| qt := zeros(n, n, n) |
| blas64.Gemm(blas.Trans, blas.NoTrans, 1, q, tmatCopy, 0, qt) |
| qtq := cloneGeneral(tmat) |
| blas64.Gemm(blas.NoTrans, blas.NoTrans, -1, qt, q, 1, qtq) |
| resid = dlange(lapack.MaxColumnSum, n, n, qtq.Data, qtq.Stride) |
| if resid > float64(n)*tol { |
| t.Errorf("%v: mismatch between Qᵀ*(initial T)*Q and (final T); resid=%v, want<=%v", |
| name, resid, float64(n)*tol) |
| } |
| } |