--- a/liboctave/CSparse.cc	Thu Mar 12 11:46:56 2009 +0100
+++ b/liboctave/CSparse.cc	Tue Mar 10 21:54:44 2009 -0400
@@ -52,6 +52,8 @@
 
 #include "Sparse-diag-op-defs.h"
 
+#include "Sparse-perm-op-defs.h"
+
 // Define whether to use a basic QR solver or one that uses a Dulmange
 // Mendelsohn factorization to seperate the problem into under-determined,
 // well-determined and over-determined parts and solves them seperately
@@ -7684,6 +7686,20 @@
   return do_sub_sm_dm<SparseComplexMatrix> (a, d);
 }
 
+// perm * sparse and sparse * perm
+
+SparseComplexMatrix
+operator * (const PermMatrix& p, const SparseComplexMatrix& a)
+{
+  return octinternal_do_mul_pm_sm (p, a);
+}
+
+SparseComplexMatrix
+operator * (const SparseComplexMatrix& a, const PermMatrix& p)
+{
+  return octinternal_do_mul_sm_pm (a, p);
+}
+
 // FIXME -- it would be nice to share code among the min/max
 // functions below.
 

--- a/liboctave/CSparse.h	Thu Mar 12 11:46:56 2009 +0100
+++ b/liboctave/CSparse.h	Tue Mar 10 21:54:44 2009 -0400
@@ -37,6 +37,7 @@
 #include "Sparse-op-defs.h"
 #include "MatrixType.h"
 
+class PermMatrix;
 class DiagMatrix;
 class ComplexDiagMatrix;
 class SparseMatrix;
@@ -498,6 +499,11 @@
 extern OCTAVE_API SparseComplexMatrix operator - (const SparseComplexMatrix&, const DiagMatrix&);
 extern OCTAVE_API SparseComplexMatrix operator - (const SparseComplexMatrix&, const ComplexDiagMatrix&);
 
