lapacke/internal: add conv package with DpbToColMajor and DpbToRowMajor

lapack/lapacke/internal/conv/conv.go

@@ -0,0 +1,19 @@
+// 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
+}

lapack/lapacke/internal/conv/pb.go

@@ -0,0 +1,81 @@
+// 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
+
+// DpbToColMajor converts a symmetric band matrix A in CBLAS row-major layout
+// to FORTRAN column-major layout and stores the result in B.
+//
+// For example, when n = 6, kd = 2 and uplo == 'U', DpbToColMajor
+// converts
+//  A = a00  a01  a02
+//      a11  a12  a13
+//      a22  a23  a24
+//      a33  a34  a35
+//      a44  a45   *
+//      a55   *    *
+// stored in a slice as
+//  a = [a00 a01 a02 a11 a12 a13 a22 a23 a24 a33 a34 a35 a44 a45 * a55 * *]
+// to
+//  B =  *   *  a02 a13 a24 a35
+//       *  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]
+//
+// When n = 6, kd = 2 and uplo == 'L', DpbToColMajor converts
+//  A =  *    *   a00
+//       *   a10  a11
+//      a20  a21  a22
+//      a31  a32  a33
+//      a42  a43  a44
+//      a53  a54  a55
+// stored in a slice as
+//  a = [* * a00 * a10 a11 a20 a21 a22 a31 a32 a33 a42 a43 a44 a53 a54 a55]
+// to
+//  B = a00 a11 a22 a33 a44 a55
+//      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 * *]
+//
+// 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' {
+		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]
+			}
+		}
+	} 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]
+			}
+		}
+	}
+}
+
+// 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' {
+		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]
+			}
+		}
+	}
+}

lapack/lapacke/internal/conv/pb_test.go

@@ -0,0 +1,123 @@
+// 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)
+					}
+				}
+			}
+		}
+	}
+}