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

     1  // Copyright ©2021 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  	"github.com/gopherd/gonum/floats"
    17  	"github.com/gopherd/gonum/lapack"
    18  )
    19  
    20  type Dorgr2er interface {
    21  	Dorgr2(m, n, k int, a []float64, lda int, tau []float64, work []float64)
    22  
    23  	Dgerqfer
    24  }
    25  
    26  func Dorgr2Test(t *testing.T, impl Dorgr2er) {
    27  	rnd := rand.New(rand.NewSource(1))
    28  	for _, k := range []int{0, 1, 2, 5} {
    29  		for _, m := range []int{k, k + 1, k + 2, k + 4} {
    30  			for _, n := range []int{m, m + 1, m + 2, m + 4, m + 7} {
    31  				for _, lda := range []int{max(1, n), n + 5} {
    32  					dorgr2Test(t, impl, rnd, m, n, k, lda)
    33  				}
    34  			}
    35  		}
    36  	}
    37  }
    38  
    39  func dorgr2Test(t *testing.T, impl Dorgr2er, rnd *rand.Rand, m, n, k, lda int) {
    40  	const tol = 1e-14
    41  
    42  	name := fmt.Sprintf("m=%v,n=%v,k=%v,lda=%v", m, n, k, lda)
    43  
    44  	// Generate a random m×n matrix A.
    45  	a := randomGeneral(m, n, lda, rnd)
    46  
    47  	// Compute the RQ decomposition of A.
    48  	rq := cloneGeneral(a)
    49  	tau := make([]float64, m)
    50  	work := make([]float64, 1)
    51  	impl.Dgerqf(m, n, rq.Data, rq.Stride, tau, work, -1)
    52  	work = make([]float64, int(work[0]))
    53  	impl.Dgerqf(m, n, rq.Data, rq.Stride, tau, work, len(work))
    54  
    55  	tauCopy := make([]float64, len(tau))
    56  	copy(tauCopy, tau)
    57  
    58  	// Compute the matrix Q using Dorg2r.
    59  	q := cloneGeneral(rq)
    60  	impl.Dorgr2(m, n, k, q.Data, q.Stride, tau[m-k:m], work)
    61  
    62  	if m == 0 {
    63  		return
    64  	}
    65  
    66  	// Check that tau hasn't been modified.
    67  	if !floats.Equal(tau, tauCopy) {
    68  		t.Errorf("%v: unexpected modification in tau", name)
    69  	}
    70  
    71  	// Check that Q has orthonormal rows.
    72  	res := residualOrthogonal(q, true)
    73  	if res > tol || math.IsNaN(res) {
    74  		t.Errorf("%v: residual |I - Q*Qᵀ| too large, got %v, want <= %v", name, res, tol)
    75  	}
    76  
    77  	if k == 0 {
    78  		return
    79  	}
    80  
    81  	// Extract the k×m upper triangular matrix R from RQ[m-k:m,n-k:n].
    82  	r := zeros(k, m, m)
    83  	for i := 0; i < k; i++ {
    84  		for j := 0; j < k; j++ {
    85  			ii := rq.Rows - k + i
    86  			jj := rq.Cols - k + j
    87  			jr := r.Cols - k + j
    88  			if i <= j {
    89  				r.Data[i*r.Stride+jr] = rq.Data[ii*rq.Stride+jj]
    90  			}
    91  		}
    92  	}
    93  
    94  	// Construct a view A[m-k:m,0:n] of the last k rows of A.
    95  	aRec := blas64.General{
    96  		Rows:   k,
    97  		Cols:   n,
    98  		Data:   a.Data[(m-k)*a.Stride:],
    99  		Stride: a.Stride,
   100  	}
   101  	// Compute A - R*Q.
   102  	blas64.Gemm(blas.NoTrans, blas.NoTrans, -1, r, q, 1, aRec)
   103  	// Check that |A - R*Q| is small.
   104  	res = dlange(lapack.MaxColumnSum, aRec.Rows, aRec.Cols, aRec.Data, aRec.Stride)
   105  	if res > tol || math.IsNaN(res) {
   106  		t.Errorf("%v: residual |A - R*Q| too large, got %v, want <= %v", name, res, tol)
   107  	}
   108  }