lapack/netlib: add Dpbtrf and Dpbtrs

lapack/lapacke/internal/conv/pb.go → 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]
 			}
 		}
 	}

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)
+					}
+				}
+			}
+		}
+	}
+}