Mercurial > octave-nkf

changeset 5282:5bdc3f24cd5f

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2005-04-14 22:17:26 by dbateman]
author dbateman
date Thu, 14 Apr 2005 22:17:27 +0000
parents f3266e7dbb99
children bb224e6c26a7
files PROJECTS doc/ChangeLog doc/interpreter/sparse.txi liboctave/ChangeLog liboctave/SparseCmplxLU.cc liboctave/SparseCmplxLU.h liboctave/SparseType.cc liboctave/SparseType.h liboctave/SparsedbleLU.cc liboctave/SparsedbleLU.h src/ChangeLog src/DLD-FUNCTIONS/luinc.cc src/DLD-FUNCTIONS/sparse.cc src/DLD-FUNCTIONS/splu.cc src/Makefile.in src/ov-bool-sparse.h src/ov-mapper.cc src/ov-re-sparse.cc src/ov-re-sparse.h
diffstat 19 files changed, 1020 insertions(+), 306 deletions(-) [+]
line wrap: on
line diff
--- a/PROJECTS	Thu Apr 14 19:35:20 2005 +0000
+++ b/PROJECTS	Thu Apr 14 22:17:27 2005 +0000
@@ -116,22 +116,13 @@
   * Once dmperm is implemented, use the technique to detect permuted
     triangular matrices. Test the permuted triangular matrix solver code
 
-  * Accelerate the copying of the data from a sparse matrix to a banded matrix
-    in the solvers, that takes a significant portion of the computation time 
-    for narrow matrices. This is not obvious, due to the treatment of zero
-    elements in the matrix. Maybe current solution is optimal.
-
-  * Perhaps split the overly long ::solve functions up, either by the type 
-    of solver, or seperate factorization functions, so that they can be 
-    reused in each of 4 different ::solve functions.
-
   * Sparse inverse function, based on Andy's code from octave-forge.
 
   * Implement fourth argument to the sprand and sprandn that the leading
     brand implements.
 
-  * Mapper functions such as real, imag, abs, etc need to be treated, either
-    with a dispatch or by changing the mapper function code.
+  * Sparse logical indexing in idx_vector class so that something like
+    "a=sprandn(1e6,1e6,1e-6); a(a<1) = 0" won't cause a memory overflow.
 
   * Write the rest of the sparse docs
 
--- a/doc/ChangeLog	Thu Apr 14 19:35:20 2005 +0000
+++ b/doc/ChangeLog	Thu Apr 14 22:17:27 2005 +0000
@@ -1,3 +1,7 @@
+2005-03-14  David Bateman  <dbateman@free.fr>
+
+	* interpreter/sparse.txi: Add luinc function.
+
 2005-03-09  John W. Eaton  <jwe@octave.org>
 
 	* Makefile.in (bin-dist): Delete target.
--- a/doc/interpreter/sparse.txi	Thu Apr 14 19:35:20 2005 +0000
+++ b/doc/interpreter/sparse.txi	Thu Apr 14 22:17:27 2005 +0000
@@ -872,7 +872,7 @@
 
 @table @asis
 @item License
-UMFPACK Version 4.3 (Jan. 16, 2004), Copyright @copyright{} 2004 by
+UMFPACK Version 4.4 (Jan. 28, 2005), Copyright @copyright{} 2005 by
 Timothy A. Davis.  All Rights Reserved.
 
 Your use or distribution of UMFPACK or any modified version of
@@ -1002,6 +1002,8 @@
 @emph{Not implemented}
 @item gmres
 @emph{Not implemented}
+@item luinc
+Produce the incomplete LU factorization of the sparse matrix A.
 @item lsqr
 @emph{Not implemented}
 @item minres
@@ -1062,6 +1064,8 @@
 		sparse
 * issparse::	Return 1 if the value of the expression EXPR is a sparse
 		matrix.
+* luinc::	Produce the incomplete LU factorization of the sparse 
+		A.
 * nnz:: 	returns number of non zero elements in SM See also: sparse
 * nonzeros::	Returns a vector of the non-zero values of the sparse
 		matrix S
@@ -1138,12 +1142,17 @@
 
 @DOCSTRING(full)
 
-@node issparse, nnz, full, Function Reference
+@node issparse, luinc, full, Function Reference
 @subsubsection issparse
 
 @DOCSTRING(issparse)
 
-@node nnz, nonzeros, issparse, Function Reference
+@node luinc, nnz, issparse, Function Reference
+@subsubsection luinc
+
+@DOCSTRING(luinc)
+
+@node nnz, nonzeros, luinc, Function Reference
 @subsubsection nnz
 
 @DOCSTRING(nnz)
--- a/liboctave/ChangeLog	Thu Apr 14 19:35:20 2005 +0000
+++ b/liboctave/ChangeLog	Thu Apr 14 22:17:27 2005 +0000
@@ -1,3 +1,13 @@
+2005-04-14  David Bateman  <dbateman@free.fr>
+
+	* SparseCmplxLU.cc: Add flags for incomplete factorization.
+	* SparsedbleLU.cc: Ditto.
+	* SparseCmplxLU.h: Definition.
+	* SparsedbleLU.h: ditto.
+
+	* SparseType.cc (transpose): New function.
+	* SparseType.h (transpose): Definition.
+	
 2005-04-11  John W. Eaton  <jwe@octave.org>
 
 	* lo-specfun.cc: Use F77_XFCN instead of F77_FUNC for calls to