+extern OCTAVE_API SparseComplexMatrix operator * (const PermMatrix&,
+						  const SparseComplexMatrix&);
+extern OCTAVE_API SparseComplexMatrix operator * (const SparseComplexMatrix&,
+						  const PermMatrix&);
+
 extern OCTAVE_API SparseComplexMatrix min (const Complex& c, 
 				const SparseComplexMatrix& m);
 extern OCTAVE_API SparseComplexMatrix min (const SparseComplexMatrix& m, 

--- a/liboctave/ChangeLog	Thu Mar 12 11:46:56 2009 +0100
+++ b/liboctave/ChangeLog	Tue Mar 10 21:54:44 2009 -0400
@@ -89,6 +89,32 @@
 
 2009-03-10  Jason Riedy  <jason@acm.org>
 
+	* Sparse-perm-op-defs.h (octinternal_do_mul_colpm_sm): New
+	template function.  Logic for the column permutation * sparse
+	matrix operator.
+	(octinternal_do_mul_pm_sm): New template function.  Logic for the
+	permutation matrix * sparse matrix operator.  Note that there is
+	no special row perm * sparse routine; the permutation is inverted
+	and the col perm routine is called.
+	(octinternal_do_mul_sm_rowpm): New template function.  Logic for
+	the sparse matrix * row permutation operator.
+	(octinternal_do_mul_sm_colpm): New template function.  Logic for
+	the sparse matrix * column permutation operator.
+	(octinternal_do_mul_sm_pm): New template function.  Logic for the
+	sparse matrix * permutation matrix operator.
+
+	* dSparse.h (operator *): Declare sparse * permutation and
+	permutation * sparse.
+	* dSparse.cc (operator *): Define sparse * permutation and
+	permutation * sparse.
+
+	* CSparse.h (operator *): Declare sparse * permutation and
+	permutation * sparse.
+	* CSparse.cc (operator *): Define sparse * permutation and
+	permutation * sparse.
+
+2009-03-10  Jason Riedy  <jason@acm.org>
+
 	* sparse-base-lu.cc (Pc_vec): The column permutation should be
 	Ufact.cols ()-long, not Lfact.rows ()-long.
 

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/Sparse-perm-op-defs.h	Tue Mar 10 21:54:44 2009 -0400
@@ -0,0 +1,164 @@
+/* -*- C++ -*-
+
+Copyright (C) 2009 Jason Riedy
+
+This file is part of Octave.
+
+Octave is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+Octave is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, see
+<http://www.gnu.org/licenses/>.
+
+*/
+
+#if !defined (octave_sparse_perm_op_defs_h)
+#define octave_sparse_perm_op_defs_h 1
+
+// Matrix multiplication
+
+template <typename SM>
+SM octinternal_do_mul_colpm_sm (const octave_idx_type *pcol, const SM& a)
+// Relabel the rows according to pcol.
+{
+  const octave_idx_type nr = a.rows ();
+  const octave_idx_type nc = a.cols ();
+  const octave_idx_type nent = a.nnz ();
+  SM r (nr, nc, nent);
+
+  for (octave_idx_type k = 0; k < nent; ++k)
+    {
+      OCTAVE_QUIT;
+      r.xridx (k) = pcol[a.ridx (k)];
+      r.xdata (k) = a.data (k);
+    }
+  for (octave_idx_type j = 0; j <= nc; ++j)
+    r.xcidx (j) = a.cidx (j);
+
+  r.maybe_compress (false);
+  return r;
+}
+
+template <typename SM>
+SM octinternal_do_mul_pm_sm (const PermMatrix& p, const SM& a)
+{
+  const octave_idx_type nr = a.rows ();
+  if (p.cols () != nr)
+    {
+      gripe_nonconformant ("operator *", p.rows (), p.cols (), a.rows (), a.cols ());
+      return SM ();
+    }
+
+  if (p.is_row_perm ())
+    {
+      // Form the column permutation and then call the colpm_sm routine.
+      const octave_idx_type *prow = p.pvec ().data ();
+      OCTAVE_LOCAL_BUFFER(octave_idx_type, pcol, nr);
+      for (octave_idx_type i = 0; i < nr; ++i)
+	pcol[prow[i]] = i;
+      return octinternal_do_mul_colpm_sm (pcol, a);
+    }
+  else
+    return octinternal_do_mul_colpm_sm (p.pvec ().data (), a);
+}
+
+template <typename SM>
+SM octinternal_do_mul_sm_rowpm (const SM& a, const octave_idx_type *prow)
+// For a row permutation, iterate across the source a and stuff the
+// results into the correct destination column in r.
+{
+  const octave_idx_type nr = a.rows ();
+  const octave_idx_type nc = a.cols ();
+  const octave_idx_type nent = a.nnz ();
+  SM r (nr, nc, nent);
+
+  for (octave_idx_type j_src = 0; j_src < nc; ++j_src)
+    r.xcidx (prow[j_src]) = a.cidx (j_src+1) - a.cidx (j_src);
+  octave_idx_type k = 0;
+  for (octave_idx_type j = 0; j < nc; ++j)
+    {
+      const octave_idx_type tmp = r.xcidx (j);
+      r.xcidx (j) = k;
+      k += tmp;
+    }
+  r.xcidx (nc) = nent;
+
+  octave_idx_type k_src = 0;
+  for (octave_idx_type j_src = 0; j_src < nc; ++j_src)
+    {
+      OCTAVE_QUIT;
+      const octave_idx_type j = prow[j_src];
+      const octave_idx_type kend_src = a.cidx (j_src + 1);
+      for (k = r.xcidx (j); k_src < kend_src; ++k, ++k_src)
+	{
+	  r.xridx (k) = a.ridx (k_src);
+	  r.xdata (k) = a.data (k_src);
+	}
+    }
+  assert (k_src == nent);
+
+  r.maybe_compress (false);
+  return r;
+}
+
+template <typename SM>
+SM octinternal_do_mul_sm_colpm (const SM& a, const octave_idx_type *pcol)
+// For a column permutation, iterate across the destination r and pull
+// data from the correct column of a.
+{
+  const octave_idx_type nr = a.rows ();
+  const octave_idx_type nc = a.cols ();
+  const octave_idx_type nent = a.nnz ();
+  SM r (nr, nc, nent);
+
+  for (octave_idx_type j = 0; j < nc; ++j)
+    {
+      const octave_idx_type j_src = pcol[j];
+      r.xcidx (j+1) = r.xcidx (j) + (a.cidx (j_src+1) - a.cidx (j_src));
+    }
+  assert (r.xcidx (nc) == nent);
+
+  octave_idx_type k = 0;
+  for (octave_idx_type j = 0; j < nc; ++j)
+    {
+      OCTAVE_QUIT;
+      const octave_idx_type j_src = pcol[j];
+      octave_idx_type k_src;
+      const octave_idx_type kend_src = a.cidx (j_src + 1);
+      for (k_src = a.cidx (j_src); k_src < kend_src; ++k_src, ++k)
+	{
+	  r.xridx (k) = a.ridx (k_src);
+	  r.xdata (k) = a.data (k_src);
+	}
+    }
+  assert (k == nent);
+
+  r.maybe_compress (false);
+  return r;
+}
+
+template <typename SM>
+SM octinternal_do_mul_sm_pm (const SM& a, const PermMatrix& p)
+{
+  const octave_idx_type nc = a.cols ();
+  if (p.rows () != nc)
+    {
+      gripe_nonconformant ("operator *", a.rows (), a.cols (), p.rows (), p.cols ());
+      return SM ();
+    }
+
+  if (p.is_row_perm ())
+    return octinternal_do_mul_sm_rowpm (a, p.pvec ().data ());
+  else
+    return octinternal_do_mul_sm_colpm (a, p.pvec ().data ());
+}
+
+#endif // octave_sparse_perm_op_defs_h

--- a/liboctave/dSparse.cc	Thu Mar 12 11:46:56 2009 +0100
+++ b/liboctave/dSparse.cc	Tue Mar 10 21:54:44 2009 -0400
@@ -53,6 +53,8 @@
 
 #include "Sparse-diag-op-defs.h"
 
+#include "Sparse-perm-op-defs.h"
+
 // Define whether to use a basic QR solver or one that uses a Dulmange
 // Mendelsohn factorization to seperate the problem into under-determined,
 // well-determined and over-determined parts and solves them seperately
@@ -7740,6 +7742,20 @@
   return do_sub_sm_dm<SparseMatrix> (a, d);
 }
 
