mat: add Dense.PermuteRows and PermuteCols

lapack/lapack64/lapack64.go

@@ -649,21 +649,21 @@ func Lantb(norm lapack.MatrixNorm, a blas64.TriangularBand, work []float64) floa
 //
 //	X[i,0:n] is moved to X[k[i],0:n] for i=0,1,...,m-1.
 //
-// k must have length m, otherwise Lapmr will panic.
+// k must have length m, otherwise Lapmr will panic. k is zero-indexed.
 func Lapmr(forward bool, x blas64.General, k []int) {
 	lapack64.Dlapmr(forward, x.Rows, x.Cols, x.Data, max(1, x.Stride), k)
 }
 
 // Lapmt rearranges the columns of the m×n matrix X as specified by the
-// permutation k_0, k_1, ..., k_{n-1} of the integers 0, ..., n-1.
+// permutation k[0],k[1],...,k[n-1] of the integers 0,...,n-1.
 //
-// If forward is true a forward permutation is performed:
+// If forward is true, a forward permutation is applied:
 //
-//	X[0:m, k[j]] is moved to X[0:m, j] for j = 0, 1, ..., n-1.
+//	X[0:m,k[j]] is moved to X[0:m,j] for j=0,1,...,n-1.
 //
-// otherwise a backward permutation is performed:
+// If forward is false, a backward permutation is applied:
 //
-//	X[0:m, j] is moved to X[0:m, k[j]] for j = 0, 1, ..., n-1.
+//	X[0:m,j] is moved to X[0:m,k[j]] for j=0,1,...,n-1.
 //
 // k must have length n, otherwise Lapmt will panic. k is zero-indexed.
 func Lapmt(forward bool, x blas64.General, k []int) {

mat/dense.go

@@ -613,3 +613,58 @@ func (m *Dense) Norm(norm float64) float64 {
 	}
 	return lapack64.Lange(lnorm, m.mat, nil)
 }
+
+// Permutation constructs an n×n permutation matrix P from the given
+// row permutation such that the nonzero entries are P[i,p[i]] = 1.
+func (m *Dense) Permutation(n int, p []int) {
+	if len(p) != n {
+		panic(badSliceLength)
+	}
+	m.reuseAsZeroed(n, n)
+	for i, v := range p {
+		if v < 0 || v >= n {
+			panic(ErrRowAccess)
+		}
+		m.mat.Data[i*m.mat.Stride+v] = 1
+	}
+}
+
+// PermuteRows rearranges the rows of the m×n matrix A in the receiver as
+// specified by the permutation p[0],p[1],...,p[m-1] of the integers 0,...,m-1.
+//
+// If inverse is false, the given permutation is applied:
+//
+//	A[p[i],0:n] is moved to A[i,0:n] for i=0,1,...,m-1.
+//
+// If inverse is true, the inverse permutation is applied:
+//
+//	A[i,0:n] is moved to A[p[i],0:n] for i=0,1,...,m-1.
+//
+// p must have length m, otherwise PermuteRows will panic.
+func (m *Dense) PermuteRows(p []int, inverse bool) {
+	r, _ := m.Dims()
+	if len(p) != r {
+		panic(badSliceLength)
+	}
+	lapack64.Lapmr(!inverse, m.mat, p)
+}
+
+// PermuteCols rearranges the columns of the m×n matrix A in the reciever as
+// specified by the permutation p[0],p[1],...,p[n-1] of the integers 0,...,n-1.
+//
+// If inverse is false, the given permutation is applied:
+//
+//	A[0:m,p[j]] is moved to A[0:m,j] for j = 0, 1, ..., n-1.
+//
+// If inverse is true, the inverse permutation is applied:
+//
+//	A[0:m,j] is moved to A[0:m,p[j]] for j = 0, 1, ..., n-1.
+//
+// p must have length n, otherwise PermuteCols will panic.
+func (m *Dense) PermuteCols(p []int, inverse bool) {
+	_, c := m.Dims()
+	if len(p) != c {
+		panic(badSliceLength)
+	}
+	lapack64.Lapmt(!inverse, m.mat, p)
+}

mat/dense_test.go

@@ -15,6 +15,7 @@ import (
 
 	"gonum.org/v1/gonum/blas/blas64"
 	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/stat/combin"
 )
 
 func TestNewDense(t *testing.T) {
@@ -2135,6 +2136,136 @@ func TestDenseInverse(t *testing.T) {
 	}
 }
 
