Added dlarf, tests, some auxiliary routines, and missing license headers

Added dlarf, tests, some auxiliary routines, and missing license headers

PR comments

lapack.go

@@ -1,3 +1,7 @@
+// Copyright ©2015 The gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
 package lapack
 
 import "github.com/gonum/blas"

native/dlarf.go

@@ -0,0 +1,70 @@
+// Copyright ©2015 The gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package native
+
+import (
+	"github.com/gonum/blas"
+	"github.com/gonum/blas/blas64"
+)
+
+// Dlarf applies an elementary reflector to a general rectangular matrix c.
+// This computes
+//  c = h * c if side == Left
+//  c = c * h if side == right
+// where
+//  h = 1 - tau * v * v^T
+// and c is an m * n matrix.
+//
+// Work is pre-allocated memory of size at least m if side == Left and at least
+// n if side == Right. Dlarf will panic if this length requirement is not met.
+func (impl Implementation) Dlarf(side blas.Side, m, n int, v []float64, incv int, tau float64, c []float64, ldc int, work []float64) {
+	applyleft := side == blas.Left
+	if (applyleft && len(work) < n) || (!applyleft && len(work) < m) {
+		panic("dlarf: insufficient work length")
+	}
+	lastv := 0 // last non-zero element of v
+	lastc := 0 // last non-zero row/column of c
+	if tau != 0 {
+		var i int
+		if applyleft {
+			lastv = m
+		} else {
+			lastv = n
+		}
+		if incv > 0 {
+			i = (lastv - 1) * incv
+		}
+		// Look for the last non-zero row in v.
+		for lastv >= 0 && v[i] == 0 {
+			lastv--
+			i -= incv
+		}
+		if applyleft {
+			// Scan for the last non-zero column in C[0:lastv, :]
+			lastc = impl.Iladlc(lastv, n, c, ldc)
+		} else {
+			// Scan for the last non-zero row in C[:, 0:lastv]
+			lastc = impl.Iladlr(m, lastv, c, ldc)
+		}
+	}
+	if lastv == -1 || lastc == -1 {
+		return
+	}
+	// Sometimes 1-indexing is nicer ...
+	bi := blas64.Implementation()
+	if applyleft {
+		// Form H * C
+		// w[0:lastc+1] = c[1:lastv+1, 1:lastc+1]^T * v[1:lastv+1,1]
+		bi.Dgemv(blas.Trans, lastv+1, lastc+1, 1, c, ldc, v, incv, 0, work, 1)
+		// c[0: lastv, 0: lastc] = c[...] - w[0:lastv, 1] * v[1:lastc, 1]^T
+		bi.Dger(lastv+1, lastc+1, -tau, v, incv, work, 1, c, ldc)
+		return
+	}
+	// Form C*H
+	// w[0:lastc+1,1] := c[0:lastc+1,0:lastv+1] * v[0:lastv+1,1]
+	bi.Dgemv(blas.NoTrans, lastc+1, lastv+1, 1, c, ldc, v, incv, 0, work, 1)
+	// c[0:lastc+1,0:lastv+1] = c[...] - w[0:lastc+1,0] * v[0:lastv+1,0]^T
+	bi.Dger(lastc+1, lastv+1, -tau, work, 1, v, incv, c, ldc)
+}

native/dpotf2.go

