Skip to content

Commit

Permalink
Merge 8c9bfb0 into 17a33f2
Browse files Browse the repository at this point in the history
  • Loading branch information
@vladimir-ch
vladimir-ch committed Nov 23, 2018
2 parents 17a33f2 + 8c9bfb0 commit 9f498c8
Show file tree
Hide file tree
Showing 4 changed files with 100 additions and 41 deletions.
78 changes: 52 additions & 26 deletions lapack/testlapack/dlatrd.go
View file
Expand Up @@ -42,11 +42,15 @@ func DlatrdTest(t *testing.T, impl Dlatrder) {
ldw = nb
}

// Allocate n×n matrix A and fill it with random numbers.
a := make([]float64, n*lda)
for i := range a {
a[i] = rnd.NormFloat64()
}

// Allocate output slices and matrix W and fill them
// with NaN. All their elements should be overwritten by
// Dlatrd.
e := make([]float64, n-1)
for i := range e {
e[i] = math.NaN()
Expand All @@ -63,30 +67,30 @@ func DlatrdTest(t *testing.T, impl Dlatrder) {
aCopy := make([]float64, len(a))
copy(aCopy, a)

// Reduce nb rows and columns of the symmetric matrix A
// defined by uplo triangle to symmetric tridiagonal
// form.
impl.Dlatrd(uplo, n, nb, a, lda, e, tau, w, ldw)

// Construct Q.
ldq := n
// Construct Q from elementary reflectors stored in
// columns of A.
q := blas64.General{
Rows: n,
Cols: n,
Stride: ldq,
Data: make([]float64, n*ldq),
Stride: n,
Data: make([]float64, n*n),
}
// Initialize Q to the identity matrix.
for i := 0; i < n; i++ {
q.Data[i*ldq+i] = 1
q.Data[i*q.Stride+i] = 1
}
if uplo == blas.Upper {
for i := n - 1; i >= n-nb; i-- {
if i == 0 {
continue
}
h := blas64.General{
Rows: n, Cols: n, Stride: n, Data: make([]float64, n*n),
}
for j := 0; j < n; j++ {
h.Data[j*n+j] = 1
}

// Extract the elementary reflector v from A.
v := blas64.Vector{
Inc: 1,
Data: make([]float64, n),
Expand All @@ -96,8 +100,16 @@ func DlatrdTest(t *testing.T, impl Dlatrder) {
}
v.Data[i-1] = 1

// Compute H = I - tau[i-1] * v * v^T.
h := blas64.General{
Rows: n, Cols: n, Stride: n, Data: make([]float64, n*n),
}
for j := 0; j < n; j++ {
h.Data[j*n+j] = 1
}
blas64.Ger(-tau[i-1], v, v, h)

// Update Q <- Q * H.
qTmp := blas64.General{
Rows: n, Cols: n, Stride: n, Data: make([]float64, n*n),
}
Expand All @@ -109,12 +121,8 @@ func DlatrdTest(t *testing.T, impl Dlatrder) {
if i == n-1 {
continue
}
h := blas64.General{
Rows: n, Cols: n, Stride: n, Data: make([]float64, n*n),
}
for j := 0; j < n; j++ {
h.Data[j*n+j] = 1
}

// Extract the elementary reflector v from A.
v := blas64.Vector{
Inc: 1,
Data: make([]float64, n),
Expand All @@ -123,8 +131,17 @@ func DlatrdTest(t *testing.T, impl Dlatrder) {
for j := i + 2; j < n; j++ {
v.Data[j] = a[j*lda+i]
}

// Compute H = I - tau[i] * v * v^T.
h := blas64.General{
Rows: n, Cols: n, Stride: n, Data: make([]float64, n*n),
}
for j := 0; j < n; j++ {
h.Data[j*n+j] = 1
}
blas64.Ger(-tau[i], v, v, h)

// Update Q <- Q * H.
qTmp := blas64.General{
Rows: n, Cols: n, Stride: n, Data: make([]float64, n*n),
}
Expand All @@ -134,56 +151,61 @@ func DlatrdTest(t *testing.T, impl Dlatrder) {
}
errStr := fmt.Sprintf("isUpper = %v, n = %v, nb = %v", uplo == blas.Upper, n, nb)
if !isOrthogonal(q) {
t.Errorf("Q not orthogonal. %s", errStr)
t.Errorf("Case %v: Q not orthogonal", errStr)
}
aGen := genFromSym(blas64.Symmetric{N: n, Stride: lda, Uplo: uplo, Data: aCopy})
if !dlatrdCheckDecomposition(t, uplo, n, nb, e, tau, a, lda, aGen, q) {
t.Errorf("Decomposition mismatch. %s", errStr)
if !dlatrdCheckDecomposition(t, uplo, n, nb, e, a, lda, aGen, q) {
t.Errorf("Case %v: Decomposition mismatch", errStr)
}
}
}
}

// dlatrdCheckDecomposition checks that the first nb rows have been successfully
// reduced.
func dlatrdCheckDecomposition(t *testing.T, uplo blas.Uplo, n, nb int, e, tau, a []float64, lda int, aGen, q blas64.General) bool {
// Compute Q^T * A * Q.
func dlatrdCheckDecomposition(t *testing.T, uplo blas.Uplo, n, nb int, e, a []float64, lda int, aGen, q blas64.General) bool {
// Compute ans = Q^T * A * Q.
// ans should be a tridiagonal matrix in the first or last nb rows and
// columns, depending on uplo.
tmp := blas64.General{
Rows: n,
Cols: n,
Stride: n,
Data: make([]float64, n*n),
}

ans := blas64.General{
Rows: n,
Cols: n,
Stride: n,
Data: make([]float64, n*n),
}

blas64.Gemm(blas.Trans, blas.NoTrans, 1, q, aGen, 0, tmp)
blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, tmp, q, 0, ans)

// Compare with T.
// Compare the output of Dlatrd (stored in a and e) with the explicit
// reduction to tridiagonal matrix Q^T * A * Q (stored in ans).
if uplo == blas.Upper {
for i := n - 1; i >= n-nb; i-- {
for i := n - nb; i < n; i++ {
for j := 0; j < n; j++ {
v := ans.Data[i*ans.Stride+j]
switch {
case i == j:
// Diagonal elements of a and ans should match.
if math.Abs(v-a[i*lda+j]) > 1e-10 {
return false
}
case i == j-1:
// Superdiagonal elements in a should be 1.
if math.Abs(a[i*lda+j]-1) > 1e-10 {
return false
}
// Superdiagonal elements of ans should match e.
if math.Abs(v-e[i]) > 1e-10 {
return false
}
case i == j+1:
default:
// All other elements should be 0.
if math.Abs(v) > 1e-10 {
return false
}
Expand All @@ -196,18 +218,22 @@ func dlatrdCheckDecomposition(t *testing.T, uplo blas.Uplo, n, nb int, e, tau, a
v := ans.Data[i*ans.Stride+j]
switch {
case i == j:
// Diagonal elements of a and ans should match.
if math.Abs(v-a[i*lda+j]) > 1e-10 {
return false
}
case i == j-1:
case i == j+1:
// Subdiagonal elements in a should be 1.
if math.Abs(a[i*lda+j]-1) > 1e-10 {
return false
}
// Subdiagonal elements of ans should match e.
if math.Abs(v-e[i-1]) > 1e-10 {
return false
}
default:
// All other elements should be 0.
if math.Abs(v) > 1e-10 {
return false
}
Expand Down
25 changes: 12 additions & 13 deletions lapack/testlapack/dorg2r.go
View file
Expand Up @@ -5,7 +5,6 @@
package testlapack

import (
"fmt"
"testing"

"golang.org/x/exp/rand"
Expand All @@ -20,7 +19,7 @@ type Dorg2rer interface {

func Dorg2rTest(t *testing.T, impl Dorg2rer) {
rnd := rand.New(rand.NewSource(1))
for _, test := range []struct {
for ti, test := range []struct {
m, n, k, lda int
}{
{3, 3, 0, 0},
Expand All @@ -39,42 +38,42 @@ func Dorg2rTest(t *testing.T, impl Dorg2rer) {
if lda == 0 {
lda = test.n
}
// Allocate m×n matrix A and fill it with random numbers.
a := make([]float64, m*lda)
for i := range a {
a[i] = rnd.NormFloat64()
}
k := min(m, n)
tau := make([]float64, k)

// Compute the QR decomposition of A.
tau := make([]float64, min(m, n))
work := make([]float64, 1)
impl.Dgeqrf(m, n, a, lda, tau, work, -1)
work = make([]float64, int(work[0]))
impl.Dgeqrf(m, n, a, lda, tau, work, len(work))

k = test.k
// Compute the matrix Q explicitly using the first k elementary reflectors.
k := test.k
if k == 0 {
k = n
}
q := constructQK("QR", m, n, k, a, lda, tau)

// Compute the matrix Q using Dorg2r.
impl.Dorg2r(m, n, k, a, lda, tau, work)

// Check that the first n columns match.
// Check that the first n columns of both results match.
same := true
loop:
for i := 0; i < m; i++ {
for j := 0; j < n; j++ {
if !floats.EqualWithinAbsOrRel(q.Data[i*q.Stride+j], a[i*lda+j], 1e-12, 1e-12) {
same = false
break
break loop
}
}
}
if !same {
fmt.Println()
fmt.Println("a =")
printRowise(a, m, n, lda, false)
fmt.Println("q =")
printRowise(q.Data, q.Rows, q.Cols, q.Stride, false)
t.Errorf("Q mismatch")
t.Errorf("Case %v: Q mismatch", ti)
}
}
}
28 changes: 26 additions & 2 deletions lapack/testlapack/dorgtr.go
View file
Expand Up @@ -45,21 +45,32 @@ func DorgtrTest(t *testing.T, impl Dorgtrer) {
if lda == 0 {
lda = n
}
// Allocate n×n matrix A and fill it with random numbers.
a := make([]float64, n*lda)
for i := range a {
a[i] = rnd.NormFloat64()
}
aCopy := make([]float64, len(a))
copy(aCopy, a)

// Allocate slices for the main diagonal and the
// first off-diagonal of the tri-diagonal matrix.
d := make([]float64, n)
e := make([]float64, n-1)
// Allocate slice for elementary reflector scales.
tau := make([]float64, n-1)

// Compute optimum workspace size for Dorgtr call.
work := make([]float64, 1)
impl.Dsytrd(uplo, n, a, lda, d, e, tau, work, -1)
work = make([]float64, int(work[0]))

// Compute elementary reflectors that reduce the
// symmetric matrix defined by the uplo triangle
// of A to a tridiagonal matrix.
impl.Dsytrd(uplo, n, a, lda, d, e, tau, work, len(work))

// Compute workspace size for Dorgtr call.
var lwork int
switch wl {
case minimumWork:
Expand All @@ -76,14 +87,23 @@ func DorgtrTest(t *testing.T, impl Dorgtrer) {
}
work = nanSlice(lwork)

// Generate an orthogonal matrix Q that reduces
// the uplo triangle of A to a tridiagonal matrix.
impl.Dorgtr(uplo, n, a, lda, tau, work, len(work))

q := blas64.General{
Rows: n,
Cols: n,
Stride: lda,
Data: a,
}

if !isOrthogonal(q) {
t.Errorf("Case uplo=%v,n=%v: Q is not orthogonal", uplo, n)
continue
}

// Create the tridiagonal matrix explicitly in
// dense representation from the diagonals d and e.
tri := blas64.General{
Rows: n,
Cols: n,
Expand All @@ -98,6 +118,8 @@ func DorgtrTest(t *testing.T, impl Dorgtrer) {
}
}

// Create the symmetric matrix A from the uplo
// triangle of aCopy, storing it explicitly in dense form.
aMat := blas64.General{
Rows: n,
Cols: n,
Expand All @@ -122,12 +144,14 @@ func DorgtrTest(t *testing.T, impl Dorgtrer) {
}
}

// Compute Q^T * A * Q and store the result in ans.
tmp := blas64.General{Rows: n, Cols: n, Stride: n, Data: make([]float64, n*n)}
blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, aMat, q, 0, tmp)

ans := blas64.General{Rows: n, Cols: n, Stride: n, Data: make([]float64, n*n)}
blas64.Gemm(blas.Trans, blas.NoTrans, 1, q, tmp, 0, ans)

// Compare the tridiagonal matrix tri from
// Dorgtr with the explicit computation ans.
if !floats.EqualApprox(ans.Data, tri.Data, 1e-13) {
t.Errorf("Recombination mismatch. n = %v, isUpper = %v", n, uplo == blas.Upper)
}
Expand Down
10 changes: 10 additions & 0 deletions lapack/testlapack/dtrti2.go
View file
Expand Up @@ -109,13 +109,20 @@ func Dtrti2Test(t *testing.T, impl Dtrti2er) {
if lda == 0 {
lda = n
}
// Allocate n×n matrix A and fill it with random numbers.
a := make([]float64, n*lda)
for i := range a {
a[i] = rnd.Float64()
}
for i := 0; i < n; i++ {
// This keeps the matrices well conditioned.
a[i*lda+i] += float64(n)
}
aCopy := make([]float64, len(a))
copy(aCopy, a)
// Compute the inverse of the uplo triangle.
impl.Dtrti2(uplo, diag, n, a, lda)
// Zero out the opposite triangle.
if uplo == blas.Upper {
for i := 1; i < n; i++ {
for j := 0; j < i; j++ {
Expand All @@ -132,13 +139,16 @@ func Dtrti2Test(t *testing.T, impl Dtrti2er) {
}
}
if diag == blas.Unit {
// Set the diagonal of A^{-1} and A explicitly to 1.
for i := 0; i < n; i++ {
a[i*lda+i] = 1
aCopy[i*lda+i] = 1
}
}
// Compute A^{-1} * A and store the result in ans.
ans := make([]float64, len(a))
bi.Dgemm(blas.NoTrans, blas.NoTrans, n, n, n, 1, a, lda, aCopy, lda, 0, ans, lda)
// Check that ans is the identity matrix.
iseye := true
for i := 0; i < n; i++ {
for j := 0; j < n; j++ {
Expand Down

0 comments on commit 9f498c8

Please sign in to comment.