+// perm * sparse and sparse * perm
+
+SparseMatrix
+operator * (const PermMatrix& p, const SparseMatrix& a)
+{
+  return octinternal_do_mul_pm_sm (p, a);
+}
+
+SparseMatrix
+operator * (const SparseMatrix& a, const PermMatrix& p)
+{
+  return octinternal_do_mul_sm_pm (a, p);
+}
+
 // FIXME -- it would be nice to share code among the min/max
 // functions below.
 

--- a/liboctave/dSparse.h	Thu Mar 12 11:46:56 2009 +0100
+++ b/liboctave/dSparse.h	Tue Mar 10 21:54:44 2009 -0400
@@ -36,6 +36,7 @@
 #include "Sparse-op-defs.h"
 #include "MatrixType.h"
 
+class PermMatrix;
 class DiagMatrix;
 class SparseComplexMatrix;
 class SparseBoolMatrix;
@@ -459,6 +460,9 @@
 extern OCTAVE_API SparseMatrix operator - (const DiagMatrix&, const SparseMatrix&);
 extern OCTAVE_API SparseMatrix operator - (const SparseMatrix&, const DiagMatrix&);
 
+extern OCTAVE_API SparseMatrix operator * (const PermMatrix&, const SparseMatrix&);
+extern OCTAVE_API SparseMatrix operator * (const SparseMatrix&, const PermMatrix&);
+
 extern OCTAVE_API SparseMatrix min (double d, const SparseMatrix& m);
 extern OCTAVE_API SparseMatrix min (const SparseMatrix& m, double d);
 extern OCTAVE_API SparseMatrix min (const SparseMatrix& a, const SparseMatrix& b);

