github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/lapack/testlapack/dlasy2.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  	"math"
    10  	"testing"
    11  
    12  	"golang.org/x/exp/rand"
    13  
    14  	"github.com/jingcheng-WU/gonum/blas"
    15  	"github.com/jingcheng-WU/gonum/blas/blas64"
    16  )
    17  
    18  type Dlasy2er interface {
    19  	Dlasy2(tranl, tranr bool, isgn, n1, n2 int, tl []float64, ldtl int, tr []float64, ldtr int, b []float64, ldb int, x []float64, ldx int) (scale, xnorm float64, ok bool)
    20  }
    21  
    22  func Dlasy2Test(t *testing.T, impl Dlasy2er) {
    23  	rnd := rand.New(rand.NewSource(1))
    24  	for _, tranl := range []bool{true, false} {
    25  		for _, tranr := range []bool{true, false} {
    26  			for _, isgn := range []int{1, -1} {
    27  				for _, n1 := range []int{0, 1, 2} {
    28  					for _, n2 := range []int{0, 1, 2} {
    29  						for _, extra := range []int{0, 1, 2, 13} {
    30  							for cas := 0; cas < 1000; cas++ {
    31  								testDlasy2(t, impl, tranl, tranr, isgn, n1, n2, extra, rnd)
    32  							}
    33  						}
    34  					}
    35  				}
    36  			}
    37  		}
    38  	}
    39  }
    40  
    41  func testDlasy2(t *testing.T, impl Dlasy2er, tranl, tranr bool, isgn, n1, n2, extra int, rnd *rand.Rand) {
    42  	const tol = 1e-10
    43  
    44  	tl := randomGeneral(n1, n1, n1+extra, rnd)
    45  	tr := randomGeneral(n2, n2, n2+extra, rnd)
    46  	b := randomGeneral(n1, n2, n2+extra, rnd)
    47  	x := randomGeneral(n1, n2, n2+extra, rnd)
    48  
    49  	scale, xnorm, ok := impl.Dlasy2(tranl, tranr, isgn, n1, n2, tl.Data, tl.Stride, tr.Data, tr.Stride, b.Data, b.Stride, x.Data, x.Stride)
    50  	if scale > 1 {
    51  		t.Errorf("invalid value of scale, want <= 1, got %v", scale)
    52  	}
    53  	if n1 == 0 || n2 == 0 {
    54  		return
    55  	}
    56  
    57  	prefix := fmt.Sprintf("Case n1=%v, n2=%v, isgn=%v", n1, n2, isgn)
    58  
    59  	// Check any invalid modifications of x.
    60  	if !generalOutsideAllNaN(x) {
    61  		t.Errorf("%v: out-of-range write to x\n%v", prefix, x.Data)
    62  	}
    63  
    64  	var xnormWant float64
    65  	for i := 0; i < n1; i++ {
    66  		var rowsum float64
    67  		for j := 0; j < n2; j++ {
    68  			rowsum += math.Abs(x.Data[i*x.Stride+j])
    69  		}
    70  		if rowsum > xnormWant {
    71  			xnormWant = rowsum
    72  		}
    73  	}
    74  	if xnormWant != xnorm {
    75  		t.Errorf("%v: unexpected xnorm: want %v, got %v", prefix, xnormWant, xnorm)
    76  	}
    77  
    78  	// Multiply b by scale to get the wanted right-hand side.
    79  	for i := 0; i < n1; i++ {
    80  		for j := 0; j < n2; j++ {
    81  			b.Data[i*b.Stride+j] *= scale
    82  		}
    83  	}
    84  	// Compute the wanted left-hand side.
    85  	lhsWant := randomGeneral(n1, n2, n2, rnd)
    86  	if tranl {
    87  		blas64.Gemm(blas.Trans, blas.NoTrans, 1, tl, x, 0, lhsWant)
    88  	} else {
    89  		blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, tl, x, 0, lhsWant)
    90  	}
    91  	if tranr {
    92  		blas64.Gemm(blas.NoTrans, blas.Trans, float64(isgn), x, tr, 1, lhsWant)
    93  	} else {
    94  		blas64.Gemm(blas.NoTrans, blas.NoTrans, float64(isgn), x, tr, 1, lhsWant)
    95  	}
    96  	// Compare them.
    97  	for i := 0; i < n1; i++ {
    98  		for j := 0; j < n2; j++ {
    99  			diff := lhsWant.Data[i*lhsWant.Stride+j] - b.Data[i*b.Stride+j]
   100  			if math.Abs(diff) > tol && ok {
   101  				t.Errorf("%v: unexpected result, diff[%v,%v]=%v", prefix, i, j, diff)
   102  			}
   103  		}
   104  	}
   105  }