github.com/gopherd/gonum@v0.0.4/lapack/testlapack/dlantr.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"
    10  	"testing"
    11  
    12  	"math/rand"
    13  
    14  	"github.com/gopherd/gonum/blas"
    15  	"github.com/gopherd/gonum/floats"
    16  	"github.com/gopherd/gonum/lapack"
    17  )
    18  
    19  type Dlantrer interface {
    20  	Dlanger
    21  	Dlantr(norm lapack.MatrixNorm, uplo blas.Uplo, diag blas.Diag, m, n int, a []float64, lda int, work []float64) float64
    22  }
    23  
    24  func DlantrTest(t *testing.T, impl Dlantrer) {
    25  	rnd := rand.New(rand.NewSource(1))
    26  	for _, m := range []int{0, 1, 2, 3, 4, 5, 10} {
    27  		for _, n := range []int{0, 1, 2, 3, 4, 5, 10} {
    28  			for _, uplo := range []blas.Uplo{blas.Lower, blas.Upper} {
    29  				if uplo == blas.Upper && m > n {
    30  					continue
    31  				}
    32  				if uplo == blas.Lower && n > m {
    33  					continue
    34  				}
    35  				for _, diag := range []blas.Diag{blas.NonUnit, blas.Unit} {
    36  					for _, lda := range []int{max(1, n), n + 3} {
    37  						dlantrTest(t, impl, rnd, uplo, diag, m, n, lda)
    38  					}
    39  				}
    40  			}
    41  		}
    42  	}
    43  }
    44  
    45  func dlantrTest(t *testing.T, impl Dlantrer, rnd *rand.Rand, uplo blas.Uplo, diag blas.Diag, m, n, lda int) {
    46  	const tol = 1e-14
    47  
    48  	// Generate a random triangular matrix. If the matrix has unit diagonal,
    49  	// don't set the diagonal elements to 1.
    50  	a := make([]float64, max(0, (m-1)*lda+n))
    51  	for i := range a {
    52  		a[i] = rnd.NormFloat64()
    53  	}
    54  	rowsum := make([]float64, m)
    55  	colsum := make([]float64, n)
    56  	var frobWant, maxabsWant float64
    57  	if diag == blas.Unit {
    58  		// Account for the unit diagonal.
    59  		for i := 0; i < min(m, n); i++ {
    60  			rowsum[i] = 1
    61  			colsum[i] = 1
    62  		}
    63  		frobWant = float64(min(m, n))
    64  		if min(m, n) > 0 {
    65  			maxabsWant = 1
    66  		}
    67  	}
    68  	if uplo == blas.Upper {
    69  		for i := 0; i < min(m, n); i++ {
    70  			start := i
    71  			if diag == blas.Unit {
    72  				start = i + 1
    73  			}
    74  			for j := start; j < n; j++ {
    75  				aij := 2*rnd.Float64() - 1
    76  				a[i*lda+j] = aij
    77  				rowsum[i] += math.Abs(aij)
    78  				colsum[j] += math.Abs(aij)
    79  				maxabsWant = math.Max(maxabsWant, math.Abs(aij))
    80  				frobWant += aij * aij
    81  			}
    82  		}
    83  	} else {
    84  		for i := 0; i < m; i++ {
    85  			end := i
    86  			if diag == blas.Unit {
    87  				end = i - 1
    88  			}
    89  			for j := 0; j <= min(end, n-1); j++ {
    90  				aij := 2*rnd.Float64() - 1
    91  				a[i*lda+j] = aij
    92  				rowsum[i] += math.Abs(aij)
    93  				colsum[j] += math.Abs(aij)
    94  				maxabsWant = math.Max(maxabsWant, math.Abs(aij))
    95  				frobWant += aij * aij
    96  			}
    97  		}
    98  	}
    99  	frobWant = math.Sqrt(frobWant)
   100  	var maxcolsumWant, maxrowsumWant float64
   101  	if n > 0 {
   102  		maxcolsumWant = floats.Max(colsum)
   103  	}
   104  	if m > 0 {
   105  		maxrowsumWant = floats.Max(rowsum)
   106  	}
   107  
   108  	aCopy := make([]float64, len(a))
   109  	copy(aCopy, a)
   110  
   111  	for _, norm := range []lapack.MatrixNorm{lapack.MaxAbs, lapack.MaxColumnSum, lapack.MaxRowSum, lapack.Frobenius} {
   112  		name := fmt.Sprintf("norm=%v,uplo=%v,diag=%v,m=%v,n=%v,lda=%v", string(norm), string(uplo), string(diag), m, n, lda)
   113  
   114  		var work []float64
   115  		if norm == lapack.MaxColumnSum {
   116  			work = make([]float64, n)
   117  		}
   118  		normGot := impl.Dlantr(norm, uplo, diag, m, n, a, lda, work)
   119  
   120  		if !floats.Equal(a, aCopy) {
   121  			t.Fatalf("%v: unexpected modification of a", name)
   122  		}
   123  
   124  		if norm == lapack.MaxAbs {
   125  			// MaxAbs norm involves no floating-point computation,
   126  			// so we expect exact equality here.
   127  			if normGot != maxabsWant {
   128  				t.Errorf("%v: unexpected result; got %v, want %v", name, normGot, maxabsWant)
   129  			}
   130  			continue
   131  		}
   132  
   133  		var normWant float64
   134  		switch norm {
   135  		case lapack.MaxColumnSum:
   136  			normWant = maxcolsumWant
   137  		case lapack.MaxRowSum:
   138  			normWant = maxrowsumWant
   139  		case lapack.Frobenius:
   140  			normWant = frobWant
   141  		}
   142  		if math.Abs(normGot-normWant) >= tol {
   143  			t.Errorf("%v: unexpected result; got %v, want %v", name, normGot, normWant)
   144  		}
   145  	}
   146  }