# HG changeset patch
# User Jason Riedy <jason@acm.org>
# Date 1236736489 14400
# Node ID 3ecbc236e2e0dc9c435e2b76667715c998fc5714
# Parent  91d53dc37f793cbe9f252d2ae85f169e89d7294d
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 -r 91d53dc37f79 -r 3ecbc236e2e0 liboctave/ChangeLog
--- a/liboctave/ChangeLog	Tue Mar 10 21:54:44 2009 -0400
+++ b/liboctave/ChangeLog	Tue Mar 10 21:54:49 2009 -0400
@@ -1,3 +1,13 @@
+2009-03-10  Jason Riedy  <jason@acm.org>
+
+	* 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-09  Jason Riedy  <jason@acm.org>
 
 	* Sparse-diag-op-defs.h (octave_impl::inner_do_add_sm_dm): New
diff -r 91d53dc37f79 -r 3ecbc236e2e0 liboctave/sparse-base-lu.cc
--- a/liboctave/sparse-base-lu.cc	Tue Mar 10 21:54:44 2009 -0400
+++ b/liboctave/sparse-base-lu.cc	Tue Mar 10 21:54:49 2009 -0400
@@ -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 -r 91d53dc37f79 -r 3ecbc236e2e0 liboctave/sparse-base-lu.h
--- a/liboctave/sparse-base-lu.h	Tue Mar 10 21:54:44 2009 -0400
+++ b/liboctave/sparse-base-lu.h	Tue Mar 10 21:54:49 2009 -0400
@@ -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 -r 91d53dc37f79 -r 3ecbc236e2e0 src/ChangeLog
--- a/src/ChangeLog	Tue Mar 10 21:54:44 2009 -0400
+++ b/src/ChangeLog	Tue Mar 10 21:54:49 2009 -0400
@@ -1,3 +1,8 @@
+2009-03-10  Jason Riedy  <jason@acm.org>
+
+	* DLD-FUNCTIONS/lu.cc (lu): Call fact.Pr_mat () and fact.Pc_mat ()
+	to return permutation matrices in the sparse case.
+
 2009-03-12  John W. Eaton  <jwe@octave.org>
 
 	* ls-mat-ascii.cc (get_mat_data_input_line): If we are looking at
diff -r 91d53dc37f79 -r 3ecbc236e2e0 src/DLD-FUNCTIONS/lu.cc
--- a/src/DLD-FUNCTIONS/lu.cc	Tue Mar 10 21:54:44 2009 -0400
+++ b/src/DLD-FUNCTIONS/lu.cc	Tue Mar 10 21:54:49 2009 -0400
@@ -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));