@@ -1,3 +1,7 @@
+// Copyright ©2015 The gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
 package native
 
 import (

native/dpotrf.go

@@ -1,3 +1,7 @@
+// Copyright ©2015 The gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
 package native
 
 import (

native/general.go

@@ -1,3 +1,7 @@
+// Copyright ©2015 The gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
 package native
 
 import "github.com/gonum/lapack"

native/iladlc.go

@@ -0,0 +1,29 @@
+// Copyright ©2015 The gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package native
+
+// Iladlc scans a matrix for its last non-zero column. Returns -1 if the matrix
+// is all zeros.
+func (Implementation) Iladlc(m, n int, a []float64, lda int) int {
+	if n == 0 {
+		return n - 1
+	}
+	// Test common case where corner is non-zero.
+	if a[n-1] != 0 || a[(m-1)*lda+(n-1)] != 0 {
+		return n - 1
+	}
+
+	// Scan each row tracking the highest column seen.
+	highest := -1
+	for i := 0; i < m; i++ {
+		for j := n - 1; j >= 0; j-- {
+			if a[i*lda+j] != 0 {
+				highest = max(highest, j)
+				break
+			}
+		}
+	}
+	return highest
+}

native/iladlr.go

@@ -0,0 +1,25 @@
+// Copyright ©2015 The gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package native
+
+// Iladlr scans a matrix for its last non-zero row. Returns -1 if the matrix
+// is all zeros.
+func (Implementation) Iladlr(m, n int, a []float64, lda int) int {
+	if m == 0 {
+		return m - 1
+	}
+	// Check the common case where the corner is non-zero
+	if a[(m-1)*lda] != 0 || a[(m-1)*lda+n-1] != 0 {
+		return m - 1
+	}
+	for i := m - 1; i >= 0; i-- {
+		for j := 0; j < n; j++ {
+			if a[i*lda+j] != 0 {
+				return i
+			}
+		}
+	}
+	return -1
+}

native/lapack_test.go

@@ -1,3 +1,7 @@
+// Copyright ©2015 The gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
 package native
 
 import (
@@ -8,6 +12,10 @@ import (
 
 var impl = Implementation{}
 
+func TestDlarf(t *testing.T) {
+	testlapack.DlarfTest(t, impl)
+}
+
 func TestDpotf2(t *testing.T) {
 	testlapack.Dpotf2Test(t, impl)
 }
@@ -15,3 +23,11 @@ func TestDpotf2(t *testing.T) {
 func TestDpotrf(t *testing.T) {
 	testlapack.DpotrfTest(t, impl)
 }
+
+func TestIladlc(t *testing.T) {
+	testlapack.IladlcTest(t, impl)
+}
+
+func TestIladlr(t *testing.T) {
+	testlapack.IladlrTest(t, impl)
+}

testlapack/dlarf.go

@@ -0,0 +1,159 @@
+// Copyright ©2015 The gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package testlapack
+
+import (
+	"math/rand"
+	"testing"
+
+	"github.com/gonum/blas"
+	"github.com/gonum/blas/blas64"
+	"github.com/gonum/floats"
+)
+
+type Dlarfer interface {
+	Dlarf(side blas.Side, m, n int, v []float64, incv int, tau float64, c []float64, ldc int, work []float64)
+}
+
+func DlarfTest(t *testing.T, impl Dlarfer) {
+	for i, test := range []struct {
+		m, n, ldc    int
+		incv, lastv  int
+		lastr, lastc int
+		tau          float64
+	}{
+		{
+			m:   3,
+			n:   2,
+			ldc: 2,
+
+			incv:  4,
+			lastv: 1,
+
+			lastr: 2,
+			lastc: 1,
+
+			tau: 2,
+		},
+		{
+			m:   2,
+			n:   3,
+			ldc: 3,
+
+			incv:  4,
+			lastv: 1,
+
+			lastr: 1,
+			lastc: 2,
+
+			tau: 2,
+		},
+		{
+			m:   2,
+			n:   3,
+			ldc: 3,
+
+			incv:  4,
+			lastv: 1,
+
+			lastr: 0,
+			lastc: 1,
+
+			tau: 2,
+		},
+		{
+			m:   10,
+			n:   10,
+			ldc: 10,
+
+			incv:  4,
+			lastv: 6,
+
+			lastr: 9,
+			lastc: 8,
+
+			tau: 2,
+		},
+	} {
+		// Construct a random matrix.
+		c := make([]float64, test.ldc*test.m)
+		for i := 0; i <= test.lastr; i++ {
+			for j := 0; j <= test.lastc; j++ {
+				c[i*test.ldc+j] = rand.Float64()
+			}
+		}
+		cCopy := make([]float64, len(c))
+		copy(cCopy, c)
+		cCopy2 := make([]float64, len(c))
+		copy(cCopy2, c)
+
+		// Test with side right.
+		sz := max(test.m, test.n) // so v works for both right and left side.
+		v := make([]float64, test.incv*sz+1)
+		// Fill with nonzero entries up until lastv.
+		for i := 0; i < test.lastv; i++ {
+			v[i*test.incv] = rand.Float64()
+		}
+		// Construct h explicitly to compare.
+		h := make([]float64, test.n*test.n)
+		for i := 0; i < test.n; i++ {
+			h[i*test.n+i] = 1
+		}
+		hMat := blas64.General{
+			Rows:   test.n,
+			Cols:   test.n,
+			Stride: test.n,
+			Data:   h,
+		}
+		vVec := blas64.Vector{
+			Inc:  test.incv,
+			Data: v,
+		}
+		blas64.Ger(-test.tau, vVec, vVec, hMat)
+
+		// Apply multiplication (2nd copy is to avoid aliasing).
+		cMat := blas64.General{
+			Rows:   test.m,
+			Cols:   test.n,
+			Stride: test.ldc,
+			Data:   cCopy,
+		}
+		cMat2 := blas64.General{
+			Rows:   test.m,
+			Cols:   test.n,
+			Stride: test.ldc,
+			Data:   cCopy2,
+		}
+		blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, cMat2, hMat, 0, cMat)
+
+		// cMat now stores the true answer. Compare with the function call.
+		work := make([]float64, sz)
+		impl.Dlarf(blas.Right, test.m, test.n, v, test.incv, test.tau, c, test.ldc, work)
+		if !floats.EqualApprox(c, cMat.Data, 1e-14) {
+			t.Errorf("Dlarf mismatch case %v. Want %v, got %v", i, cMat.Data, c)
+		}
+
+		// Test on the left side.
+		copy(c, cCopy2)
+		copy(cCopy, c)
+		// Construct h.
+		h = make([]float64, test.m*test.m)
+		for i := 0; i < test.m; i++ {
+			h[i*test.m+i] = 1
+		}
+		hMat = blas64.General{
+			Rows:   test.m,
+			Cols:   test.m,
+			Stride: test.m,
+			Data:   h,
+		}
+		blas64.Ger(-test.tau, vVec, vVec, hMat)
+		blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, hMat, cMat2, 0, cMat)
+		impl.Dlarf(blas.Left, test.m, test.n, v, test.incv, test.tau, c, test.ldc, work)
+		if !floats.EqualApprox(c, cMat.Data, 1e-14) {
+			t.Errorf("Dlarf mismatch case %v. Want %v, got %v", i, cMat.Data, c)
+		}
+	}
+}

testlapack/dpotf2.go

@@ -1,3 +1,7 @@
+// Copyright ©2015 The gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
 package testlapack
 
 import (
@@ -33,6 +37,17 @@ func Dpotf2Test(t *testing.T, impl Dpotf2er) {
 				{0, 0, 0, 3.401680257083044},
 			},
 		},
+		{
+			a: [][]float64{
+				{8, 2},
+				{2, 4},
+			},
+			pos: true,
+			U: [][]float64{
+				{2.82842712474619, 0.707106781186547},
+				{0, 1.870828693386971},
+			},
+		},
 	} {
 		testDpotf2(t, impl, test.pos, test.a, test.U, len(test.a[0]), blas.Upper)
 		testDpotf2(t, impl, test.pos, test.a, test.U, len(test.a[0])+5, blas.Upper)

testlapack/general.go

@@ -0,0 +1,12 @@
+// Copyright ©2015 The gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package testlapack
+
+func max(a, b int) int {
+	if a > b {
+		return a
+	}
+	return b
+}

testlapack/iladlc.go

@@ -0,0 +1,83 @@
+// Copyright ©2015 The gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package testlapack
+
+import "testing"
+
+type Iladlcer interface {
+	Iladlc(m, n int, a []float64, lda int) int
+}
+
+func IladlcTest(t *testing.T, impl Iladlcer) {
+	for i, test := range []struct {
+		a         []float64
+		m, n, lda int
+		ans       int
+	}{
+		{
+			a:   []float64{0, 0, 0, 0},
+			m:   1,
+			n:   1,
+			lda: 2,
+			ans: -1,
+		},
+		{
+			a:   []float64{0, 0, 0, 0},
+			m:   2,
+			n:   2,
+			lda: 2,
+			ans: -1,
+		},
+		{
+			a:   []float64{0, 0, 0, 0},
+			m:   4,
+			n:   1,
+			lda: 1,
+			ans: -1,
+		},
+		{
+			a:   []float64{0, 0, 0, 0},
+			m:   1,
+			n:   4,
+			lda: 4,
+			ans: -1,
+		},
+		{
+			a: []float64{
+				1, 2, 3, 4,
+				5, 6, 7, 8,
+			},
+			m:   2,
+			n:   4,
+			lda: 4,
+			ans: 3,
+		},
+		{
+			a: []float64{
+				1, 2, 3, 0,
+				0, 0, 0, 0,
+			},
+			m:   2,
+			n:   4,
+			lda: 4,
+			ans: 2,
+		},
+		{
+			a: []float64{
+				0, 0, 3, 4,
+				0, 0, 0, 0,
+			},
+			m:   2,
+			n:   2,
+			lda: 4,
+			ans: -1,
+		},
+	} {
+		ans := impl.Iladlc(test.m, test.n, test.a, test.lda)
+		if ans != test.ans {
+			t.Errorf("Column mismatch case %v. Want: %v, got: %v", i, test.ans, ans)
+		}
+	}
+}

testlapack/iladlr.go

@@ -0,0 +1,83 @@
+// Copyright ©2015 The gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package testlapack
+
+import "testing"
+
+type Iladlrer interface {
+	Iladlr(m, n int, a []float64, lda int) int
+}
+
+func IladlrTest(t *testing.T, impl Iladlrer) {
+	for i, test := range []struct {
+		a         []float64
+		m, n, lda int
+		ans       int
+	}{
+		{
+			a:   []float64{0, 0, 0, 0},
+			m:   1,
+			n:   1,
+			lda: 2,
+			ans: -1,
+		},
+		{
+			a:   []float64{0, 0, 0, 0},
+			m:   2,
+			n:   2,
+			lda: 2,
+			ans: -1,
+		},
+		{
+			a:   []float64{0, 0, 0, 0},
+			m:   4,
+			n:   1,
+			lda: 1,
+			ans: -1,
+		},
+		{
+			a:   []float64{0, 0, 0, 0},
+			m:   1,
+			n:   4,
+			lda: 4,
+			ans: -1,
+		},
+		{
+			a: []float64{
+				1, 2, 3, 4,
+				5, 6, 7, 8,
+			},
+			m:   2,
+			n:   4,
+			lda: 4,
+			ans: 1,
+		},
+		{
+			a: []float64{
+				1, 2, 3, 0,
+				0, 0, 0, 0,
+			},
+			m:   2,
+			n:   4,
+			lda: 4,
+			ans: 0,
+		},
+		{
+			a: []float64{
+				0, 0, 3, 4,
+				0, 0, 0, 0,
+			},
+			m:   2,
+			n:   2,
+			lda: 4,
+			ans: -1,
+		},
+	} {
+		ans := impl.Iladlr(test.m, test.n, test.a, test.lda)
+		if ans != test.ans {
+			t.Errorf("Column mismatch case %v. Want: %v, got: %v", i, test.ans, ans)
+		}
+	}
+}