--- a/liboctave/SparseCmplxLU.cc	Thu Apr 14 19:35:20 2005 +0000
+++ b/liboctave/SparseCmplxLU.cc	Thu Apr 14 22:17:27 2005 +0000
@@ -223,117 +223,113 @@
 
 SparseComplexLU::SparseComplexLU (const SparseComplexMatrix& a, 
 				  const ColumnVector& Qinit, 
-				  double piv_thres, bool FixedQ)
+				  double piv_thres, bool FixedQ,
+				  double droptol, bool milu, bool udiag)
 {
 #ifdef HAVE_UMFPACK
-  octave_idx_type nr = a.rows ();
-  octave_idx_type nc = a.cols ();
-
-  // Setup the control parameters
-  Matrix Control (UMFPACK_CONTROL, 1);
-  double *control = Control.fortran_vec ();
-  umfpack_zi_defaults (control);
-
-  double tmp = Voctave_sparse_controls.get_key ("spumoni");
-  if (!xisnan (tmp))
-    Control (UMFPACK_PRL) = tmp;
-  if (piv_thres >= 0.)
-    {
-      piv_thres = (piv_thres > 1. ? 1. : piv_thres);
-      Control (UMFPACK_SYM_PIVOT_TOLERANCE) = piv_thres;
-      Control (UMFPACK_PIVOT_TOLERANCE) = piv_thres;
-    }
-  else
-    {
-      tmp = Voctave_sparse_controls.get_key ("piv_tol");
-      if (!xisnan (tmp))
-	{
-	  Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
-	  Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
-	}
-    }
-
-  // Set whether we are allowed to modify Q or not
-  if (FixedQ)
-    Control (UMFPACK_FIXQ) = 1.0;
+  if (milu)
+    (*current_liboctave_error_handler) 
+      ("Modified incomplete LU not implemented");   
   else
     {
-      tmp = Voctave_sparse_controls.get_key ("autoamd");
-      if (!xisnan (tmp))
-	Control (UMFPACK_FIXQ) = tmp;
-    }
+      octave_idx_type nr = a.rows ();
+      octave_idx_type nc = a.cols ();
 
-  // Turn-off UMFPACK scaling for LU 
-  Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
-
-  umfpack_zi_report_control (control);
+      // Setup the control parameters
+      Matrix Control (UMFPACK_CONTROL, 1);
+      double *control = Control.fortran_vec ();
+      umfpack_zi_defaults (control);
 
-  const octave_idx_type *Ap = a.cidx ();
-  const octave_idx_type *Ai = a.ridx ();
-  const Complex *Ax = a.data ();
-
-  umfpack_zi_report_matrix (nr, nc, Ap, Ai, X_CAST (const double *, Ax), NULL,
-			    1, control);
+      double tmp = Voctave_sparse_controls.get_key ("spumoni");
+      if (!xisnan (tmp))
+	Control (UMFPACK_PRL) = tmp;
+      if (piv_thres >= 0.)
+	{
+	  piv_thres = (piv_thres > 1. ? 1. : piv_thres);
+	  Control (UMFPACK_SYM_PIVOT_TOLERANCE) = piv_thres;
+	  Control (UMFPACK_PIVOT_TOLERANCE) = piv_thres;
+	}
+      else
+	{
+	  tmp = Voctave_sparse_controls.get_key ("piv_tol");
+	  if (!xisnan (tmp))
+	    {
+	      Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
+	      Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
+	    }
+	}
 
-  void *Symbolic;
-  Matrix Info (1, UMFPACK_INFO);
-  double *info = Info.fortran_vec ();
-  int status;
-
-  // Null loop so that qinit is imediately deallocated when not needed
-  do {
-    OCTAVE_LOCAL_BUFFER (int, qinit, nc);
+      if (droptol >= 0.)
+	Control (UMFPACK_DROPTOL) = droptol;
 
-    for (int i = 0; i < nc; i++)
-      qinit [i] = static_cast<int> (Qinit (i));
+      // Set whether we are allowed to modify Q or not
+      if (FixedQ)
+	Control (UMFPACK_FIXQ) = 1.0;
+      else
+	{
+	  tmp = Voctave_sparse_controls.get_key ("autoamd");
+	  if (!xisnan (tmp))
+	    Control (UMFPACK_FIXQ) = tmp;
+	}
 
-    status = umfpack_zi_qsymbolic (nr, nc, Ap, Ai, X_CAST (const double *, Ax),
-				   NULL, qinit, &Symbolic, control, info);
-  } while (0);
+      // Turn-off UMFPACK scaling for LU 
+      Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
 
-  if (status < 0)
-    {
-      (*current_liboctave_error_handler) 
-	    ("SparseComplexLU::SparseComplexLU symbolic factorization failed");
+      umfpack_zi_report_control (control);
+
+      const octave_idx_type *Ap = a.cidx ();
+      const octave_idx_type *Ai = a.ridx ();
+      const Complex *Ax = a.data ();
 
-      umfpack_zi_report_status (control, status);
-      umfpack_zi_report_info (control, info);
+      umfpack_zi_report_matrix (nr, nc, Ap, Ai, 
+				X_CAST (const double *, Ax), NULL,
+				1, control);
+
+      void *Symbolic;
+      Matrix Info (1, UMFPACK_INFO);
+      double *info = Info.fortran_vec ();
+      int status;
 
-      umfpack_zi_free_symbolic (&Symbolic) ;
-    }
-  else
-    {
-      umfpack_zi_report_symbolic (Symbolic, control);
+      // Null loop so that qinit is imediately deallocated when not
+      // needed
+      do {
+	OCTAVE_LOCAL_BUFFER (int, qinit, nc);
 
-      void *Numeric;
-      status = umfpack_zi_numeric (Ap, Ai, X_CAST (const double *, Ax), NULL,
-				   Symbolic, &Numeric, control, info) ;
-      umfpack_zi_free_symbolic (&Symbolic) ;
+	for (int i = 0; i < nc; i++)
+	  qinit [i] = static_cast<int> (Qinit (i));
 
-      cond = Info (UMFPACK_RCOND);
+	status = umfpack_zi_qsymbolic (nr, nc, Ap, Ai, 
+				       X_CAST (const double *, Ax),
+				       NULL, qinit, &Symbolic, control, 
+				       info);
+      } while (0);
 
       if (status < 0)
 	{
 	  (*current_liboctave_error_handler) 
-	    ("SparseComplexLU::SparseComplexLU numeric factorization failed");
+	    ("SparseComplexLU::SparseComplexLU symbolic factorization failed");
 
 	  umfpack_zi_report_status (control, status);
 	  umfpack_zi_report_info (control, info);
 
-	  umfpack_zi_free_numeric (&Numeric);
+	  umfpack_zi_free_symbolic (&Symbolic) ;
 	}
       else
 	{
-	  umfpack_zi_report_numeric (Numeric, control);
+	  umfpack_zi_report_symbolic (Symbolic, control);
 
-	  int lnz, unz, ignore1, ignore2, ignore3;
-	  status = umfpack_zi_get_lunz (&lnz, &unz, &ignore1, &ignore2,
-					&ignore3, Numeric) ;
-	  
+	  void *Numeric;
+	  status = umfpack_zi_numeric (Ap, Ai, 
+				       X_CAST (const double *, Ax), NULL,
+				       Symbolic, &Numeric, control, info) ;
+	  umfpack_zi_free_symbolic (&Symbolic) ;
+
+	  cond = Info (UMFPACK_RCOND);
+
 	  if (status < 0)
 	    {
 	      (*current_liboctave_error_handler) 
-		("SparseComplexLU::SparseComplexLU extracting LU factors failed");
+		("SparseComplexLU::SparseComplexLU numeric factorization failed");
 
 	      umfpack_zi_report_status (control, status);
 	      umfpack_zi_report_info (control, info);
@@ -342,71 +338,102 @@
 	    }
 	  else
 	    {
-	      int n_inner = (nr < nc ? nr : nc);
-
-	      if (lnz < 1)
-		Lfact = SparseComplexMatrix (static_cast<octave_idx_type> (n_inner), nr,
-					     static_cast<octave_idx_type> (1));
-	      else
-		Lfact = SparseComplexMatrix (static_cast<octave_idx_type> (n_inner), nr,
-					     static_cast<octave_idx_type> (lnz));
-
-	      octave_idx_type *Ltp = Lfact.cidx ();
-	      octave_idx_type *Ltj = Lfact.ridx ();
-	      Complex *Ltx = Lfact.data ();
-
-	      if (unz < 1)
-		Ufact = SparseComplexMatrix (static_cast<octave_idx_type> (n_inner), nc,
-					     static_cast<octave_idx_type> (1));
-	      else
-		Ufact = SparseComplexMatrix (static_cast<octave_idx_type> (n_inner), nc, unz);
+	      umfpack_zi_report_numeric (Numeric, control);
 
-	      octave_idx_type *Up = Ufact.cidx ();
-	      octave_idx_type *Uj = Ufact.ridx ();
-	      Complex *Ux = Ufact.data ();
-	      
-	      P.resize (nr);
-	      int *p = P.fortran_vec ();
-
-	      Q.resize (nc);
-	      int *q = Q.fortran_vec ();
-
-	      int do_recip;
-	      status = umfpack_zi_get_numeric (Ltp, Ltj, X_CAST (double *, Ltx),
-					       NULL, Up, Uj,
-					       X_CAST (double *, Ux), NULL, p, 
-					       q, NULL, NULL, &do_recip,
-					       NULL, Numeric) ;
-
-	      umfpack_zi_free_numeric (&Numeric) ;
-
-	      if (status < 0 || do_recip)
+	      int lnz, unz, ignore1, ignore2, ignore3;
+	      status = umfpack_zi_get_lunz (&lnz, &unz, &ignore1, 
+					    &ignore2, &ignore3, Numeric);
+	  
+	      if (status < 0)
 		{
 		  (*current_liboctave_error_handler) 
 		    ("SparseComplexLU::SparseComplexLU extracting LU factors failed");
 
 		  umfpack_zi_report_status (control, status);
+		  umfpack_zi_report_info (control, info);
+
+		  umfpack_zi_free_numeric (&Numeric);
 		}
 	      else
 		{
-		  Lfact = Lfact.transpose ();
+		  int n_inner = (nr < nc ? nr : nc);
+
+		  if (lnz < 1)
+		    Lfact = SparseComplexMatrix 
+		      (static_cast<octave_idx_type> (n_inner), nr,
+		       static_cast<octave_idx_type> (1));
+		  else
+		    Lfact = SparseComplexMatrix 
+		      (static_cast<octave_idx_type> (n_inner), nr,
+		       static_cast<octave_idx_type> (lnz));
+
+		  octave_idx_type *Ltp = Lfact.cidx ();
+		  octave_idx_type *Ltj = Lfact.ridx ();
+		  Complex *Ltx = Lfact.data ();
 
-		  umfpack_zi_report_matrix (nr, n_inner, Lfact.cidx (), 
-					    Lfact.ridx (), 
-					    X_CAST (double *, Lfact.data()), 
-					    NULL, 1, control);
+		  if (unz < 1)
+		    Ufact = SparseComplexMatrix 
+		      (static_cast<octave_idx_type> (n_inner), nc,
+		       static_cast<octave_idx_type> (1));
+		  else
+		    Ufact = SparseComplexMatrix 
+		      (static_cast<octave_idx_type> (n_inner), nc, unz);
+
+		  octave_idx_type *Up = Ufact.cidx ();
+		  octave_idx_type *Uj = Ufact.ridx ();
+		  Complex *Ux = Ufact.data ();
+	      
+		  P.resize (nr);
+		  int *p = P.fortran_vec ();
+
+		  Q.resize (nc);
+		  int *q = Q.fortran_vec ();
 
-		  umfpack_zi_report_matrix (n_inner, nc, Ufact.cidx (), 
-					    Ufact.ridx (), 
-					    X_CAST (double *, Ufact.data()), 
-					    NULL, 1, control);
-		  umfpack_zi_report_perm (nr, p, control);
-		  umfpack_zi_report_perm (nc, q, control);
+		  int do_recip;
+		  status = 
+		    umfpack_zi_get_numeric (Ltp, Ltj, 
+					    X_CAST (double *, Ltx),
+					    NULL, Up, Uj,
+					    X_CAST (double *, Ux), 
+					    NULL, p, q, NULL, NULL, 
+					    &do_recip, NULL, Numeric) ;
+
+		  umfpack_zi_free_numeric (&Numeric) ;
+
+		  if (status < 0 || do_recip)
+		    {
+		      (*current_liboctave_error_handler) 
+			("SparseComplexLU::SparseComplexLU extracting LU factors failed");
+
+		      umfpack_zi_report_status (control, status);
+		    }
+		  else
+		    {
+		      Lfact = Lfact.transpose ();
+
+		      umfpack_zi_report_matrix (nr, n_inner, 
+						Lfact.cidx (), 
+						Lfact.ridx (), 
+						X_CAST (double *, Lfact.data()), 
+						NULL, 1, control);
+
+		      umfpack_zi_report_matrix (n_inner, nc, 
+						Ufact.cidx (), 
+						Ufact.ridx (), 
+						X_CAST (double *, Ufact.data()), 
+						NULL, 1, control);
+		      umfpack_zi_report_perm (nr, p, control);
+		      umfpack_zi_report_perm (nc, q, control);
+		    }
+
+		  umfpack_zi_report_info (control, info);
 		}
-
-	      umfpack_zi_report_info (control, info);
 	    }
 	}
+
+      if (udiag)
+	(*current_liboctave_error_handler) 
+	  ("Option udiag of incomplete LU not implemented");   
     }
 #else
   (*current_liboctave_error_handler) ("UMFPACK not installed");
--- a/liboctave/SparseCmplxLU.h	Thu Apr 14 19:35:20 2005 +0000
+++ b/liboctave/SparseCmplxLU.h	Thu Apr 14 22:17:27 2005 +0000
@@ -38,7 +38,9 @@
   SparseComplexLU (const SparseComplexMatrix& a, double piv_thres = -1);
 
   SparseComplexLU (const SparseComplexMatrix& a, const ColumnVector& Qinit,
-		   double piv_thres = -1, bool FixedQ = false);
+		   double piv_thres = -1, bool FixedQ = false,
+		   double droptol = -1., bool milu = false,
+		   bool udiag = false);
 
   SparseComplexLU (const SparseComplexLU& a) 
     : sparse_base_lu <SparseComplexMatrix, Complex, SparseMatrix, double> (a) { }
