lapack/netlib: add Dpbtrf and Dpbtrs

Also move the banded matrix conversion code to the lapack/netlib package
because that's where any matrix conversion should be done. The lapacke
package should accept exactly what LAPACKE accepts which means that
unfortunately for banded matrices there will be the inevitable overhead
of two conversions: one from BLAS (Gonum) row-major format to LAPACKE
row-major format for banded matrices and one inside LAPACKE to the
FORTRAN column-major banded format. In case of Dpbtrf the inverse
conversion must be performed for the factored matrix.

lapack/netlib: add Dpbtrs

conv

Author: vladimir-ch

lapack/lapacke/internal/conv/conv.go

@@ -2,13 +2,32 @@
 // 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"
+
+func DpbToRowMajor(uplo byte, n, kd int, a []float64, lda int, b []float64, ldb int) {
+	if uplo == 'U' {
+		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]
+			}
+		}
+	} 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]
+			}
+		}
+	}
+}
+
+// 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 +41,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 +57,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 DpbToColMajor.
+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]
 			}
 		}
 	}

lapack/lapacke/internal/conv/pb_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)
+					}
+				}
+			}
+		}
+	}
+}