Small cleanups

native/general.go

@@ -90,13 +90,10 @@ var (
 
 	// dlamchP is 2 * eps
 	dlamchP = math.Float64frombits(0x3cb0000000000000)
-	//dlamchP = 2.2204460492503131E-016
 
 	// dlamchS is the "safe min", that is, the lowest number such that 1/sfmin does
 	// not overflow. The Netlib code for calculating this number is not correct --
-	// it overflows. Found by trial and error, it is equal to (1/math.MaxFloat64) * (1+ 6*eps)
-	//dlamchS = math.Float64frombits(0x4000000000001)
-
+	// it overflows.
 	// Found from printing out from FORTRAN
 	dlamchS = 2.2250738585072014E-308
 

testlapack/dlasq1.go

@@ -40,7 +40,9 @@ type Dlasq1er interface {
 
 func Dlasq1Test(t *testing.T, impl Dlasq1er) {
 	bi := blas64.Implementation()
-	for _, n := range []int{100} {
+	// TODO(btracey): Increase the size of this test when we have a more numerically
+	// stable way to test the singular values.
+	for _, n := range []int{1, 2, 5, 8} {
 		work := make([]float64, 4*n)
 		d := make([]float64, n)
 		e := make([]float64, n-1)
@@ -53,15 +55,13 @@ func Dlasq1Test(t *testing.T, impl Dlasq1er) {
 			}
 			for i := range e {
 				e[i] = rand.NormFloat64()
-				//e[i] = 0
 			}
 			ldm := n
 			m := make([]float64, n*ldm)
 			// Set up the matrix
 			for i := 0; i < n; i++ {
 				m[i*ldm+i] = d[i]
 				if i != n-1 {
-					//m[i*ldm+i+1] = e[i]
 					m[(i+1)*ldm+i] = e[i]
 				}
 			}
@@ -70,26 +70,11 @@ func Dlasq1Test(t *testing.T, impl Dlasq1er) {
 			mm := make([]float64, n*ldmm)
 			bi.Dgemm(blas.Trans, blas.NoTrans, n, n, n, 1, m, ldm, m, ldm, 0, mm, ldmm)
 
-			//fmt.Println("m = ")
-			//printRowise(m, n, n, ldm, false)
-			//fmt.Println("mm = ")
-			//printRowise(mm, n, n, ldm, false)
-
-			// fmt.Println("d = ", d)
-			// fmt.Println("e = ", e)
-			// fmt.Println("m = ", m)
-			// Compute the singular values
-
 			//printDlasq1FortranInput(d, e, work, n)
 			//os.Exit(1)
 
 			impl.Dlasq1(n, d, e, work)
 
-			fmt.Println("m = ", m)
-
-			// Check that the singular values are sorted
-			fmt.Println("singular values", d[0:n])
-
 			// Check that they are singular values. The
 			// singular values are the square roots of the
 			// eigenvalues of X^T * X
@@ -105,14 +90,16 @@ func Dlasq1Test(t *testing.T, impl Dlasq1er) {
 
 				// Compute LU.
 				ok := impl.Dgetrf(n, n, mm, ldmm, ipiv)
-				fmt.Println("ok", ok)
+				if !ok {
+					// Definitely singular.
+					continue
+				}
 				// Compute determinant
 				var logdet float64
 				for i := 0; i < n; i++ {
 					v := mm[i*ldm+i]
 					logdet += math.Log(math.Abs(v))
 				}
-				//fmt.Println(math.Exp(logdet))
 				if math.Exp(logdet) > 2 {
 					t.Errorf("Incorrect singular value. n = %d, cas = %d, elem = %d, det = %v", n, cas, elem, math.Exp(logdet))
 				}

testlapack/dlasq2.go

@@ -6,8 +6,6 @@ import (
 	"math/rand"
 	"testing"
 
-	"github.com/gonum/blas"
-	"github.com/gonum/blas/blas64"
 	"github.com/gonum/floats"
 )
 
@@ -643,111 +641,66 @@ func Dlasq2Test(t *testing.T, impl Dlasq2er) {
 
 		info := impl.Dlasq2(test.n, z)
 		if !floats.EqualApprox(test.zOut, z, dTol) {
-			fmt.Println("not approx")
-			//t.Errorf("Case %v: Z mismatch. Want %v, got %v", c, test.zOut, test.z)
 			diff := make([]float64, len(z))
 			floats.SubTo(diff, z, test.zOut)
 			for i := range diff {
 				diff[i] = math.Abs(diff[i])
 			}
-			fmt.Println("max diff = ", floats.Max(diff))
-			maxidx := floats.MaxIdx(diff)
-			fmt.Println("max idx = ", maxidx)
-			fmt.Println("got", z[maxidx], "want", test.zOut[maxidx])
 			t.Errorf("Case %v, Z Mismatch", c)
 		}
 		if test.info != info {
 			t.Errorf("Info mismatch. Want %v, got %v", test.info, info)
 		}
 	}
 
-	// This still gets stuck in a for loop. Figure out which one does and
-	// test that case.
-	// Woo! No longer stuck.
-	bi := blas64.Implementation()
-	//for _, n := range []int{5, 8, 20, 21} {
-	for _, n := range []int{21} {
-		for k := 0; k < 1; k++ {
+	// Perform a bunch of random tests to check for access out of bounds or
+	// infinite loops.
+	// TODO(btracey): Implement direct tests.
+	// bi := blas64.Implementation()
+	for _, n := range []int{5, 8, 20, 25} {
+		for k := 0; k < 10; k++ {
 			z := make([]float64, 4*n)
 			for i := range z {
 				z[i] = rand.Float64()
 			}
 			zCopy := make([]float64, len(z))
 			copy(zCopy, z)
 
-			z = []float64{0.7940281584071446E+00, 0.8540600349699839E+00, 0.8158431165852809E-01, 0.5431841788581357E+00, 0.3696613346727944E+00, 0.2343742079469738E+00, 0.4891909888056500E-01, 0.6769876984160987E+00, 0.4777767465052760E+00, 0.1867381312399053E+00, 0.2018744873845245E+00, 0.5511201479607295E+00, 0.6938788283912793E+00, 0.8167542438070282E+00, 0.7904606414789531E+00, 0.9443564310071292E+00, 0.7287247677237652E-01, 0.8645122013586991E+00, 0.1884651475116826E+00, 0.3844755283611681E+00, 0.9959264361467982E+00, 0.6424370932833342E+00, 0.1972122925077952E+00, 0.2842024247377670E+00, 0.9819646913482807E+00, 0.9118347224008859E+00, 0.8184691845197246E+00, 0.7051587281589254E+00, 0.7604703230109544E+00, 0.6312964755149379E+00, 0.5240863862347888E+00, 0.3442050916384676E-01, 0.2415614308212055E+00, 0.2814868323669945E+00, 0.6529284673126197E+00, 0.3727305084153835E+00, 0.5033733868757848E+00, 0.2317122058804952E+00, 0.7555584130128312E+00, 0.5854566742645219E+00, 0.5481204696337160E+00, 0.8479425268049923E+00, 0.2310874615764000E+00, 0.1250993726775007E-01, 0.6243285982203539E-01, 0.8533587246073391E+00, 0.9203815588639257E+00, 0.9256849509751471E+00, 0.6691405057262187E+00, 0.8847091531299658E+00, 0.6783572983386376E+00, 0.4701257141291857E+00, 0.8976078424378102E+00, 0.8575018884445876E+00, 0.4119363561363949E+00, 0.2824477027676924E+00, 0.2787507690368071E+00, 0.7994878185780909E+00, 0.6141832897278305E+00, 0.6772728066124333E+00, 0.1568652581579784E+00, 0.8025492691231176E+00, 0.2609459151100056E+00, 0.4956700691019098E+00, 0.1008839464621498E+00, 0.6129709499983976E+00, 0.4551038858718992E-02, 0.8382785474023564E+00, 0.9327452694814308E+00, 0.9710431593941808E+00, 0.3785578217695214E+00, 0.9620839159000718E+00, 0.3183561960196257E-01, 0.9167635157854341E+00, 0.8989971039988554E+00, 0.2723769512210017E-01, 0.4176537489735596E+00, 0.9619881273217982E+00, 0.8761769579995293E+00, 0.6385245520487358E+00, 0.6821739872929905E+00, 0.3927943300877799E+00, 0.3299501391296433E-01, 0.6026481165267817E+00}
-
-			// Construct the full tridiagonal matrices.
-			ldl := n
-			ldu := n
-			u := make([]float64, n*n)
-			for i := 0; i < n; i++ {
-				u[i*ldu+i] = z[2*i]
-				if i != n-1 {
-					u[i*ldu+i+1] = 1
-				}
-			}
-			l := make([]float64, n*n)
-			for i := 0; i < n; i++ {
-				l[i*ldl+i] = 1
-				if i != n-1 {
-					l[(i+1)*ldl+i] = z[2*i+1]
-				}
-			}
-			fmt.Println("z[0:4]", z[0:4])
-
-			fmt.Println("u = ")
-			printRowise(u, n, n, ldu, false)
-
-			fmt.Println("l = ")
-			printRowise(l, n, n, ldl, false)
-
-			ldTriDi := n
-			triDi := make([]float64, n*n)
-			bi.Dgemm(blas.NoTrans, blas.NoTrans, n, n, n, 1, l, ldl, u, ldu, 0, triDi, ldTriDi)
-
-			fmt.Println("u = ")
-			printRowise(u, n, n, ldu, false)
-
-			fmt.Println("tridi = ")
-
-			printRowise(triDi, n, n, n, false)
-
-			//triDi2 := make([]float64, n*n)
-
 			// Compute the eigenvalues
 			impl.Dlasq2(n, z)
 
-			tridi2 := make([]float64, n*n)
-			bi.Dgemm(blas.Trans, blas.NoTrans, n, n, n, 1, triDi, n, triDi, n, 0, tridi2, n)
+			// Below is the code to test the eigenvalues. Eventually implement
+			// real tests.
+			// The code below is missing the transformation from L and U into
+			// the symmetric tridiagonal matrix.
+			// See discussion http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=5&t=4839
+			// for format.
 
-			// Check that they are singular values. The
-			// singular values are the square roots of the
-			// eigenvalues of X^T * X
-			triDiCopy := make([]float64, len(tridi2))
-			copy(triDiCopy, tridi2)
-			ipiv := make([]int, n)
-			for _, sv := range z[0:n] {
-				copy(tridi2, triDiCopy)
-				//lambda := sv * sv
-				lambda := sv
+			/*
+				ldl := n
+				ldu := n
+				u := make([]float64, n*n)
 				for i := 0; i < n; i++ {
-					tridi2[i*n+i] -= lambda
+					u[i*ldu+i] = zCopy[2*i]
+					if i != n-1 {
+						u[i*ldu+i+1] = 1
+					}
 				}
-				// Compute LU.
-				impl.Dgetrf(n, n, tridi2, n, ipiv)
-				// Compute determinant
-				var logdet float64
+				l := make([]float64, n*n)
 				for i := 0; i < n; i++ {
-					v := tridi2[i*n+i]
-					logdet += math.Log(math.Abs(v))
-				}
-				if math.Exp(logdet) > 1 {
-					t.Errorf("Incorrect singular value. n = %d, k = %d, det = %v", n, k, math.Exp(logdet))
+					l[i*ldl+i] = 1
+					if i != n-1 {
+						l[(i+1)*ldl+i] = zCopy[2*i+1]
+					}
 				}
-			}
 
-			/*
+				ldTriDi := n
+				triDi := make([]float64, n*n)
+				bi.Dgemm(blas.NoTrans, blas.NoTrans, n, n, n, 1, l, ldl, u, ldu, 0, triDi, ldTriDi)
+
+				tridi2 := make([]float64, n*n)
+				bi.Dgemm(blas.Trans, blas.NoTrans, n, n, n, 1, triDi, n, triDi, n, 0, tridi2, n)
+
 				// Eigenvalues have the property that det(A - lambda I ) = 0
 				triDiCopy := make([]float64, len(triDi))
 				copy(triDiCopy, triDi)

testlapack/dlasq3.go

@@ -1,7 +1,6 @@
 package testlapack
 
 import (
-	"fmt"
 	"math"
 	"testing"
 
@@ -2641,13 +2640,6 @@ func Dlasq3Test(t *testing.T, impl Dlasq3er) {
 			impl.Dlasq3(i0, n0, z, test.pp, test.dmin, test.sigma, test.desig, test.qmax, test.nFail, test.iter, test.nDiv, test.ttype, test.dmin1, test.dmin2, test.dn, test.dn1, test.dn2, test.g, test.tau)
 
 		if !floats.EqualApprox(z, test.zOut, dTol) {
-			//fmt.Println("z want")
-			//fmt.Printf("%0.16v\n", test.zOut)
-			//fmt.Println("z got")
-			//fmt.Printf("%0.16v\n", z)
-			floats.Sub(z, test.zOut)
-			//fmt.Println("diff = ", z)
-			fmt.Println("max z diff = ", floats.Max(z))
 			t.Errorf("Z mismatch")
 		}
 		if i0Out != test.i0Out-1 {

testlapack/dlasq4.go

@@ -3107,7 +3107,7 @@ func Dlasq4Test(t *testing.T, impl Dlasq4er) {
 		copy(z, test.z)
 
 		// Print for fortran input
-		printDlasq4FortranInput(test)
+		//printDlasq4FortranInput(test)
 
 		i0 := test.i0 - 1 // zero index
 		n0 := test.n0 - 1 // zero index