github.com/gopherd/gonum@v0.0.4/lapack/testlapack/dorgql.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  	"testing"
    10  
    11  	"math/rand"
    12  
    13  	"github.com/gopherd/gonum/blas"
    14  	"github.com/gopherd/gonum/blas/blas64"
    15  )
    16  
    17  type Dorgqler interface {
    18  	Dorgql(m, n, k int, a []float64, lda int, tau, work []float64, lwork int)
    19  
    20  	Dlarfger
    21  }
    22  
    23  func DorgqlTest(t *testing.T, impl Dorgqler) {
    24  	const tol = 1e-14
    25  
    26  	type Dorg2ler interface {
    27  		Dorg2l(m, n, k int, a []float64, lda int, tau, work []float64)
    28  	}
    29  	dorg2ler, hasDorg2l := impl.(Dorg2ler)
    30  
    31  	rnd := rand.New(rand.NewSource(1))
    32  	for _, m := range []int{0, 1, 2, 3, 4, 5, 7, 10, 15, 30, 50, 150} {
    33  		for _, extra := range []int{0, 11} {
    34  			for _, wl := range []worklen{minimumWork, mediumWork, optimumWork} {
    35  				var k int
    36  				if m >= 129 {
    37  					// For large matrices make sure that k
    38  					// is large enough to trigger blocked
    39  					// path.
    40  					k = 129 + rnd.Intn(m-129+1)
    41  				} else {
    42  					k = rnd.Intn(m + 1)
    43  				}
    44  				n := k + rnd.Intn(m-k+1)
    45  				if m == 0 || n == 0 {
    46  					m = 0
    47  					n = 0
    48  					k = 0
    49  				}
    50  
    51  				// Generate k elementary reflectors in the last
    52  				// k columns of A.
    53  				a := nanGeneral(m, n, n+extra)
    54  				tau := make([]float64, k)
    55  				for l := 0; l < k; l++ {
    56  					jj := m - k + l
    57  					v := randomSlice(jj, rnd)
    58  					_, tau[l] = impl.Dlarfg(len(v)+1, rnd.NormFloat64(), v, 1)
    59  					j := n - k + l
    60  					for i := 0; i < jj; i++ {
    61  						a.Data[i*a.Stride+j] = v[i]
    62  					}
    63  				}
    64  				aCopy := cloneGeneral(a)
    65  
    66  				// Compute the full matrix Q by forming the
    67  				// Householder reflectors explicitly.
    68  				q := eye(m, m)
    69  				qCopy := eye(m, m)
    70  				for l := 0; l < k; l++ {
    71  					h := eye(m, m)
    72  					jj := m - k + l
    73  					j := n - k + l
    74  					v := blas64.Vector{Data: make([]float64, m), Inc: 1}
    75  					for i := 0; i < jj; i++ {
    76  						v.Data[i] = a.Data[i*a.Stride+j]
    77  					}
    78  					v.Data[jj] = 1
    79  					blas64.Ger(-tau[l], v, v, h)
    80  					copy(qCopy.Data, q.Data)
    81  					blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, h, qCopy, 0, q)
    82  				}
    83  				// View the last n columns of Q as 'want'.
    84  				want := blas64.General{
    85  					Rows:   m,
    86  					Cols:   n,
    87  					Stride: q.Stride,
    88  					Data:   q.Data[m-n:],
    89  				}
    90  
    91  				var lwork int
    92  				switch wl {
    93  				case minimumWork:
    94  					lwork = max(1, n)
    95  				case mediumWork:
    96  					work := make([]float64, 1)
    97  					impl.Dorgql(m, n, k, a.Data, a.Stride, tau, work, -1)
    98  					lwork = (int(work[0]) + n) / 2
    99  					lwork = max(1, lwork)
   100  				case optimumWork:
   101  					work := make([]float64, 1)
   102  					impl.Dorgql(m, n, k, a.Data, a.Stride, tau, work, -1)
   103  					lwork = int(work[0])
   104  				}
   105  				work := make([]float64, lwork)
   106  
   107  				// Compute the last n columns of Q by a call to
   108  				// Dorgql.
   109  				impl.Dorgql(m, n, k, a.Data, a.Stride, tau, work, len(work))
   110  
   111  				prefix := fmt.Sprintf("Case m=%v,n=%v,k=%v,wl=%v", m, n, k, wl)
   112  				if !generalOutsideAllNaN(a) {
   113  					t.Errorf("%v: out-of-range write to A", prefix)
   114  				}
   115  				if !equalApproxGeneral(want, a, tol) {
   116  					t.Errorf("%v: unexpected Q", prefix)
   117  				}
   118  
   119  				// Compute the last n columns of Q by a call to
   120  				// Dorg2l and check that we get the same result.
   121  				if !hasDorg2l {
   122  					continue
   123  				}
   124  				dorg2ler.Dorg2l(m, n, k, aCopy.Data, aCopy.Stride, tau, work)
   125  				if !equalApproxGeneral(aCopy, a, tol) {
   126  					t.Errorf("%v: mismatch between Dorgql and Dorg2l", prefix)
   127  				}
   128  			}
   129  		}
   130  	}
   131  }