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

     1  // Copyright ©2017 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  	"golang.org/x/exp/rand"
    12  
    13  	"github.com/jingcheng-WU/gonum/blas"
    14  	"github.com/jingcheng-WU/gonum/blas/blas64"
    15  	"github.com/jingcheng-WU/gonum/lapack"
    16  )
    17  
    18  type Dlaqpser interface {
    19  	Dlapmter
    20  	Dlaqps(m, n, offset, nb int, a []float64, lda int, jpvt []int, tau, vn1, vn2, auxv, f []float64, ldf int) (kb int)
    21  }
    22  
    23  func DlaqpsTest(t *testing.T, impl Dlaqpser) {
    24  	const tol = 1e-14
    25  
    26  	rnd := rand.New(rand.NewSource(1))
    27  	for ti, test := range []struct {
    28  		m, n, nb, offset int
    29  	}{
    30  		{m: 4, n: 3, nb: 2, offset: 0},
    31  		{m: 4, n: 3, nb: 1, offset: 2},
    32  		{m: 3, n: 4, nb: 2, offset: 0},
    33  		{m: 3, n: 4, nb: 1, offset: 2},
    34  		{m: 8, n: 3, nb: 2, offset: 0},
    35  		{m: 8, n: 3, nb: 1, offset: 4},
    36  		{m: 3, n: 8, nb: 2, offset: 0},
    37  		{m: 3, n: 8, nb: 1, offset: 1},
    38  		{m: 10, n: 10, nb: 3, offset: 0},
    39  		{m: 10, n: 10, nb: 2, offset: 5},
    40  	} {
    41  		m := test.m
    42  		n := test.n
    43  		jpiv := make([]int, n)
    44  
    45  		for _, extra := range []int{0, 11} {
    46  			a := randomGeneral(m, n, n+extra, rnd)
    47  			aCopy := cloneGeneral(a)
    48  
    49  			for j := range jpiv {
    50  				jpiv[j] = j
    51  			}
    52  
    53  			tau := make([]float64, n)
    54  			vn1 := columnNorms(m, n, a.Data, a.Stride)
    55  			vn2 := columnNorms(m, n, a.Data, a.Stride)
    56  			auxv := make([]float64, test.nb)
    57  			f := zeros(test.n, test.nb, n)
    58  
    59  			kb := impl.Dlaqps(m, n, test.offset, test.nb, a.Data, a.Stride, jpiv, tau, vn1, vn2, auxv, f.Data, f.Stride)
    60  
    61  			prefix := fmt.Sprintf("Case %v (m=%v,n=%v,offset=%d,nb=%v,extra=%v)", ti, m, n, test.offset, test.nb, extra)
    62  
    63  			if !generalOutsideAllNaN(a) {
    64  				t.Errorf("%v: out-of-range write to A", prefix)
    65  			}
    66  			if kb != test.nb {
    67  				t.Logf("%v: %v columns out of %v factorized", prefix, kb, test.nb)
    68  			}
    69  
    70  			mo := m - test.offset
    71  			if mo == 0 {
    72  				continue
    73  			}
    74  
    75  			q := constructQ("QR", mo, kb, a.Data[test.offset*a.Stride:], a.Stride, tau)
    76  
    77  			// Check that Q is orthogonal.
    78  			if resid := residualOrthogonal(q, false); resid > tol {
    79  				t.Errorf("%v: Q is not orthogonal; resid=%v, want<=%v", prefix, resid, tol)
    80  			}
    81  
    82  			// Check that |A*P - Q*R| is small.
    83  			impl.Dlapmt(true, aCopy.Rows, aCopy.Cols, aCopy.Data, aCopy.Stride, jpiv)
    84  			qrap := blas64.General{
    85  				Rows:   mo,
    86  				Cols:   kb,
    87  				Stride: aCopy.Stride,
    88  				Data:   aCopy.Data[test.offset*aCopy.Stride:],
    89  			}
    90  			r := zeros(mo, kb, n)
    91  			for i := 0; i < mo; i++ {
    92  				for j := i; j < kb; j++ {
    93  					r.Data[i*r.Stride+j] = a.Data[(test.offset+i)*a.Stride+j]
    94  				}
    95  			}
    96  			blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, q, r, -1, qrap)
    97  			resid := dlange(lapack.MaxColumnSum, qrap.Rows, qrap.Cols, qrap.Data, qrap.Stride)
    98  			if resid > tol {
    99  				t.Errorf("%v: |Q*R - A*P|=%v, want<=%v", prefix, resid, tol)
   100  			}
   101  		}
   102  	}
   103  }