github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/lapack/testlapack/dgelq2.go (about)

     1  // Copyright ©2015 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  	"testing"
     9  
    10  	"golang.org/x/exp/rand"
    11  
    12  	"github.com/jingcheng-WU/gonum/blas"
    13  	"github.com/jingcheng-WU/gonum/blas/blas64"
    14  	"github.com/jingcheng-WU/gonum/floats"
    15  )
    16  
    17  type Dgelq2er interface {
    18  	Dgelq2(m, n int, a []float64, lda int, tau, work []float64)
    19  }
    20  
    21  func Dgelq2Test(t *testing.T, impl Dgelq2er) {
    22  	const tol = 1e-14
    23  
    24  	rnd := rand.New(rand.NewSource(1))
    25  	for c, test := range []struct {
    26  		m, n, lda int
    27  	}{
    28  		{1, 1, 0},
    29  		{2, 2, 0},
    30  		{3, 2, 0},
    31  		{2, 3, 0},
    32  		{1, 12, 0},
    33  		{2, 6, 0},
    34  		{3, 4, 0},
    35  		{4, 3, 0},
    36  		{6, 2, 0},
    37  		{1, 12, 0},
    38  		{1, 1, 20},
    39  		{2, 2, 20},
    40  		{3, 2, 20},
    41  		{2, 3, 20},
    42  		{1, 12, 20},
    43  		{2, 6, 20},
    44  		{3, 4, 20},
    45  		{4, 3, 20},
    46  		{6, 2, 20},
    47  		{1, 12, 20},
    48  	} {
    49  		n := test.n
    50  		m := test.m
    51  		lda := test.lda
    52  		if lda == 0 {
    53  			lda = test.n
    54  		}
    55  		k := min(m, n)
    56  		tau := make([]float64, k)
    57  		for i := range tau {
    58  			tau[i] = rnd.Float64()
    59  		}
    60  		work := make([]float64, m)
    61  		for i := range work {
    62  			work[i] = rnd.Float64()
    63  		}
    64  		a := make([]float64, m*lda)
    65  		for i := 0; i < m*lda; i++ {
    66  			a[i] = rnd.Float64()
    67  		}
    68  		aCopy := make([]float64, len(a))
    69  		copy(aCopy, a)
    70  		impl.Dgelq2(m, n, a, lda, tau, work)
    71  
    72  		Q := constructQ("LQ", m, n, a, lda, tau)
    73  
    74  		// Check that Q is orthogonal.
    75  		if resid := residualOrthogonal(Q, false); resid > tol {
    76  			t.Errorf("Case %v: Q not orthogonal; resid=%v, want<=%v", c, resid, tol)
    77  		}
    78  
    79  		L := blas64.General{
    80  			Rows:   m,
    81  			Cols:   n,
    82  			Stride: n,
    83  			Data:   make([]float64, m*n),
    84  		}
    85  		for i := 0; i < m; i++ {
    86  			for j := 0; j <= min(i, n-1); j++ {
    87  				L.Data[i*L.Stride+j] = a[i*lda+j]
    88  			}
    89  		}
    90  
    91  		ans := blas64.General{
    92  			Rows:   m,
    93  			Cols:   n,
    94  			Stride: lda,
    95  			Data:   make([]float64, m*lda),
    96  		}
    97  		copy(ans.Data, aCopy)
    98  		blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, L, Q, 0, ans)
    99  		if !floats.EqualApprox(aCopy, ans.Data, tol) {
   100  			t.Errorf("Case %v, LQ mismatch. Want %v, got %v.", c, aCopy, ans.Data)
   101  		}
   102  	}
   103  }