github.com/gopherd/gonum@v0.0.4/lapack/testlapack/dorghr.go (about)

     1  // Copyright ©2016 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 testlapack
     6  
     7  import (
     8  	"fmt"
     9  	"math"
    10  	"testing"
    11  
    12  	"math/rand"
    13  
    14  	"github.com/gopherd/gonum/blas"
    15  	"github.com/gopherd/gonum/blas/blas64"
    16  )
    17  
    18  type Dorghrer interface {
    19  	Dorghr(n, ilo, ihi int, a []float64, lda int, tau, work []float64, lwork int)
    20  
    21  	Dgehrder
    22  }
    23  
    24  func DorghrTest(t *testing.T, impl Dorghrer) {
    25  	rnd := rand.New(rand.NewSource(1))
    26  
    27  	for _, n := range []int{1, 2, 3, 4, 5, 6, 7, 8, 23, 34} {
    28  		for _, extra := range []int{0, 1, 13} {
    29  			for _, optwork := range []bool{true, false} {
    30  				for cas := 0; cas < 100; cas++ {
    31  					ilo := rnd.Intn(n)
    32  					ihi := rnd.Intn(n)
    33  					if ilo > ihi {
    34  						ilo, ihi = ihi, ilo
    35  					}
    36  					testDorghr(t, impl, n, ilo, ihi, extra, optwork, rnd)
    37  				}
    38  			}
    39  		}
    40  	}
    41  	testDorghr(t, impl, 0, 0, -1, 0, false, rnd)
    42  	testDorghr(t, impl, 0, 0, -1, 0, true, rnd)
    43  }
    44  
    45  func testDorghr(t *testing.T, impl Dorghrer, n, ilo, ihi, extra int, optwork bool, rnd *rand.Rand) {
    46  	const tol = 1e-14
    47  
    48  	// Construct the matrix A with elementary reflectors and scalar factors tau.
    49  	a := randomGeneral(n, n, n+extra, rnd)
    50  	var tau []float64
    51  	if n > 1 {
    52  		tau = nanSlice(n - 1)
    53  	}
    54  	work := nanSlice(max(1, n)) // Minimum work for Dgehrd.
    55  	impl.Dgehrd(n, ilo, ihi, a.Data, a.Stride, tau, work, len(work))
    56  
    57  	// Extract Q for later comparison.
    58  	q := eye(n, n)
    59  	qCopy := cloneGeneral(q)
    60  	for j := ilo; j < ihi; j++ {
    61  		h := eye(n, n)
    62  		v := blas64.Vector{
    63  			Inc:  1,
    64  			Data: make([]float64, n),
    65  		}
    66  		v.Data[j+1] = 1
    67  		for i := j + 2; i < ihi+1; i++ {
    68  			v.Data[i] = a.Data[i*a.Stride+j]
    69  		}
    70  		blas64.Ger(-tau[j], v, v, h)
    71  		copy(qCopy.Data, q.Data)
    72  		blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, qCopy, h, 0, q)
    73  	}
    74  
    75  	if optwork {
    76  		work = nanSlice(1)
    77  		impl.Dorghr(n, ilo, ihi, a.Data, a.Stride, tau, work, -1)
    78  		work = nanSlice(int(work[0]))
    79  	} else {
    80  		work = nanSlice(max(1, ihi-ilo))
    81  	}
    82  	impl.Dorghr(n, ilo, ihi, a.Data, a.Stride, tau, work, len(work))
    83  
    84  	prefix := fmt.Sprintf("Case n=%v, ilo=%v, ihi=%v, extra=%v, optwork=%v", n, ilo, ihi, extra, optwork)
    85  	if !generalOutsideAllNaN(a) {
    86  		t.Errorf("%v: out-of-range write to A\n%v", prefix, a.Data)
    87  	}
    88  	if resid := residualOrthogonal(a, false); resid > tol {
    89  		t.Errorf("%v: A is not orthogonal; resid=%v, want<=%v", prefix, resid, tol)
    90  	}
    91  	for i := 0; i < n; i++ {
    92  		for j := 0; j < n; j++ {
    93  			aij := a.Data[i*a.Stride+j]
    94  			qij := q.Data[i*q.Stride+j]
    95  			if math.Abs(aij-qij) > tol {
    96  				t.Errorf("%v: unexpected value of A[%v,%v]. want %v, got %v", prefix, i, j, qij, aij)
    97  			}
    98  		}
    99  	}
   100  }