diff --git a/lapack/lapacke/internal/conv/conv.go b/lapack/lapacke/internal/conv/conv.go new file mode 100644 index 00000000..33247990 --- /dev/null +++ b/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 +} diff --git a/lapack/lapacke/internal/conv/pb.go b/lapack/lapacke/internal/conv/pb.go new file mode 100644 index 00000000..16fd13be --- /dev/null +++ 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] + } + } + } +} diff --git a/lapack/lapacke/internal/conv/pb_test.go b/lapack/lapacke/internal/conv/pb_test.go new file mode 100644 index 00000000..a6d1c2d8 --- /dev/null +++ b/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) + } + } + } + } + } +}