--- a/src/ChangeLog	Thu Mar 12 11:46:56 2009 +0100
+++ b/src/ChangeLog	Tue Mar 10 21:54:44 2009 -0400
@@ -47,6 +47,25 @@
 
 2009-03-10  Jason Riedy  <jason@acm.org>
 
+	* OPERATORS/op-pm-sm.cc (mul_pm_sm): New Octave binding for
+	perm * sparse.
+	(ldiv_pm_sm): New Octave binding for perm \ sparse.
+	(mul_sm_pm): New Octave binding for sparse * perm.
+	(div_sm_pm): New Octave binding for sparse / perm.
+	(install_pm_sm_ops): Install the above bindings.
+
+	* OPERATORS/op-pm-scm.cc (mul_pm_scm): New Octave binding for
+	perm * sparse complex.
+	(ldiv_pm_scm): New Octave binding for perm \ sparse complex.
+	(mul_scm_pm): New Octave binding for sparse complex * perm.
+	(div_scm_pm): New Octave binding for sparse complex / perm.
+	(install_pm_scm_ops): Install the above bindings.
+
+	* Makefile.in (PERM_OP_XSRC): Add op-pm-sm.cc and op-pm-scm.cc for
+	operations between permutations and sparse matrices.
+
+2009-03-10  Jason Riedy  <jason@acm.org>
+
 	* DLD-FUNCTIONS/find.cc (find_nonzero_elem_idx): New override
 	for find on PermMatrix.
 	(find): Add a branch testing arg.is_perm_matrix () and calling the

--- a/src/Makefile.in	Thu Mar 12 11:46:56 2009 +0100
+++ b/src/Makefile.in	Tue Mar 10 21:54:44 2009 -0400
@@ -174,7 +174,8 @@
 	op-fdm-fm.cc op-fdm-fs.cc op-fm-fcdm.cc op-fm-fdm.cc
 
 PERM_OP_XSRC := op-cm-pm.cc op-fcm-pm.cc op-fm-pm.cc op-pm-fcm.cc \