+func TestDensePermutation(t *testing.T) {
+	for n := 1; n <= 6; n++ {
+		for idx, perm := range combin.Permutations(n, n) {
+			want := NewDense(n, n, nil)
+			for i := 0; i < n; i++ {
+				want.Set(i, perm[i], 1)
+			}
+
+			var got Dense
+			got.Permutation(n, perm)
+			if !Equal(&got, want) {
+				t.Errorf("n=%d,idx=%d: unexpected permutation matrix\n got=%v\nwant=%v",
+					n, idx, Formatted(&got, Prefix("     ")), Formatted(want, Prefix("     ")))
+			}
+		}
+	}
+}
+
+func TestDensePermuteRows(t *testing.T) {
+	rnd := rand.New(rand.NewSource(1))
+	for m := 1; m <= 5; m++ {
+		for idx, perm := range combin.Permutations(m, m) {
+			// Construct a permutation matrix P from perm.
+			p := NewDense(m, m, nil)
+			for i := 0; i < m; i++ {
+				p.Set(i, perm[i], 1)
+			}
+
+			for n := 1; n <= 5; n++ {
+				name := fmt.Sprintf("m=%d,n=%d,idx=%d", m, n, idx)
+
+				// Generate a random m×n matrix A.
+				a := NewDense(m, n, nil)
+				for i := 0; i < m; i++ {
+					for j := 0; j < n; j++ {
+						a.Set(i, j, rnd.Float64())
+					}
+				}
+
+				// Compute P*A.
+				var want Dense
+				want.Mul(p, a)
+				// Permute rows of A by applying perm.
+				var got Dense
+				got.CloneFrom(a)
+				got.PermuteRows(perm, false)
+				if !Equal(&got, &want) {
+					t.Errorf("%s: unexpected result\n got=%v\nwant=%v",
+						name, Formatted(&got, Prefix("     ")), Formatted(&want, Prefix("     ")))
+				}
+
+				// Compute Pᵀ*A.
+				want.Mul(p.T(), a)
+				// Permute rows of A by applying inverse perm.
+				got.Copy(a)
+				got.PermuteRows(perm, true)
+				if !Equal(&got, &want) {
+					t.Errorf("%s: unexpected result with inverse permutation\n got=%v\nwant=%v",
+						name, Formatted(&got, Prefix("     ")), Formatted(&want, Prefix("     ")))
+				}
+
+				// Permute rows of A by applying perm and inverse perm.
+				got.Copy(a)
+				got.PermuteRows(perm, false)
+				got.PermuteRows(perm, true)
+				if !Equal(&got, a) {
+					t.Errorf("%s: original A not recovered\n got=%v\nwant=%v",
+						name, Formatted(&got, Prefix("     ")), Formatted(a, Prefix("     ")))
+				}
+			}
+		}
+	}
+}
+
+func TestDensePermuteCols(t *testing.T) {
+	rnd := rand.New(rand.NewSource(1))
+	for n := 1; n <= 5; n++ {
+		for idx, perm := range combin.Permutations(n, n) {
+			// Construct a permutation matrix P from perm.
+			p := NewDense(n, n, nil)
+			for j := 0; j < n; j++ {
+				p.Set(perm[j], j, 1)
+			}
+
+			for m := 1; m <= 5; m++ {
+				name := fmt.Sprintf("m=%d,n=%d,idx=%d", m, n, idx)
+
+				// Generate a random m×n matrix A.
+				a := NewDense(m, n, nil)
+				for i := 0; i < m; i++ {
+					for j := 0; j < n; j++ {
+						a.Set(i, j, rnd.Float64())
+					}
+				}
+
+				// Compute A*P.
+				var want Dense
+				want.Mul(a, p)
+				// Permute columns of A by applying perm.
+				var got Dense
+				got.CloneFrom(a)
+				got.PermuteCols(perm, false)
+				if !Equal(&got, &want) {
+					t.Errorf("%s: unexpected result\n got=%v\nwant=%v",
+						name, Formatted(&got, Prefix("     ")), Formatted(&want, Prefix("     ")))
+				}
+
+				// Compute A*Pᵀ.
+				want.Mul(a, p.T())
+				// Permute columns of A by applying inverse perm.
+				got.Copy(a)
+				got.PermuteCols(perm, true)
+				if !Equal(&got, &want) {
+					t.Errorf("%s: unexpected result with inverse permutation\n got=%v\nwant=%v",
+						name, Formatted(&got, Prefix("     ")), Formatted(&want, Prefix("     ")))
+				}
+
+				// Permute columns of A by applying perm and inverse perm.
+				got.Copy(a)
+				got.PermuteCols(perm, false)
+				got.PermuteCols(perm, true)
+				if !Equal(&got, a) {
+					t.Errorf("%s: original A not recovered\n got=%v\nwant=%v",
+						name, Formatted(&got, Prefix("     ")), Formatted(a, Prefix("     ")))
+				}
+			}
+		}
+	}
+}
+
 var (
 	wd *Dense
 )

mat/lu.go

@@ -314,22 +314,6 @@ func (lu *LU) UTo(dst *TriDense) {
 	}
 }
 
-// Permutation constructs an r×r permutation matrix with the given row swaps.
-// A permutation matrix has exactly one element equal to one in each row and column
-// and all other elements equal to zero. swaps[i] specifies the row with which
-// i will be swapped, which is equivalent to the non-zero column of row i.
-func (m *Dense) Permutation(r int, swaps []int) {
-	m.reuseAsNonZeroed(r, r)
-	for i := 0; i < r; i++ {
-		zero(m.mat.Data[i*m.mat.Stride : i*m.mat.Stride+r])
-		v := swaps[i]
-		if v < 0 || v >= r {
-			panic(ErrRowAccess)
-		}
-		m.mat.Data[i*m.mat.Stride+v] = 1
-	}
-}
-
 // SolveTo solves a system of linear equations using the LU decomposition of a matrix.
 // It computes
 //