changeset 8969:3ecbc236e2e0

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.
author Jason Riedy <jason@acm.org>
date Tue, 10 Mar 2009 21:54:49 -0400
parents 91d53dc37f79
children b37a6c27c23f
files liboctave/ChangeLog liboctave/sparse-base-lu.cc liboctave/sparse-base-lu.h src/ChangeLog src/DLD-FUNCTIONS/lu.cc
diffstat 5 files changed, 41 insertions(+), 8 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -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
--- 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++ ***
--- 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 (); }
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -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
--- 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));