github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/lapack/testlapack/dlansy.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  	"math"
     9  	"testing"
    10  
    11  	"golang.org/x/exp/rand"
    12  
    13  	"github.com/jingcheng-WU/gonum/blas"
    14  	"github.com/jingcheng-WU/gonum/lapack"
    15  )
    16  
    17  type Dlansyer interface {
    18  	Dlanger
    19  	Dlansy(norm lapack.MatrixNorm, uplo blas.Uplo, n int, a []float64, lda int, work []float64) float64
    20  }
    21  
    22  func DlansyTest(t *testing.T, impl Dlansyer) {
    23  	rnd := rand.New(rand.NewSource(1))
    24  	for _, norm := range []lapack.MatrixNorm{lapack.MaxAbs, lapack.MaxColumnSum, lapack.MaxRowSum, lapack.Frobenius} {
    25  		for _, uplo := range []blas.Uplo{blas.Lower, blas.Upper} {
    26  			for _, test := range []struct {
    27  				n, lda int
    28  			}{
    29  				{1, 0},
    30  				{3, 0},
    31  
    32  				{1, 10},
    33  				{3, 10},
    34  			} {
    35  				for trial := 0; trial < 100; trial++ {
    36  					n := test.n
    37  					lda := test.lda
    38  					if lda == 0 {
    39  						lda = n
    40  					}
    41  					// Allocate n×n matrix A and fill it.
    42  					// Only the uplo triangle of A will be used below
    43  					// to represent a symmetric matrix.
    44  					a := make([]float64, lda*n)
    45  					if trial == 0 {
    46  						// In the first trial fill the matrix
    47  						// with predictable integers.
    48  						for i := range a {
    49  							a[i] = float64(i)
    50  						}
    51  					} else {
    52  						// Otherwise fill it with random numbers.
    53  						for i := range a {
    54  							a[i] = rnd.NormFloat64()
    55  						}
    56  					}
    57  
    58  					// Create a dense representation of the symmetric matrix
    59  					// stored in the uplo triangle of A.
    60  					aDense := make([]float64, n*n)
    61  					if uplo == blas.Upper {
    62  						for i := 0; i < n; i++ {
    63  							for j := i; j < n; j++ {
    64  								v := a[i*lda+j]
    65  								aDense[i*n+j] = v
    66  								aDense[j*n+i] = v
    67  							}
    68  						}
    69  					} else {
    70  						for i := 0; i < n; i++ {
    71  							for j := 0; j <= i; j++ {
    72  								v := a[i*lda+j]
    73  								aDense[i*n+j] = v
    74  								aDense[j*n+i] = v
    75  							}
    76  						}
    77  					}
    78  
    79  					work := make([]float64, n)
    80  					// Compute the norm of the symmetric matrix A.
    81  					got := impl.Dlansy(norm, uplo, n, a, lda, work)
    82  					// Compute the reference norm value using Dlange
    83  					// and the dense representation of A.
    84  					want := impl.Dlange(norm, n, n, aDense, n, work)
    85  					if math.Abs(want-got) > 1e-14 {
    86  						t.Errorf("Norm mismatch. norm = %c, upper = %v, n = %v, lda = %v, want %v, got %v.",
    87  							norm, uplo == blas.Upper, n, lda, got, want)
    88  					}
    89  				}
    90  			}
    91  		}
    92  	}
    93  }