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