--- a/liboctave/SparseType.cc	Thu Apr 14 19:35:20 2005 +0000
+++ b/liboctave/SparseType.cc	Thu Apr 14 22:17:27 2005 +0000
@@ -688,6 +688,37 @@
     typ = SparseType::Lower;
 }
 
+SparseType
+SparseType::transpose (void) const
+{
+  SparseType retval (*this);
+  if (typ == SparseType::Upper)
+    retval.typ = Lower;
+  else if (typ == SparseType::Permuted_Upper)
+    {
+      int *tmp = retval.row_perm;
+      retval.row_perm = retval.col_perm;
+      retval.col_perm = tmp;
+      retval.typ = Lower;
+    }
+  else if (typ == SparseType::Lower)
+    retval.typ = Upper;
+  else if (typ == SparseType::Permuted_Upper)
+    {
+      int *tmp = retval.row_perm;
+      retval.row_perm = retval.col_perm;
+      retval.col_perm = tmp;
+      retval.typ = Upper;
+    }
+  else if (typ == SparseType::Banded)
+    {
+      retval.upper_band = lower_band;
+      retval.lower_band = upper_band;
+    }
+
+  return retval;
+}
+
 /*
 ;;; Local Variables: ***
 ;;; mode: C++ ***
--- a/liboctave/SparseType.h	Thu Apr 14 19:35:20 2005 +0000
+++ b/liboctave/SparseType.h	Thu Apr 14 22:17:27 2005 +0000
@@ -132,6 +132,8 @@
 
   void mark_as_unpermuted (void);
 
+  SparseType transpose (void) const;
+
 private:
   void type (int new_typ) { typ = static_cast<matrix_type>(new_typ); }
 
--- a/liboctave/SparsedbleLU.cc	Thu Apr 14 19:35:20 2005 +0000
+++ b/liboctave/SparsedbleLU.cc	Thu Apr 14 22:17:27 2005 +0000
@@ -217,116 +217,108 @@
 }
 
 SparseLU::SparseLU (const SparseMatrix& a, const ColumnVector& Qinit,
-		    double piv_thres, bool FixedQ)
+		    double piv_thres, bool FixedQ, double droptol,
+		    bool milu, bool udiag)
 {
 #ifdef HAVE_UMFPACK
-  octave_idx_type nr = a.rows ();
-  octave_idx_type nc = a.cols ();
-
-  // Setup the control parameters
-  Matrix Control (UMFPACK_CONTROL, 1);
-  double *control = Control.fortran_vec ();
-  umfpack_di_defaults (control);
-
-  double tmp = Voctave_sparse_controls.get_key ("spumoni");
-  if (!xisnan (tmp))
-    Control (UMFPACK_PRL) = tmp;
-  if (piv_thres >= 0.)
-    {
-      piv_thres = (piv_thres > 1. ? 1. : piv_thres);
-      Control (UMFPACK_SYM_PIVOT_TOLERANCE) = piv_thres;
-      Control (UMFPACK_PIVOT_TOLERANCE) = piv_thres;
-    }
-  else
-    {
-      tmp = Voctave_sparse_controls.get_key ("piv_tol");
-      if (!xisnan (tmp))
-	{
-	  Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
-	  Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
-	}
-    }
-
-  // Set whether we are allowed to modify Q or not
-  if (FixedQ)
-    Control (UMFPACK_FIXQ) = 1.0;
+  if (milu)
+    (*current_liboctave_error_handler) 
+      ("Modified incomplete LU not implemented");   
   else
     {
-      tmp = Voctave_sparse_controls.get_key ("autoamd");
-      if (!xisnan (tmp))
-	Control (UMFPACK_FIXQ) = tmp;
-    }
+      octave_idx_type nr = a.rows ();
+      octave_idx_type nc = a.cols ();
 
-  // Turn-off UMFPACK scaling for LU 
-  Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
-
-  umfpack_di_report_control (control);
+      // Setup the control parameters
+      Matrix Control (UMFPACK_CONTROL, 1);
+      double *control = Control.fortran_vec ();
+      umfpack_di_defaults (control);
 
-  const octave_idx_type *Ap = a.cidx ();
-  const octave_idx_type *Ai = a.ridx ();
-  const double *Ax = a.data ();
-
-  umfpack_di_report_matrix (nr, nc, Ap, Ai, Ax, 1, control);
+      double tmp = Voctave_sparse_controls.get_key ("spumoni");
+      if (!xisnan (tmp))
+	Control (UMFPACK_PRL) = tmp;
+      if (piv_thres >= 0.)
+	{
+	  piv_thres = (piv_thres > 1. ? 1. : piv_thres);
+	  Control (UMFPACK_SYM_PIVOT_TOLERANCE) = piv_thres;
+	  Control (UMFPACK_PIVOT_TOLERANCE) = piv_thres;
+	}
+      else
+	{
+	  tmp = Voctave_sparse_controls.get_key ("piv_tol");
+	  if (!xisnan (tmp))
+	    {
+	      Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
+	      Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
+	    }
+	}
 
-  void *Symbolic;
-  Matrix Info (1, UMFPACK_INFO);
-  double *info = Info.fortran_vec ();
-  int status;
+      if (droptol >= 0.)
+	Control (UMFPACK_DROPTOL) = droptol;
 
-  // Null loop so that qinit is imediately deallocated when not needed
-  do {
-    OCTAVE_LOCAL_BUFFER (int, qinit, nc);
 
-    for (int i = 0; i < nc; i++)
-      qinit [i] = static_cast<int> (Qinit (i));
+      // Set whether we are allowed to modify Q or not
+      if (FixedQ)
+	Control (UMFPACK_FIXQ) = 1.0;
+      else
+	{
+	  tmp = Voctave_sparse_controls.get_key ("autoamd");
+	  if (!xisnan (tmp))
+	    Control (UMFPACK_FIXQ) = tmp;
+	}
 
-    status = umfpack_di_qsymbolic (nr, nc, Ap, Ai, Ax, qinit,
-				   &Symbolic, control, info);
-  } while (0);
+      // Turn-off UMFPACK scaling for LU 
+      Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
 
-  if (status < 0)
-    {
-      (*current_liboctave_error_handler) 
-	    ("SparseLU::SparseLU symbolic factorization failed");
+      umfpack_di_report_control (control);
 
-      umfpack_di_report_status (control, status);
-      umfpack_di_report_info (control, info);
+      const octave_idx_type *Ap = a.cidx ();
+      const octave_idx_type *Ai = a.ridx ();
+      const double *Ax = a.data ();
+
+      umfpack_di_report_matrix (nr, nc, Ap, Ai, Ax, 1, control);
+
+      void *Symbolic;
+      Matrix Info (1, UMFPACK_INFO);
+      double *info = Info.fortran_vec ();
+      int status;
 
-      umfpack_di_free_symbolic (&Symbolic) ;
-    }
-  else
-    {
-      umfpack_di_report_symbolic (Symbolic, control);
+      // Null loop so that qinit is imediately deallocated when not needed
+      do {
+	OCTAVE_LOCAL_BUFFER (int, qinit, nc);
 
-      void *Numeric;
-      status = umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric,
-				   control, info) ;
-      umfpack_di_free_symbolic (&Symbolic) ;
+	for (int i = 0; i < nc; i++)
+	  qinit [i] = static_cast<int> (Qinit (i));
 
-      cond = Info (UMFPACK_RCOND);
+	status = umfpack_di_qsymbolic (nr, nc, Ap, Ai, Ax, qinit,
+				       &Symbolic, control, info);
+      } while (0);
 
       if (status < 0)
 	{
 	  (*current_liboctave_error_handler) 
-	    ("SparseLU::SparseLU numeric factorization failed");
+	    ("SparseLU::SparseLU symbolic factorization failed");
 
 	  umfpack_di_report_status (control, status);
 	  umfpack_di_report_info (control, info);
 
-	  umfpack_di_free_numeric (&Numeric);
+	  umfpack_di_free_symbolic (&Symbolic) ;
 	}
       else
 	{
-	  umfpack_di_report_numeric (Numeric, control);
+	  umfpack_di_report_symbolic (Symbolic, control);
 
-	  int lnz, unz, ignore1, ignore2, ignore3;
-	  status = umfpack_di_get_lunz (&lnz, &unz, &ignore1, &ignore2,
-					&ignore3, Numeric) ;
-	  
+	  void *Numeric;
+	  status = umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric,
+				       control, info) ;
+	  umfpack_di_free_symbolic (&Symbolic) ;
+
+	  cond = Info (UMFPACK_RCOND);
+
 	  if (status < 0)
 	    {
 	      (*current_liboctave_error_handler) 
-		("SparseLU::SparseLU extracting LU factors failed");
+		("SparseLU::SparseLU numeric factorization failed");
 
 	      umfpack_di_report_status (control, status);
 	      umfpack_di_report_info (control, info);
@@ -335,67 +327,100 @@
 	    }
 	  else
 	    {
-	      int n_inner = (nr < nc ? nr : nc);
-
-	      if (lnz < 1)
-		Lfact = SparseMatrix (static_cast<octave_idx_type> (n_inner), nr,
-				      static_cast<octave_idx_type> (1));
-	      else
-		Lfact = SparseMatrix (static_cast<octave_idx_type> (n_inner), nr,
-				      static_cast<octave_idx_type> (lnz));
-
-	      octave_idx_type *Ltp = Lfact.cidx ();
-	      octave_idx_type *Ltj = Lfact.ridx ();
-	      double *Ltx = Lfact.data ();
-
-	      if (unz < 1)
-		Ufact = SparseMatrix (static_cast<octave_idx_type> (n_inner), nc,
-				      static_cast<octave_idx_type> (1));
-	      else
-		Ufact = SparseMatrix (static_cast<octave_idx_type> (n_inner), nc,
-				      static_cast<octave_idx_type> (unz));
+	      umfpack_di_report_numeric (Numeric, control);
 
-	      octave_idx_type *Up = Ufact.cidx ();
-	      octave_idx_type *Uj = Ufact.ridx ();
-	      double *Ux = Ufact.data ();
-
-	      P.resize (nr);
-	      int *p = P.fortran_vec ();
-
-	      Q.resize (nc);
-	      int *q = Q.fortran_vec ();
-
-	      int do_recip;
-	      status = umfpack_di_get_numeric (Ltp, Ltj, Ltx, Up, Uj,
-					       Ux, p, q, (double *) NULL,
-					       &do_recip, (double *) NULL, 
-					       Numeric) ;
-
-	      umfpack_di_free_numeric (&Numeric) ;
-
-	      if (status < 0 || do_recip)
+	      int lnz, unz, ignore1, ignore2, ignore3;
+	      status = umfpack_di_get_lunz (&lnz, &unz, &ignore1, &ignore2,
+					    &ignore3, Numeric) ;
+	  
+	      if (status < 0)
 		{
 		  (*current_liboctave_error_handler) 
 		    ("SparseLU::SparseLU extracting LU factors failed");
 
 		  umfpack_di_report_status (control, status);
+		  umfpack_di_report_info (control, info);
+
+		  umfpack_di_free_numeric (&Numeric);
 		}
 	      else
 		{
-		  Lfact = Lfact.transpose ();
-		  umfpack_di_report_matrix (nr, n_inner, Lfact.cidx (), 
-					    Lfact.ridx (), Lfact.data (),
-					    1, control);
-		  umfpack_di_report_matrix (n_inner, nc, Ufact.cidx (), 
-					    Ufact.ridx (), Ufact.data (),
-					    1, control);
-		  umfpack_di_report_perm (nr, p, control);
-		  umfpack_di_report_perm (nc, q, control);
+		  int n_inner = (nr < nc ? nr : nc);
+
+		  if (lnz < 1)
+		    Lfact = SparseMatrix 
+		      (static_cast<octave_idx_type> (n_inner), nr,
+		       static_cast<octave_idx_type> (1));
+		  else
+		    Lfact = SparseMatrix 
+		      (static_cast<octave_idx_type> (n_inner), nr,
+		       static_cast<octave_idx_type> (lnz));
+
+		  octave_idx_type *Ltp = Lfact.cidx ();
+		  octave_idx_type *Ltj = Lfact.ridx ();
+		  double *Ltx = Lfact.data ();
+
+		  if (unz < 1)
+		    Ufact = SparseMatrix 
+		      (static_cast<octave_idx_type> (n_inner), nc,
+		       static_cast<octave_idx_type> (1));
+		  else
+		    Ufact = SparseMatrix 
+		      (static_cast<octave_idx_type> (n_inner), nc,
+		       static_cast<octave_idx_type> (unz));
+
+		  octave_idx_type *Up = Ufact.cidx ();
+		  octave_idx_type *Uj = Ufact.ridx ();
+		  double *Ux = Ufact.data ();
+
+		  P.resize (nr);
+		  int *p = P.fortran_vec ();
+
+		  Q.resize (nc);
+		  int *q = Q.fortran_vec ();
+
+		  int do_recip;
+		  status = umfpack_di_get_numeric (Ltp, Ltj, Ltx, Up, Uj,
+						   Ux, p, q, 
+						   (double *) NULL,
+						   &do_recip, 
+						   (double *) NULL, 
+						   Numeric) ;
+
+		  umfpack_di_free_numeric (&Numeric) ;
+
+		  if (status < 0 || do_recip)
+		    {
+		      (*current_liboctave_error_handler) 
+			("SparseLU::SparseLU extracting LU factors failed");
+
+		      umfpack_di_report_status (control, status);
+		    }
+		  else
+		    {
+		      Lfact = Lfact.transpose ();
+		      umfpack_di_report_matrix (nr, n_inner, 
+						Lfact.cidx (), 
+						Lfact.ridx (), 
+						Lfact.data (),
+						1, control);
+		      umfpack_di_report_matrix (n_inner, nc, 
+						Ufact.cidx (), 
+						Ufact.ridx (), 
+						Ufact.data (),
+						1, control);
+		      umfpack_di_report_perm (nr, p, control);
+		      umfpack_di_report_perm (nc, q, control);
+		    }
+
+		  umfpack_di_report_info (control, info);
 		}
-
-	      umfpack_di_report_info (control, info);
 	    }
 	}
+
+      if (udiag)
+	(*current_liboctave_error_handler) 
+	  ("Option udiag of incomplete LU not implemented");   
     }
 #else
   (*current_liboctave_error_handler) ("UMFPACK not installed");
--- a/liboctave/SparsedbleLU.h	Thu Apr 14 19:35:20 2005 +0000
+++ b/liboctave/SparsedbleLU.h	Thu Apr 14 22:17:27 2005 +0000
@@ -36,7 +36,8 @@
   SparseLU (const SparseMatrix& a, double piv_thres = -1.0);
 
   SparseLU (const SparseMatrix& a, const ColumnVector& Qinit, 
-	    double piv_thres = -1.0, bool FixedQ = false);
+	    double piv_thres = -1.0, bool FixedQ = false,
+	    double droptol = -1., bool milu = false, bool udiag = false);
 
   SparseLU (const SparseLU& a) 
     : sparse_base_lu <SparseMatrix, double, SparseMatrix, double> (a) { }
--- a/src/ChangeLog	Thu Apr 14 19:35:20 2005 +0000
+++ b/src/ChangeLog	Thu Apr 14 22:17:27 2005 +0000
@@ -1,3 +1,27 @@
+2005-04-14  David Bateman  <dbateman@free.fr>
+
+	* Makefile.in: Add luinc.cc to DLD_XSRC.
+	* DLD-FUNCTIONS/luinc.cc: New file for incomplete LU factorization.
+
+	* ov-bool-sparse.h (index_vector): New function.
+	* ov-re-sparse.cc (index_vector): Ditto.
+	* ov-re-sparse.h (index_vector): Definition.
+
+	* ov.mapper (any_element_less_than, any_element_greater_than):
+	New versions for SparseMatrix
+	(SPARSE_MAPPER_LOOP_2, SPARSE_MAPPER_LOOP_1, SPARSE_MAPPER_LOOP):
+	New macros.
+	(apply): Add special cases for sparse arguments to the mapper
+	functions
+
+	* ov-re-sparse.cc (streamoff_array_value): Use octave_idx_type.
+	(convert_to_str_internal): New function.
+	* ov-re-sparse.h (convert_to_str_internal): Definition.
+	
+	* DLD-FUNCTIONS/sparse.cc (Fsparse): More care for nargin=2 case.
+
+	* DLD-FUNCTIONS/splu.cc (Fsplu): Use octave_idx_type.
+	
 2005-04-14  John W. Eaton  <jwe@octave.org>
 
 	* strfns.cc (Fchar): If arg is a dq string, return a dq string.
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/DLD-FUNCTIONS/luinc.cc	Thu Apr 14 22:17:27 2005 +0000
@@ -0,0 +1,364 @@
+/*
+
+Copyright (C) 2005 David Bateman
+
+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 2, 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 this program; see the file COPYING.  If not, write to the Free
+Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "defun-dld.h"
+#include "error.h"
+#include "gripes.h"
+#include "oct-obj.h"
+#include "utils.h"
+#include "oct-map.h"
+
+#include "SparseCmplxLU.h"
+#include "SparsedbleLU.h"
+#include "ov-re-sparse.h"
+#include "ov-cx-sparse.h"
+
+DEFUN_DLD (luinc, args, nargout,
+  "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {[@var{l}, @var{u}, @var{p}, @var{q}] =} luinc (@var{a}, '0')\n\
+@deftypefnx {Loadable Function} {[@var{l}, @var{u}, @var{p}, @var{q}] =} luinc (@var{a}, @var{droptol})\n\
+@deftypefnx {Loadable Function} {[@var{l}, @var{u}, @var{p}, @var{q}] =} luinc (@var{a}, @var{opts})\n\
+@cindex LU decomposition\n\
+Produce the incomplete LU factorization of the sparse matrix @var{a}.\n\
+Two types of incomplete factorization are possible, and the type\n\
+is determined by the second argument to @dfn{luinc}.\n\
+\n\
+Called with a second argument of '0', the zero-level incomplete\n\
+LU factorization is produced. This creates a factorization of @var{a}\n\
+where the position of the non-zero arguments correspond to the same\n\
+positions as in the matrix @var{a}.\n\
+\n\
+Alternatively, the fill-in of the incomplete LU factorization can\n\
+be controlled through the variable @var{droptol} or the structure\n\
+@var{opts}. The UMFPACK multifrontal factorization code by Tim A.\n\
+Davis is used for the incomplete LU factorication, (availability\n\
+@url{http://www.cise.ufl.edu/research/sparse/umfpack/})\n\
+\n\
+@var{droptol} determines the values below which the values in the LU\n\
+factorization are dropped and replaced by zero. It must be a positive\n\
+scalar, and any values in the factorization whose absolute value are\n\
+less than this value are dropped, expect if leaving them increase the\n\
+sparsity of the matrix. Setting @var{droptol} to zero results in a\n\
+complete LU factorization which is the default.\n\
+\n\
+@var{opts} is a structure containing one or more of the fields\n\
+\n\
+@table @code\n\
+@item droptol\n\
+The drop tolerance as above. If @var{opts} only contains @code{droptol}\n\
+then this is equivalent to using the variable @var{droptol}.\n\
+\n\
+@item milu\n\
+A logical variable flagging whether to use the modified incomplete LU\n\
+factorization. In the case that @code{milu} is true, the dropped values\n\
+are subtract from the diagonal of the matrix U of the factorization.\n\
+The default is @code{false}.\n\
+\n\
+@item udiag\n\
+A logical variable that flags whether zero elements on the diagonal of U\n\
+should be replaced with @var{droptol} to attempt to avoid singular\n\
+factors. The default is @code{false}.\n\
+\n\
+@item thresh\n\
+Defines the pivot threshold in the interval [0,1]. Values outside that\n\
+range are ignored.\n\
+@end table\n\
+\n\
+All other fields in @var{opts} are ignored. The outputs from @dfn{luinc}\n\
+are the same as for @dfn{lu}.\n\
+@end deftypefn\n\
+@seealso{sparse, lu, cholinc}")
+{
+  int nargin = args.length ();
+  octave_value_list retval;
+
+  if (nargin == 0)
+    print_usage ("luinc");
+  else if (nargin != 2)
+    error ("luinc: incorrect number of arguments");
+  else
+    {
+      bool zero_level = false;
+      bool milu = false;
+      bool udiag = false;
+      bool thresh = -1;
+      double droptol = -1.;
+
+      if (args(1).is_string ())
+	{
+	  if (args(1).string_value () == "0")
+	    zero_level = true;
+	  else
+	    error ("luinc: unrecognized string argument");
+	}
+      else if (args(1).is_map ())
+	{
+	  Octave_map map = args(1).map_value ();
+
+	  if (map.contains ("droptol"))
+	    droptol = map.contents ("droptol")(0).double_value ();
+
+	  if (map.contains ("milu"))
+	    {
+	      double tmp = map.contents ("milu")(0).double_value ();
+
+	      milu = (tmp == 0. ? false : true);
+	    }
+
+	  if (map.contains ("udiag"))
+	    {
+	      double tmp = map.contents ("udiag")(0).double_value ();
+
+	      udiag = (tmp == 0. ? false : true);
+	    }
+
+	  if (map.contains ("thresh"))
+	    thresh = map.contents ("thresh")(0).double_value ();
+	}
+      else
+	droptol = args(1).double_value ();
+
+      // XXX FIXME XXX Add code for zero-level factorization
+      if (zero_level)
+	error ("luinc: zero-level factorization not implemented");
+
+      if (!error_state)
+	{
+	  if (args(0).class_name () == "sparse") 
+	    {
+	      SparseMatrix sm = args(0).sparse_matrix_value ();
+	      octave_idx_type sm_nc = sm.cols ();
+	      ColumnVector Qinit (sm_nc);
+
+	      for (octave_idx_type i = 0; i < sm_nc; i++)
+		Qinit (i) = i;
+
+	      if (! error_state)
+		{
+		  switch (nargout)
+		    {
+		    case 0:
+		    case 1:
+		    case 2:
+		      {
+			SparseLU fact (sm, Qinit, thresh, true, droptol,
+				       milu, udiag);
+
+			SparseMatrix P = fact.Pr ();
+			SparseMatrix L = P.transpose () * fact.L ();
+			retval(1) = fact.U ();
+			retval(0) = L;
+		      }
+		      break;
+
+		    case 3:
+		      {
+			SparseLU fact (sm, Qinit, thresh, true, droptol,
+				       milu, udiag);
+
+			retval(2) = fact.Pr ();
+			retval(1) = fact.U ();
+			retval(0) = fact.L ();
+		      }
+		      break;
+
+		    case 4:
+		    default:
+		      {
+			SparseLU fact (sm, Qinit, thresh, false, droptol,
+				       milu, udiag);
+
+			retval(3) = fact.Pc ();
+			retval(2) = fact.Pr ();
+			retval(1) = fact.U ();
+			retval(0) = fact.L ();
+		      }
+		      break;
+		    }
+		}
+	    }
+	  else if (args(0).class_name () == "sparse complex") 
+	    {
+	      SparseComplexMatrix sm = 
+		args(0).sparse_complex_matrix_value ();
+	      octave_idx_type sm_nc = sm.cols ();
+	      ColumnVector Qinit (sm_nc);
+
+	      for (octave_idx_type i = 0; i < sm_nc; i++)
+		Qinit (i) = i;
+
+	      if (! error_state)
+		{
+		  switch (nargout)
+		    {
+		    case 0:
+		    case 1:
+		    case 2:
+		      {
+			SparseComplexLU fact (sm, Qinit, thresh, true, 
+					      droptol, milu, udiag);
+
+			SparseMatrix P = fact.Pr ();
+			SparseComplexMatrix L = P.transpose () * fact.L ();
+			retval(1) = fact.U ();
+			retval(0) = L;
+		      }
+		      break;
+
+		    case 3:
+		      {
+			SparseComplexLU fact (sm, Qinit, thresh, true,
+					      droptol, milu, udiag);
+
+			retval(2) = fact.Pr ();
+			retval(1) = fact.U ();
+			retval(0) = fact.L ();
+		      }
+		      break;
+
+		    case 4:
+		    default:
+		      {
+			SparseComplexLU fact (sm, Qinit, thresh, false,
+					      droptol, milu, udiag);
+
+			retval(3) = fact.Pc ();
+			retval(2) = fact.Pr ();
+			retval(1) = fact.U ();
+			retval(0) = fact.L ();
+		      }
+		      break;
+		    }
+		}
+	    }
+	  else
+	    error ("luinc: first argument must be sparse");
+	}
+    }
+
+  return retval;
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
+
+/* 
+ LUINC  Sparse Incomplete LU factorization.
+    LUINC produces two different kinds of incomplete LU factorizations -- the
+    drop tolerance and the 0 level of fill-in factorizations.  These factors
+    may be useful as preconditioners for a system of linear equations being
+    solved by an iterative method such as BICG (BiConjugate Gradients).
+ 
+    LUINC(X,DROPTOL) performs the incomplete LU factorization of
+    X with drop tolerance DROPTOL.
+ 
+    LUINC(X,OPTS) allows additional options to the incomplete LU
+    factorization.  OPTS is a structure with up to four fields:
+        droptol --- the drop tolerance of incomplete LU
+        milu    --- modified incomplete LU
+        udiag   --- replace zeros on the diagonal of U
+        thresh  --- the pivot threshold (see also LU)
+ 
+    Only the fields of interest need to be set.
+ 
+    droptol is a non-negative scalar used as the drop
+    tolerance for the incomplete LU factorization.  This factorization
+    is computed in the same (column-oriented) manner as the
+    LU factorization except after each column of L and U has
+    been calculated, all entries in that column which are smaller
+    in magnitude than the local drop tolerance, which is 
+    droptol * NORM of the column of X, are "dropped" from L or U.
+    The only exception to this dropping rule is the diagonal of the
+    upper triangular factor U which remains even if it is too small.
+    Note that entries of the lower triangular factor L are tested
+    before being scaled by the pivot.  Setting droptol = 0
+    produces the complete LU factorization, which is the default.
+ 
+    milu stands for modified incomplete LU factorization.  Its
+    value is either 0 (unmodified, the default) or 1 (modified).
+    Instead of discarding those entries from the newly-formed
+    column of the factors, they are subtracted from the diagonal
+    of the upper triangular factor, U.
+ 
+    udiag is either 0 or 1.  If it is 1, any zero diagonal entries
+    of the upper triangular factor U are replaced by the local drop
+    tolerance in an attempt to avoid a singular factor.  The default
+    is 0.
+ 
+    thresh is a pivot threshold in [0,1].  Pivoting occurs
+    when the diagonal entry in a column has magnitude less
+    than thresh times the magnitude of any sub-diagonal entry in
+    that column.  thresh = 0 forces diagonal pivoting.  thresh = 1 is
+    the default.
+ 
+    Example:
+ 
+       load west0479
+       A = west0479;
+       nnz(A)
+       nnz(lu(A))
+       nnz(luinc(A,1e-6))
+ 
+       This shows that A has 1887 nonzeros, its complete LU factorization
+       has 16777 nonzeros, and its incomplete LU factorization with a
+       drop tolerance of 1e-6 has 10311 nonzeros.
+ 
+ 
+    [L,U,P] = LUINC(X,'0') produces the incomplete LU factors of a sparse
+    matrix with 0 level of fill-in (i.e. no fill-in).  L is unit lower
+    trianglar, U is upper triangular and P is a permutation matrix.  U has the
+    same sparsity pattern as triu(P*X).  L has the same sparsity pattern as
+    tril(P*X), except for 1's on the diagonal of L where P*X may be zero.  Both
+    L and U may have a zero because of cancellation where P*X is nonzero.  L*U
+    differs from P*X only outside of the sparsity pattern of P*X.
+ 
+    [L,U] = LUINC(X,'0') produces upper triangular U and L is a permutation of
+    unit lower triangular matrix.  Thus no comparison can be made between the
+    sparsity patterns of L,U and X, although nnz(L) + nnz(U) = nnz(X) + n.  L*U
+    differs from X only outside of its sparsity pattern.
+ 
+    LU = LUINC(X,'0') returns "L+U-I", where L is unit lower triangular, U is
+    upper triangular and the permutation information is lost.
+ 
+    Example:
+ 
+       load west0479
+       A = west0479;
+       [L,U,P] = luinc(A,'0');
+       isequal(spones(U),spones(triu(P*A)))
+       spones(L) ~= spones(tril(P*A))
+       D = (L*U) .* spones(P*A) - P*A
+ 
+       spones(L) differs from spones(tril(P*A)) at some positions on the
+       diagonal and at one position in L where cancellation zeroed out a
+       nonzero element of P*A.  The entries of D are of the order of eps.
+ 
+    LUINC works only for sparse matrices.
+ 
+    See also LU, CHOLINC, BICG.
+
+*/
--- a/src/DLD-FUNCTIONS/sparse.cc	Thu Apr 14 19:35:20 2005 +0000
+++ b/src/DLD-FUNCTIONS/sparse.cc	Thu Apr 14 22:17:27 2005 +0000
@@ -139,7 +139,8 @@
    if (nargin == 2 && ! args(0).is_scalar_type() && args(1).is_scalar_type())
        mutate = (args(1).double_value() != 0.);
 
