package/gonum
1
0
Fork 0
mirror of https://github.com/gonum/gonum.git synced 2025-10-22 22:59:24 +08:00
Code Issues Actions Packages Projects Releases Wiki Activity

native: implement dlaqp2 and test

Browse Source
This commit is contained in:
kortschak
2017-01-20 12:25:51 +10:30
parent 238ef89b13
commit fadf409885
6 changed files with 250 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
cgo/lapack.go
View File

@@ -16,6 +16,7 @@ import (
// Copied from lapack/native. Keep in sync.
const (
absIncNotOne = "lapack: increment not one or negative one"
badAuxv = "lapack: auxv has insufficient length"
badD = "lapack: d has insufficient length"
badDecompUpdate = "lapack: bad decomp update"
badDiag = "lapack: bad diag"
@@ -34,6 +35,7 @@ const (
badK2 = "lapack: k2 out of range"
badKperm = "lapack: incorrect permutation length"
badLdA = "lapack: index of a out of range"
badNb = "lapack: nb out of range"
badNorm = "lapack: bad norm"
badPivot = "lapack: bad pivot"
badS = "lapack: s has insufficient length"
@@ -46,6 +48,8 @@ const (
badTauQ = "lapack: tauQ has insufficient length"
badTauP = "lapack: tauP has insufficient length"
badTrans = "lapack: bad trans"
badVn1 = "lapack: vn1 has insufficient length"
badVn2 = "lapack: vn2 has insufficient length"
badUplo = "lapack: illegal triangle"
badWork = "lapack: insufficient working memory"
badWorkStride = "lapack: insufficient working array stride"
@@ -60,6 +64,7 @@ const (
negZ = "lapack: negative z value"
nLT0 = "lapack: n < 0"
nLTM = "lapack: n < m"
offsetGTM = "lapack: offset > m"
shortWork = "lapack: working array shorter than declared"
zeroDiv = "lapack: zero divisor"
)

111
native/dlaqp2.go Normal file
View File

@@ -0,0 +1,111 @@
// Copyright ©2017 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 (
"math"
"github.com/gonum/blas"
"github.com/gonum/blas/blas64"
)
// Dlaqp2 computes a QR factorization with column pivoting of the block A[offset:m, 0:n]
// of the m×n matrix A. The block A[0:offset, 0:n] is accordingly pivoted, but not factorized.
//
// On exit, the upper triangle of block A[offset:m, 0:n] is the triangular factor obtained.
// The elements in block A[offset:m, 0:n] below the diagonal, together with tau, represent
// the orthogonal matrix Q as a product of elementary reflectors.
//
// offset is number of rows of the matrix A that must be pivoted but not factorized.
// offset must not be negative otherwise Dlaqp2 will panic.
//
// On exit, jpvt holds the permutation that was applied; the jth column of A*P was the
// jpvt[j] column of A. jpvt must have length n, otherwise Dlaqp2 will panic.
//
// On exit tau holds the scalar factors of the elementary reflectors. It must have length
// at least min(m-offset, n) otherwise Dlaqp2 will panic.
//
// vn1 and vn2 hold the partial and complete column norms respectively. They must have length n,
// otherwise Dlaqp2 will panic.
//
// work must have length n, otherwise Dlaqp2 will panic.
//
// Dlaqp2 is an internal routine. It is exported for testing purposes.
func (impl Implementation) Dlaqp2(m, n, offset int, a []float64, lda int, jpvt []int, tau, vn1, vn2, work []float64) {
checkMatrix(m, n, a, lda)
if len(jpvt) != n {
panic(badIpiv)
}
mn := min(m-offset, n)
if len(tau) < mn {
panic(badTau)
}
if len(vn1) < n {
panic(badVn1)
}
if len(vn2) < n {
panic(badVn2)
}
if len(work) < n {
panic(badWork)
}
tol3z := math.Sqrt(dlamchE)
bi := blas64.Implementation()
// Compute factorization.
for i := 0; i < mn; i++ {
offpi := offset + i
// Determine ith pivot column and swap if necessary.
p := i + bi.Idamax(n-i, vn1[i:], 1)
if p != i {
bi.Dswap(m, a[p:], lda, a[i:], lda)
jpvt[p], jpvt[i] = jpvt[i], jpvt[p]
vn1[p] = vn1[i]
vn2[p] = vn2[i]
}
// Generate elementary reflector H_i.
if offpi < m-1 {
a[offpi*lda+i], tau[i] = impl.Dlarfg(m-offpi, a[offpi*lda+i], a[(offpi+1)*lda+i:], lda)
} else {
tau[i] = 0
}
if i < n-1 {
// Apply H_i^T to A[offset+i:m, i:n] from the left.
aii := a[offpi*lda+i]
a[offpi*lda+i] = 1
impl.Dlarf(blas.Left, m-offpi, n-i-1, a[offpi*lda+i:], lda, tau[i], a[offpi*lda+i+1:], lda, work)
a[offpi*lda+i] = aii
}
// Update partial column norms.
for j := i + 1; j < n; j++ {
if vn1[j] == 0 {
continue
}
// The following marked lines follow from the
// analysis in Lapack Working Note 176.
r := math.Abs(a[offpi*lda+j]) / vn1[j] // *
temp := math.Max(0, 1-r*r) // *
r = vn1[j] / vn2[j] // *
temp2 := temp * r * r // *
if temp2 < tol3z {
var v float64
if offpi < m-1 {
v = bi.Dnrm2(m-offpi-1, a[(offpi+1)*lda+j:], lda)
}
vn1[j] = v
vn2[j] = v
} else {
vn1[j] *= math.Sqrt(temp) // *
}
}
}
}

5
native/general.go
View File

@@ -20,6 +20,7 @@ var _ lapack.Float64 = Implementation{}
// This list is duplicated in lapack/cgo. Keep in sync.
const (
absIncNotOne = "lapack: increment not one or negative one"
badAuxv = "lapack: auxv has insufficient length"
badD = "lapack: d has insufficient length"
badDecompUpdate = "lapack: bad decomp update"
badDiag = "lapack: bad diag"
@@ -38,6 +39,7 @@ const (
badK2 = "lapack: k2 out of range"
badKperm = "lapack: incorrect permutation length"
badLdA = "lapack: index of a out of range"
badNb = "lapack: nb out of range"
badNorm = "lapack: bad norm"
badPivot = "lapack: bad pivot"
badS = "lapack: s has insufficient length"
@@ -50,6 +52,8 @@ const (
badTauQ = "lapack: tauQ has insufficient length"
badTauP = "lapack: tauP has insufficient length"
badTrans = "lapack: bad trans"
badVn1 = "lapack: vn1 has insufficient length"
badVn2 = "lapack: vn2 has insufficient length"
badUplo = "lapack: illegal triangle"
badWork = "lapack: insufficient working memory"
badWorkStride = "lapack: insufficient working array stride"
@@ -64,6 +68,7 @@ const (
negZ = "lapack: negative z value"
nLT0 = "lapack: n < 0"
nLTM = "lapack: n < m"
offsetGTM = "lapack: offset > m"
shortWork = "lapack: working array shorter than declared"
zeroDiv = "lapack: zero divisor"
)

4
native/lapack_test.go
View File

@@ -192,6 +192,10 @@ func TestDlaqr04(t *testing.T) {
testlapack.Dlaqr04Test(t, impl)
}
func TestDlaqp2(t *testing.T) {
testlapack.Dlaqp2Test(t, impl)
}
func TestDlaqr1(t *testing.T) {
testlapack.Dlaqr1Test(t, impl)
}

115
testlapack/dlaqp2.go Normal file
View File

@@ -0,0 +1,115 @@
// Copyright ©2017 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 (
"fmt"
"math"
"testing"
"github.com/gonum/blas"
"github.com/gonum/blas/blas64"
)
type Dlaqp2er interface {
Dlapmter
Dlaqp2(m, n, offset int, a []float64, lda int, jpvt []int, tau, vn1, vn2, work []float64)
}
func Dlaqp2Test(t *testing.T, impl Dlaqp2er) {
for ti, test := range []struct {
m, n, offset int
}{
{m: 4, n: 3, offset: 0},
{m: 4, n: 3, offset: 2},
{m: 4, n: 3, offset: 4},
{m: 3, n: 4, offset: 0},
{m: 3, n: 4, offset: 1},
{m: 3, n: 4, offset: 2},
{m: 8, n: 3, offset: 0},
{m: 8, n: 3, offset: 4},
{m: 8, n: 3, offset: 8},
{m: 3, n: 8, offset: 0},
{m: 3, n: 8, offset: 1},
{m: 3, n: 8, offset: 2},
{m: 10, n: 10, offset: 0},
{m: 10, n: 10, offset: 5},
{m: 10, n: 10, offset: 10},
} {
m := test.m
n := test.n
jpiv := make([]int, n)
for _, extra := range []int{0, 11} {
a := zeros(m, n, n+extra)
c := 1
for i := 0; i < m; i++ {
for j := 0; j < n; j++ {
a.Data[i*a.Stride+j] = float64(c)
c++
}
}
aCopy := cloneGeneral(a)
for j := range jpiv {
jpiv[j] = j
}
tau := make([]float64, n)
vn1 := columnNorms(m, n, a.Data, a.Stride)
vn2 := columnNorms(m, n, a.Data, a.Stride)
work := make([]float64, n)
impl.Dlaqp2(m, n, test.offset, a.Data, a.Stride, jpiv, tau, vn1, vn2, work)
prefix := fmt.Sprintf("Case %v (offset=%t,m=%v,n=%v,extra=%v)", ti, test.offset, m, n, extra)
if !generalOutsideAllNaN(a) {
t.Errorf("%v: out-of-range write to A", prefix)
}
if test.offset == m {
continue
}
mo := m - test.offset
q := constructQ("QR", mo, n, a.Data[test.offset*a.Stride:], a.Stride, tau)
// Check that q is orthonormal
for i := 0; i < mo; i++ {
nrm := blas64.Nrm2(mo, blas64.Vector{Inc: 1, Data: q.Data[i*mo:]})
if math.Abs(nrm-1) > 1e-13 {
t.Errorf("Case %v, q not normal", ti)
}
for j := 0; j < i; j++ {
dot := blas64.Dot(mo, blas64.Vector{Inc: 1, Data: q.Data[i*mo:]}, blas64.Vector{Inc: 1, Data: q.Data[j*mo:]})
if math.Abs(dot) > 1e-14 {
t.Errorf("Case %v, q not orthogonal", ti)
}
}
}
// Check that A * P = Q * R
r := blas64.General{
Rows: mo,
Cols: n,
Stride: n,
Data: make([]float64, mo*n),
}
for i := 0; i < mo; i++ {
for j := i; j < n; j++ {
r.Data[i*n+j] = a.Data[(test.offset+i)*a.Stride+j]
}
}
got := nanGeneral(mo, n, n)
blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, q, r, 0, got)
want := aCopy
impl.Dlapmt(true, want.Rows, want.Cols, want.Data, want.Stride, jpiv)
want.Rows = mo
want.Data = want.Data[test.offset*want.Stride:]
if !equalApproxGeneral(got, want, 1e-12) {
t.Errorf("Case %v, Q*R != A*P\nQ*R=%v\nA*P=%v", ti, got, want)
}
}
}
}

10
testlapack/general.go
View File

@@ -262,6 +262,16 @@ func transposeGeneral(a blas64.General) blas64.General {
return ans
}
// columnNorms returns the column norms of a.
func columnNorms(m, n int, a []float64, lda int) []float64 {
bi := blas64.Implementation()
norms := make([]float64, n)
for j := 0; j < n; j++ {
norms[j] = bi.Dnrm2(m, a[j:], lda)
}
return norms
}
// extractVMat collects the single reflectors from a into a matrix.
func extractVMat(m, n int, a []float64, lda int, direct lapack.Direct, store lapack.StoreV) blas64.General {
k := min(m, n)