github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/testlapack/dormqr.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  	"fmt"
     9  	"math/rand"
    10  	"testing"
    11  
    12  	"github.com/gonum/blas"
    13  	"github.com/gonum/floats"
    14  )
    15  
    16  type Dormqrer interface {
    17  	Dorm2rer
    18  	Dormqr(side blas.Side, trans blas.Transpose, m, n, k int, a []float64, lda int, tau, c []float64, ldc int, work []float64, lwork int)
    19  }
    20  
    21  func DormqrTest(t *testing.T, impl Dormqrer) {
    22  	rnd := rand.New(rand.NewSource(1))
    23  	for _, side := range []blas.Side{blas.Left, blas.Right} {
    24  		for _, trans := range []blas.Transpose{blas.NoTrans, blas.Trans} {
    25  			for _, test := range []struct {
    26  				common, adim, cdim, lda, ldc int
    27  			}{
    28  				{6, 7, 8, 0, 0},
    29  				{6, 8, 7, 0, 0},
    30  				{7, 6, 8, 0, 0},
    31  				{7, 8, 6, 0, 0},
    32  				{8, 6, 7, 0, 0},
    33  				{8, 7, 6, 0, 0},
    34  				{100, 200, 300, 0, 0},
    35  				{100, 300, 200, 0, 0},
    36  				{200, 100, 300, 0, 0},
    37  				{200, 300, 100, 0, 0},
    38  				{300, 100, 200, 0, 0},
    39  				{300, 200, 100, 0, 0},
    40  				{100, 200, 300, 400, 500},
    41  				{100, 300, 200, 400, 500},
    42  				{200, 100, 300, 400, 500},
    43  				{200, 300, 100, 400, 500},
    44  				{300, 100, 200, 400, 500},
    45  				{300, 200, 100, 400, 500},
    46  				{100, 200, 300, 500, 400},
    47  				{100, 300, 200, 500, 400},
    48  				{200, 100, 300, 500, 400},
    49  				{200, 300, 100, 500, 400},
    50  				{300, 100, 200, 500, 400},
    51  				{300, 200, 100, 500, 400},
    52  			} {
    53  				var ma, na, mc, nc int
    54  				if side == blas.Left {
    55  					ma = test.common
    56  					na = test.adim
    57  					mc = test.common
    58  					nc = test.cdim
    59  				} else {
    60  					ma = test.common
    61  					na = test.adim
    62  					mc = test.cdim
    63  					nc = test.common
    64  				}
    65  				// Generate a random matrix
    66  				lda := test.lda
    67  				if lda == 0 {
    68  					lda = na
    69  				}
    70  				a := make([]float64, ma*lda)
    71  				for i := range a {
    72  					a[i] = rnd.Float64()
    73  				}
    74  				// Compute random C matrix
    75  				ldc := test.ldc
    76  				if ldc == 0 {
    77  					ldc = nc
    78  				}
    79  				c := make([]float64, mc*ldc)
    80  				for i := range c {
    81  					c[i] = rnd.Float64()
    82  				}
    83  
    84  				// Compute QR
    85  				k := min(ma, na)
    86  				tau := make([]float64, k)
    87  				work := make([]float64, 1)
    88  				impl.Dgeqrf(ma, na, a, lda, tau, work, -1)
    89  				work = make([]float64, int(work[0]))
    90  				impl.Dgeqrf(ma, na, a, lda, tau, work, len(work))
    91  
    92  				cCopy := make([]float64, len(c))
    93  				copy(cCopy, c)
    94  				ans := make([]float64, len(c))
    95  				copy(ans, cCopy)
    96  
    97  				if side == blas.Left {
    98  					work = make([]float64, nc)
    99  				} else {
   100  					work = make([]float64, mc)
   101  				}
   102  				impl.Dorm2r(side, trans, mc, nc, k, a, lda, tau, ans, ldc, work)
   103  
   104  				// Make sure Dorm2r and Dormqr match with small work
   105  				for i := range work {
   106  					work[i] = rnd.Float64()
   107  				}
   108  				copy(c, cCopy)
   109  				impl.Dormqr(side, trans, mc, nc, k, a, lda, tau, c, ldc, work, len(work))
   110  				if !floats.EqualApprox(c, ans, 1e-12) {
   111  					t.Errorf("Dormqr and Dorm2r mismatch for small work")
   112  				}
   113  
   114  				// Try with the optimum amount of work
   115  				copy(c, cCopy)
   116  				impl.Dormqr(side, trans, mc, nc, k, nil, lda, nil, nil, ldc, work, -1)
   117  				work = make([]float64, int(work[0]))
   118  				for i := range work {
   119  					work[i] = rnd.Float64()
   120  				}
   121  				impl.Dormqr(side, trans, mc, nc, k, a, lda, tau, c, ldc, work, len(work))
   122  				if !floats.EqualApprox(c, ans, 1e-12) {
   123  					t.Errorf("Dormqr and Dorm2r mismatch for full work")
   124  					fmt.Println("ccopy")
   125  					for i := 0; i < mc; i++ {
   126  						fmt.Println(cCopy[i*ldc : (i+1)*ldc])
   127  					}
   128  					fmt.Println("ans =")
   129  					for i := 0; i < mc; i++ {
   130  						fmt.Println(ans[i*ldc : (i+1)*ldc])
   131  					}
   132  					fmt.Println("c =")
   133  					for i := 0; i < mc; i++ {
   134  						fmt.Println(c[i*ldc : (i+1)*ldc])
   135  					}
   136  				}
   137  
   138  				// Try with amount of work that is less than
   139  				// optimal but still long enough to use the
   140  				// blocked code.
   141  				copy(c, cCopy)
   142  				if side == blas.Left {
   143  					work = make([]float64, 3*nc)
   144  				} else {
   145  					work = make([]float64, 3*mc)
   146  				}
   147  				impl.Dormqr(side, trans, mc, nc, k, a, lda, tau, c, ldc, work, len(work))
   148  				if !floats.EqualApprox(c, ans, 1e-12) {
   149  					t.Errorf("Dormqr and Dorm2r mismatch for medium work")
   150  				}
   151  			}
   152  		}
   153  	}
   154  }