github.com/gopherd/gonum@v0.0.4/lapack/testlapack/dlacn2.go (about)

     1  // Copyright ©2016 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  	"math/rand"
    12  
    13  	"github.com/gopherd/gonum/blas"
    14  	"github.com/gopherd/gonum/blas/blas64"
    15  )
    16  
    17  type Dlacn2er interface {
    18  	Dlacn2(n int, v, x []float64, isgn []int, est float64, kase int, isave *[3]int) (float64, int)
    19  }
    20  
    21  func Dlacn2Test(t *testing.T, impl Dlacn2er) {
    22  	rnd := rand.New(rand.NewSource(1))
    23  	for _, n := range []int{1, 2, 3, 4, 5, 7, 10, 15, 20, 100} {
    24  		for cas := 0; cas < 10; cas++ {
    25  			a := randomGeneral(n, n, n, rnd)
    26  
    27  			// Compute the 1-norm of A explicitly.
    28  			var norm1 float64
    29  			for j := 0; j < n; j++ {
    30  				var sum float64
    31  				for i := 0; i < n; i++ {
    32  					sum += math.Abs(a.Data[i*a.Stride+j])
    33  				}
    34  				if sum > norm1 {
    35  					norm1 = sum
    36  				}
    37  			}
    38  
    39  			// Compute the estimate of 1-norm using Dlanc2.
    40  			x := make([]float64, n)
    41  			work := make([]float64, n)
    42  			v := make([]float64, n)
    43  			isgn := make([]int, n)
    44  			var (
    45  				kase  int
    46  				isave [3]int
    47  				got   float64
    48  			)
    49  		loop:
    50  			for {
    51  				got, kase = impl.Dlacn2(n, v, x, isgn, got, kase, &isave)
    52  				switch kase {
    53  				default:
    54  					panic("Dlacn2 returned invalid value of kase")
    55  				case 0:
    56  					break loop
    57  				case 1:
    58  					blas64.Gemv(blas.NoTrans, 1, a, blas64.Vector{Data: x, Inc: 1}, 0, blas64.Vector{Data: work, Inc: 1})
    59  					copy(x, work)
    60  				case 2:
    61  					blas64.Gemv(blas.Trans, 1, a, blas64.Vector{Data: x, Inc: 1}, 0, blas64.Vector{Data: work, Inc: 1})
    62  					copy(x, work)
    63  				}
    64  			}
    65  
    66  			// Check that got is either accurate enough or a
    67  			// lower estimate of the 1-norm of A.
    68  			if math.Abs(got-norm1) > 1e-8 && got > norm1 {
    69  				t.Errorf("Case n=%v: not lower estimate. 1-norm %v, estimate %v", n, norm1, got)
    70  			}
    71  		}
    72  	}
    73  }