diff --git a/lapack/lapacke/internal/conv/conv.go b/lapack/lapacke/internal/conv/conv.go deleted file mode 100644 index 33247990..00000000 --- a/lapack/lapacke/internal/conv/conv.go +++ /dev/null @@ -1,19 +0,0 @@ -// Copyright ©2019 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 conv - -func min(a, b int) int { - if a < b { - return a - } - return b -} - -func max(a, b int) int { - if a > b { - return a - } - return b -} diff --git a/lapack/lapacke/internal/conv/pb_test.go b/lapack/lapacke/internal/conv/pb_test.go deleted file mode 100644 index a6d1c2d8..00000000 --- a/lapack/lapacke/internal/conv/pb_test.go +++ /dev/null @@ -1,123 +0,0 @@ -// Copyright ©2019 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 conv - -import ( - "testing" - - "golang.org/x/exp/rand" - - "gonum.org/v1/gonum/floats" -) - -func TestDpb(t *testing.T) { - for ti, test := range []struct { - uplo byte - n, kd int - a, b []float64 - }{ - { - uplo: 'U', - n: 6, - kd: 2, - a: []float64{ - 1, 2, 3, // 1. row - 4, 5, 6, - 7, 8, 9, - 10, 11, 12, - 13, 14, -1, - 15, -1, -1, // 6. row - }, - b: []float64{ - -1, -1, 1, // 1. column - -1, 2, 4, - 3, 5, 7, - 6, 8, 10, - 9, 11, 13, - 12, 14, 15, // 6. column - }, - }, - { - uplo: 'L', - n: 6, - kd: 2, - a: []float64{ - -1, -1, 1, // 1. row - -1, 2, 3, - 4, 5, 6, - 7, 8, 9, - 10, 11, 12, - 13, 14, 15, // 6. row - }, - b: []float64{ - 1, 2, 4, // 1. column - 3, 5, 7, - 6, 8, 10, - 9, 11, 13, - 12, 14, -1, - 15, -1, -1, // 6. column - }, - }, - } { - uplo := test.uplo - n := test.n - kd := test.kd - lda := kd + 1 - - a := make([]float64, len(test.a)) - copy(a, test.a) - - b := make([]float64, len(test.b)) - copy(b, test.b) - - got := make([]float64, len(a)) - for i := range got { - got[i] = -1 - } - DpbToColMajor(uplo, n, kd, a, lda, got, lda) - if !floats.Equal(test.b, got) { - t.Errorf("Case %v (uplo=%v,n=%v,kd=%v): unexpected conversion to column-major;\ngot %v\nwant %v", - ti, string(uplo), n, kd, got, test.b) - } - - for i := range got { - got[i] = -1 - } - DpbToRowMajor(uplo, n, kd, b, lda, got, lda) - if !floats.Equal(test.a, got) { - t.Errorf("Case %v (uplo=%v,n=%v,kd=%v): unexpected conversion to row-major;\ngot %v\nwant %v", - ti, string(uplo), n, kd, got, test.b) - } - } - - rnd := rand.New(rand.NewSource(1)) - for _, n := range []int{0, 1, 2, 3, 4, 5, 10} { - for _, kd := range []int{0, (n + 1) / 4, (3*n - 1) / 4, (5*n + 1) / 4} { - for _, uplo := range []byte{'U', 'L'} { - for _, lda := range []int{kd + 1, kd + 1 + 7} { - a := make([]float64, n*lda) - for i := range a { - a[i] = rnd.NormFloat64() - } - aCopy := make([]float64, len(a)) - copy(aCopy, a) - - ldb := lda - b := make([]float64, ldb*n) - for i := range b { - b[i] = rnd.NormFloat64() - } - - DpbToColMajor(uplo, n, kd, a, lda, b, ldb) - DpbToRowMajor(uplo, n, kd, b, ldb, a, lda) - - if !floats.Equal(a, aCopy) { - t.Errorf("uplo=%v,n=%v,kd=%v: conversion does not roundtrip", string(uplo), n, kd) - } - } - } - } - } -} diff --git a/lapack/lapacke/internal/conv/pb.go b/lapack/netlib/conv.go similarity index 58% rename from lapack/lapacke/internal/conv/pb.go rename to lapack/netlib/conv.go index 16fd13be..303e938f 100644 --- a/lapack/lapacke/internal/conv/pb.go +++ b/lapack/netlib/conv.go @@ -2,13 +2,14 @@ // Use of this source code is governed by a BSD-style // license that can be found in the LICENSE file. -package conv +package netlib -// DpbToColMajor converts a symmetric band matrix A in CBLAS row-major layout -// to FORTRAN column-major layout and stores the result in B. +import "gonum.org/v1/gonum/blas" + +// convDpbToLapacke converts a symmetric band matrix A in CBLAS row-major layout +// to LAPACKE row-major layout and stores the result in B. // -// For example, when n = 6, kd = 2 and uplo == 'U', DpbToColMajor -// converts +// For example, when n = 6, kd = 2 and uplo == 'U', convDpbToLapacke converts // A = a00 a01 a02 // a11 a12 a13 // a22 a23 a24 @@ -22,9 +23,9 @@ package conv // * a01 a12 a23 a34 a45 // a00 a11 a22 a33 a44 a55 // stored in a slice as -// b = [* * a00 * a01 a11 a02 a12 a22 a13 a23 a33 a24 a34 a44 a35 a45 a55] +// b = [* * a02 a13 a24 a35 * a01 a12 a23 a34 a45 a00 a11 a22 a33 a44 a55] // -// When n = 6, kd = 2 and uplo == 'L', DpbToColMajor converts +// When n = 6, kd = 2 and uplo == 'L', convDpbToLapacke converts // A = * * a00 // * a10 a11 // a20 a21 a22 @@ -38,43 +39,43 @@ package conv // a10 a21 a32 a43 a54 * // a20 a31 a42 a53 * * // stored in a slice as -// b = [a00 a10 a20 a11 a21 a31 a22 a32 a42 a33 a43 a53 a44 a54 * a55 * *] +// b = [a00 a11 a22 a33 a44 a55 a10 a21 a32 a43 a54 * a20 a31 a42 a53 * * ] // // In these example elements marked as * are not referenced. -func DpbToColMajor(uplo byte, n, kd int, a []float64, lda int, b []float64, ldb int) { - if uplo == 'U' { +func convDpbToLapacke(uplo blas.Uplo, n, kd int, a []float64, lda int, b []float64, ldb int) { + if uplo == blas.Upper { for i := 0; i < n; i++ { for jb := 0; jb < min(n-i, kd+1); jb++ { j := i + jb // Column index in the full matrix - b[kd-jb+j*ldb] = a[i*lda+jb] + b[(kd-jb)*ldb+j] = a[i*lda+jb] } } } else { for i := 0; i < n; i++ { for jb := max(0, kd-i); jb < kd+1; jb++ { j := i - kd + jb // Column index in the full matrix - b[kd-jb+j*ldb] = a[i*lda+jb] + b[(kd-jb)*ldb+j] = a[i*lda+jb] } } } } -// DpbToRowMajor converts a symmetric band matrix A in FORTRAN column-major -// layout to CBLAS row-major layout and stores the result in B. In other words, -// it performs the inverse conversion to DpbToColMajor. -func DpbToRowMajor(uplo byte, n, kd int, a []float64, lda int, b []float64, ldb int) { - if uplo == 'U' { +// convDpbToGonum converts a symmetric band matrix A in LAPACKE row-major layout +// to CBLAS row-major layout and stores the result in B. In other words, it +// performs the inverse conversion to convDpbToLapacke. +func convDpbToGonum(uplo blas.Uplo, n, kd int, a []float64, lda int, b []float64, ldb int) { + if uplo == blas.Upper { for j := 0; j < n; j++ { for ib := max(0, kd-j); ib < kd+1; ib++ { i := j - kd + ib // Row index in the full matrix - b[i*ldb+kd-ib] = a[ib+j*lda] + b[i*ldb+kd-ib] = a[ib*lda+j] } } } else { for j := 0; j < n; j++ { for ib := 0; ib < min(n-j, kd+1); ib++ { i := j + ib // Row index in the full matrix - b[i*ldb+kd-ib] = a[ib+j*lda] + b[i*ldb+kd-ib] = a[ib*lda+j] } } } diff --git a/lapack/netlib/conv_test.go b/lapack/netlib/conv_test.go new file mode 100644 index 00000000..d50fa71d --- /dev/null +++ b/lapack/netlib/conv_test.go @@ -0,0 +1,131 @@ +// Copyright ©2019 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 netlib + +import ( + "fmt" + "testing" + + "golang.org/x/exp/rand" + + "gonum.org/v1/gonum/blas" + "gonum.org/v1/gonum/floats" +) + +func TestConvDpb(t *testing.T) { + for ti, test := range []struct { + uplo blas.Uplo + n, kd int + a, b []float64 + }{ + { + uplo: blas.Upper, + n: 6, + kd: 2, + a: []float64{ + 1, 2, 3, // 1. row + 4, 5, 6, + 7, 8, 9, + 10, 11, 12, + 13, 14, -1, + 15, -1, -1, // 6. row + }, + b: []float64{ + -1, -1, 3, 6, 9, 12, // 2. super-diagonal + -1, 2, 5, 8, 11, 14, + 1, 4, 7, 10, 13, 15, // main diagonal + }, + }, + { + uplo: blas.Lower, + n: 6, + kd: 2, + a: []float64{ + -1, -1, 1, // 1. row + -1, 2, 3, + 4, 5, 6, + 7, 8, 9, + 10, 11, 12, + 13, 14, 15, // 6. row + }, + b: []float64{ + 1, 3, 6, 9, 12, 15, // main diagonal + 2, 5, 8, 11, 14, -1, + 4, 7, 10, 13, -1, -1, // 2. sub-diagonal + }, + }, + } { + uplo := test.uplo + n := test.n + kd := test.kd + name := fmt.Sprintf("Case %v (uplo=%c,n=%v,kd=%v)", ti, uplo, n, kd) + + a := make([]float64, len(test.a)) + copy(a, test.a) + lda := kd + 1 + + got := make([]float64, len(test.b)) + for i := range got { + got[i] = -1 + } + ldb := max(1, n) + + convDpbToLapacke(uplo, n, kd, a, lda, got, ldb) + if !floats.Equal(test.a, a) { + t.Errorf("%v: unexpected modification of A in conversion to LAPACKE row-major", name) + } + if !floats.Equal(test.b, got) { + t.Errorf("%v: unexpected conversion to LAPACKE row-major;\ngot %v\nwant %v", name, got, test.b) + } + + b := make([]float64, len(test.b)) + copy(b, test.b) + + got = make([]float64, len(test.a)) + for i := range got { + got[i] = -1 + } + + convDpbToGonum(uplo, n, kd, b, ldb, got, lda) + if !floats.Equal(test.b, b) { + t.Errorf("%v: unexpected modification of B in conversion to Gonum row-major", name) + } + if !floats.Equal(test.a, got) { + t.Errorf("%v: unexpected conversion to Gonum row-major;\ngot %v\nwant %v", name, got, test.b) + } + } + + rnd := rand.New(rand.NewSource(1)) + for _, n := range []int{0, 1, 2, 3, 4, 5, 10} { + for _, kd := range []int{0, (n + 1) / 4, (3*n - 1) / 4, (5*n + 1) / 4} { + for _, uplo := range []blas.Uplo{blas.Upper, blas.Lower} { + for _, ldextra := range []int{0, 3} { + name := fmt.Sprintf("uplo=%c,n=%v,kd=%v", uplo, n, kd) + + lda := kd + 1 + ldextra + a := make([]float64, n*lda) + for i := range a { + a[i] = rnd.NormFloat64() + } + aCopy := make([]float64, len(a)) + copy(aCopy, a) + + ldb := max(1, n) + ldextra + b := make([]float64, (kd+1)*ldb) + for i := range b { + b[i] = rnd.NormFloat64() + } + + convDpbToLapacke(uplo, n, kd, a, lda, b, ldb) + convDpbToGonum(uplo, n, kd, b, ldb, a, lda) + + if !floats.Equal(a, aCopy) { + t.Errorf("%v: conversion does not roundtrip", name) + } + } + } + } + } +} diff --git a/lapack/netlib/lapack.go b/lapack/netlib/lapack.go index e513a836..769799fe 100644 --- a/lapack/netlib/lapack.go +++ b/lapack/netlib/lapack.go @@ -794,6 +794,105 @@ func (impl Implementation) Dlaswp(n int, a []float64, lda, k1, k2 int, ipiv []in lapacke.Dlaswp(n, a, lda, k1+1, k2+1, ipiv32, incX) } +// Dpbtrf computes the Cholesky factorization of an n×n symmetric positive +// definite band matrix +// A = U^T * U if uplo == blas.Upper +// A = L * L^T if uplo == blas.Lower +// where U is an upper triangular band matrix and L is lower triangular. kd is +// the number of super- or sub-diagonals of A. +// +// The band storage scheme is illustrated below when n = 6 and kd = 2. Elements +// marked * are not used by the function. +// +// uplo == blas.Upper +// On entry: On return: +// a00 a01 a02 u00 u01 u02 +// a11 a12 a13 u11 u12 u13 +// a22 a23 a24 u22 u23 u24 +// a33 a34 a35 u33 u34 u35 +// a44 a45 * u44 u45 * +// a55 * * u55 * * +// +// uplo == blas.Lower +// On entry: On return: +// * * a00 * * l00 +// * a10 a11 * l10 l11 +// a20 a21 a22 l20 l21 l22 +// a31 a32 a33 l31 l32 l33 +// a42 a43 a44 l42 l43 l44 +// a53 a54 a55 l53 l54 l55 +func (impl Implementation) Dpbtrf(uplo blas.Uplo, n, kd int, ab []float64, ldab int) (ok bool) { + switch { + case uplo != blas.Upper && uplo != blas.Lower: + panic(badUplo) + case n < 0: + panic(nLT0) + case kd < 0: + panic(kdLT0) + case ldab < kd+1: + panic(badLdA) + } + + // Quick return if possible. + if n == 0 { + return true + } + + if len(ab) < (n-1)*ldab+kd+1 { + panic(shortAB) + } + + ldabConv := n + abConv := make([]float64, (kd+1)*ldabConv) + convDpbToLapacke(uplo, n, kd, ab, ldab, abConv, ldabConv) + info := lapacke.Dpbtrf(byte(uplo), n, kd, abConv, ldabConv) + convDpbToGonum(uplo, n, kd, abConv, ldabConv, ab, ldab) + return info +} + +// Dpbtrs solves a system of linear equations A*X = B with an n×n symmetric +// positive definite band matrix A using the Cholesky factorization +// A = U^T * U if uplo == blas.Upper +// A = L * L^T if uplo == blas.Lower +// computed by Dpbtrf. kd is the number of super- or sub-diagonals of A. See the +// documentation for Dpbtrf for a description of the band storage format of A. +// +// On entry, b contains the n×nrhs right hand side matrix B. On return, it is +// overwritten with the solution matrix X. +func (Implementation) Dpbtrs(uplo blas.Uplo, n, kd, nrhs int, ab []float64, ldab int, b []float64, ldb int) { + switch { + case uplo != blas.Upper && uplo != blas.Lower: + panic(badUplo) + case n < 0: + panic(nLT0) + case kd < 0: + panic(kdLT0) + case nrhs < 0: + panic(nrhsLT0) + case ldab < kd+1: + panic(badLdA) + case ldb < max(1, nrhs): + panic(badLdB) + } + + // Quick return if possible. + if n == 0 || nrhs == 0 { + return + } + + if len(ab) < (n-1)*ldab+kd { + panic(shortAB) + } + if len(b) < (n-1)*ldb+nrhs { + panic(shortB) + } + + ldabConv := n + abConv := make([]float64, (kd+1)*ldabConv) + convDpbToLapacke(uplo, n, kd, ab, ldab, abConv, ldabConv) + lapacke.Dpbtrs(byte(uplo), n, kd, nrhs, abConv, ldabConv, b, ldb) +} + // Dpotrf computes the Cholesky decomposition of the symmetric positive definite // matrix a. If ul == blas.Upper, then a is stored as an upper-triangular matrix, // and a = U U^T is stored in place into a. If ul == blas.Lower, then a = L L^T diff --git a/lapack/netlib/lapack_test.go b/lapack/netlib/lapack_test.go index 5738eb2a..7fad08b4 100644 --- a/lapack/netlib/lapack_test.go +++ b/lapack/netlib/lapack_test.go @@ -104,6 +104,14 @@ func TestDlaswp(t *testing.T) { testlapack.DlaswpTest(t, impl) } +func TestDpbtrf(t *testing.T) { + testlapack.DpbtrfTest(t, impl) +} + +func TestDpbtrs(t *testing.T) { + testlapack.DpbtrsTest(t, impl) +} + func TestDpotrf(t *testing.T) { testlapack.DpotrfTest(t, impl) }