github.com/gopherd/gonum@v0.0.4/lapack/testlapack/dgesc2.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  package testlapack
     5  
     6  import (
     7  	"fmt"
     8  	"math"
     9  	"testing"
    10  
    11  	"math/rand"
    12  
    13  	"github.com/gopherd/gonum/blas"
    14  	"github.com/gopherd/gonum/blas/blas64"
    15  	"github.com/gopherd/gonum/floats"
    16  	"github.com/gopherd/gonum/lapack"
    17  )
    18  
    19  type Dgesc2er interface {
    20  	Dgesc2(n int, a []float64, lda int, rhs []float64, ipiv, jpiv []int) (scale float64)
    21  
    22  	Dgetc2er
    23  }
    24  
    25  func Dgesc2Test(t *testing.T, impl Dgesc2er) {
    26  	rnd := rand.New(rand.NewSource(1))
    27  	for _, n := range []int{0, 1, 2, 3, 4, 5, 10, 20, 50} {
    28  		for _, lda := range []int{n, n + 3} {
    29  			testDgesc2(t, impl, rnd, n, lda, false)
    30  			testDgesc2(t, impl, rnd, n, lda, true)
    31  		}
    32  	}
    33  }
    34  
    35  func testDgesc2(t *testing.T, impl Dgesc2er, rnd *rand.Rand, n, lda int, big bool) {
    36  	const tol = 1e-14
    37  
    38  	name := fmt.Sprintf("n=%v,lda=%v,big=%v", n, lda, big)
    39  
    40  	// Generate random general matrix.
    41  	a := randomGeneral(n, n, max(1, lda), rnd)
    42  
    43  	// Generate a random right hand side vector.
    44  	b := randomGeneral(n, 1, 1, rnd)
    45  	if big {
    46  		for i := 0; i < n; i++ {
    47  			b.Data[i] *= bignum
    48  		}
    49  	}
    50  
    51  	// Compute the LU factorization of A with full pivoting.
    52  	lu := cloneGeneral(a)
    53  	ipiv := make([]int, n)
    54  	jpiv := make([]int, n)
    55  	impl.Dgetc2(n, lu.Data, lu.Stride, ipiv, jpiv)
    56  
    57  	// Make copies of const input to Dgesc2.
    58  	luCopy := cloneGeneral(lu)
    59  	ipivCopy := make([]int, len(ipiv))
    60  	copy(ipivCopy, ipiv)
    61  	jpivCopy := make([]int, len(jpiv))
    62  	copy(jpivCopy, jpiv)
    63  
    64  	// Call Dgesc2 to solve A*x = scale*b.
    65  	x := cloneGeneral(b)
    66  	scale := impl.Dgesc2(n, lu.Data, lu.Stride, x.Data, ipiv, jpiv)
    67  
    68  	if n == 0 {
    69  		return
    70  	}
    71  
    72  	// Check that const input to Dgesc2 hasn't been modified.
    73  	if !floats.Same(lu.Data, luCopy.Data) {
    74  		t.Errorf("%v: unexpected modification in lu", name)
    75  	}
    76  	if !intsEqual(ipiv, ipivCopy) {
    77  		t.Errorf("%v: unexpected modification in ipiv", name)
    78  	}
    79  	if !intsEqual(jpiv, jpivCopy) {
    80  		t.Errorf("%v: unexpected modification in jpiv", name)
    81  	}
    82  
    83  	if scale <= 0 || 1 < scale {
    84  		t.Errorf("%v: scale %v out of bounds (0,1]", name, scale)
    85  	}
    86  	if !big && scale != 1 {
    87  		t.Errorf("%v: unexpected scaling, scale=%v", name, scale)
    88  	}
    89  
    90  	// Compute the difference rhs := A*x - scale*b.
    91  	diff := b
    92  	for i := 0; i < n; i++ {
    93  		diff.Data[i] *= scale
    94  	}
    95  	blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, a, x, -1, diff)
    96  
    97  	// Compute the residual |A*x - scale*b| / |x|.
    98  	xnorm := dlange(lapack.MaxColumnSum, n, 1, x.Data, 1)
    99  	resid := dlange(lapack.MaxColumnSum, n, 1, diff.Data, 1) / xnorm
   100  	if resid > tol || math.IsNaN(resid) {
   101  		t.Errorf("%v: unexpected result; resid=%v, want<=%v", name, resid, tol)
   102  	}
   103  }