github.com/gopherd/gonum@v0.0.4/lapack/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  	"sort"
    11  	"testing"
    12  
    13  	"math/rand"
    14  
    15  	"github.com/gopherd/gonum/floats"
    16  )
    17  
    18  type Dlasq1er interface {
    19  	Dlasq1(n int, d, e, work []float64) int
    20  
    21  	Dgebrd(m, n int, a []float64, lda int, d, e, tauQ, tauP, work []float64, lwork int)
    22  }
    23  
    24  func Dlasq1Test(t *testing.T, impl Dlasq1er) {
    25  	const tol = 1e-14
    26  	rnd := rand.New(rand.NewSource(1))
    27  	for _, n := range []int{0, 1, 2, 3, 4, 5, 8, 10, 30, 50} {
    28  		for typ := 0; typ <= 7; typ++ {
    29  			name := fmt.Sprintf("n=%v,typ=%v", n, typ)
    30  
    31  			// Generate a diagonal matrix D with positive entries.
    32  			d := make([]float64, n)
    33  			switch typ {
    34  			case 0:
    35  				// The zero matrix.
    36  			case 1:
    37  				// The identity matrix.
    38  				for i := range d {
    39  					d[i] = 1
    40  				}
    41  			case 2:
    42  				// A diagonal matrix with evenly spaced entries 1, ..., eps.
    43  				for i := 0; i < n; i++ {
    44  					if i == 0 {
    45  						d[0] = 1
    46  					} else {
    47  						d[i] = 1 - (1-dlamchE)*float64(i)/float64(n-1)
    48  					}
    49  				}
    50  			case 3, 4, 5:
    51  				// A diagonal matrix with geometrically spaced entries 1, ..., eps.
    52  				for i := 0; i < n; i++ {
    53  					if i == 0 {
    54  						d[0] = 1
    55  					} else {
    56  						d[i] = math.Pow(dlamchE, float64(i)/float64(n-1))
    57  					}
    58  				}
    59  				switch typ {
    60  				case 4:
    61  					// Multiply by SQRT(overflow threshold).
    62  					floats.Scale(math.Sqrt(1/dlamchS), d)
    63  				case 5:
    64  					// Multiply by SQRT(underflow threshold).
    65  					floats.Scale(math.Sqrt(dlamchS), d)
    66  				}
    67  			case 6:
    68  				// A diagonal matrix with "clustered" entries 1, eps, ..., eps.
    69  				for i := range d {
    70  					if i == 0 {
    71  						d[i] = 1
    72  					} else {
    73  						d[i] = dlamchE
    74  					}
    75  				}
    76  			case 7:
    77  				// Diagonal matrix with random entries.
    78  				for i := range d {
    79  					d[i] = math.Abs(rnd.NormFloat64())
    80  				}
    81  			}
    82  
    83  			dWant := make([]float64, n)
    84  			copy(dWant, d)
    85  			sort.Sort(sort.Reverse(sort.Float64Slice(dWant)))
    86  
    87  			// Allocate work slice to the maximum length needed below.
    88  			work := make([]float64, max(1, 4*n))
    89  
    90  			// Generate an n×n matrix A by pre- and post-multiplying D with
    91  			// random orthogonal matrices:
    92  			//  A = U*D*V.
    93  			lda := max(1, n)
    94  			a := make([]float64, n*lda)
    95  			Dlagge(n, n, 0, 0, d, a, lda, rnd, work)
    96  
    97  			// Reduce A to bidiagonal form B represented by the diagonal d and
    98  			// off-diagonal e.
    99  			tauQ := make([]float64, n)
   100  			tauP := make([]float64, n)
   101  			e := make([]float64, max(0, n-1))
   102  			impl.Dgebrd(n, n, a, lda, d, e, tauQ, tauP, work, len(work))
   103  
   104  			// Compute the singular values of B.
   105  			for i := range work {
   106  				work[i] = math.NaN()
   107  			}
   108  			info := impl.Dlasq1(n, d, e, work)
   109  			if info != 0 {
   110  				t.Fatalf("%v: Dlasq1 returned non-zero info=%v", name, info)
   111  			}
   112  
   113  			if n == 0 {
   114  				continue
   115  			}
   116  
   117  			if !sort.IsSorted(sort.Reverse(sort.Float64Slice(d))) {
   118  				t.Errorf("%v: singular values not sorted", name)
   119  			}
   120  
   121  			diff := floats.Distance(d, dWant, math.Inf(1))
   122  			if diff > tol*floats.Max(dWant) {
   123  				t.Errorf("%v: unexpected result; diff=%v", name, diff)
   124  			}
   125  		}
   126  	}
   127  }