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