# HG changeset patch
# User Jason Riedy <address@hidden>
# Date 1236736489 14400
# Node ID 55d40cb27682ddf21a4c2f8d5d23f4a382956c4a
# Parent  b12a3be04411cb55054fe556b7b20ed6969faac8
Have sparse LU return permutation matrices rather than sparse matrices.

This could well impact user code.  It'd be interesting to see if there
is any actual fall-out...  Quite often, the permutation matrices are
applied to *dense* vectors.  Returning permutation matrices rather
than sparse matrices is a slight performance enhancement, but likely
lost in the noise.

diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,13 @@
+2009-03-10  Jason Riedy  <address@hidden>
+
+	* sparse-base-lu.cc (Pr_mat): New member function.  Return the row
+	permutation as a PermMatrix.
+	(Pc_mat): New member function.  Return the col permutation as a
+	PermMatrix.
+
+	* sparse-base-lu.h (sparse_base_lu): Declare Pc_mat and Pr_mat
+	member functions.
+
 2009-03-10  Jason Riedy  <address@hidden>
 
 	* Sparse-perm-op-defs.h (octinternal_do_mul_colpm_sm): New
diff --git a/liboctave/sparse-base-lu.cc b/liboctave/sparse-base-lu.cc
--- a/liboctave/sparse-base-lu.cc
+++ b/liboctave/sparse-base-lu.cc
@@ -98,6 +98,13 @@
 }
 
 template <class lu_type, class lu_elt_type, class p_type, class p_elt_type>
+PermMatrix
+sparse_base_lu <lu_type, lu_elt_type, p_type, p_elt_type> :: Pr_mat (void) const
+{
+  return PermMatrix (P, false);
+}
+
+template <class lu_type, class lu_elt_type, class p_type, class p_elt_type>
 p_type
 sparse_base_lu <lu_type, lu_elt_type, p_type, p_elt_type> :: Pc (void) const
 {
@@ -131,6 +138,13 @@
   return Pout;
 }
 
+template <class lu_type, class lu_elt_type, class p_type, class p_elt_type>
+PermMatrix
+sparse_base_lu <lu_type, lu_elt_type, p_type, p_elt_type> :: Pc_mat (void) const
+{
+  return PermMatrix (Q, true);
+}
+
 /*
 ;;; Local Variables: ***
 ;;; mode: C++ ***
diff --git a/liboctave/sparse-base-lu.h b/liboctave/sparse-base-lu.h
--- a/liboctave/sparse-base-lu.h
+++ b/liboctave/sparse-base-lu.h
@@ -70,6 +70,10 @@
 
   ColumnVector Pr_vec (void) const;
 
+  PermMatrix Pc_mat (void) const;
+
+  PermMatrix Pr_mat (void) const;
+
   const octave_idx_type * row_perm (void) const { return P.fortran_vec (); }
 
   const octave_idx_type * col_perm (void) const { return Q.fortran_vec (); }
diff --git a/src/ChangeLog b/src/ChangeLog
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,3 +1,8 @@
+2009-03-10  Jason Riedy  <address@hidden>
+
+	* DLD-FUNCTIONS/lu.cc (lu): Call fact.Pr_mat () and fact.Pc_mat ()
+	to return permutation matrices in the sparse case.
+
 2009-03-10  Jason Riedy  <address@hidden>
 
 	* OPERATORS/op-pm-sm.cc (mul_pm_sm): New Octave binding for
diff --git a/src/DLD-FUNCTIONS/lu.cc b/src/DLD-FUNCTIONS/lu.cc
--- a/src/DLD-FUNCTIONS/lu.cc
+++ b/src/DLD-FUNCTIONS/lu.cc
@@ -226,7 +226,7 @@
 		  retval (0) = fact.Y ();
 		else
 		  {
-		    SparseMatrix P = fact.Pr ();
+		    PermMatrix P = fact.Pr_mat ();
 		    SparseMatrix L = P.transpose () * fact.L ();
 		    retval(1) = octave_value (fact.U (), 
 					      MatrixType (MatrixType::Upper));
@@ -245,7 +245,7 @@
 		if (vecout)
 		  retval (2) = fact.Pr_vec ();
 		else
-		  retval(2) = fact.Pr ();
+		  retval(2) = fact.Pr_mat ();
 
 		retval(1) = octave_value (fact.U (), 
 					  MatrixType (MatrixType::Upper));
@@ -269,8 +269,8 @@
 		  }
 		else
 		  {
-		    retval(3) = fact.Pc ();
-		    retval(2) = fact.Pr ();
+		    retval(3) = fact.Pc_mat ();
+		    retval(2) = fact.Pr_mat ();
 		  }
 		retval(1) = octave_value (fact.U (), 
 					  MatrixType (MatrixType::Upper));
@@ -296,7 +296,7 @@
 		  retval (0) = fact.Y ();
 		else
 		  {
-		    SparseMatrix P = fact.Pr ();
+		    PermMatrix P = fact.Pr_mat ();
 		    SparseComplexMatrix L = P.transpose () * fact.L ();
 		    retval(1) = octave_value (fact.U (), 
 					      MatrixType (MatrixType::Upper));
@@ -315,7 +315,7 @@
 		if (vecout)
 		  retval (2) = fact.Pr_vec ();
 		else
-		  retval(2) = fact.Pr ();
+		  retval(2) = fact.Pr_mat ();
 
 		retval(1) = octave_value (fact.U (), 
 					  MatrixType (MatrixType::Upper));
@@ -339,8 +339,8 @@
 		  }
 		else
 		  {
-		    retval(3) = fact.Pc ();
-		    retval(2) = fact.Pr ();
+		    retval(3) = fact.Pc_mat ();
+		    retval(2) = fact.Pr_mat ();
 		  }
 		retval(1) = octave_value (fact.U (), 
 					  MatrixType (MatrixType::Upper));