github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/testlapack/dlasq1.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  	"math/rand"
    11  	"testing"
    12  
    13  	"github.com/gonum/blas"
    14  	"github.com/gonum/blas/blas64"
    15  )
    16  
    17  func printDlasq1FortranInput(d, e, work []float64, n int) {
    18  	printFortranArray(d, "d")
    19  	printFortranArray(e, "e")
    20  	printFortranArray(work, "work")
    21  	fmt.Println("n = ", n)
    22  	fmt.Println("info = 0")
    23  }
    24  
    25  type Dlasq1er interface {
    26  	Dlasq1(n int, d, e, work []float64) int
    27  	Dgetrfer
    28  }
    29  
    30  func Dlasq1Test(t *testing.T, impl Dlasq1er) {
    31  	rnd := rand.New(rand.NewSource(1))
    32  	bi := blas64.Implementation()
    33  	// TODO(btracey): Increase the size of this test when we have a more numerically
    34  	// stable way to test the singular values.
    35  	for _, n := range []int{1, 2, 5, 8} {
    36  		work := make([]float64, 4*n)
    37  		d := make([]float64, n)
    38  		e := make([]float64, n-1)
    39  		for cas := 0; cas < 1; cas++ {
    40  			for i := range work {
    41  				work[i] = rnd.Float64()
    42  			}
    43  			for i := range d {
    44  				d[i] = rnd.NormFloat64() + 10
    45  			}
    46  			for i := range e {
    47  				e[i] = rnd.NormFloat64()
    48  			}
    49  			ldm := n
    50  			m := make([]float64, n*ldm)
    51  			// Set up the matrix
    52  			for i := 0; i < n; i++ {
    53  				m[i*ldm+i] = d[i]
    54  				if i != n-1 {
    55  					m[(i+1)*ldm+i] = e[i]
    56  				}
    57  			}
    58  
    59  			ldmm := n
    60  			mm := make([]float64, n*ldmm)
    61  			bi.Dgemm(blas.Trans, blas.NoTrans, n, n, n, 1, m, ldm, m, ldm, 0, mm, ldmm)
    62  
    63  			impl.Dlasq1(n, d, e, work)
    64  
    65  			// Check that they are singular values. The
    66  			// singular values are the square roots of the
    67  			// eigenvalues of X^T * X
    68  			mmCopy := make([]float64, len(mm))
    69  			copy(mmCopy, mm)
    70  			ipiv := make([]int, n)
    71  			for elem, sv := range d[0:n] {
    72  				copy(mm, mmCopy)
    73  				lambda := sv * sv
    74  				for i := 0; i < n; i++ {
    75  					mm[i*ldm+i] -= lambda
    76  				}
    77  
    78  				// Compute LU.
    79  				ok := impl.Dgetrf(n, n, mm, ldmm, ipiv)
    80  				if !ok {
    81  					// Definitely singular.
    82  					continue
    83  				}
    84  				// Compute determinant
    85  				var logdet float64
    86  				for i := 0; i < n; i++ {
    87  					v := mm[i*ldm+i]
    88  					logdet += math.Log(math.Abs(v))
    89  				}
    90  				if math.Exp(logdet) > 2 {
    91  					t.Errorf("Incorrect singular value. n = %d, cas = %d, elem = %d, det = %v", n, cas, elem, math.Exp(logdet))
    92  				}
    93  			}
    94  		}
    95  	}
    96  }