blob: 002be5e7028da10e7ab4f0d8602aed0495b729ca [file] [log] [blame]
 // 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) } }