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

     1  // Copyright ©2019 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/floats"
    16  	"github.com/jingcheng-WU/gonum/lapack"
    17  )
    18  
    19  type Dlansber interface {
    20  	Dlansb(norm lapack.MatrixNorm, uplo blas.Uplo, n, kd int, ab []float64, ldab int, work []float64) float64
    21  }
    22  
    23  func DlansbTest(t *testing.T, impl Dlansber) {
    24  	rnd := rand.New(rand.NewSource(1))
    25  	for _, n := range []int{0, 1, 2, 3, 4, 5, 10} {
    26  		for _, kd := range []int{0, (n + 1) / 4, (3*n - 1) / 4, (5*n + 1) / 4} {
    27  			for _, uplo := range []blas.Uplo{blas.Upper, blas.Lower} {
    28  				for _, ldab := range []int{kd + 1, kd + 1 + 7} {
    29  					dlansbTest(t, impl, rnd, uplo, n, kd, ldab)
    30  				}
    31  			}
    32  		}
    33  	}
    34  }
    35  
    36  func dlansbTest(t *testing.T, impl Dlansber, rnd *rand.Rand, uplo blas.Uplo, n, kd int, ldab int) {
    37  	const tol = 1e-15
    38  
    39  	// Generate a random symmetric band matrix and compute all its norms.
    40  	ab := make([]float64, max(0, (n-1)*ldab+kd+1))
    41  	rowsum := make([]float64, n)
    42  	colsum := make([]float64, n)
    43  	var frobWant, maxabsWant float64
    44  	if uplo == blas.Upper {
    45  		for i := 0; i < n; i++ {
    46  			for jb := 0; jb < min(n-i, kd+1); jb++ {
    47  				aij := 2*rnd.Float64() - 1
    48  				ab[i*ldab+jb] = aij
    49  
    50  				j := jb + i
    51  				colsum[j] += math.Abs(aij)
    52  				rowsum[i] += math.Abs(aij)
    53  				maxabsWant = math.Max(maxabsWant, math.Abs(aij))
    54  				frobWant += aij * aij
    55  				if i != j {
    56  					// Take into account the symmetric elements.
    57  					colsum[i] += math.Abs(aij)
    58  					rowsum[j] += math.Abs(aij)
    59  					frobWant += aij * aij
    60  				}
    61  			}
    62  		}
    63  	} else {
    64  		for i := 0; i < n; i++ {
    65  			for jb := max(0, kd-i); jb < kd+1; jb++ {
    66  				aij := 2*rnd.Float64() - 1
    67  				ab[i*ldab+jb] = aij
    68  
    69  				j := jb - kd + i
    70  				colsum[j] += math.Abs(aij)
    71  				rowsum[i] += math.Abs(aij)
    72  				maxabsWant = math.Max(maxabsWant, math.Abs(aij))
    73  				frobWant += aij * aij
    74  				if i != j {
    75  					// Take into account the symmetric elements.
    76  					colsum[i] += math.Abs(aij)
    77  					rowsum[j] += math.Abs(aij)
    78  					frobWant += aij * aij
    79  				}
    80  			}
    81  		}
    82  	}
    83  	frobWant = math.Sqrt(frobWant)
    84  	var maxcolsumWant, maxrowsumWant float64
    85  	if n > 0 {
    86  		maxcolsumWant = floats.Max(colsum)
    87  		maxrowsumWant = floats.Max(rowsum)
    88  	}
    89  
    90  	abCopy := make([]float64, len(ab))
    91  	copy(abCopy, ab)
    92  
    93  	work := make([]float64, n)
    94  	var maxcolsumGot, maxrowsumGot float64
    95  	for _, norm := range []lapack.MatrixNorm{lapack.MaxAbs, lapack.MaxColumnSum, lapack.MaxRowSum, lapack.Frobenius} {
    96  		name := fmt.Sprintf("norm=%v,uplo=%v,n=%v,kd=%v,ldab=%v", string(norm), string(uplo), n, kd, ldab)
    97  
    98  		normGot := impl.Dlansb(norm, uplo, n, kd, ab, ldab, work)
    99  
   100  		if !floats.Equal(ab, abCopy) {
   101  			t.Fatalf("%v: unexpected modification of ab", name)
   102  		}
   103  
   104  		if norm == lapack.MaxAbs {
   105  			// MaxAbs norm involves no computation, so we expect
   106  			// exact equality here.
   107  			if normGot != maxabsWant {
   108  				t.Errorf("%v: unexpected result; got %v, want %v", name, normGot, maxabsWant)
   109  			}
   110  			continue
   111  		}
   112  
   113  		var normWant float64
   114  		switch norm {
   115  		case lapack.MaxColumnSum:
   116  			normWant = maxcolsumWant
   117  			maxcolsumGot = normGot
   118  		case lapack.MaxRowSum:
   119  			normWant = maxrowsumWant
   120  			maxrowsumGot = normGot
   121  		case lapack.Frobenius:
   122  			normWant = frobWant
   123  		}
   124  		if math.Abs(normGot-normWant) > tol*float64(n) {
   125  			t.Errorf("%v: unexpected result; got %v, want %v", name, normGot, normWant)
   126  		}
   127  	}
   128  	// MaxColSum and MaxRowSum norms should be exactly equal because the
   129  	// matrix is symmetric.
   130  	if maxcolsumGot != maxrowsumGot {
   131  		name := fmt.Sprintf("uplo=%v,n=%v,kd=%v,ldab=%v", string(uplo), n, kd, ldab)
   132  		t.Errorf("%v: unexpected mismatch between MaxColSum and MaxRowSum norms of A; MaxColSum %v, MaxRowSum %v", name, maxcolsumGot, maxrowsumGot)
   133  	}
   134  }