-   if (nargin == 1 || (nargin == 2 && mutate)) 
+   if (nargin == 1 || (nargin == 2 && ! args(0).is_scalar_type() && 
+		       args(1).is_scalar_type()))
      {
        octave_value arg = args (0);
 
--- a/src/DLD-FUNCTIONS/splu.cc	Thu Apr 14 19:35:20 2005 +0000
+++ b/src/DLD-FUNCTIONS/splu.cc	Thu Apr 14 22:17:27 2005 +0000
@@ -82,8 +82,8 @@
 
   octave_value arg = args(0);
 
-  int nr = arg.rows ();
-  int nc = arg.columns ();
+  octave_idx_type nr = arg.rows ();
+  octave_idx_type nc = arg.columns ();
 
   int arg_is_empty = empty_arg ("splu", nr, nc);
 
@@ -114,21 +114,21 @@
 	    thres = tmp (0);
 	  else if (dv(0) == 1 || dv(1) == 1)
 	    {
-	      int nel = tmp.numel ();
+	      octave_idx_type nel = tmp.numel ();
 	      Qinit.resize (nel);
-	      for (int i = 0; i < nel; i++)
+	      for (octave_idx_type i = 0; i < nel; i++)
 		Qinit (i) = tmp (i) - 1;
 	      have_Qinit = true;
 	    }
 	  else
 	    {
-	      int t_nc = tmp.cols ();
+	      octave_idx_type t_nc = tmp.cols ();
 	      
 	      if (tmp.nnz () != t_nc)
 		error ("splu: Not a valid permutation matrix");
 	      else
 		{
-		  for (int i = 0; i < t_nc + 1; i++)
+		  for (octave_idx_type i = 0; i < t_nc + 1; i++)
 		    if (tmp.cidx(i) != i)
 		      {
 			error ("splu: Not a valid permutation matrix");
@@ -138,7 +138,7 @@
 		  
 	      if (!error_state)
 		{
-		  for (int i = 0; i < t_nc; i++)
+		  for (octave_idx_type i = 0; i < t_nc; i++)
 		    if (tmp.data (i) != 1.)
 		      {
 			error ("splu: Not a valid permutation matrix");
@@ -168,9 +168,9 @@
 	    thres = tmp (0);
 	  else if (dv(0) == 1 || dv(1) == 1)
 	    {
-	      int nel = tmp.numel ();
+	      octave_idx_type nel = tmp.numel ();
 	      Qinit.resize (nel);
-	      for (int i = 0; i < nel; i++)
+	      for (octave_idx_type i = 0; i < nel; i++)
 		Qinit (i) = tmp (i) - 1;
 	      have_Qinit = true;
 	    }
@@ -178,13 +178,13 @@
 	    {
 	      SparseMatrix tmp2 (tmp);
 
-	      int t_nc = tmp2.cols ();
+	      octave_idx_type t_nc = tmp2.cols ();
 	      
 	      if (tmp2.nnz () != t_nc)
 		error ("splu: Not a valid permutation matrix");
 	      else
 		{
-		  for (int i = 0; i < t_nc + 1; i++)
+		  for (octave_idx_type i = 0; i < t_nc + 1; i++)
 		    if (tmp2.cidx(i) != i)
 		      {
 			error ("splu: Not a valid permutation matrix");
@@ -194,7 +194,7 @@
 		  
 	      if (!error_state)
 		{
-		  for (int i = 0; i < t_nc; i++)
+		  for (octave_idx_type i = 0; i < t_nc; i++)
 		    if (tmp2.data (i) != 1.)
 		      {
 			error ("splu: Not a valid permutation matrix");
@@ -219,9 +219,9 @@
 
       if (nargout < 4 && ! have_Qinit)
 	{
-	  int m_nc = m.cols ();
+	  octave_idx_type m_nc = m.cols ();
 	  Qinit.resize (m_nc);
-	  for (int i = 0; i < m_nc; i++)
+	  for (octave_idx_type i = 0; i < m_nc; i++)
 	    Qinit (i) = i;
 	}
 
@@ -291,9 +291,9 @@
 
       if (nargout < 4 && ! have_Qinit)
 	{
-	  int m_nc = m.cols ();
+	  octave_idx_type m_nc = m.cols ();
 	  Qinit.resize (m_nc);
-	  for (int i = 0; i < m_nc; i++)
+	  for (octave_idx_type i = 0; i < m_nc; i++)
 	    Qinit (i) = i;
 	}
 
--- a/src/Makefile.in	Thu Apr 14 19:35:20 2005 +0000
+++ b/src/Makefile.in	Thu Apr 14 22:17:27 2005 +0000
@@ -50,7 +50,7 @@
 	lpsolve.cc lsode.cc lu.cc minmax.cc pinv.cc qr.cc \
 	quad.cc qz.cc rand.cc schur.cc sort.cc sparse.cc spdet.cc \
 	splu.cc spparms.cc sqrtm.cc svd.cc syl.cc time.cc gplot.l \
-	__glpk__.cc
+	__glpk__.cc luinc.cc
 
 DLD_SRC := $(addprefix DLD-FUNCTIONS/, $(DLD_XSRC))
 
--- a/src/ov-bool-sparse.h	Thu Apr 14 19:35:20 2005 +0000
+++ b/src/ov-bool-sparse.h	Thu Apr 14 22:17:27 2005 +0000
@@ -82,9 +82,9 @@
 
   octave_value *try_narrowing_conversion (void);
 
-#if 0
-  idx_vector index_vector (void) const { return idx_vector (matrix); }
-#endif
+  // XXX FIXME XXX Adapt idx_vector to allow sparse logical indexing!!
+  idx_vector index_vector (void) const 
+    { return idx_vector (bool_array_value ()); }
 
   bool is_bool_matrix (void) const { return true; }
 
--- a/src/ov-mapper.cc	Thu Apr 14 19:35:20 2005 +0000
+++ b/src/ov-mapper.cc	Thu Apr 14 22:17:27 2005 +0000
@@ -57,6 +57,25 @@
 }
 
 static bool
+any_element_less_than (const SparseMatrix& a, double val)
+{
+  octave_idx_type len = a.nonzero ();
+
+  if (val > 0. && len != a.numel ())
+    return true;
+
+  for (octave_idx_type i = 0; i < len; i++)
+    {
+      OCTAVE_QUIT;
+
+      if (a.data(i) < val)
+	return true;
+    }
+
+  return false;
+}
+
+static bool
 any_element_greater_than (const NDArray& a, double val)
 {
   octave_idx_type len = a.length ();
@@ -72,6 +91,25 @@
   return false;
 }
 
+static bool
+any_element_greater_than (const SparseMatrix& a, double val)
+{
+  octave_idx_type len = a.nonzero ();
+
+  if (val < 0. && len != a.numel ())
+    return true;
+
+  for (octave_idx_type i = 0; i < len; i++)
+    {
+      OCTAVE_QUIT;
+
+      if (a.data(i) > val)
+	return true;
+    }
+
+  return false;
+}
+
 // In most cases, we could use the map member function from the NDArray
 // classes, but as currently implemented, they don't allow us to
 // detect errors and abort properly.  So use these macros to do the
@@ -104,6 +142,72 @@
 #define MAPPER_LOOP(T, F, M) \
   MAPPER_LOOP_1 (T, F, M, )
 
+#define SPARSE_MAPPER_LOOP_2(T, ET, F, M, CONV, R) \
+  do \
+    { \
+      ET f_zero = CONV (F (0.)); \
+      \
+      if (f_zero != 0.) \
+	{ \
+	  octave_idx_type nr = M.rows (); \
+	  octave_idx_type nc = M.cols (); \
+	  \
+	  T result (nr, nc, f_zero); \
+	  \
+	  for (octave_idx_type j = 0; j < nc; j++) \
+	    for (octave_idx_type i = M.cidx(j); i < M.cidx (j+1); i++) \
+	      { \
+		OCTAVE_QUIT; \
+	        result.elem (M.ridx (i), j) = CONV (F (M.data(i))); \
+		\
+		if (error_state) \
+		  return retval; \
+	      } \
+	  \
+	  result.maybe_compress (true);	\
+          retval = R; \
+	} \
+      else \
+	{ \
+	  octave_idx_type nnz = M.nonzero (); \
+	  octave_idx_type nr = M.rows (); \
+	  octave_idx_type nc = M.cols (); \
+	  \
+	  T result (nr, nc, nnz); \
+	  ET zero = ET (0.); \
+	  octave_idx_type ii = 0; \
+	  result.cidx (ii) = 0; \
+	  \
+	  for (octave_idx_type j = 0; j < nc; j++) \
+	    { \
+	      for (octave_idx_type i = M.cidx(j); i < M.cidx (j+1); i++) \
+		{ \
+		  ET val = CONV (F (M.data (i))); \
+		  if (val != zero) \
+		    { \
+		      result.data (ii) = val; \
+		      result.ridx (ii++) = M.ridx (i); \
+		    } \
+		  OCTAVE_QUIT; \
+		  \
+		  if (error_state) \
+		    return retval; \
+		} \
+	      result.cidx (j+1) = ii; \
+	    } \
+	  \
+	  result.maybe_compress (false); \
+          retval = R; \
+	} \
+    } \
+  while (0)
+
+#define SPARSE_MAPPER_LOOP_1(T, ET, F, M, CONV)	\
+  SPARSE_MAPPER_LOOP_2 (T, ET, F, M, CONV, result)
+
+#define SPARSE_MAPPER_LOOP(T, ET, F, M) \
+  SPARSE_MAPPER_LOOP_1 (T, ET, F, M, )
+
 octave_value
 octave_mapper::apply (const octave_value& arg) const
 {
@@ -136,6 +240,32 @@
 	    error ("%s: unable to handle real arguments",
 		   name().c_str ());
 	}
+      else if (arg.class_name () == "sparse")
+	{
+	  const SparseMatrix m = arg.sparse_matrix_value ();
+
+	  if (error_state)
+	    return retval;
+
+	  if (can_ret_cmplx_for_real
+	      && (any_element_less_than (m, lower_limit)
+		  || any_element_greater_than (m, upper_limit)))
+	    {
+	      if (c_c_map_fcn)
+		SPARSE_MAPPER_LOOP (SparseComplexMatrix, Complex, 
+				    c_c_map_fcn, m);
+	      else
+		error ("%s: unable to handle real arguments",
+		       name().c_str ());
+	    }
+	  else if (d_d_map_fcn)
+	    SPARSE_MAPPER_LOOP (SparseMatrix, double, d_d_map_fcn, m);
+	  else if (d_b_map_fcn)
+	    SPARSE_MAPPER_LOOP (SparseBoolMatrix, bool, d_b_map_fcn, m);
+	  else
+	    error ("%s: unable to handle real arguments",
+		   name().c_str ());
+	}
       else
 	{
 	  NDArray m = arg.array_value ();
@@ -178,6 +308,25 @@
 	    error ("%s: unable to handle complex arguments",
 		   name().c_str ());
 	}
+      else if (arg.class_name () == "sparse")
+	{
+	  SparseComplexMatrix cm = arg.sparse_complex_matrix_value ();
+
+	  if (error_state)
+	    return retval;
+
+	  if (d_c_map_fcn)
+	    SPARSE_MAPPER_LOOP (SparseMatrix, double, d_c_map_fcn, cm);
+	  else if (c_c_map_fcn)
+	    SPARSE_MAPPER_LOOP (SparseComplexMatrix, Complex, 
+				c_c_map_fcn, cm);
+	  else if (c_b_map_fcn)
+	    SPARSE_MAPPER_LOOP (SparseBoolMatrix, bool, 
+				c_b_map_fcn, cm);
+	  else
+	    error ("%s: unable to handle complex arguments",
+		   name().c_str ());
+	}
       else
 	{
 	  ComplexNDArray cm = arg.complex_array_value ();
--- a/src/ov-re-sparse.cc	Thu Apr 14 19:35:20 2005 +0000
+++ b/src/ov-re-sparse.cc	Thu Apr 14 22:17:27 2005 +0000
@@ -47,6 +47,19 @@
 
 DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA (octave_sparse_matrix, "sparse matrix", "sparse");
 
+idx_vector
+octave_sparse_matrix::index_vector (void) const
+{
+  if (matrix.numel () == matrix.nonzero ())
+    return idx_vector (array_value ());
+  else
+    {
+      std::string nm = type_name ();
+      error ("%s type invalid as index value", nm.c_str ());
+      return idx_vector ();
+    }
+}
+
 octave_value *
 octave_sparse_matrix::try_narrowing_conversion (void)
 {
@@ -146,11 +159,11 @@
 octave_sparse_matrix::streamoff_array_value (void) const
 {
   streamoff_array retval (dims ());
-  int nc = matrix.cols ();
-  int nr = matrix.rows ();
+  octave_idx_type nc = matrix.cols ();
+  octave_idx_type nr = matrix.rows ();
 
-  for (int j = 0; j < nc; j++)
-    for (int i = matrix.cidx(i); i < matrix.cidx(i+1); i++)
+  for (octave_idx_type j = 0; j < nc; j++)
+    for (octave_idx_type i = matrix.cidx(j); i < matrix.cidx(j+1); i++)
       {
 	double d = matrix.data(i);
 
@@ -169,6 +182,67 @@
   return retval;
 }
 
+octave_value
+octave_sparse_matrix::convert_to_str_internal (bool, bool) const
+{
+  octave_value retval;
+  dim_vector dv = dims ();
+  octave_idx_type nel = dv.numel ();
+
+  if (nel == 0)
+    {
+      char s = '\0';
+      retval = octave_value (&s);
+    }
+  else
+    {
+      octave_idx_type nr = matrix.rows ();
+      octave_idx_type nc = matrix.cols ();
+      charNDArray chm (dv, static_cast<char> (0));
+	  
+      bool warned = false;
+
+      for (octave_idx_type j = 0; j < nc; j++)
+	for (octave_idx_type i = matrix.cidx(j); 
+	     i < matrix.cidx(j+1); i++)
+	  {
+	    OCTAVE_QUIT;
+
+	    double d = matrix.data (i);
+
+	      if (xisnan (d))
+		{
+		  ::error ("invalid conversion from NaN to character");
+		  return retval;
+		}
+	      else
+		{
+		  int ival = NINT (d);
+
+		  if (ival < 0 || ival > UCHAR_MAX)
+		    {
+		      // XXX FIXME XXX -- is there something
+		      // better we could do?
+
+		      ival = 0;
+
+		      if (! warned)
+			{
+			  ::warning ("range error for conversion to character value");
+			  warned = true;
+			}
+		    }
+
+		  chm (matrix.ridx(i) + j * nr) = 
+		    static_cast<char> (ival);
+		}
+	  }
+      retval = octave_value (chm, 1);
+    }
+
+  return retval;
+}
+
 bool 
 octave_sparse_matrix::save_binary (std::ostream& os, bool&save_as_floats)
 {
--- a/src/ov-re-sparse.h	Thu Apr 14 19:35:20 2005 +0000
+++ b/src/ov-re-sparse.h	Thu Apr 14 22:17:27 2005 +0000
@@ -81,9 +81,7 @@
 
   octave_value *try_narrowing_conversion (void);
 
-#if 0
-  idx_vector index_vector (void) const { return idx_vector (matrix); }
-#endif
+  idx_vector index_vector (void) const;
 
   bool is_real_matrix (void) const { return true; }
 
@@ -114,6 +112,8 @@
 
   streamoff_array streamoff_array_value (void) const;
 
+  octave_value convert_to_str_internal (bool pad, bool force) const;
+
 #if 0
   int write (octave_stream& os, int block_size,
 	     oct_data_conv::data_type output_type, int skip,