-	op-pm-fm.cc op-m-pm.cc op-pm-cm.cc op-pm-m.cc op-pm-pm.cc
+	op-pm-fm.cc op-m-pm.cc op-pm-cm.cc op-pm-m.cc op-pm-pm.cc \
+	op-pm-sm.cc op-pm-scm.cc
 
 OP_XSRC :=  op-b-b.cc op-b-bm.cc op-bm-b.cc op-bm-bm.cc op-cell.cc \
 	op-chm.cc op-class.cc op-list.cc op-range.cc op-str-m.cc \

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/OPERATORS/op-pm-scm.cc	Tue Mar 10 21:54:44 2009 -0400
@@ -0,0 +1,97 @@
+/*
+
+Copyright (C) 2009 Jason Riedy
+
+This file is part of Octave.
+
+Octave is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+Octave is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, see
+<http://www.gnu.org/licenses/>.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "gripes.h"
+#include "oct-obj.h"
+#include "ov.h"
+#include "ov-typeinfo.h"
+#include "ops.h"
+
+#include "ov-perm.h"
+#include "ov-cx-sparse.h"
+
+// permutation matrix by sparse matrix ops
+
+DEFBINOP (mul_pm_scm, perm_matrix, sparse_complex_matrix)
+{
+  CAST_BINOP_ARGS (const octave_perm_matrix&, const octave_sparse_complex_matrix&);
+
+  if (v2.rows() == 1 && v2.columns() == 1)
+    {
+      std::complex<double> d = v2.complex_value ();
+
+      return octave_value (v1.sparse_matrix_value () * d);
+    }
+  else if (v1.rows() == 1 && v1.columns() == 1)
+    return octave_value (v2.sparse_complex_matrix_value ());
+  else
+    return v1.perm_matrix_value  () * v2.sparse_complex_matrix_value ();
+}
+
+DEFBINOP (ldiv_pm_scm, perm_matrix, sparse_complex_matrix)
+{
+  CAST_BINOP_ARGS (const octave_perm_matrix&, const octave_sparse_complex_matrix&);
+
+  return v1.perm_matrix_value ().inverse () * v2.sparse_complex_matrix_value (); 
+}
+
+// sparse matrix by diagonal matrix ops
+
+DEFBINOP (mul_scm_pm, sparse_complex_matrix, perm_matrix)
+{
+  CAST_BINOP_ARGS (const octave_sparse_complex_matrix&, const octave_perm_matrix&);
+
+  if (v1.rows() == 1 && v1.columns() == 1)
+    {
+      std::complex<double> d = v1.scalar_value ();
+
+      return octave_value (d * v2.sparse_matrix_value ());
+    }
+  else if (v2.rows() == 1 && v2.columns() == 1)
+    return octave_value (v1.sparse_complex_matrix_value ());
+  else
+    return v1.sparse_complex_matrix_value  () * v2.perm_matrix_value ();
+}
+
+DEFBINOP (div_scm_pm, sparse_complex_matrix, perm_matrix)
+{
+  CAST_BINOP_ARGS (const octave_sparse_complex_matrix&, const octave_perm_matrix&);
+
+  return v1.sparse_complex_matrix_value () * v2.perm_matrix_value ().inverse ();
+}
+
+void
+install_pm_scm_ops (void)
+{
+  INSTALL_BINOP (op_mul, octave_perm_matrix, octave_sparse_complex_matrix,
+		 mul_pm_scm);
+  INSTALL_BINOP (op_ldiv, octave_perm_matrix, octave_sparse_complex_matrix,
+		 ldiv_pm_scm);
+  INSTALL_BINOP (op_mul, octave_sparse_complex_matrix, octave_perm_matrix,
+		 mul_scm_pm);
+  INSTALL_BINOP (op_div, octave_sparse_complex_matrix, octave_perm_matrix,
+		 div_scm_pm);
+}

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/OPERATORS/op-pm-sm.cc	Tue Mar 10 21:54:44 2009 -0400
@@ -0,0 +1,97 @@
+/*
+
+Copyright (C) 2009 Jason Riedy
+
+This file is part of Octave.
+
+Octave is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+Octave is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, see
+<http://www.gnu.org/licenses/>.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "gripes.h"
+#include "oct-obj.h"
+#include "ov.h"
+#include "ov-typeinfo.h"
+#include "ops.h"
+
+#include "ov-perm.h"
+#include "ov-re-sparse.h"
+
+// permutation matrix by sparse matrix ops
+
+DEFBINOP (mul_pm_sm, perm_matrix, sparse_matrix)
+{
+  CAST_BINOP_ARGS (const octave_perm_matrix&, const octave_sparse_matrix&);
+
+  if (v2.rows() == 1 && v2.columns() == 1)
+    {
+      double d = v2.scalar_value ();
+
+      return octave_value (v1.sparse_matrix_value () * d);
+    }
+  else if (v1.rows() == 1 && v1.columns() == 1)
+    return octave_value (v2.sparse_matrix_value ());
+  else
+    return v1.perm_matrix_value  () * v2.sparse_matrix_value ();
+}
+
+DEFBINOP (ldiv_pm_sm, perm_matrix, sparse_matrix)
+{
+  CAST_BINOP_ARGS (const octave_perm_matrix&, const octave_sparse_matrix&);
+
+  return v1.perm_matrix_value ().inverse () * v2.sparse_matrix_value (); 
+}
+
+// sparse matrix by diagonal matrix ops
+
+DEFBINOP (mul_sm_pm, sparse_matrix, perm_matrix)
+{
+  CAST_BINOP_ARGS (const octave_sparse_matrix&, const octave_perm_matrix&);
+
+  if (v1.rows() == 1 && v1.columns() == 1)
+    {
+      double d = v1.scalar_value ();
+
+      return octave_value (d * v2.sparse_matrix_value ());
+    }
+  else if (v2.rows() == 1 && v2.columns() == 1)
+    return octave_value (v1.sparse_matrix_value ());
+  else
+    return v1.sparse_matrix_value  () * v2.perm_matrix_value ();
+}
+
+DEFBINOP (div_sm_pm, sparse_matrix, perm_matrix)
+{
+  CAST_BINOP_ARGS (const octave_sparse_matrix&, const octave_perm_matrix&);
+
+  return v1.sparse_matrix_value () * v2.perm_matrix_value ().inverse ();
+}
+
+void
+install_pm_sm_ops (void)
+{
+  INSTALL_BINOP (op_mul, octave_perm_matrix, octave_sparse_matrix,
+		 mul_pm_sm);
+  INSTALL_BINOP (op_ldiv, octave_perm_matrix, octave_sparse_matrix,
+		 ldiv_pm_sm);
+  INSTALL_BINOP (op_mul, octave_sparse_matrix, octave_perm_matrix,
+		 mul_sm_pm);
+  INSTALL_BINOP (op_div, octave_sparse_matrix, octave_perm_matrix,
+		 div_sm_pm);
+}

--- a/test/ChangeLog	Thu Mar 12 11:46:56 2009 +0100
+++ b/test/ChangeLog	Tue Mar 10 21:54:44 2009 -0400
@@ -1,3 +1,9 @@
+2009-03-10  Jason Riedy  <jason@acm.org>
+
+	* test_diag_perm.m: Add tests for permuting sparse matrices and
+	for the correct types from interactions between
+	pseudo-scalars (1x1 matrices).
+
 2009-03-10  Jason Riedy  <jason@acm.org>
 
 	* build_sparse_tests.sh: Add LU tests to the rectangular tests.

--- a/test/test_diag_perm.m	Thu Mar 12 11:46:56 2009 +0100
+++ b/test/test_diag_perm.m	Tue Mar 10 21:54:44 2009 -0400
@@ -57,6 +57,52 @@
 %! assert (typeinfo (rand (n) * Pc), "matrix");
 %! assert (typeinfo (Pr * rand (n)), "matrix");
 
+## preserve sparse matrix structure
+%!test
+%! n = 7;
+%! Pc = eye (n) (:, randperm (n));
+%! Ac = sprand (n-3, n, .5) + I () * sprand (n-3, n, .5);
+%! Pr = eye (n) (randperm (n), :);
+%! Ar = sprand (n, n+2, .5);
+%! assert (typeinfo (Ac * Pc), "sparse complex matrix");
+%! assert (full (Ac * Pc), full (Ac) * Pc);
+%! assert (full (Ac / Pc), full (Ac) / Pc);
+%! assert (typeinfo (Pr * Ar), "sparse matrix");
+%! assert (full (Pr * Ar), Pr * full (Ar));
+%! assert (full (Pr \ Ar), Pr \ full (Ar));
+
+## structure rules for 1x1 dense / scalar and 1x1 perm
+%!test
+%! n = 7;
+%! P1 = eye (1) (:, [1]);
+%! A1 = 1;
+%! P = eye (n) (:, randperm (n));
+%! A = rand (n-3, n, .5);
+%! assert (typeinfo (A * P1), "matrix");
+%! assert (full (A * P1), full (A) * P1);
+%! assert (typeinfo (P1 * A), "matrix");
+%! assert (full (P1 * A), P1 * full (A));
+%! assert (typeinfo (A1 * P), "matrix");
+%! assert (full (A1 * P), full (A1) * P);
+%! assert (typeinfo (P * A1), "matrix");
+%! assert (full (P * A1), P * full (A1));
+
+## structure rules for 1x1 sparse and 1x1 perm
+%!test
+%! n = 7;
+%! P1 = eye (1) (:, [1]);
+%! A1 = sparse (1, 1, 2);
+%! P = eye (n) (:, randperm (n));
+%! A = sprand (n-3, n, .5);
+%! assert (typeinfo (A * P1), "sparse matrix");
+%! assert (full (A * P1), full (A) * P1);
+%! assert (typeinfo (P1 * A), "sparse matrix");
+%! assert (full (P1 * A), P1 * full (A));
+%! assert (typeinfo (A1 * P), "sparse matrix");
+%! assert (full (A1 * P), full (A1) * P);
+%! assert (typeinfo (P * A1), "sparse matrix");
+%! assert (full (P * A1), P * full (A1));
+
 ## permuting a matrix with exceptional values does not introduce new ones.
 %!test
 %! n = 5;