--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libinterp/corefcn/__dsearchn__.cc	Fri Feb 27 19:44:28 2015 -0500
@@ -0,0 +1,115 @@
+/*
+
+Copyright (C) 2007-2015 David Bateman
+
+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 <iostream>
+#include <fstream>
+#include <string>
+
+#include "lo-math.h"
+
+#include "defun.h"
+#include "error.h"
+#include "oct-obj.h"
+
+DEFUN (__dsearchn__, args, ,
+       "-*- texinfo -*-\n\
+@deftypefn {Built-in Function} {[@var{idx}, @var{d}] =} dsearch (@var{x}, @var{xi})\n\
+Undocumented internal function.\n\
+@end deftypefn")
+{
+  int nargin = args.length ();
+  octave_value_list retval;
+
+  if (nargin != 2)
+    {
+      print_usage ();
+      return retval;
+    }
+
+  Matrix x = args(0).matrix_value ().transpose ();
+  Matrix xi = args(1).matrix_value ().transpose ();
+
+  if (! error_state)
+    {
+      if (x.rows () != xi.rows () || x.columns () < 1)
+        error ("__dsearch__: number of rows of X and XI must match");
+      else
+        {
+          octave_idx_type n = x.rows ();
+          octave_idx_type nx = x.columns ();
+          octave_idx_type nxi = xi.columns ();
+
+          ColumnVector idx (nxi);
+          double *pidx = idx.fortran_vec ();
+          ColumnVector dist (nxi);
+          double *pdist = dist.fortran_vec ();
+
+#define DIST(dd, y, yi, m) \
+  dd = 0.; \
+  for (octave_idx_type k = 0; k < m; k++) \
+   { \
+     double yd = y[k] - yi[k]; \
+     dd += yd * yd; \
+   } \
+  dd = sqrt (dd);
+
+          const double *pxi = xi.fortran_vec ();
+          for (octave_idx_type i = 0; i < nxi; i++)
+            {
+              double d0;
+              const double *px = x.fortran_vec ();
+              DIST(d0, px, pxi, n);
+              *pidx = 1.;
+              for (octave_idx_type j = 1; j < nx; j++)
+                {
+                  px += n;
+                  double d;
+                  DIST (d, px, pxi, n);
+                  if (d < d0)
+                    {
+                      d0 = d;
+                      *pidx = static_cast<double>(j + 1);
+                    }
+                  OCTAVE_QUIT;
+                }
+
+              *pdist++ = d0;
+              pidx++;
+              pxi += n;
+            }
+
+          retval(1) = dist;
+          retval(0) = idx;
+        }
+    }
+
+  return retval;
+}
+
+/*
+## No test needed for internal helper function.
+%!assert (1)
+*/

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libinterp/corefcn/__ichol__.cc	Fri Feb 27 19:44:28 2015 -0500
@@ -0,0 +1,519 @@
+/*
+
+Copyright (C) 2014-2015 Eduardo Ramos Fernández <eduradical951@gmail.com>
+Copyright (C) 2013-2015 Kai T. Ohlhus <k.ohlhus@gmail.com>
+
+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 "oct-locbuf.h"
+
+#include "defun.h"
+#include "error.h"
+#include "parse.h"
+
+// Secondary functions for complex and real case used in ichol algorithms.
+Complex ichol_mult_complex (Complex a, Complex b)
+{
+#if defined (HAVE_CXX_COMPLEX_SETTERS)
+  b.imag (-std::imag (b));
+#elif defined (HAVE_CXX_COMPLEX_REFERENCE_ACCESSORS)
+  b.imag () = -std::imag (b);
+#else
+  b = std::conj (b);
+#endif
+  return a * b;
+}
+
+double ichol_mult_real (double a, double b)
+{
+  return a * b;
+}
+
+bool ichol_checkpivot_complex (Complex pivot)
+{
+  if (pivot.imag () != 0)
+    {
+      error ("ichol: non-real pivot encountered.  The matrix must be hermitian.");
+      return false;
+    }
+  else if (pivot.real () < 0)
+    {
+      error ("ichol: negative pivot encountered");
+      return false;
+    }
+  return true;
+}
+
+bool ichol_checkpivot_real (double pivot)
+{
+  if (pivot < 0)
+    {
+      error ("ichol: negative pivot encountered");
+      return false;
+    }
+  return true;
+}
+
+template <typename octave_matrix_t, typename T, T (*ichol_mult) (T, T),
+          bool (*ichol_checkpivot) (T)>
+void ichol_0 (octave_matrix_t& sm, const std::string michol = "off")
+{
+
+  const octave_idx_type n = sm.cols ();
+  octave_idx_type j1, jend, j2, jrow, jjrow, j, jw, i, k, jj, r;
+  T tl;
+  char opt;
+  enum {OFF, ON};
+  if (michol == "on")
+    opt = ON;
+  else
+    opt = OFF;
+
+  // Input matrix pointers
+  octave_idx_type* cidx = sm.cidx ();
+  octave_idx_type* ridx = sm.ridx ();
+  T* data = sm.data ();
+
+  // Working arrays
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, Lfirst, n);
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, Llist, n);
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, iw, n);
+  OCTAVE_LOCAL_BUFFER (T, dropsums, n);
+
+  // Initialize working arrays
+  for (i = 0; i < n; i++)
+    {
+      iw[i] = -1;
+      Llist[i] = -1;
+      Lfirst[i] = -1;
+      dropsums[i] = 0;
+    }
+
+  // Main loop
+  for (k = 0; k < n; k++)
+    {
+      j1 = cidx[k];
+      j2 = cidx[k+1];
+      for (j = j1; j < j2; j++)
+        iw[ridx[j]] = j;
+
+      jrow = Llist [k];
+      // Iterate over each non-zero element in the actual row.
+      while (jrow != -1)
+        {
+          jjrow = Lfirst[jrow];
+          jend = cidx[jrow+1];
+          for (jj = jjrow; jj < jend; jj++)
+            {
+              r = ridx[jj];
+              jw = iw[r];
+              tl = ichol_mult (data[jj], data[jjrow]);
+              if (jw != -1)
+                data[jw] -= tl;
+              else
+                // Because of the symmetry of the matrix, we know
+                // the drops in the column r are also in the column k.
+                if (opt == ON)
+                  {
+                    dropsums[r] -= tl;
+                    dropsums[k] -= tl;
+                  }
+            }
+          // Update the linked list and the first entry of the actual column.
+          if ((jjrow + 1) < jend)
+            {
+              Lfirst[jrow]++;
+              j = jrow;
+              jrow = Llist[jrow];
+              Llist[j] = Llist[ridx[Lfirst[j]]];
+              Llist[ridx[Lfirst[j]]] = j;
+            }
+          else
+            jrow = Llist[jrow];
+        }
+
+      if (opt == ON)
+        data[j1] += dropsums[k];
+
+      if (ridx[j1] != k)
+        {
+          error ("ichol: encountered a pivot equal to 0");
+          break;
+        }
+
+      if (! ichol_checkpivot (data[j1]))
+        break;
+
+      data[cidx[k]] = std::sqrt (data[j1]);
+
+      // Update Llist and Lfirst with the k-column information.  Also,
+      // scale the column elements by the pivot and reset the working array iw.
+      if (k < (n - 1))
+        {
+          iw[ridx[j1]] = -1;
+          for (i = j1 + 1; i < j2; i++)
+            {
+              iw[ridx[i]] = -1;
+              data[i] /= data[j1];
+            }
+          Lfirst[k] = j1;
+          if ((Lfirst[k] + 1) < j2)
+            {
+              Lfirst[k]++;
+              jjrow = ridx[Lfirst[k]];
+              Llist[k] = Llist[jjrow];
+              Llist[jjrow] = k;
+            }
+        }
+    }
+}
+
+DEFUN (__ichol0__, args, nargout,
+       "-*- texinfo -*-\n\
+@deftypefn  {Built-in Function} {@var{L} =} __ichol0__ (@var{A})\n\
+@deftypefnx {Built-in Function} {@var{L} =} __ichol0__ (@var{A}, @var{michol})\n\
+Undocumented internal function.\n\
+@end deftypefn")
+
+{
+  octave_value_list retval;
+
+  int nargin = args.length ();
+  std::string michol = "off";
+
+  if (nargout > 1 || nargin < 1 || nargin > 2)
+    {
+      print_usage ();
+      return retval;
+    }
+
+  if (nargin == 2)
+    michol = args(1).string_value ();
+
+  // In ICHOL0 algorithm the zero-pattern of the input matrix is preserved
+  // so it's structure does not change during the algorithm.  The same input
+  // matrix is used to build the output matrix due to that fact.
+  octave_value_list param_list;
+  if (!args(0).is_complex_type ())
+    {
+      SparseMatrix sm = args(0).sparse_matrix_value ();
+      param_list.append (sm);
+      sm = feval ("tril", param_list)(0).sparse_matrix_value ();
+      ichol_0 <SparseMatrix, double, ichol_mult_real,
+               ichol_checkpivot_real> (sm, michol);
+      if (! error_state)
+        retval(0) = sm;
+    }
+  else
+    {
+      SparseComplexMatrix sm = args(0).sparse_complex_matrix_value ();
+      param_list.append (sm);
+      sm = feval ("tril", param_list)(0).sparse_complex_matrix_value ();
+      ichol_0 <SparseComplexMatrix, Complex, ichol_mult_complex,
+               ichol_checkpivot_complex> (sm, michol);
+      if (! error_state)
+        retval(0) = sm;
+    }
+
+  return retval;
+}
+
+template <typename octave_matrix_t, typename T,  T (*ichol_mult) (T, T),
+          bool (*ichol_checkpivot) (T)>
+void ichol_t (const octave_matrix_t& sm, octave_matrix_t& L, const T* cols_norm,
+              const T droptol, const std::string michol = "off")
+
+{
+
+  const octave_idx_type n = sm.cols ();
+  octave_idx_type j, jrow, jend, jjrow, i, k, jj, total_len,
+                  w_len, max_len, ind;
+  char opt;
+  enum {OFF, ON};
+  if (michol == "on")
+    opt = ON;
+  else
+    opt = OFF;
+
+  // Input matrix pointers
+  octave_idx_type* cidx = sm.cidx ();
+  octave_idx_type* ridx = sm.ridx ();
+  T* data = sm.data ();
+
+  // Output matrix data structures.  Because the final zero pattern pattern of
+  // the output matrix is not known due to fill-in elements, a heuristic
+  // approach has been adopted for memory allocation.  The size of ridx_out_l
+  // and data_out_l is incremented 10% of their actual size (nnz (A) in the
+  // beginning).  If that amount is less than n, their size is just incremented
+  // in n elements.  This way the number of reallocations decreases throughout
+  // the process, obtaining a good performance.
+  max_len = sm.nnz ();
+  max_len += (0.1 * max_len) > n ? 0.1 * max_len : n;
+  Array <octave_idx_type> cidx_out_l (dim_vector (n + 1, 1));
+  octave_idx_type* cidx_l = cidx_out_l.fortran_vec ();
+  Array <octave_idx_type> ridx_out_l (dim_vector (max_len ,1));
+  octave_idx_type* ridx_l = ridx_out_l.fortran_vec ();
+  Array <T> data_out_l (dim_vector (max_len, 1));
+  T* data_l = data_out_l.fortran_vec ();
+
+  // Working arrays
+  OCTAVE_LOCAL_BUFFER (T, w_data, n);
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, Lfirst, n);
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, Llist, n);
+  OCTAVE_LOCAL_BUFFER (T, col_drops, n);
+  std::vector <octave_idx_type> vec;
+  vec.resize (n);
+
+  T zero = T (0);
+  cidx_l[0] = cidx[0];
+  for (i = 0; i < n; i++)
+    {
+      Llist[i] = -1;
+      Lfirst[i] = -1;
+      w_data[i] = 0;
+      col_drops[i] = zero;
+      vec[i] = 0;
+    }
+
+  total_len = 0;
+  for (k = 0; k < n; k++)
+    {
+      ind = 0;
+      for (j = cidx[k]; j < cidx[k+1]; j++)
+        {
+          w_data[ridx[j]] = data[j];
+          if (ridx[j] != k)
+            {
+              vec[ind] = ridx[j];
+              ind++;
+            }
+        }
+      jrow = Llist[k];
+      while (jrow != -1)
+        {
+          jjrow = Lfirst[jrow];
+          jend = cidx_l[jrow+1];
+          for (jj = jjrow; jj < jend; jj++)
+            {
+              j = ridx_l[jj];
+              // If the element in the j position of the row is zero,
+              // then it will become non-zero, so we add it to the
+              // vector that tracks non-zero elements in the working row.
+              if (w_data[j] == zero)
+                {
+                  vec[ind] = j;
+                  ind++;
+                }
+              w_data[j] -=  ichol_mult (data_l[jj], data_l[jjrow]);
+            }
+          // Update the actual column first element and
+          // update the linked list of the jrow row.
+          if ((jjrow + 1) < jend)
+            {
+              Lfirst[jrow]++;
+              j = jrow;
+              jrow = Llist[jrow];
+              Llist[j] = Llist[ridx_l[Lfirst[j]]];
+              Llist[ridx_l[Lfirst[j]]] = j;
+            }
+          else
+            jrow = Llist[jrow];
+        }
+
+      // Resizing output arrays
+      if ((max_len - total_len) < n)
+        {
+          max_len += (0.1 * max_len) > n ? 0.1 * max_len : n;
+          data_out_l.resize (dim_vector (max_len, 1));
+          data_l = data_out_l.fortran_vec ();
+          ridx_out_l.resize (dim_vector (max_len, 1));
+          ridx_l = ridx_out_l.fortran_vec ();
+        }
+
+      // The sorting of the non-zero elements of the working column can be
+      // handled in a couple of ways.  The most efficient two I found, are
+      // keeping the elements in an ordered binary search tree dynamically or
+      // keep them unsorted in a vector and at the end of the outer iteration
+      // order them.  The last approach exhibits lower execution times.
+      std::sort (vec.begin (), vec.begin () + ind);
+
+      data_l[total_len] = w_data[k];
+      ridx_l[total_len] = k;
+      w_len = 1;
+
+      // Extract the non-zero elements of working column and
+      // drop the elements that are lower than droptol * cols_norm[k].
+      for (i = 0; i < ind ; i++)
+        {
+          jrow = vec[i];
+          if (w_data[jrow] != zero)
+            {
+              if (std::abs (w_data[jrow]) < (droptol * cols_norm[k]))
+                {
+                  if (opt == ON)
+                    {
+                      col_drops[k] += w_data[jrow];
+                      col_drops[jrow] += w_data[jrow];
+                    }
+                }
+              else
+                {
+                  data_l[total_len + w_len] = w_data[jrow];
+                  ridx_l[total_len + w_len] = jrow;
+                  w_len++;
+                }
+              vec[i] = 0;
+            }
+          w_data[jrow] = zero;
+        }
+
+      // Compensate column sums --> michol option
+      if (opt == ON)
+        data_l[total_len] += col_drops[k];
+
+      if (data_l[total_len] == zero)
+        {
+          error ("ichol: encountered a pivot equal to 0");
+          break;
+        }
+      else if (! ichol_checkpivot (data_l[total_len]))
+        break;
+
+      // Once elements are dropped and compensation of column sums are done,
+      // scale the elements by the pivot.
+      data_l[total_len] = std::sqrt (data_l[total_len]);
+      for (jj = total_len + 1; jj < (total_len + w_len); jj++)
+        data_l[jj] /=  data_l[total_len];
+      total_len += w_len;
+      // Check if there are too many elements to be indexed with
+      // octave_idx_type type due to fill-in during the process.
+      if (total_len < 0)
+        {
+          error ("ichol: integer overflow.  Too many fill-in elements in L");
+          break;
+        }
+      cidx_l[k+1] = cidx_l[k] - cidx_l[0] + w_len;
+
+      // Update Llist and Lfirst with the k-column information.
+      if (k < (n - 1))
+        {
+          Lfirst[k] = cidx_l[k];
+          if ((Lfirst[k] + 1) < cidx_l[k+1])
+            {
+              Lfirst[k]++;
+              jjrow = ridx_l[Lfirst[k]];
+              Llist[k] = Llist[jjrow];
+              Llist[jjrow] = k;
+            }
+        }
+    }
+
+  if (! error_state)
+    {
+      // Build the output matrices
+      L = octave_matrix_t (n, n, total_len);
+      for (i = 0; i <= n; i++)
+        L.cidx (i) = cidx_l[i];
+      for (i = 0; i < total_len; i++)
+        {
+          L.ridx (i) = ridx_l[i];
+          L.data (i) = data_l[i];
+        }
+    }
+}
+
+DEFUN (__icholt__, args, nargout,
+       "-*- texinfo -*-\n\
+@deftypefn  {Built-in Function} {@var{L} =} __icholt__ (@var{A})\n\
+@deftypefnx {Built-in Function} {@var{L} =} __icholt__ (@var{A}, @var{droptol})\n\
+@deftypefnx {Built-in Function} {@var{L} =} __icholt__ (@var{A}, @var{droptol}, @var{michol})\n\
+Undocumented internal function.\n\
+@end deftypefn")
+{
+  octave_value_list retval;
+  int nargin = args.length ();
+  // Default values of parameters
+  std::string michol = "off";
+  double droptol = 0;
+
+  if (nargout > 1 || nargin < 1 || nargin > 3)
+    {
+      print_usage ();
+      return retval;
+    }
+
+  // Don't repeat input validation of arguments done in ichol.m
+
+  if (nargin >= 2)
+    droptol = args(1).double_value ();
+
+  if (nargin == 3)
+    michol = args(2).string_value ();
+
+  octave_value_list param_list;
+  if (! args(0).is_complex_type ())
+    {
+      Array <double> cols_norm;
+      SparseMatrix L;
+      param_list.append (args(0).sparse_matrix_value ());
+      SparseMatrix sm_l =
+        feval ("tril", param_list)(0).sparse_matrix_value ();
+      param_list(0) = sm_l;
+      param_list(1) = 1;
+      param_list(2) = "cols";
+      cols_norm = feval ("norm", param_list)(0).vector_value ();
+      param_list.clear ();
+      ichol_t <SparseMatrix,
+               double, ichol_mult_real, ichol_checkpivot_real>
+               (sm_l, L, cols_norm.fortran_vec (), droptol, michol);
+      if (! error_state)
+        retval(0) = L;
+    }
+  else
+    {
+      Array <Complex> cols_norm;
+      SparseComplexMatrix L;
+      param_list.append (args(0).sparse_complex_matrix_value ());
+      SparseComplexMatrix sm_l =
+        feval ("tril", param_list)(0).sparse_complex_matrix_value ();
+      param_list(0) = sm_l;
+      param_list(1) = 1;
+      param_list(2) = "cols";
+      cols_norm = feval ("norm", param_list)(0).complex_vector_value ();
+      param_list.clear ();
+      ichol_t <SparseComplexMatrix,
+               Complex, ichol_mult_complex, ichol_checkpivot_complex>
+               (sm_l, L, cols_norm.fortran_vec (),
+                Complex (droptol), michol);
+      if (! error_state)
+        retval(0) = L;
+    }
+
+  return retval;
+}
+
+/*
+## No test needed for internal helper function.
+%!assert (1)
+*/
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libinterp/corefcn/__ilu__.cc	Fri Feb 27 19:44:28 2015 -0500
@@ -0,0 +1,1108 @@
+/*
+
+Copyright (C) 2014-2015 Eduardo Ramos Fernández <eduradical951@gmail.com>
+Copyright (C) 2013-2015 Kai T. Ohlhus <k.ohlhus@gmail.com>
+
+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 "oct-locbuf.h"
+
+#include "defun.h"
+#include "error.h"
+#include "parse.h"
+
+// That function implements the IKJ and JKI variants of Gaussian elimination to
+// perform the ILUTP decomposition.  The behaviour is controlled by milu
+// parameter.  If milu = ['off'|'col'] the JKI version is performed taking
+// advantage of CCS format of the input matrix.  If milu = 'row' the input
+// matrix has to be transposed to obtain the equivalent CRS structure so we can
+// work efficiently with rows.  In this case IKJ version is used.
+template <typename octave_matrix_t, typename T>
+void ilu_0 (octave_matrix_t& sm, const std::string milu = "off")
+{
+
+  const octave_idx_type n = sm.cols ();
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, iw, n);
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, uptr, n);
+  octave_idx_type j1, j2, jrow, jw, i, k, jj;
+  T tl, r;
+
+  enum {OFF, ROW, COL};
+  char opt;
+  if (milu == "row")
+    {
+      opt = ROW;
+      sm = sm.transpose ();
+    }
+  else if (milu == "col")
+    opt = COL;
+  else
+    opt = OFF;
+
+  octave_idx_type* cidx = sm.cidx ();
+  octave_idx_type* ridx = sm.ridx ();
+  T* data = sm.data ();
+  for (i = 0; i < n; i++)
+    iw[i] = -1;
+  for (k = 0; k < n; k++)
+    {
+      j1 = cidx[k];
+      j2 = cidx[k+1] - 1;
+      octave_idx_type j;
+      for (j = j1; j <= j2; j++)
+        {
+          iw[ridx[j]] = j;
+        }
+      r = 0;
+      j = j1;
+      jrow = ridx[j];
+      while ((jrow < k) && (j <= j2))
+        {
+          if (opt == ROW)
+            {
+              tl = data[j] / data[uptr[jrow]];
+              data[j] = tl;
+            }
+          for (jj = uptr[jrow] + 1; jj < cidx[jrow+1]; jj++)
+            {
+              jw = iw[ridx[jj]];
+              if (jw != -1)
+                if (opt == ROW)
+                  data[jw] -= tl * data[jj];
+                else
+                  data[jw] -= data[j] * data[jj];
+
+              else
+                // That is for the milu='row'
+                if (opt == ROW)
+                  r += tl * data[jj];
+                else if (opt == COL)
+                  r += data[j] * data[jj];
+            }
+          j++;
+          jrow = ridx[j];
+        }
+      uptr[k] = j;
+      if (opt != OFF)
+        data[uptr[k]] -= r;
+      if (opt != ROW)
+        for (jj = uptr[k] + 1; jj < cidx[k+1]; jj++)
+          data[jj] /=  data[uptr[k]];
+      if (k != jrow)
+        {
+          error ("ilu: A has a zero on the diagonal");
+          break;
+        }
+
+      if (data[j] == T(0))
+        {
+          error ("ilu: encountered a pivot equal to 0");
+          break;
+        }
+      for (i = j1; i <= j2; i++)
+        iw[ridx[i]] = -1;
+    }
+  if (opt == ROW)
+    sm = sm.transpose ();
+}
+
+DEFUN (__ilu0__, args, nargout,
+       "-*- texinfo -*-\n\
+@deftypefn  {Built-in Function} {[@var{L}, @var{U}] =} __ilu0__ (@var{A})\n\
+@deftypefnx {Built-in Function} {[@var{L}, @var{U}] =} __ilu0__ (@var{A}, @var{milu})\n\
+@deftypefnx {Built-in Function} {[@var{L}, @var{U}, @var{P}] =} __ilu0__ (@var{A}, @dots{})\n\
+Undocumented internal function.\n\
+@end deftypefn")
+{
+  octave_value_list retval;
+
+  int nargin = args.length ();
+  std::string milu;
+
+  if (nargout > 2 || nargin < 1 || nargin > 2)
+    {
+      print_usage ();
+      return retval;
+    }
+
+  // In ILU0 algorithm the zero-pattern of the input matrix is preserved so
+  // it's structure does not change during the algorithm.  The same input
+  // matrix is used to build the output matrix due to that fact.
+  octave_value_list param_list;
+  if (! args(0).is_complex_type ())
+    {
+      SparseMatrix sm = args(0).sparse_matrix_value ();
+      ilu_0 <SparseMatrix, double> (sm, milu);
+      if (!error_state)
+        {
+          param_list.append (sm);
+          retval(1) = feval ("triu", param_list)(0).sparse_matrix_value ();
+          SparseMatrix eye =
+            feval ("speye", octave_value_list (
+                     octave_value (sm.cols ())))(0).sparse_matrix_value ();
+          param_list.append (-1);
+          retval(0) = eye +
+                      feval ("tril", param_list)(0).sparse_matrix_value ();
+        }
+    }
+  else
+    {
+      SparseComplexMatrix sm = args(0).sparse_complex_matrix_value ();
+      ilu_0 <SparseComplexMatrix, Complex> (sm, milu);
+      if (! error_state)
+        {
+          param_list.append (sm);
+          retval(1) =
+            feval ("triu", param_list)(0).sparse_complex_matrix_value ();
+          SparseComplexMatrix eye =
+            feval ("speye", octave_value_list (
+                     octave_value (sm.cols ())))(0).sparse_complex_matrix_value ();
+          param_list.append (-1);
+          retval(0) =
+            eye + feval ("tril", param_list)(0).sparse_complex_matrix_value ();
+        }
+    }
+
+  return retval;
+}
+
+template <typename octave_matrix_t, typename T>
+void ilu_crout (octave_matrix_t& sm_l, octave_matrix_t& sm_u,
+                octave_matrix_t& L, octave_matrix_t& U, T* cols_norm,
+                T* rows_norm, const T droptol = 0,
+                const std::string milu = "off")
+{
+
+  // Map the strings into chars for faster comparing inside loops
+  char opt;
+  enum {OFF, ROW, COL};
+  if (milu == "row")
+    opt = ROW;
+  else if (milu == "col")
+    opt = COL;
+  else
+    opt = OFF;
+
+  octave_idx_type jrow, i, j, k, jj, total_len_l, total_len_u, max_len_u,
+                  max_len_l, w_len_u, w_len_l, cols_list_len, rows_list_len;
+
+  const octave_idx_type n = sm_u.cols ();
+  sm_u = sm_u.transpose ();
+
+  max_len_u = sm_u.nnz ();
+  max_len_u += (0.1 * max_len_u) > n ? 0.1 * max_len_u : n;
+  max_len_l = sm_l.nnz ();
+  max_len_l += (0.1 * max_len_l) > n ? 0.1 * max_len_l : n;
+  // Extract pointers to the arrays for faster access inside loops
+  octave_idx_type* cidx_in_u = sm_u.cidx ();
+  octave_idx_type* ridx_in_u = sm_u.ridx ();
+  T* data_in_u = sm_u.data ();
+  octave_idx_type* cidx_in_l = sm_l.cidx ();
+  octave_idx_type* ridx_in_l = sm_l.ridx ();
+  T* data_in_l = sm_l.data ();
+
+  // L output arrays
+  Array <octave_idx_type> ridx_out_l (dim_vector (max_len_l, 1));
+  octave_idx_type* ridx_l = ridx_out_l.fortran_vec ();
+  Array <T> data_out_l (dim_vector (max_len_l, 1));
+  T* data_l = data_out_l.fortran_vec ();
+
+  // U output arrays
+  Array <octave_idx_type> ridx_out_u (dim_vector (max_len_u, 1));
+  octave_idx_type* ridx_u = ridx_out_u.fortran_vec ();
+  Array <T> data_out_u (dim_vector (max_len_u, 1));
+  T* data_u = data_out_u.fortran_vec ();
+
+  // Working arrays
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, cidx_l, n + 1);
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, cidx_u, n + 1);
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, cols_list, n);
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, rows_list, n);
+  OCTAVE_LOCAL_BUFFER (T, w_data_l, n);
+  OCTAVE_LOCAL_BUFFER (T, w_data_u, n);
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, Ufirst, n);
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, Lfirst, n);
+  OCTAVE_LOCAL_BUFFER (T, cr_sum, n);
+
+  T zero = T (0);
+
+  cidx_u[0] = cidx_in_u[0];
+  cidx_l[0] = cidx_in_l[0];
+  for (i = 0; i < n; i++)
+    {
+      w_data_u[i] = zero;
+      w_data_l[i] = zero;
+      cr_sum[i] = zero;
+    }
+
+  total_len_u = 0;
+  total_len_l = 0;
+  cols_list_len = 0;
+  rows_list_len = 0;
+
+  for (k = 0; k < n; k++)
+    {
+      // Load the working column and working row
+      for (i = cidx_in_l[k]; i < cidx_in_l[k+1]; i++)
+        w_data_l[ridx_in_l[i]] = data_in_l[i];
+
+      for (i = cidx_in_u[k]; i < cidx_in_u[k+1]; i++)
+        w_data_u[ridx_in_u[i]] = data_in_u[i];
+
+      // Update U working row
+      for (j = 0; j < rows_list_len; j++)
+        {
+          if ((Ufirst[rows_list[j]] != -1))
+            for (jj = Ufirst[rows_list[j]]; jj < cidx_u[rows_list[j]+1]; jj++)
+              {
+                jrow = ridx_u[jj];
+                w_data_u[jrow] -= data_u[jj] * data_l[Lfirst[rows_list[j]]];
+              }
+        }
+      // Update L working column
+      for (j = 0; j < cols_list_len; j++)
+        {
+          if (Lfirst[cols_list[j]] != -1)
+            for (jj = Lfirst[cols_list[j]]; jj < cidx_l[cols_list[j]+1]; jj++)
+              {
+                jrow = ridx_l[jj];
+                w_data_l[jrow] -= data_l[jj] * data_u[Ufirst[cols_list[j]]];
+              }
+        }
+
+      if ((max_len_u - total_len_u) < n)
+        {
+          max_len_u += (0.1 * max_len_u) > n ? 0.1 * max_len_u : n;
+          data_out_u.resize (dim_vector (max_len_u, 1));
+          data_u = data_out_u.fortran_vec ();
+          ridx_out_u.resize (dim_vector (max_len_u, 1));
+          ridx_u = ridx_out_u.fortran_vec ();
+        }
+
+      if ((max_len_l - total_len_l) < n)
+        {
+          max_len_l += (0.1 * max_len_l) > n ? 0.1 * max_len_l : n;
+          data_out_l.resize (dim_vector (max_len_l, 1));
+          data_l = data_out_l.fortran_vec ();
+          ridx_out_l.resize (dim_vector (max_len_l, 1));
+          ridx_l = ridx_out_l.fortran_vec ();
+        }
+
+      // Expand the working row into the U output data structures
+      w_len_l = 0;
+      data_u[total_len_u] = w_data_u[k];
+      ridx_u[total_len_u] = k;
+      w_len_u = 1;
+      for (i = k + 1; i < n; i++)
+        {
+          if (w_data_u[i] != zero)
+            {
+              if (std::abs (w_data_u[i]) < (droptol * rows_norm[k]))
+                {
+                  if (opt == ROW)
+                    cr_sum[k] += w_data_u[i];
+                  else if (opt == COL)
+                    cr_sum[i] += w_data_u[i];
+                }
+              else
+                {
+                  data_u[total_len_u + w_len_u] = w_data_u[i];
+                  ridx_u[total_len_u + w_len_u] = i;
+                  w_len_u++;
+                }
+            }
+
+          // Expand the working column into the L output data structures
+          if (w_data_l[i] != zero)
+            {
+              if (std::abs (w_data_l[i]) < (droptol * cols_norm[k]))
+                {
+                  if (opt == COL)
+                    cr_sum[k] += w_data_l[i];
+                  else if (opt == ROW)
+                    cr_sum[i] += w_data_l[i];
+                }
+              else
+                {
+                  data_l[total_len_l + w_len_l] = w_data_l[i];
+                  ridx_l[total_len_l + w_len_l] = i;
+                  w_len_l++;
+                }
+            }
+          w_data_u[i] = zero;
+          w_data_l[i] = zero;
+        }
+
+      // Compensate row and column sums --> milu option
+      if (opt == COL || opt == ROW)
+        data_u[total_len_u] += cr_sum[k];
+
+      // Check if the pivot is zero
+      if (data_u[total_len_u] == zero)
+        {
+          error ("ilu: encountered a pivot equal to 0");
+          break;
+        }
+
+      // Scale the elements in L by the pivot
+      for (i = total_len_l ; i < (total_len_l + w_len_l); i++)
+        data_l[i] /= data_u[total_len_u];
+
+
+      total_len_u += w_len_u;
+      total_len_l += w_len_l;
+      // Check if there are too many elements to be indexed with
+      // octave_idx_type type due to fill-in during the process.
+      if (total_len_l < 0 || total_len_u < 0)
+        {
+          error ("ilu: integer overflow.  Too many fill-in elements in L or U");
+          break;
+        }
+      cidx_u[k+1] = cidx_u[k] - cidx_u[0] + w_len_u;
+      cidx_l[k+1] = cidx_l[k] - cidx_l[0] + w_len_l;
+
+      // The tricky part of the algorithm.  The arrays pointing to the first
+      // working element of each column in the next iteration (Lfirst) or
+      // the first working element of each row (Ufirst) are updated.  Also the
+      // arrays working as lists cols_list and rows_list are filled with
+      // indices pointing to Ufirst and Lfirst respectively.
+      // TODO: Maybe the -1 indicating in Ufirst and Lfirst, that no elements
+      // have to be considered in a certain column or row in next iteration,
+      // can be removed.  It feels safer to me using such an indicator.
+      if (k < (n - 1))
+        {
+          if (w_len_u > 0)
+            Ufirst[k] = cidx_u[k];
+          else
+            Ufirst[k] = -1;
+          if (w_len_l > 0)
+            Lfirst[k] = cidx_l[k];
+          else
+            Lfirst[k] = -1;
+          cols_list_len = 0;
+          rows_list_len = 0;
+          for (i = 0; i <= k; i++)
+            {
+              if (Ufirst[i] != -1)
+                {
+                  jj = ridx_u[Ufirst[i]];
+                  if (jj < (k + 1))
+                    {
+                      if (Ufirst[i] < (cidx_u[i+1]))
+                        {
+                          Ufirst[i]++;
+                          if (Ufirst[i] == cidx_u[i+1])
+                            Ufirst[i] = -1;
+                          else
+                            jj = ridx_u[Ufirst[i]];
+                        }
+                    }
+                  if (jj == (k + 1))
+                    {
+                      cols_list[cols_list_len] = i;
+                      cols_list_len++;
+                    }
+                }
+
+              if (Lfirst[i] != -1)
+                {
+                  jj = ridx_l[Lfirst[i]];
+                  if (jj < (k + 1))
+                    if (Lfirst[i] < (cidx_l[i+1]))
+                      {
+                        Lfirst[i]++;
+                        if (Lfirst[i] == cidx_l[i+1])
+                          Lfirst[i] = -1;
+                        else
+                          jj = ridx_l[Lfirst[i]];
+                      }
+                  if (jj == (k + 1))
+                    {
+                      rows_list[rows_list_len] = i;
+                      rows_list_len++;
+                    }
+                }
+            }
+        }
+    }
+
+  if (! error_state)
+    {
+      // Build the output matrices
+      L = octave_matrix_t (n, n, total_len_l);
+      U = octave_matrix_t (n, n, total_len_u);
+      for (i = 0; i <= n; i++)
+        L.cidx (i) = cidx_l[i];
+      for (i = 0; i < total_len_l; i++)
+        {
+          L.ridx (i) = ridx_l[i];
+          L.data (i) = data_l[i];
+        }
+      for (i = 0; i <= n; i++)
+        U.cidx (i) = cidx_u[i];
+      for (i = 0; i < total_len_u; i++)
+        {
+          U.ridx (i) = ridx_u[i];
+          U.data (i) = data_u[i];
+        }
+      U = U.transpose ();
+    }
+}
+
+DEFUN (__iluc__, args, nargout,
+       "-*- texinfo -*-\n\
+@deftypefn  {Built-in Function} {[@var{L}, @var{U}] =} __iluc__ (@var{A})\n\
+@deftypefnx {Built-in Function} {[@var{L}, @var{U}] =} __iluc__ (@var{A}, @var{droptol}) \n\
+@deftypefnx {Built-in Function} {[@var{L}, @var{U}] =} __iluc__ (@var{A}, @var{droptol}, @var{milu})\n\
+@deftypefnx {Built-in Function} {[@var{L}, @var{U}, @var{P}] =} __iluc__ (@var{A}, @dots{})\n\
+Undocumented internal function.\n\
+@end deftypefn")
+{
+  octave_value_list retval;
+  int nargin = args.length ();
+  std::string milu = "off";
+  double droptol = 0;
+
+  if (nargout != 2 || nargin < 1 || nargin > 3)
+    {
+      print_usage ();
+      return retval;
+    }
+
+  // Don't repeat input validation of arguments done in ilu.m
+  if (nargin >= 2)
+    droptol = args(1).double_value ();
+
+  if (nargin == 3)
+    milu = args(2).string_value ();
+
+  octave_value_list param_list;
+  if (! args(0).is_complex_type ())
+    {
+      Array<double> cols_norm, rows_norm;
+      param_list.append (args(0).sparse_matrix_value ());
+      SparseMatrix sm_u = feval ("triu", param_list)(0).sparse_matrix_value ();
+      param_list.append (-1);
+      SparseMatrix sm_l = feval ("tril", param_list)(0).sparse_matrix_value ();
+      param_list(1) = "rows";
+      rows_norm = feval ("norm", param_list)(0).vector_value ();
+      param_list(1) = "cols";
+      cols_norm = feval ("norm", param_list)(0).vector_value ();
+      param_list.clear ();
+      SparseMatrix U;
+      SparseMatrix L;
+      ilu_crout <SparseMatrix, double> (sm_l, sm_u, L, U,
+                                        cols_norm.fortran_vec (),
+                                        rows_norm.fortran_vec (),
+                                        droptol, milu);
+      if (! error_state)
+        {
+          param_list.append (octave_value (L.cols ()));
+          SparseMatrix eye =
+            feval ("speye", param_list)(0).sparse_matrix_value ();
+          retval(1) = U;
+          retval(0) = L + eye;
+        }
+    }
+  else
+    {
+      Array<Complex> cols_norm, rows_norm;
+      param_list.append (args(0).sparse_complex_matrix_value ());
+      SparseComplexMatrix sm_u =
+        feval("triu", param_list)(0).sparse_complex_matrix_value ();
+      param_list.append (-1);
+      SparseComplexMatrix sm_l =
+        feval("tril", param_list)(0).sparse_complex_matrix_value ();
+      param_list(1) = "rows";
+      rows_norm = feval ("norm", param_list)(0).complex_vector_value ();
+      param_list(1) = "cols";
+      cols_norm = feval ("norm", param_list)(0).complex_vector_value ();
+      param_list.clear ();
+      SparseComplexMatrix U;
+      SparseComplexMatrix L;
+      ilu_crout < SparseComplexMatrix, Complex >
+                (sm_l, sm_u, L, U, cols_norm.fortran_vec () ,
+                 rows_norm.fortran_vec (), Complex (droptol), milu);
+      if (! error_state)
+        {
+          param_list.append (octave_value (L.cols ()));
+          SparseComplexMatrix eye =
+            feval ("speye", param_list)(0).sparse_complex_matrix_value ();
+          retval(1) = U;
+          retval(0) = L + eye;
+        }
+    }
+
+  return retval;
+}
+
+// That function implements the IKJ and JKI variants of gaussian elimination
+// to perform the ILUTP decomposition.  The behaviour is controlled by milu
+// parameter.  If milu = ['off'|'col'] the JKI version is performed taking
+// advantage of CCS format of the input matrix.  Row pivoting is performed.
+// If milu = 'row' the input matrix has to be transposed to obtain the
+// equivalent CRS structure so we can work efficiently with rows.  In that
+// case IKJ version is used and column pivoting is performed.
+
+template <typename octave_matrix_t, typename T>
+void ilu_tp (octave_matrix_t& sm, octave_matrix_t& L, octave_matrix_t& U,
+             octave_idx_type nnz_u, octave_idx_type nnz_l, T* cols_norm,
+             Array <octave_idx_type>& perm_vec, const T droptol = T(0),
+             const T thresh = T(0), const  std::string milu = "off",
+             const double udiag = 0)
+{
+  char opt;
+  enum {OFF, ROW, COL};
+  if (milu == "row")
+    opt = ROW;
+  else if (milu == "col")
+    opt = COL;
+  else
+    opt = OFF;
+
+  const octave_idx_type n = sm.cols ();
+
+  // That is necessary for the JKI (milu = "row") variant.
+  if (opt == ROW)
+    sm = sm.transpose();
+
+  // Extract pointers to the arrays for faster access inside loops
+  octave_idx_type* cidx_in = sm.cidx ();
+  octave_idx_type* ridx_in = sm.ridx ();
+  T* data_in = sm.data ();
+  octave_idx_type jrow, i, j, k, jj, c, total_len_l, total_len_u, p_perm,
+                  max_ind, max_len_l, max_len_u;
+  T zero = T(0);
+
+  T tl = zero, aux, maximum;
+
+  max_len_u = nnz_u;
+  max_len_u += (0.1 * max_len_u) > n ? 0.1 * max_len_u : n;
+  max_len_l = nnz_l;
+  max_len_l += (0.1 * max_len_l) > n ? 0.1 * max_len_l : n;
+
+  Array <octave_idx_type> cidx_out_l (dim_vector (n + 1, 1));
+  octave_idx_type* cidx_l = cidx_out_l.fortran_vec ();
+  Array <octave_idx_type> ridx_out_l (dim_vector (max_len_l, 1));
+  octave_idx_type* ridx_l = ridx_out_l.fortran_vec ();
+  Array <T> data_out_l (dim_vector (max_len_l ,1));
+  T* data_l = data_out_l.fortran_vec ();
+  // Data for U
+  Array <octave_idx_type> cidx_out_u (dim_vector (n + 1, 1));
+  octave_idx_type* cidx_u = cidx_out_u.fortran_vec ();
+  Array <octave_idx_type> ridx_out_u (dim_vector (max_len_u, 1));
+  octave_idx_type* ridx_u = ridx_out_u.fortran_vec ();
+  Array <T> data_out_u (dim_vector (max_len_u, 1));
+  T* data_u = data_out_u.fortran_vec();
+
+  // Working arrays and permutation arrays
+  octave_idx_type w_len_u, w_len_l;
+  T total_sum, partial_col_sum = zero, partial_row_sum = zero;
+  std::set <octave_idx_type> iw_l;
+  std::set <octave_idx_type> iw_u;
+  std::set <octave_idx_type>::iterator it, it2;
+  OCTAVE_LOCAL_BUFFER (T, w_data, n);
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, iperm, n);
+  octave_idx_type* perm = perm_vec.fortran_vec ();
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, uptr, n);
+
+
+  cidx_l[0] = cidx_in[0];
+  cidx_u[0] = cidx_in[0];
+  for (i = 0; i < n; i++)
+    {
+      w_data[i] = 0;
+      perm[i] = i;
+      iperm[i] = i;
+    }
+  total_len_u = 0;
+  total_len_l = 0;
+
+  for (k = 0; k < n; k++)
+    {
+
+      for (j = cidx_in[k]; j < cidx_in[k+1]; j++)
+        {
+          p_perm = iperm[ridx_in[j]];
+          w_data[iperm[ridx_in[j]]] = data_in[j];
+          if (p_perm > k)
+            iw_l.insert (ridx_in[j]);
+          else
+            iw_u.insert (p_perm);
+        }
+
+      it = iw_u.begin ();
+      jrow = *it;
+      total_sum = zero;
+      while ((jrow < k) && (it != iw_u.end ()))
+        {
+          if (opt == COL)
+            partial_col_sum = w_data[jrow];
+          if (w_data[jrow] != zero)
+            {
+              if (opt == ROW)
+                {
+                  partial_row_sum = w_data[jrow];
+                  tl = w_data[jrow] / data_u[uptr[jrow]];
+                }
+              for (jj = cidx_l[jrow]; jj < cidx_l[jrow+1]; jj++)
+                {
+                  p_perm = iperm[ridx_l[jj]];
+                  aux = w_data[p_perm];
+                  if (opt == ROW)
+                    {
+                      w_data[p_perm] -= tl * data_l[jj];
+                      partial_row_sum += tl * data_l[jj];
+                    }
+                  else
+                    {
+                      tl = data_l[jj] * w_data[jrow];
+                      w_data[p_perm] -= tl;
+                      if (opt == COL)
+                        partial_col_sum += tl;
+                    }
+
+                  if ((aux == zero) && (w_data[p_perm] != zero))
+                    {
+                      if (p_perm > k)
+                        iw_l.insert (ridx_l[jj]);
+                      else
+                        iw_u.insert (p_perm);
+                    }
+                }
+
+              // Drop element from the U part in IKJ and L part in JKI
+              // variant (milu = [col|off])
+              if ((std::abs (w_data[jrow]) < (droptol * cols_norm[k]))
+                  && (w_data[jrow] != zero))
+                {
+                  if (opt == COL)
+                    total_sum += partial_col_sum;
+                  else if (opt == ROW)
+                    total_sum += partial_row_sum;
+                  w_data[jrow] = zero;
+                  it2 = it;
+                  it++;
+                  iw_u.erase (it2);
+                  jrow = *it;
+                  continue;
+                }
+              else
+                // This is the element scaled by the pivot
+                // in the actual iteration
+                if (opt == ROW)
+                  w_data[jrow] = tl;
+            }
+          jrow = *(++it);
+        }
+
+      // Search for the pivot and update iw_l and iw_u if the pivot is not the
+      // diagonal element
+      if ((thresh > zero) && (k < (n - 1)))
+        {
+          maximum = std::abs (w_data[k]) / thresh;
+          max_ind = perm[k];
+          for (it = iw_l.begin (); it != iw_l.end (); ++it)
+            {
+              p_perm = iperm[*it];
+              if (std::abs (w_data[p_perm]) > maximum)
+                {
+                  maximum = std::abs (w_data[p_perm]);
+                  max_ind = *it;
+                  it2 = it;
+                }
+            }
+          // If the pivot is not the diagonal element update all.
+          p_perm = iperm[max_ind];
+          if (max_ind != perm[k])
+            {
+              iw_l.erase (it2);
+              if (w_data[k] != zero)
+                iw_l.insert (perm[k]);
+              else
+                iw_u.insert (k);
+              // Swap data and update permutation vectors
+              aux = w_data[k];
+              iperm[perm[p_perm]] = k;
+              iperm[perm[k]] = p_perm;
+              c = perm[k];
+              perm[k] = perm[p_perm];
+              perm[p_perm] = c;
+              w_data[k] = w_data[p_perm];
+              w_data[p_perm] = aux;
+            }
+
+        }
+
+      // Drop elements in the L part in the IKJ and from the U part in the JKI
+      // version.
+      it = iw_l.begin ();
+      while (it != iw_l.end ())
+        {
+          p_perm = iperm[*it];
+          if (droptol > zero)
+            if (std::abs (w_data[p_perm]) < (droptol * cols_norm[k]))
+              {
+                if (opt != OFF)
+                  total_sum += w_data[p_perm];
+                w_data[p_perm] = zero;
+                it2 = it;
+                it++;
+                iw_l.erase (it2);
+                continue;
+              }
+
+          it++;
+        }
+
+      // If milu == [row|col] summation is preserved.
+      // Compensate diagonal element.
+      if (opt != OFF)
+        {
+          if ((total_sum > zero) && (w_data[k] == zero))
+            iw_u.insert (k);
+          w_data[k] += total_sum;
+        }
+
+
+
+      // Check if the pivot is zero and if udiag is active.
+      // NOTE: If the pivot == 0 and udiag is active, then if milu = [col|row]
+      //       will not preserve the row sum for that column/row.
+      if (w_data[k] == zero)
+        {
+          if (udiag == 1)
+            {
+              w_data[k] = droptol;
+              iw_u.insert (k);
+            }
+          else
+            {
+              error ("ilu: encountered a pivot equal to 0");
+              break;
+            }
+        }
+
+      // Scale the elements on the L part for IKJ version (milu = [col|off])
+      if (opt != ROW)
+        for (it = iw_l.begin (); it != iw_l.end (); ++it)
+          {
+            p_perm = iperm[*it];
+            w_data[p_perm] = w_data[p_perm] / w_data[k];
+          }
+
+
+      if ((max_len_u - total_len_u) < n)
+        {
+          max_len_u += (0.1 * max_len_u) > n ? 0.1 * max_len_u : n;
+          data_out_u.resize (dim_vector (max_len_u, 1));
+          data_u = data_out_u.fortran_vec ();
+          ridx_out_u.resize (dim_vector (max_len_u, 1));
+          ridx_u = ridx_out_u.fortran_vec ();
+        }
+
+      if ((max_len_l - total_len_l) < n)
+        {
+          max_len_l += (0.1 * max_len_l) > n ? 0.1 * max_len_l : n;
+          data_out_l.resize (dim_vector (max_len_l, 1));
+          data_l = data_out_l.fortran_vec ();
+          ridx_out_l.resize (dim_vector (max_len_l, 1));
+          ridx_l = ridx_out_l.fortran_vec ();
+        }
+
+      // Expand working vector into U.
+      w_len_u = 0;
+      for (it = iw_u.begin (); it != iw_u.end (); ++it)
+        {
+          if (w_data[*it] != zero)
+            {
+              data_u[total_len_u + w_len_u] = w_data[*it];
+              ridx_u[total_len_u + w_len_u] = *it;
+              w_len_u++;
+            }
+          w_data[*it] = 0;
+        }
+      // Expand working vector into L.
+      w_len_l = 0;
+      for (it = iw_l.begin (); it != iw_l.end (); ++it)
+        {
+          p_perm = iperm[*it];
+          if (w_data[p_perm] != zero)
+            {
+              data_l[total_len_l + w_len_l] = w_data[p_perm];
+              ridx_l[total_len_l + w_len_l] = *it;
+              w_len_l++;
+            }
+          w_data[p_perm] = 0;
+        }
+      total_len_u += w_len_u;
+      total_len_l += w_len_l;
+      // Check if there are too many elements to be indexed with
+      // octave_idx_type type due to fill-in during the process.
+      if (total_len_l < 0 || total_len_u < 0)
+        {
+          error ("ilu: Integer overflow.  Too many fill-in elements in L or U");
+          break;
+        }
+      if (opt == ROW)
+        uptr[k] = total_len_u - 1;
+      cidx_u[k+1] = cidx_u[k] - cidx_u[0] + w_len_u;
+      cidx_l[k+1] = cidx_l[k] - cidx_l[0] + w_len_l;
+
+      iw_l.clear ();
+      iw_u.clear ();
+    }
+
+  if (! error_state)
+    {
+      octave_matrix_t *L_ptr;
+      octave_matrix_t *U_ptr;
+      octave_matrix_t diag (n, n, n);
+
+      // L and U are interchanged if milu = 'row'.  It is a matter
+      // of nomenclature to re-use code with both IKJ and JKI
+      // versions of the algorithm.
+      if (opt == ROW)
+        {
+          L_ptr = &U;
+          U_ptr = &L;
+          L = octave_matrix_t (n, n, total_len_u - n);
+          U = octave_matrix_t (n, n, total_len_l);
+        }
+      else
+        {
+          L_ptr = &L;
+          U_ptr = &U;
+          L = octave_matrix_t (n, n, total_len_l);
+          U = octave_matrix_t (n, n, total_len_u);
+        }
+
+      for (i = 0; i <= n; i++)
+        {
+          L_ptr->cidx (i) = cidx_l[i];
+          U_ptr->cidx (i) = cidx_u[i];
+          if (opt == ROW)
+            U_ptr->cidx (i) -= i;
+        }
+
+      for (i = 0; i < n; i++)
+        {
+          if (opt == ROW)
+            diag.elem (i,i) = data_u[uptr[i]];
+          j = cidx_l[i];
+
+          while (j < cidx_l[i+1])
+            {
+              L_ptr->ridx (j) = ridx_l[j];
+              L_ptr->data (j) = data_l[j];
+              j++;
+            }
+          j = cidx_u[i];
+
+          while (j < cidx_u[i+1])
+            {
+              c = j;
+              if (opt == ROW)
+                {
+                  // The diagonal is removed from L if milu = 'row'.
+                  // That is because is convenient to have it inside
+                  // the L part to carry out the process.
+                  if (ridx_u[j] == i)
+                    {
+                      j++;
+                      continue;
+                    }
+                  else
+                    c -= i;
+                }
+              U_ptr->data (c) = data_u[j];
+              U_ptr->ridx (c) = ridx_u[j];
+              j++;
+            }
+        }
+
+      if (opt == ROW)
+        {
+          U = U.transpose ();
+          // The diagonal, conveniently permuted is added to U
+          U += diag.index (idx_vector::colon, perm_vec);
+          L = L.transpose ();
+        }
+    }
+}
+
+DEFUN (__ilutp__, args, nargout,
+       "-*- texinfo -*-\n\
+@deftypefn  {Built-in Function} {[@var{L}, @var{U}] =} __ilutp__ (@var{A})\n\
+@deftypefnx {Built-in Function} {[@var{L}, @var{U}] =} __ilutp__ (@var{A}, @var{droptol})\n\
+@deftypefnx {Built-in Function} {[@var{L}, @var{U}] =} __ilutp__ (@var{A}, @var{droptol}, @var{thresh})\n\
+@deftypefnx {Built-in Function} {[@var{L}, @var{U}] =} __ilutp__ (@var{A}, @var{droptol}, @var{thresh}, @var{milu})\n\
+@deftypefnx {Built-in Function} {[@var{L}, @var{U}] =} __ilutp__ (@var{A}, @var{droptol}, @var{thresh}, @var{milu}, @var{udiag})\n\
+@deftypefnx {Built-in Function} {[@var{L}, @var{U}, @var{P}] =} __ilutp__ (@var{A}, @dots{})\n\
+Undocumented internal function.\n\
+@end deftypefn")
+{
+  octave_value_list retval;
+
+  int nargin = args.length ();
+  std::string milu = "";
+  double droptol = 0, thresh = 1;
+  double udiag = 0;
+
+  if (nargout < 2 || nargout > 3 || nargin < 1 || nargin > 5)
+    {
+      print_usage ();
+      return retval;
+    }
+
+  // Don't repeat input validation of arguments done in ilu.m
+  if (nargin >= 2)
+    droptol = args(1).double_value ();
+
+  if (nargin >= 3)
+    thresh = args(2).double_value ();
+
+  if (nargin >= 4)
+    milu = args(3).string_value ();
+
+  if (nargin == 5)
+    udiag = args(4).double_value ();
+
+  octave_value_list param_list;
+  octave_idx_type nnz_u, nnz_l;
+  if (! args(0).is_complex_type ())
+    {
+      Array <double> rc_norm;
+      SparseMatrix sm = args(0).sparse_matrix_value ();
+      param_list.append (sm);
+      nnz_u =  (feval ("triu", param_list)(0).sparse_matrix_value ()).nnz ();
+      param_list.append (-1);
+      nnz_l =  (feval ("tril", param_list)(0).sparse_matrix_value ()).nnz ();
+      if (milu == "row")
+        param_list (1) = "rows";
+      else
+        param_list (1) = "cols";
+      rc_norm = feval ("norm", param_list)(0).vector_value ();
+      param_list.clear ();
+      Array <octave_idx_type> perm (dim_vector (sm.cols (), 1));
+      SparseMatrix U;
+      SparseMatrix L;
+      ilu_tp <SparseMatrix, double> (sm, L, U, nnz_u, nnz_l,
+                                     rc_norm.fortran_vec (),
+                                     perm, droptol, thresh, milu, udiag);
+      if (! error_state)
+        {
+          param_list.append (octave_value (L.cols ()));
+          SparseMatrix eye =
+            feval ("speye", param_list)(0).sparse_matrix_value ();
+          if (milu == "row")
+            {
+              if (nargout == 3)
+                {
+                  retval(2) = eye.index (idx_vector::colon, perm);
+                  retval(1) = U.index (idx_vector::colon, perm);
+                }
+              else if (nargout == 2)
+                retval(1) = U;
+              retval(0) = L + eye;
+            }
+          else
+            {
+              if (nargout == 3)
+                {
+                  retval(2) = eye.index (perm, idx_vector::colon);
+                  retval(1) = U;
+                  retval(0) = L.index (perm, idx_vector::colon) + eye;
+                }
+              else
+                {
+                  retval(1) = U;
+                  retval(0) = L + eye.index (perm, idx_vector::colon);
+                }
+            }
+        }
+    }
+  else
+    {
+      Array <Complex> rc_norm;
+      SparseComplexMatrix sm = args(0).sparse_complex_matrix_value ();
+      param_list.append (sm);
+      nnz_u =
+        feval ("triu", param_list)(0).sparse_complex_matrix_value ().nnz ();
+      param_list.append (-1);
+      nnz_l =
+        feval ("tril", param_list)(0).sparse_complex_matrix_value ().nnz ();
+      if (milu == "row")
+        param_list(1) = "rows";
+      else
+        param_list(1) = "cols";
+      rc_norm = feval ("norm", param_list)(0).complex_vector_value ();
+      Array <octave_idx_type> perm (dim_vector (sm.cols (), 1));
+      param_list.clear ();
+      SparseComplexMatrix U;
+      SparseComplexMatrix L;
+      ilu_tp < SparseComplexMatrix, Complex>
+              (sm, L, U, nnz_u, nnz_l, rc_norm.fortran_vec (), perm,
+               Complex (droptol), Complex (thresh), milu, udiag);
+
+      if (! error_state)
+        {
+          param_list.append (octave_value (L.cols ()));
+          SparseComplexMatrix eye =
+            feval ("speye", param_list)(0).sparse_complex_matrix_value ();
+          if (milu == "row")
+            {
+              if (nargout == 3)
+                {
+                  retval(2) = eye.index (idx_vector::colon, perm);
+                  retval(1) = U.index (idx_vector::colon, perm);
+                }
+              else if (nargout == 2)
+                retval(1) = U;
+              retval(0) = L + eye;
+            }
+          else
+            {
+              if (nargout == 3)
+                {
+                  retval(2) = eye.index (perm, idx_vector::colon);
+                  retval(1) = U;
+                  retval(0) = L.index (perm, idx_vector::colon) + eye;
+                }
+              else
+                {
+                  retval(1) = U;
+                  retval(0) = L + eye.index (perm, idx_vector::colon);
+                }
+            }
+        }
+    }
+
+  return retval;
+}
+
+/*
+## No test needed for internal helper function.
+%!assert (1)
+*/
+

--- a/libinterp/corefcn/module.mk	Sat Feb 28 07:42:26 2015 -0800
+++ b/libinterp/corefcn/module.mk	Fri Feb 27 19:44:28 2015 -0500
@@ -138,6 +138,9 @@
   corefcn/Cell.cc \
   corefcn/__contourc__.cc \
   corefcn/__dispatch__.cc \
+  corefcn/__dsearchn__.cc \
+  corefcn/__ichol__.cc \
+  corefcn/__ilu__.cc \
   corefcn/__lin_interpn__.cc \
   corefcn/__pchip_deriv__.cc \
   corefcn/__qp__.cc \
@@ -255,8 +258,9 @@
   corefcn/time.cc \
   corefcn/toplev.cc \
   corefcn/tril.cc \
+  corefcn/tsearch.cc \
+  corefcn/txt-eng-ft.cc \
   corefcn/txt-eng.cc \
-  corefcn/txt-eng-ft.cc \
   corefcn/typecast.cc \
   corefcn/urlwrite.cc \
   corefcn/utils.cc \

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libinterp/corefcn/tsearch.cc	Fri Feb 27 19:44:28 2015 -0500
@@ -0,0 +1,183 @@
+/*
+
+Copyright (C) 2002-2015 Andreas Stahel
+
+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/>.
+
+*/
+
+// Author: Andreas Stahel <Andreas.Stahel@hta-bi.bfh.ch>
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "lo-ieee.h"
+#include "lo-math.h"
+
+#include "defun.h"
+#include "error.h"
+#include "oct-obj.h"
+
+inline double max (double a, double b, double c)
+{
+  if (a < b)
+    return (b < c ? c : b);
+  else
+    return (a < c ? c : a);
+}
+
+inline double min (double a, double b, double c)
+{
+  if (a > b)
+    return (b > c ? c : b);
+  else
+    return (a > c ? c : a);
+}
+
+#define REF(x,k,i) x(static_cast<octave_idx_type>(elem((k), (i))) - 1)
+
+// for large data set the algorithm is very slow
+// one should presort (how?) either the elements of the points of evaluation
+// to cut down the time needed to decide which triangle contains the
+// given point
+
+// e.g., build up a neighbouring triangle structure and use a simplex-like
+// method to traverse it
+
+DEFUN (tsearch, args, ,
+       "-*- texinfo -*-\n\
+@deftypefn {Built-in Function} {@var{idx} =} tsearch (@var{x}, @var{y}, @var{t}, @var{xi}, @var{yi})\n\
+Search for the enclosing Delaunay convex hull.  For @code{@var{t} =\n\
+delaunay (@var{x}, @var{y})}, finds the index in @var{t} containing the\n\
+points @code{(@var{xi}, @var{yi})}.  For points outside the convex hull,\n\
+@var{idx} is NaN.\n\
+@seealso{delaunay, delaunayn}\n\
+@end deftypefn")
+{
+  const double eps=1.0e-12;
+
+  octave_value_list retval;
+  const int nargin = args.length ();
+  if (nargin != 5)
+    {
+      print_usage ();
+      return retval;
+    }
+
+  const ColumnVector x (args(0).vector_value ());
+  const ColumnVector y (args(1).vector_value ());
+  const Matrix elem (args(2).matrix_value ());
+  const ColumnVector xi (args(3).vector_value ());
+  const ColumnVector yi (args(4).vector_value ());
+
+  if (error_state)
+    return retval;
+
+  const octave_idx_type nelem = elem.rows ();
+
+  ColumnVector minx (nelem);
+  ColumnVector maxx (nelem);
+  ColumnVector miny (nelem);
+  ColumnVector maxy (nelem);
+  for (octave_idx_type k = 0; k < nelem; k++)
+    {
+      minx(k) = min (REF (x, k, 0), REF (x, k, 1), REF (x, k, 2)) - eps;
+      maxx(k) = max (REF (x, k, 0), REF (x, k, 1), REF (x, k, 2)) + eps;
+      miny(k) = min (REF (y, k, 0), REF (y, k, 1), REF (y, k, 2)) - eps;
+      maxy(k) = max (REF (y, k, 0), REF (y, k, 1), REF (y, k, 2)) + eps;
+    }
+
+  const octave_idx_type np = xi.length ();
+  ColumnVector values (np);
+
+  double x0, y0, a11, a12, a21, a22, det;
+  x0 = y0 = 0.0;
+  a11 = a12 = a21 = a22 = 0.0;
+  det = 0.0;
+
+  octave_idx_type k = nelem; // k is a counter of elements
+  for (octave_idx_type kp = 0; kp < np; kp++)
+    {
+      const double xt = xi(kp);
+      const double yt = yi(kp);
+
+      // check if last triangle contains the next point
+      if (k < nelem)
+        {
+          const double dx1 = xt - x0;
+          const double dx2 = yt - y0;
+          const double c1 = (a22 * dx1 - a21 * dx2) / det;
+          const double c2 = (-a12 * dx1 + a11 * dx2) / det;
+          if (c1 >= -eps && c2 >= -eps && (c1 + c2) <= (1 + eps))
+            {
+              values(kp) = double(k+1);
+              continue;
+            }
+        }
+
+      // it doesn't, so go through all elements
+      for (k = 0; k < nelem; k++)
+        {
+          OCTAVE_QUIT;
+          if (xt >= minx(k) && xt <= maxx(k) && yt >= miny(k) && yt <= maxy(k))
+            {
+              // element inside the minimum rectangle: examine it closely
+              x0  = REF (x, k, 0);
+              y0  = REF (y, k, 0);
+              a11 = REF (x, k, 1) - x0;
+              a12 = REF (y, k, 1) - y0;
+              a21 = REF (x, k, 2) - x0;
+              a22 = REF (y, k, 2) - y0;
+              det = a11 * a22 - a21 * a12;
+
+              // solve the system
+              const double dx1 = xt - x0;
+              const double dx2 = yt - y0;
+              const double c1 = (a22 * dx1 - a21 * dx2) / det;
+              const double c2 = (-a12 * dx1 + a11 * dx2) / det;
+              if ((c1 >= -eps) && (c2 >= -eps) && ((c1 + c2) <= (1 + eps)))
+                {
+                  values(kp) = double(k+1);
+                  break;
+                }
+            } //endif # examine this element closely
+        } //endfor # each element
+
+      if (k == nelem)
+        values(kp) = lo_ieee_nan_value ();
+
+    } //endfor # kp
+
+  retval(0) = values;
+
+  return retval;
+}
+
+/*
+%!shared x, y, tri
+%! x = [-1;-1;1];
+%! y = [-1;1;-1];
+%! tri = [1, 2, 3];
+%!assert (tsearch (x,y,tri,-1,-1), 1)
+%!assert (tsearch (x,y,tri, 1,-1), 1)
+%!assert (tsearch (x,y,tri,-1, 1), 1)
+%!assert (tsearch (x,y,tri,-1/3, -1/3), 1)
+%!assert (tsearch (x,y,tri, 1, 1), NaN)
+
+%!error tsearch ()
+*/

--- a/libinterp/dldfcn/__dsearchn__.cc	Sat Feb 28 07:42:26 2015 -0800
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,115 +0,0 @@
-/*
-
-Copyright (C) 2007-2015 David Bateman
-
-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 <iostream>
-#include <fstream>
-#include <string>
-
-#include "lo-math.h"
-
-#include "defun-dld.h"
-#include "error.h"
-#include "oct-obj.h"
-
-DEFUN_DLD (__dsearchn__, args, ,
-           "-*- texinfo -*-\n\
-@deftypefn {Loadable Function} {[@var{idx}, @var{d}] =} dsearch (@var{x}, @var{xi})\n\
-Undocumented internal function.\n\
-@end deftypefn")
-{
-  int nargin = args.length ();
-  octave_value_list retval;
-
-  if (nargin != 2)
-    {
-      print_usage ();
-      return retval;
-    }
-
-  Matrix x = args(0).matrix_value ().transpose ();
-  Matrix xi = args(1).matrix_value ().transpose ();
-
-  if (! error_state)
-    {
-      if (x.rows () != xi.rows () || x.columns () < 1)
-        error ("__dsearch__: number of rows of X and XI must match");
-      else
-        {
-          octave_idx_type n = x.rows ();
-          octave_idx_type nx = x.columns ();
-          octave_idx_type nxi = xi.columns ();
-
-          ColumnVector idx (nxi);
-          double *pidx = idx.fortran_vec ();
-          ColumnVector dist (nxi);
-          double *pdist = dist.fortran_vec ();
-
-#define DIST(dd, y, yi, m) \
-  dd = 0.; \
-  for (octave_idx_type k = 0; k < m; k++) \
-   { \
-     double yd = y[k] - yi[k]; \
-     dd += yd * yd; \
-   } \
-  dd = sqrt (dd);
-
-          const double *pxi = xi.fortran_vec ();
-          for (octave_idx_type i = 0; i < nxi; i++)
-            {
-              double d0;
-              const double *px = x.fortran_vec ();
-              DIST(d0, px, pxi, n);
-              *pidx = 1.;
-              for (octave_idx_type j = 1; j < nx; j++)
-                {
-                  px += n;
-                  double d;
-                  DIST (d, px, pxi, n);
-                  if (d < d0)
-                    {
-                      d0 = d;
-                      *pidx = static_cast<double>(j + 1);
-                    }
-                  OCTAVE_QUIT;
-                }
-
-              *pdist++ = d0;
-              pidx++;
-              pxi += n;
-            }
-
-          retval(1) = dist;
-          retval(0) = idx;
-        }
-    }
-
-  return retval;
-}
-
-/*
-## No test needed for internal helper function.
-%!assert (1)
-*/

--- a/libinterp/dldfcn/__ichol__.cc	Sat Feb 28 07:42:26 2015 -0800
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,516 +0,0 @@
-/*
-
-Copyright (C) 2014-2015 Eduardo Ramos Fernández <eduradical951@gmail.com>
-Copyright (C) 2013-2015 Kai T. Ohlhus <k.ohlhus@gmail.com>
-
-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 "oct-locbuf.h"
-
-#include "defun-dld.h"
-#include "parse.h"
-
-// Secondary functions for complex and real case used in ichol algorithms.
-Complex ichol_mult_complex (Complex a, Complex b)
-{
-#if defined (HAVE_CXX_COMPLEX_SETTERS)
-  b.imag (-std::imag (b));
-#elif defined (HAVE_CXX_COMPLEX_REFERENCE_ACCESSORS)
-  b.imag () = -std::imag (b);
-#else
-  b = std::conj (b);
-#endif
-  return a * b;
-}
-
-double ichol_mult_real (double a, double b)
-{
-  return a * b;
-}
-
-bool ichol_checkpivot_complex (Complex pivot)
-{
-  if (pivot.imag () != 0)
-    {
-      error ("ichol: non-real pivot encountered.  The matrix must be hermitian.");
-      return false;
-    }
-  else if (pivot.real () < 0)
-    {
-      error ("ichol: negative pivot encountered");
-      return false;
-    }
-  return true;
-}
-
-bool ichol_checkpivot_real (double pivot)
-{
-  if (pivot < 0)
-    {
-      error ("ichol: negative pivot encountered");
-      return false;
-    }
-  return true;
-}
-
-template <typename octave_matrix_t, typename T, T (*ichol_mult) (T, T),
-          bool (*ichol_checkpivot) (T)>
-void ichol_0 (octave_matrix_t& sm, const std::string michol = "off")
-{
-
-  const octave_idx_type n = sm.cols ();
-  octave_idx_type j1, jend, j2, jrow, jjrow, j, jw, i, k, jj, r;
-  T tl;
-  char opt;
-  enum {OFF, ON};
-  if (michol == "on")
-    opt = ON;
-  else
-    opt = OFF;
-
-  // Input matrix pointers
-  octave_idx_type* cidx = sm.cidx ();
-  octave_idx_type* ridx = sm.ridx ();
-  T* data = sm.data ();
-
-  // Working arrays
-  OCTAVE_LOCAL_BUFFER (octave_idx_type, Lfirst, n);
-  OCTAVE_LOCAL_BUFFER (octave_idx_type, Llist, n);
-  OCTAVE_LOCAL_BUFFER (octave_idx_type, iw, n);
-  OCTAVE_LOCAL_BUFFER (T, dropsums, n);
-
-  // Initialize working arrays
-  for (i = 0; i < n; i++)
-    {
-      iw[i] = -1;
-      Llist[i] = -1;
-      Lfirst[i] = -1;
-      dropsums[i] = 0;
-    }
-
-  // Main loop
-  for (k = 0; k < n; k++)
-    {
-      j1 = cidx[k];
-      j2 = cidx[k+1];
-      for (j = j1; j < j2; j++)
-        iw[ridx[j]] = j;
-
-      jrow = Llist [k];
-      // Iterate over each non-zero element in the actual row.
-      while (jrow != -1)
-        {
-          jjrow = Lfirst[jrow];
-          jend = cidx[jrow+1];
-          for (jj = jjrow; jj < jend; jj++)
-            {
-              r = ridx[jj];
-              jw = iw[r];
-              tl = ichol_mult (data[jj], data[jjrow]);
-              if (jw != -1)
-                data[jw] -= tl;
-              else
-                // Because of the symmetry of the matrix, we know
-                // the drops in the column r are also in the column k.
-                if (opt == ON)
-                  {
-                    dropsums[r] -= tl;
-                    dropsums[k] -= tl;
-                  }
-            }
-          // Update the linked list and the first entry of the actual column.
-          if ((jjrow + 1) < jend)
-            {
-              Lfirst[jrow]++;
-              j = jrow;
-              jrow = Llist[jrow];
-              Llist[j] = Llist[ridx[Lfirst[j]]];
-              Llist[ridx[Lfirst[j]]] = j;
-            }
-          else
-            jrow = Llist[jrow];
-        }
-
-      if (opt == ON)
-        data[j1] += dropsums[k];
-
-      if (ridx[j1] != k)
-        {
-          error ("ichol: encountered a pivot equal to 0");
-          break;
-        }
-
-      if (! ichol_checkpivot (data[j1]))
-        break;
-
-      data[cidx[k]] = std::sqrt (data[j1]);
-
-      // Update Llist and Lfirst with the k-column information.  Also,
-      // scale the column elements by the pivot and reset the working array iw.
-      if (k < (n - 1))
-        {
-          iw[ridx[j1]] = -1;
-          for (i = j1 + 1; i < j2; i++)
-            {
-              iw[ridx[i]] = -1;
-              data[i] /= data[j1];
-            }
-          Lfirst[k] = j1;
-          if ((Lfirst[k] + 1) < j2)
-            {
-              Lfirst[k]++;
-              jjrow = ridx[Lfirst[k]];
-              Llist[k] = Llist[jjrow];
-              Llist[jjrow] = k;
-            }
-        }
-    }
-}
-
-DEFUN_DLD (__ichol0__, args, nargout, "-*- texinfo -*-\n\
-@deftypefn  {Loadable Function} {@var{L} =} __ichol0__ (@var{A})\n\
-@deftypefnx {Loadable Function} {@var{L} =} __ichol0__ (@var{A}, @var{michol})\n\
-Undocumented internal function.\n\
-@end deftypefn")
-
-{
-  octave_value_list retval;
-
-  int nargin = args.length ();
-  std::string michol = "off";
-
-  if (nargout > 1 || nargin < 1 || nargin > 2)
-    {
-      print_usage ();
-      return retval;
-    }
-
-  if (nargin == 2)
-    michol = args(1).string_value ();
-
-  // In ICHOL0 algorithm the zero-pattern of the input matrix is preserved
-  // so it's structure does not change during the algorithm.  The same input
-  // matrix is used to build the output matrix due to that fact.
-  octave_value_list param_list;
-  if (!args(0).is_complex_type ())
-    {
-      SparseMatrix sm = args(0).sparse_matrix_value ();
-      param_list.append (sm);
-      sm = feval ("tril", param_list)(0).sparse_matrix_value ();
-      ichol_0 <SparseMatrix, double, ichol_mult_real,
-               ichol_checkpivot_real> (sm, michol);
-      if (! error_state)
-        retval(0) = sm;
-    }
-  else
-    {
-      SparseComplexMatrix sm = args(0).sparse_complex_matrix_value ();
-      param_list.append (sm);
-      sm = feval ("tril", param_list)(0).sparse_complex_matrix_value ();
-      ichol_0 <SparseComplexMatrix, Complex, ichol_mult_complex,
-               ichol_checkpivot_complex> (sm, michol);
-      if (! error_state)
-        retval(0) = sm;
-    }
-
-  return retval;
-}
-
-template <typename octave_matrix_t, typename T,  T (*ichol_mult) (T, T),
-          bool (*ichol_checkpivot) (T)>
-void ichol_t (const octave_matrix_t& sm, octave_matrix_t& L, const T* cols_norm,
-              const T droptol, const std::string michol = "off")
-
-{
-
-  const octave_idx_type n = sm.cols ();
-  octave_idx_type j, jrow, jend, jjrow, i, k, jj, total_len,
-                  w_len, max_len, ind;
-  char opt;
-  enum {OFF, ON};
-  if (michol == "on")
-    opt = ON;
-  else
-    opt = OFF;
-
-  // Input matrix pointers
-  octave_idx_type* cidx = sm.cidx ();
-  octave_idx_type* ridx = sm.ridx ();
-  T* data = sm.data ();
-
-  // Output matrix data structures.  Because the final zero pattern pattern of
-  // the output matrix is not known due to fill-in elements, a heuristic
-  // approach has been adopted for memory allocation.  The size of ridx_out_l
-  // and data_out_l is incremented 10% of their actual size (nnz (A) in the
-  // beginning).  If that amount is less than n, their size is just incremented
-  // in n elements.  This way the number of reallocations decreases throughout
-  // the process, obtaining a good performance.
-  max_len = sm.nnz ();
-  max_len += (0.1 * max_len) > n ? 0.1 * max_len : n;
-  Array <octave_idx_type> cidx_out_l (dim_vector (n + 1, 1));
-  octave_idx_type* cidx_l = cidx_out_l.fortran_vec ();
-  Array <octave_idx_type> ridx_out_l (dim_vector (max_len ,1));
-  octave_idx_type* ridx_l = ridx_out_l.fortran_vec ();
-  Array <T> data_out_l (dim_vector (max_len, 1));
-  T* data_l = data_out_l.fortran_vec ();
-
-  // Working arrays
-  OCTAVE_LOCAL_BUFFER (T, w_data, n);
-  OCTAVE_LOCAL_BUFFER (octave_idx_type, Lfirst, n);
-  OCTAVE_LOCAL_BUFFER (octave_idx_type, Llist, n);
-  OCTAVE_LOCAL_BUFFER (T, col_drops, n);
-  std::vector <octave_idx_type> vec;
-  vec.resize (n);
-
-  T zero = T (0);
-  cidx_l[0] = cidx[0];
-  for (i = 0; i < n; i++)
-    {
-      Llist[i] = -1;
-      Lfirst[i] = -1;
-      w_data[i] = 0;
-      col_drops[i] = zero;
-      vec[i] = 0;
-    }
-
-  total_len = 0;
-  for (k = 0; k < n; k++)
-    {
-      ind = 0;
-      for (j = cidx[k]; j < cidx[k+1]; j++)
-        {
-          w_data[ridx[j]] = data[j];
-          if (ridx[j] != k)
-            {
-              vec[ind] = ridx[j];
-              ind++;
-            }
-        }
-      jrow = Llist[k];
-      while (jrow != -1)
-        {
-          jjrow = Lfirst[jrow];
-          jend = cidx_l[jrow+1];
-          for (jj = jjrow; jj < jend; jj++)
-            {
-              j = ridx_l[jj];
-              // If the element in the j position of the row is zero,
-              // then it will become non-zero, so we add it to the
-              // vector that tracks non-zero elements in the working row.
-              if (w_data[j] == zero)
-                {
-                  vec[ind] = j;
-                  ind++;
-                }
-              w_data[j] -=  ichol_mult (data_l[jj], data_l[jjrow]);
-            }
-          // Update the actual column first element and
-          // update the linked list of the jrow row.
-          if ((jjrow + 1) < jend)
-            {
-              Lfirst[jrow]++;
-              j = jrow;
-              jrow = Llist[jrow];
-              Llist[j] = Llist[ridx_l[Lfirst[j]]];
-              Llist[ridx_l[Lfirst[j]]] = j;
-            }
-          else
-            jrow = Llist[jrow];
-        }
-
-      // Resizing output arrays
-      if ((max_len - total_len) < n)
-        {
-          max_len += (0.1 * max_len) > n ? 0.1 * max_len : n;
-          data_out_l.resize (dim_vector (max_len, 1));
-          data_l = data_out_l.fortran_vec ();
-          ridx_out_l.resize (dim_vector (max_len, 1));
-          ridx_l = ridx_out_l.fortran_vec ();
-        }
-
-      // The sorting of the non-zero elements of the working column can be
-      // handled in a couple of ways.  The most efficient two I found, are
-      // keeping the elements in an ordered binary search tree dynamically or
-      // keep them unsorted in a vector and at the end of the outer iteration
-      // order them.  The last approach exhibits lower execution times.
-      std::sort (vec.begin (), vec.begin () + ind);
-
-      data_l[total_len] = w_data[k];
-      ridx_l[total_len] = k;
-      w_len = 1;
-
-      // Extract the non-zero elements of working column and
-      // drop the elements that are lower than droptol * cols_norm[k].
-      for (i = 0; i < ind ; i++)
-        {
-          jrow = vec[i];
-          if (w_data[jrow] != zero)
-            {
-              if (std::abs (w_data[jrow]) < (droptol * cols_norm[k]))
-                {
-                  if (opt == ON)
-                    {
-                      col_drops[k] += w_data[jrow];
-                      col_drops[jrow] += w_data[jrow];
-                    }
-                }
-              else
-                {
-                  data_l[total_len + w_len] = w_data[jrow];
-                  ridx_l[total_len + w_len] = jrow;
-                  w_len++;
-                }
-              vec[i] = 0;
-            }
-          w_data[jrow] = zero;
-        }
-
-      // Compensate column sums --> michol option
-      if (opt == ON)
-        data_l[total_len] += col_drops[k];
-
-      if (data_l[total_len] == zero)
-        {
-          error ("ichol: encountered a pivot equal to 0");
-          break;
-        }
-      else if (! ichol_checkpivot (data_l[total_len]))
-        break;
-
-      // Once elements are dropped and compensation of column sums are done,
-      // scale the elements by the pivot.
-      data_l[total_len] = std::sqrt (data_l[total_len]);
-      for (jj = total_len + 1; jj < (total_len + w_len); jj++)
-        data_l[jj] /=  data_l[total_len];
-      total_len += w_len;
-      // Check if there are too many elements to be indexed with
-      // octave_idx_type type due to fill-in during the process.
-      if (total_len < 0)
-        {
-          error ("ichol: integer overflow.  Too many fill-in elements in L");
-          break;
-        }
-      cidx_l[k+1] = cidx_l[k] - cidx_l[0] + w_len;
-
-      // Update Llist and Lfirst with the k-column information.
-      if (k < (n - 1))
-        {
-          Lfirst[k] = cidx_l[k];
-          if ((Lfirst[k] + 1) < cidx_l[k+1])
-            {
-              Lfirst[k]++;
-              jjrow = ridx_l[Lfirst[k]];
-              Llist[k] = Llist[jjrow];
-              Llist[jjrow] = k;
-            }
-        }
-    }
-
-  if (! error_state)
-    {
-      // Build the output matrices
-      L = octave_matrix_t (n, n, total_len);
-      for (i = 0; i <= n; i++)
-        L.cidx (i) = cidx_l[i];
-      for (i = 0; i < total_len; i++)
-        {
-          L.ridx (i) = ridx_l[i];
-          L.data (i) = data_l[i];
-        }
-    }
-}
-
-DEFUN_DLD (__icholt__, args, nargout, "-*- texinfo -*-\n\
-@deftypefn  {Loadable Function} {@var{L} =} __icholt__ (@var{A})\n\
-@deftypefnx {Loadable Function} {@var{L} =} __icholt__ (@var{A}, @var{droptol})\n\
-@deftypefnx {Loadable Function} {@var{L} =} __icholt__ (@var{A}, @var{droptol}, @var{michol})\n\
-Undocumented internal function.\n\
-@end deftypefn")
-{
-  octave_value_list retval;
-  int nargin = args.length ();
-  // Default values of parameters
-  std::string michol = "off";
-  double droptol = 0;
-
-  if (nargout > 1 || nargin < 1 || nargin > 3)
-    {
-      print_usage ();
-      return retval;
-    }
-
-  // Don't repeat input validation of arguments done in ichol.m
-
-  if (nargin >= 2)
-    droptol = args(1).double_value ();
-
-  if (nargin == 3)
-    michol = args(2).string_value ();
-
-  octave_value_list param_list;
-  if (! args(0).is_complex_type ())
-    {
-      Array <double> cols_norm;
-      SparseMatrix L;
-      param_list.append (args(0).sparse_matrix_value ());
-      SparseMatrix sm_l =
-        feval ("tril", param_list)(0).sparse_matrix_value ();
-      param_list(0) = sm_l;
-      param_list(1) = 1;
-      param_list(2) = "cols";
-      cols_norm = feval ("norm", param_list)(0).vector_value ();
-      param_list.clear ();
-      ichol_t <SparseMatrix,
-               double, ichol_mult_real, ichol_checkpivot_real>
-               (sm_l, L, cols_norm.fortran_vec (), droptol, michol);
-      if (! error_state)
-        retval(0) = L;
-    }
-  else
-    {
-      Array <Complex> cols_norm;
-      SparseComplexMatrix L;
-      param_list.append (args(0).sparse_complex_matrix_value ());
-      SparseComplexMatrix sm_l =
-        feval ("tril", param_list)(0).sparse_complex_matrix_value ();
-      param_list(0) = sm_l;
-      param_list(1) = 1;
-      param_list(2) = "cols";
-      cols_norm = feval ("norm", param_list)(0).complex_vector_value ();
-      param_list.clear ();
-      ichol_t <SparseComplexMatrix,
-               Complex, ichol_mult_complex, ichol_checkpivot_complex>
-               (sm_l, L, cols_norm.fortran_vec (),
-                Complex (droptol), michol);
-      if (! error_state)
-        retval(0) = L;
-    }
-
-  return retval;
-}
-
-/*
-## No test needed for internal helper function.
-%!assert (1)
-*/
-

--- a/libinterp/dldfcn/__ilu__.cc	Sat Feb 28 07:42:26 2015 -0800
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,1104 +0,0 @@
-/*
-
-Copyright (C) 2014-2015 Eduardo Ramos Fernández <eduradical951@gmail.com>
-Copyright (C) 2013-2015 Kai T. Ohlhus <k.ohlhus@gmail.com>
-
-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 "oct-locbuf.h"
-
-#include "defun-dld.h"
-#include "parse.h"
-
-// That function implements the IKJ and JKI variants of Gaussian elimination to
-// perform the ILUTP decomposition.  The behaviour is controlled by milu
-// parameter.  If milu = ['off'|'col'] the JKI version is performed taking
-// advantage of CCS format of the input matrix.  If milu = 'row' the input
-// matrix has to be transposed to obtain the equivalent CRS structure so we can
-// work efficiently with rows.  In this case IKJ version is used.
-template <typename octave_matrix_t, typename T>
-void ilu_0 (octave_matrix_t& sm, const std::string milu = "off")
-{
-
-  const octave_idx_type n = sm.cols ();
-  OCTAVE_LOCAL_BUFFER (octave_idx_type, iw, n);
-  OCTAVE_LOCAL_BUFFER (octave_idx_type, uptr, n);
-  octave_idx_type j1, j2, jrow, jw, i, k, jj;
-  T tl, r;
-
-  enum {OFF, ROW, COL};
-  char opt;
-  if (milu == "row")
-    {
-      opt = ROW;
-      sm = sm.transpose ();
-    }
-  else if (milu == "col")
-    opt = COL;
-  else
-    opt = OFF;
-
-  octave_idx_type* cidx = sm.cidx ();
-  octave_idx_type* ridx = sm.ridx ();
-  T* data = sm.data ();
-  for (i = 0; i < n; i++)
-    iw[i] = -1;
-  for (k = 0; k < n; k++)
-    {
-      j1 = cidx[k];
-      j2 = cidx[k+1] - 1;
-      octave_idx_type j;
-      for (j = j1; j <= j2; j++)
-        {
-          iw[ridx[j]] = j;
-        }
-      r = 0;
-      j = j1;
-      jrow = ridx[j];
-      while ((jrow < k) && (j <= j2))
-        {
-          if (opt == ROW)
-            {
-              tl = data[j] / data[uptr[jrow]];
-              data[j] = tl;
-            }
-          for (jj = uptr[jrow] + 1; jj < cidx[jrow+1]; jj++)
-            {
-              jw = iw[ridx[jj]];
-              if (jw != -1)
-                if (opt == ROW)
-                  data[jw] -= tl * data[jj];
-                else
-                  data[jw] -= data[j] * data[jj];
-
-              else
-                // That is for the milu='row'
-                if (opt == ROW)
-                  r += tl * data[jj];
-                else if (opt == COL)
-                  r += data[j] * data[jj];
-            }
-          j++;
-          jrow = ridx[j];
-        }
-      uptr[k] = j;
-      if (opt != OFF)
-        data[uptr[k]] -= r;
-      if (opt != ROW)
-        for (jj = uptr[k] + 1; jj < cidx[k+1]; jj++)
-          data[jj] /=  data[uptr[k]];
-      if (k != jrow)
-        {
-          error ("ilu: A has a zero on the diagonal");
-          break;
-        }
-
-      if (data[j] == T(0))
-        {
-          error ("ilu: encountered a pivot equal to 0");
-          break;
-        }
-      for (i = j1; i <= j2; i++)
-        iw[ridx[i]] = -1;
-    }
-  if (opt == ROW)
-    sm = sm.transpose ();
-}
-
-DEFUN_DLD (__ilu0__, args, nargout, "-*- texinfo -*-\n\
-@deftypefn  {Loadable Function} {[@var{L}, @var{U}] =} __ilu0__ (@var{A})\n\
-@deftypefnx {Loadable Function} {[@var{L}, @var{U}] =} __ilu0__ (@var{A}, @var{milu})\n\
-@deftypefnx {Loadable Function} {[@var{L}, @var{U}, @var{P}] =} __ilu0__ (@var{A}, @dots{})\n\
-Undocumented internal function.\n\
-@end deftypefn")
-{
-  octave_value_list retval;
-
-  int nargin = args.length ();
-  std::string milu;
-
-  if (nargout > 2 || nargin < 1 || nargin > 2)
-    {
-      print_usage ();
-      return retval;
-    }
-
-  // In ILU0 algorithm the zero-pattern of the input matrix is preserved so
-  // it's structure does not change during the algorithm.  The same input
-  // matrix is used to build the output matrix due to that fact.
-  octave_value_list param_list;
-  if (! args(0).is_complex_type ())
-    {
-      SparseMatrix sm = args(0).sparse_matrix_value ();
-      ilu_0 <SparseMatrix, double> (sm, milu);
-      if (!error_state)
-        {
-          param_list.append (sm);
-          retval(1) = feval ("triu", param_list)(0).sparse_matrix_value ();
-          SparseMatrix eye =
-            feval ("speye", octave_value_list (
-                     octave_value (sm.cols ())))(0).sparse_matrix_value ();
-          param_list.append (-1);
-          retval(0) = eye +
-                      feval ("tril", param_list)(0).sparse_matrix_value ();
-        }
-    }
-  else
-    {
-      SparseComplexMatrix sm = args(0).sparse_complex_matrix_value ();
-      ilu_0 <SparseComplexMatrix, Complex> (sm, milu);
-      if (! error_state)
-        {
-          param_list.append (sm);
-          retval(1) =
-            feval ("triu", param_list)(0).sparse_complex_matrix_value ();
-          SparseComplexMatrix eye =
-            feval ("speye", octave_value_list (
-                     octave_value (sm.cols ())))(0).sparse_complex_matrix_value ();
-          param_list.append (-1);
-          retval(0) =
-            eye + feval ("tril", param_list)(0).sparse_complex_matrix_value ();
-        }
-    }
-
-  return retval;
-}
-
-template <typename octave_matrix_t, typename T>
-void ilu_crout (octave_matrix_t& sm_l, octave_matrix_t& sm_u,
-                octave_matrix_t& L, octave_matrix_t& U, T* cols_norm,
-                T* rows_norm, const T droptol = 0,
-                const std::string milu = "off")
-{
-
-  // Map the strings into chars for faster comparing inside loops
-  char opt;
-  enum {OFF, ROW, COL};
-  if (milu == "row")
-    opt = ROW;
-  else if (milu == "col")
-    opt = COL;
-  else
-    opt = OFF;
-
-  octave_idx_type jrow, i, j, k, jj, total_len_l, total_len_u, max_len_u,
-                  max_len_l, w_len_u, w_len_l, cols_list_len, rows_list_len;
-
-  const octave_idx_type n = sm_u.cols ();
-  sm_u = sm_u.transpose ();
-
-  max_len_u = sm_u.nnz ();
-  max_len_u += (0.1 * max_len_u) > n ? 0.1 * max_len_u : n;
-  max_len_l = sm_l.nnz ();
-  max_len_l += (0.1 * max_len_l) > n ? 0.1 * max_len_l : n;
-  // Extract pointers to the arrays for faster access inside loops
-  octave_idx_type* cidx_in_u = sm_u.cidx ();
-  octave_idx_type* ridx_in_u = sm_u.ridx ();
-  T* data_in_u = sm_u.data ();
-  octave_idx_type* cidx_in_l = sm_l.cidx ();
-  octave_idx_type* ridx_in_l = sm_l.ridx ();
-  T* data_in_l = sm_l.data ();
-
-  // L output arrays
-  Array <octave_idx_type> ridx_out_l (dim_vector (max_len_l, 1));
-  octave_idx_type* ridx_l = ridx_out_l.fortran_vec ();
-  Array <T> data_out_l (dim_vector (max_len_l, 1));
-  T* data_l = data_out_l.fortran_vec ();
-
-  // U output arrays
-  Array <octave_idx_type> ridx_out_u (dim_vector (max_len_u, 1));
-  octave_idx_type* ridx_u = ridx_out_u.fortran_vec ();
-  Array <T> data_out_u (dim_vector (max_len_u, 1));
-  T* data_u = data_out_u.fortran_vec ();
-
-  // Working arrays
-  OCTAVE_LOCAL_BUFFER (octave_idx_type, cidx_l, n + 1);
-  OCTAVE_LOCAL_BUFFER (octave_idx_type, cidx_u, n + 1);
-  OCTAVE_LOCAL_BUFFER (octave_idx_type, cols_list, n);
-  OCTAVE_LOCAL_BUFFER (octave_idx_type, rows_list, n);
-  OCTAVE_LOCAL_BUFFER (T, w_data_l, n);
-  OCTAVE_LOCAL_BUFFER (T, w_data_u, n);
-  OCTAVE_LOCAL_BUFFER (octave_idx_type, Ufirst, n);
-  OCTAVE_LOCAL_BUFFER (octave_idx_type, Lfirst, n);
-  OCTAVE_LOCAL_BUFFER (T, cr_sum, n);
-
-  T zero = T (0);
-
-  cidx_u[0] = cidx_in_u[0];
-  cidx_l[0] = cidx_in_l[0];
-  for (i = 0; i < n; i++)
-    {
-      w_data_u[i] = zero;
-      w_data_l[i] = zero;
-      cr_sum[i] = zero;
-    }
-
-  total_len_u = 0;
-  total_len_l = 0;
-  cols_list_len = 0;
-  rows_list_len = 0;
-
-  for (k = 0; k < n; k++)
-    {
-      // Load the working column and working row
-      for (i = cidx_in_l[k]; i < cidx_in_l[k+1]; i++)
-        w_data_l[ridx_in_l[i]] = data_in_l[i];
-
-      for (i = cidx_in_u[k]; i < cidx_in_u[k+1]; i++)
-        w_data_u[ridx_in_u[i]] = data_in_u[i];
-
-      // Update U working row
-      for (j = 0; j < rows_list_len; j++)
-        {
-          if ((Ufirst[rows_list[j]] != -1))
-            for (jj = Ufirst[rows_list[j]]; jj < cidx_u[rows_list[j]+1]; jj++)
-              {
-                jrow = ridx_u[jj];
-                w_data_u[jrow] -= data_u[jj] * data_l[Lfirst[rows_list[j]]];
-              }
-        }
-      // Update L working column
-      for (j = 0; j < cols_list_len; j++)
-        {
-          if (Lfirst[cols_list[j]] != -1)
-            for (jj = Lfirst[cols_list[j]]; jj < cidx_l[cols_list[j]+1]; jj++)
-              {
-                jrow = ridx_l[jj];
-                w_data_l[jrow] -= data_l[jj] * data_u[Ufirst[cols_list[j]]];
-              }
-        }
-
-      if ((max_len_u - total_len_u) < n)
-        {
-          max_len_u += (0.1 * max_len_u) > n ? 0.1 * max_len_u : n;
-          data_out_u.resize (dim_vector (max_len_u, 1));
-          data_u = data_out_u.fortran_vec ();
-          ridx_out_u.resize (dim_vector (max_len_u, 1));
-          ridx_u = ridx_out_u.fortran_vec ();
-        }
-
-      if ((max_len_l - total_len_l) < n)
-        {
-          max_len_l += (0.1 * max_len_l) > n ? 0.1 * max_len_l : n;
-          data_out_l.resize (dim_vector (max_len_l, 1));
-          data_l = data_out_l.fortran_vec ();
-          ridx_out_l.resize (dim_vector (max_len_l, 1));
-          ridx_l = ridx_out_l.fortran_vec ();
-        }
-
-      // Expand the working row into the U output data structures
-      w_len_l = 0;
-      data_u[total_len_u] = w_data_u[k];
-      ridx_u[total_len_u] = k;
-      w_len_u = 1;
-      for (i = k + 1; i < n; i++)
-        {
-          if (w_data_u[i] != zero)
-            {
-              if (std::abs (w_data_u[i]) < (droptol * rows_norm[k]))
-                {
-                  if (opt == ROW)
-                    cr_sum[k] += w_data_u[i];
-                  else if (opt == COL)
-                    cr_sum[i] += w_data_u[i];
-                }
-              else
-                {
-                  data_u[total_len_u + w_len_u] = w_data_u[i];
-                  ridx_u[total_len_u + w_len_u] = i;
-                  w_len_u++;
-                }
-            }
-
-          // Expand the working column into the L output data structures
-          if (w_data_l[i] != zero)
-            {
-              if (std::abs (w_data_l[i]) < (droptol * cols_norm[k]))
-                {
-                  if (opt == COL)
-                    cr_sum[k] += w_data_l[i];
-                  else if (opt == ROW)
-                    cr_sum[i] += w_data_l[i];
-                }
-              else
-                {
-                  data_l[total_len_l + w_len_l] = w_data_l[i];
-                  ridx_l[total_len_l + w_len_l] = i;
-                  w_len_l++;
-                }
-            }
-          w_data_u[i] = zero;
-          w_data_l[i] = zero;
-        }
-
-      // Compensate row and column sums --> milu option
-      if (opt == COL || opt == ROW)
-        data_u[total_len_u] += cr_sum[k];
-
-      // Check if the pivot is zero
-      if (data_u[total_len_u] == zero)
-        {
-          error ("ilu: encountered a pivot equal to 0");
-          break;
-        }
-
-      // Scale the elements in L by the pivot
-      for (i = total_len_l ; i < (total_len_l + w_len_l); i++)
-        data_l[i] /= data_u[total_len_u];
-
-
-      total_len_u += w_len_u;
-      total_len_l += w_len_l;
-      // Check if there are too many elements to be indexed with
-      // octave_idx_type type due to fill-in during the process.
-      if (total_len_l < 0 || total_len_u < 0)
-        {
-          error ("ilu: integer overflow.  Too many fill-in elements in L or U");
-          break;
-        }
-      cidx_u[k+1] = cidx_u[k] - cidx_u[0] + w_len_u;
-      cidx_l[k+1] = cidx_l[k] - cidx_l[0] + w_len_l;
-
-      // The tricky part of the algorithm.  The arrays pointing to the first
-      // working element of each column in the next iteration (Lfirst) or
-      // the first working element of each row (Ufirst) are updated.  Also the
-      // arrays working as lists cols_list and rows_list are filled with
-      // indices pointing to Ufirst and Lfirst respectively.
-      // TODO: Maybe the -1 indicating in Ufirst and Lfirst, that no elements
-      // have to be considered in a certain column or row in next iteration,
-      // can be removed.  It feels safer to me using such an indicator.
-      if (k < (n - 1))
-        {
-          if (w_len_u > 0)
-            Ufirst[k] = cidx_u[k];
-          else
-            Ufirst[k] = -1;
-          if (w_len_l > 0)
-            Lfirst[k] = cidx_l[k];
-          else
-            Lfirst[k] = -1;
-          cols_list_len = 0;
-          rows_list_len = 0;
-          for (i = 0; i <= k; i++)
-            {
-              if (Ufirst[i] != -1)
-                {
-                  jj = ridx_u[Ufirst[i]];
-                  if (jj < (k + 1))
-                    {
-                      if (Ufirst[i] < (cidx_u[i+1]))
-                        {
-                          Ufirst[i]++;
-                          if (Ufirst[i] == cidx_u[i+1])
-                            Ufirst[i] = -1;
-                          else
-                            jj = ridx_u[Ufirst[i]];
-                        }
-                    }
-                  if (jj == (k + 1))
-                    {
-                      cols_list[cols_list_len] = i;
-                      cols_list_len++;
-                    }
-                }
-
-              if (Lfirst[i] != -1)
-                {
-                  jj = ridx_l[Lfirst[i]];
-                  if (jj < (k + 1))
-                    if (Lfirst[i] < (cidx_l[i+1]))
-                      {
-                        Lfirst[i]++;
-                        if (Lfirst[i] == cidx_l[i+1])
-                          Lfirst[i] = -1;
-                        else
-                          jj = ridx_l[Lfirst[i]];
-                      }
-                  if (jj == (k + 1))
-                    {
-                      rows_list[rows_list_len] = i;
-                      rows_list_len++;
-                    }
-                }
-            }
-        }
-    }
-
-  if (! error_state)
-    {
-      // Build the output matrices
-      L = octave_matrix_t (n, n, total_len_l);
-      U = octave_matrix_t (n, n, total_len_u);
-      for (i = 0; i <= n; i++)
-        L.cidx (i) = cidx_l[i];
-      for (i = 0; i < total_len_l; i++)
-        {
-          L.ridx (i) = ridx_l[i];
-          L.data (i) = data_l[i];
-        }
-      for (i = 0; i <= n; i++)
-        U.cidx (i) = cidx_u[i];
-      for (i = 0; i < total_len_u; i++)
-        {
-          U.ridx (i) = ridx_u[i];
-          U.data (i) = data_u[i];
-        }
-      U = U.transpose ();
-    }
-}
-
-DEFUN_DLD (__iluc__, args, nargout, "-*- texinfo -*-\n\
-@deftypefn  {Loadable Function} {[@var{L}, @var{U}] =} __iluc__ (@var{A})\n\
-@deftypefnx {Loadable Function} {[@var{L}, @var{U}] =} __iluc__ (@var{A}, @var{droptol}) \n\
-@deftypefnx {Loadable Function} {[@var{L}, @var{U}] =} __iluc__ (@var{A}, @var{droptol}, @var{milu})\n\
-@deftypefnx {Loadable Function} {[@var{L}, @var{U}, @var{P}] =} __iluc__ (@var{A}, @dots{})\n\
-Undocumented internal function.\n\
-@end deftypefn")
-{
-  octave_value_list retval;
-  int nargin = args.length ();
-  std::string milu = "off";
-  double droptol = 0;
-
-  if (nargout != 2 || nargin < 1 || nargin > 3)
-    {
-      print_usage ();
-      return retval;
-    }
-
-  // Don't repeat input validation of arguments done in ilu.m
-  if (nargin >= 2)
-    droptol = args(1).double_value ();
-
-  if (nargin == 3)
-    milu = args(2).string_value ();
-
-  octave_value_list param_list;
-  if (! args(0).is_complex_type ())
-    {
-      Array<double> cols_norm, rows_norm;
-      param_list.append (args(0).sparse_matrix_value ());
-      SparseMatrix sm_u = feval ("triu", param_list)(0).sparse_matrix_value ();
-      param_list.append (-1);
-      SparseMatrix sm_l = feval ("tril", param_list)(0).sparse_matrix_value ();
-      param_list(1) = "rows";
-      rows_norm = feval ("norm", param_list)(0).vector_value ();
-      param_list(1) = "cols";
-      cols_norm = feval ("norm", param_list)(0).vector_value ();
-      param_list.clear ();
-      SparseMatrix U;
-      SparseMatrix L;
-      ilu_crout <SparseMatrix, double> (sm_l, sm_u, L, U,
-                                        cols_norm.fortran_vec (),
-                                        rows_norm.fortran_vec (),
-                                        droptol, milu);
-      if (! error_state)
-        {
-          param_list.append (octave_value (L.cols ()));
-          SparseMatrix eye =
-            feval ("speye", param_list)(0).sparse_matrix_value ();
-          retval(1) = U;
-          retval(0) = L + eye;
-        }
-    }
-  else
-    {
-      Array<Complex> cols_norm, rows_norm;
-      param_list.append (args(0).sparse_complex_matrix_value ());
-      SparseComplexMatrix sm_u =
-        feval("triu", param_list)(0).sparse_complex_matrix_value ();
-      param_list.append (-1);
-      SparseComplexMatrix sm_l =
-        feval("tril", param_list)(0).sparse_complex_matrix_value ();
-      param_list(1) = "rows";
-      rows_norm = feval ("norm", param_list)(0).complex_vector_value ();
-      param_list(1) = "cols";
-      cols_norm = feval ("norm", param_list)(0).complex_vector_value ();
-      param_list.clear ();
-      SparseComplexMatrix U;
-      SparseComplexMatrix L;
-      ilu_crout < SparseComplexMatrix, Complex >
-                (sm_l, sm_u, L, U, cols_norm.fortran_vec () ,
-                 rows_norm.fortran_vec (), Complex (droptol), milu);
-      if (! error_state)
-        {
-          param_list.append (octave_value (L.cols ()));
-          SparseComplexMatrix eye =
-            feval ("speye", param_list)(0).sparse_complex_matrix_value ();
-          retval(1) = U;
-          retval(0) = L + eye;
-        }
-    }
-
-  return retval;
-}
-
-// That function implements the IKJ and JKI variants of gaussian elimination
-// to perform the ILUTP decomposition.  The behaviour is controlled by milu
-// parameter.  If milu = ['off'|'col'] the JKI version is performed taking
-// advantage of CCS format of the input matrix.  Row pivoting is performed.
-// If milu = 'row' the input matrix has to be transposed to obtain the
-// equivalent CRS structure so we can work efficiently with rows.  In that
-// case IKJ version is used and column pivoting is performed.
-
-template <typename octave_matrix_t, typename T>
-void ilu_tp (octave_matrix_t& sm, octave_matrix_t& L, octave_matrix_t& U,
-             octave_idx_type nnz_u, octave_idx_type nnz_l, T* cols_norm,
-             Array <octave_idx_type>& perm_vec, const T droptol = T(0),
-             const T thresh = T(0), const  std::string milu = "off",
-             const double udiag = 0)
-{
-  char opt;
-  enum {OFF, ROW, COL};
-  if (milu == "row")
-    opt = ROW;
-  else if (milu == "col")
-    opt = COL;
-  else
-    opt = OFF;
-
-  const octave_idx_type n = sm.cols ();
-
-  // That is necessary for the JKI (milu = "row") variant.
-  if (opt == ROW)
-    sm = sm.transpose();
-
-  // Extract pointers to the arrays for faster access inside loops
-  octave_idx_type* cidx_in = sm.cidx ();
-  octave_idx_type* ridx_in = sm.ridx ();
-  T* data_in = sm.data ();
-  octave_idx_type jrow, i, j, k, jj, c, total_len_l, total_len_u, p_perm,
-                  max_ind, max_len_l, max_len_u;
-  T zero = T(0);
-
-  T tl = zero, aux, maximum;
-
-  max_len_u = nnz_u;
-  max_len_u += (0.1 * max_len_u) > n ? 0.1 * max_len_u : n;
-  max_len_l = nnz_l;
-  max_len_l += (0.1 * max_len_l) > n ? 0.1 * max_len_l : n;
-
-  Array <octave_idx_type> cidx_out_l (dim_vector (n + 1, 1));
-  octave_idx_type* cidx_l = cidx_out_l.fortran_vec ();
-  Array <octave_idx_type> ridx_out_l (dim_vector (max_len_l, 1));
-  octave_idx_type* ridx_l = ridx_out_l.fortran_vec ();
-  Array <T> data_out_l (dim_vector (max_len_l ,1));
-  T* data_l = data_out_l.fortran_vec ();
-  // Data for U
-  Array <octave_idx_type> cidx_out_u (dim_vector (n + 1, 1));
-  octave_idx_type* cidx_u = cidx_out_u.fortran_vec ();
-  Array <octave_idx_type> ridx_out_u (dim_vector (max_len_u, 1));
-  octave_idx_type* ridx_u = ridx_out_u.fortran_vec ();
-  Array <T> data_out_u (dim_vector (max_len_u, 1));
-  T* data_u = data_out_u.fortran_vec();
-
-  // Working arrays and permutation arrays
-  octave_idx_type w_len_u, w_len_l;
-  T total_sum, partial_col_sum = zero, partial_row_sum = zero;
-  std::set <octave_idx_type> iw_l;
-  std::set <octave_idx_type> iw_u;
-  std::set <octave_idx_type>::iterator it, it2;
-  OCTAVE_LOCAL_BUFFER (T, w_data, n);
-  OCTAVE_LOCAL_BUFFER (octave_idx_type, iperm, n);
-  octave_idx_type* perm = perm_vec.fortran_vec ();
-  OCTAVE_LOCAL_BUFFER (octave_idx_type, uptr, n);
-
-
-  cidx_l[0] = cidx_in[0];
-  cidx_u[0] = cidx_in[0];
-  for (i = 0; i < n; i++)
-    {
-      w_data[i] = 0;
-      perm[i] = i;
-      iperm[i] = i;
-    }
-  total_len_u = 0;
-  total_len_l = 0;
-
-  for (k = 0; k < n; k++)
-    {
-
-      for (j = cidx_in[k]; j < cidx_in[k+1]; j++)
-        {
-          p_perm = iperm[ridx_in[j]];
-          w_data[iperm[ridx_in[j]]] = data_in[j];
-          if (p_perm > k)
-            iw_l.insert (ridx_in[j]);
-          else
-            iw_u.insert (p_perm);
-        }
-
-      it = iw_u.begin ();
-      jrow = *it;
-      total_sum = zero;
-      while ((jrow < k) && (it != iw_u.end ()))
-        {
-          if (opt == COL)
-            partial_col_sum = w_data[jrow];
-          if (w_data[jrow] != zero)
-            {
-              if (opt == ROW)
-                {
-                  partial_row_sum = w_data[jrow];
-                  tl = w_data[jrow] / data_u[uptr[jrow]];
-                }
-              for (jj = cidx_l[jrow]; jj < cidx_l[jrow+1]; jj++)
-                {
-                  p_perm = iperm[ridx_l[jj]];
-                  aux = w_data[p_perm];
-                  if (opt == ROW)
-                    {
-                      w_data[p_perm] -= tl * data_l[jj];
-                      partial_row_sum += tl * data_l[jj];
-                    }
-                  else
-                    {
-                      tl = data_l[jj] * w_data[jrow];
-                      w_data[p_perm] -= tl;
-                      if (opt == COL)
-                        partial_col_sum += tl;
-                    }
-
-                  if ((aux == zero) && (w_data[p_perm] != zero))
-                    {
-                      if (p_perm > k)
-                        iw_l.insert (ridx_l[jj]);
-                      else
-                        iw_u.insert (p_perm);
-                    }
-                }
-
-              // Drop element from the U part in IKJ and L part in JKI
-              // variant (milu = [col|off])
-              if ((std::abs (w_data[jrow]) < (droptol * cols_norm[k]))
-                  && (w_data[jrow] != zero))
-                {
-                  if (opt == COL)
-                    total_sum += partial_col_sum;
-                  else if (opt == ROW)
-                    total_sum += partial_row_sum;
-                  w_data[jrow] = zero;
-                  it2 = it;
-                  it++;
-                  iw_u.erase (it2);
-                  jrow = *it;
-                  continue;
-                }
-              else
-                // This is the element scaled by the pivot
-                // in the actual iteration
-                if (opt == ROW)
-                  w_data[jrow] = tl;
-            }
-          jrow = *(++it);
-        }
-
-      // Search for the pivot and update iw_l and iw_u if the pivot is not the
-      // diagonal element
-      if ((thresh > zero) && (k < (n - 1)))
-        {
-          maximum = std::abs (w_data[k]) / thresh;
-          max_ind = perm[k];
-          for (it = iw_l.begin (); it != iw_l.end (); ++it)
-            {
-              p_perm = iperm[*it];
-              if (std::abs (w_data[p_perm]) > maximum)
-                {
-                  maximum = std::abs (w_data[p_perm]);
-                  max_ind = *it;
-                  it2 = it;
-                }
-            }
-          // If the pivot is not the diagonal element update all.
-          p_perm = iperm[max_ind];
-          if (max_ind != perm[k])
-            {
-              iw_l.erase (it2);
-              if (w_data[k] != zero)
-                iw_l.insert (perm[k]);
-              else
-                iw_u.insert (k);
-              // Swap data and update permutation vectors
-              aux = w_data[k];
-              iperm[perm[p_perm]] = k;
-              iperm[perm[k]] = p_perm;
-              c = perm[k];
-              perm[k] = perm[p_perm];
-              perm[p_perm] = c;
-              w_data[k] = w_data[p_perm];
-              w_data[p_perm] = aux;
-            }
-
-        }
-
-      // Drop elements in the L part in the IKJ and from the U part in the JKI
-      // version.
-      it = iw_l.begin ();
-      while (it != iw_l.end ())
-        {
-          p_perm = iperm[*it];
-          if (droptol > zero)
-            if (std::abs (w_data[p_perm]) < (droptol * cols_norm[k]))
-              {
-                if (opt != OFF)
-                  total_sum += w_data[p_perm];
-                w_data[p_perm] = zero;
-                it2 = it;
-                it++;
-                iw_l.erase (it2);
-                continue;
-              }
-
-          it++;
-        }
-
-      // If milu == [row|col] summation is preserved.
-      // Compensate diagonal element.
-      if (opt != OFF)
-        {
-          if ((total_sum > zero) && (w_data[k] == zero))
-            iw_u.insert (k);
-          w_data[k] += total_sum;
-        }
-
-
-
-      // Check if the pivot is zero and if udiag is active.
-      // NOTE: If the pivot == 0 and udiag is active, then if milu = [col|row]
-      //       will not preserve the row sum for that column/row.
-      if (w_data[k] == zero)
-        {
-          if (udiag == 1)
-            {
-              w_data[k] = droptol;
-              iw_u.insert (k);
-            }
-          else
-            {
-              error ("ilu: encountered a pivot equal to 0");
-              break;
-            }
-        }
-
-      // Scale the elements on the L part for IKJ version (milu = [col|off])
-      if (opt != ROW)
-        for (it = iw_l.begin (); it != iw_l.end (); ++it)
-          {
-            p_perm = iperm[*it];
-            w_data[p_perm] = w_data[p_perm] / w_data[k];
-          }
-
-
-      if ((max_len_u - total_len_u) < n)
-        {
-          max_len_u += (0.1 * max_len_u) > n ? 0.1 * max_len_u : n;
-          data_out_u.resize (dim_vector (max_len_u, 1));
-          data_u = data_out_u.fortran_vec ();
-          ridx_out_u.resize (dim_vector (max_len_u, 1));
-          ridx_u = ridx_out_u.fortran_vec ();
-        }
-
-      if ((max_len_l - total_len_l) < n)
-        {
-          max_len_l += (0.1 * max_len_l) > n ? 0.1 * max_len_l : n;
-          data_out_l.resize (dim_vector (max_len_l, 1));
-          data_l = data_out_l.fortran_vec ();
-          ridx_out_l.resize (dim_vector (max_len_l, 1));
-          ridx_l = ridx_out_l.fortran_vec ();
-        }
-
-      // Expand working vector into U.
-      w_len_u = 0;
-      for (it = iw_u.begin (); it != iw_u.end (); ++it)
-        {
-          if (w_data[*it] != zero)
-            {
-              data_u[total_len_u + w_len_u] = w_data[*it];
-              ridx_u[total_len_u + w_len_u] = *it;
-              w_len_u++;
-            }
-          w_data[*it] = 0;
-        }
-      // Expand working vector into L.
-      w_len_l = 0;
-      for (it = iw_l.begin (); it != iw_l.end (); ++it)
-        {
-          p_perm = iperm[*it];
-          if (w_data[p_perm] != zero)
-            {
-              data_l[total_len_l + w_len_l] = w_data[p_perm];
-              ridx_l[total_len_l + w_len_l] = *it;
-              w_len_l++;
-            }
-          w_data[p_perm] = 0;
-        }
-      total_len_u += w_len_u;
-      total_len_l += w_len_l;
-      // Check if there are too many elements to be indexed with
-      // octave_idx_type type due to fill-in during the process.
-      if (total_len_l < 0 || total_len_u < 0)
-        {
-          error ("ilu: Integer overflow.  Too many fill-in elements in L or U");
-          break;
-        }
-      if (opt == ROW)
-        uptr[k] = total_len_u - 1;
-      cidx_u[k+1] = cidx_u[k] - cidx_u[0] + w_len_u;
-      cidx_l[k+1] = cidx_l[k] - cidx_l[0] + w_len_l;
-
-      iw_l.clear ();
-      iw_u.clear ();
-    }
-
-  if (! error_state)
-    {
-      octave_matrix_t *L_ptr;
-      octave_matrix_t *U_ptr;
-      octave_matrix_t diag (n, n, n);
-
-      // L and U are interchanged if milu = 'row'.  It is a matter
-      // of nomenclature to re-use code with both IKJ and JKI
-      // versions of the algorithm.
-      if (opt == ROW)
-        {
-          L_ptr = &U;
-          U_ptr = &L;
-          L = octave_matrix_t (n, n, total_len_u - n);
-          U = octave_matrix_t (n, n, total_len_l);
-        }
-      else
-        {
-          L_ptr = &L;
-          U_ptr = &U;
-          L = octave_matrix_t (n, n, total_len_l);
-          U = octave_matrix_t (n, n, total_len_u);
-        }
-
-      for (i = 0; i <= n; i++)
-        {
-          L_ptr->cidx (i) = cidx_l[i];
-          U_ptr->cidx (i) = cidx_u[i];
-          if (opt == ROW)
-            U_ptr->cidx (i) -= i;
-        }
-
-      for (i = 0; i < n; i++)
-        {
-          if (opt == ROW)
-            diag.elem (i,i) = data_u[uptr[i]];
-          j = cidx_l[i];
-
-          while (j < cidx_l[i+1])
-            {
-              L_ptr->ridx (j) = ridx_l[j];
-              L_ptr->data (j) = data_l[j];
-              j++;
-            }
-          j = cidx_u[i];
-
-          while (j < cidx_u[i+1])
-            {
-              c = j;
-              if (opt == ROW)
-                {
-                  // The diagonal is removed from L if milu = 'row'.
-                  // That is because is convenient to have it inside
-                  // the L part to carry out the process.
-                  if (ridx_u[j] == i)
-                    {
-                      j++;
-                      continue;
-                    }
-                  else
-                    c -= i;
-                }
-              U_ptr->data (c) = data_u[j];
-              U_ptr->ridx (c) = ridx_u[j];
-              j++;
-            }
-        }
-
-      if (opt == ROW)
-        {
-          U = U.transpose ();
-          // The diagonal, conveniently permuted is added to U
-          U += diag.index (idx_vector::colon, perm_vec);
-          L = L.transpose ();
-        }
-    }
-}
-
-DEFUN_DLD (__ilutp__, args, nargout, "-*- texinfo -*-\n\
-@deftypefn  {Loadable Function} {[@var{L}, @var{U}] =} __ilutp__ (@var{A})\n\
-@deftypefnx {Loadable Function} {[@var{L}, @var{U}] =} __ilutp__ (@var{A}, @var{droptol})\n\
-@deftypefnx {Loadable Function} {[@var{L}, @var{U}] =} __ilutp__ (@var{A}, @var{droptol}, @var{thresh})\n\
-@deftypefnx {Loadable Function} {[@var{L}, @var{U}] =} __ilutp__ (@var{A}, @var{droptol}, @var{thresh}, @var{milu})\n\
-@deftypefnx {Loadable Function} {[@var{L}, @var{U}] =} __ilutp__ (@var{A}, @var{droptol}, @var{thresh}, @var{milu}, @var{udiag})\n\
-@deftypefnx {Loadable Function} {[@var{L}, @var{U}, @var{P}] =} __ilutp__ (@var{A}, @dots{})\n\
-Undocumented internal function.\n\
-@end deftypefn")
-{
-  octave_value_list retval;
-
-  int nargin = args.length ();
-  std::string milu = "";
-  double droptol = 0, thresh = 1;
-  double udiag = 0;
-
-  if (nargout < 2 || nargout > 3 || nargin < 1 || nargin > 5)
-    {
-      print_usage ();
-      return retval;
-    }
-
-  // Don't repeat input validation of arguments done in ilu.m
-  if (nargin >= 2)
-    droptol = args(1).double_value ();
-
-  if (nargin >= 3)
-    thresh = args(2).double_value ();
-
-  if (nargin >= 4)
-    milu = args(3).string_value ();
-
-  if (nargin == 5)
-    udiag = args(4).double_value ();
-
-  octave_value_list param_list;
-  octave_idx_type nnz_u, nnz_l;
-  if (! args(0).is_complex_type ())
-    {
-      Array <double> rc_norm;
-      SparseMatrix sm = args(0).sparse_matrix_value ();
-      param_list.append (sm);
-      nnz_u =  (feval ("triu", param_list)(0).sparse_matrix_value ()).nnz ();
-      param_list.append (-1);
-      nnz_l =  (feval ("tril", param_list)(0).sparse_matrix_value ()).nnz ();
-      if (milu == "row")
-        param_list (1) = "rows";
-      else
-        param_list (1) = "cols";
-      rc_norm = feval ("norm", param_list)(0).vector_value ();
-      param_list.clear ();
-      Array <octave_idx_type> perm (dim_vector (sm.cols (), 1));
-      SparseMatrix U;
-      SparseMatrix L;
-      ilu_tp <SparseMatrix, double> (sm, L, U, nnz_u, nnz_l,
-                                     rc_norm.fortran_vec (),
-                                     perm, droptol, thresh, milu, udiag);
-      if (! error_state)
-        {
-          param_list.append (octave_value (L.cols ()));
-          SparseMatrix eye =
-            feval ("speye", param_list)(0).sparse_matrix_value ();
-          if (milu == "row")
-            {
-              if (nargout == 3)
-                {
-                  retval(2) = eye.index (idx_vector::colon, perm);
-                  retval(1) = U.index (idx_vector::colon, perm);
-                }
-              else if (nargout == 2)
-                retval(1) = U;
-              retval(0) = L + eye;
-            }
-          else
-            {
-              if (nargout == 3)
-                {
-                  retval(2) = eye.index (perm, idx_vector::colon);
-                  retval(1) = U;
-                  retval(0) = L.index (perm, idx_vector::colon) + eye;
-                }
-              else
-                {
-                  retval(1) = U;
-                  retval(0) = L + eye.index (perm, idx_vector::colon);
-                }
-            }
-        }
-    }
-  else
-    {
-      Array <Complex> rc_norm;
-      SparseComplexMatrix sm = args(0).sparse_complex_matrix_value ();
-      param_list.append (sm);
-      nnz_u =
-        feval ("triu", param_list)(0).sparse_complex_matrix_value ().nnz ();
-      param_list.append (-1);
-      nnz_l =
-        feval ("tril", param_list)(0).sparse_complex_matrix_value ().nnz ();
-      if (milu == "row")
-        param_list(1) = "rows";
-      else
-        param_list(1) = "cols";
-      rc_norm = feval ("norm", param_list)(0).complex_vector_value ();
-      Array <octave_idx_type> perm (dim_vector (sm.cols (), 1));
-      param_list.clear ();
-      SparseComplexMatrix U;
-      SparseComplexMatrix L;
-      ilu_tp < SparseComplexMatrix, Complex>
-              (sm, L, U, nnz_u, nnz_l, rc_norm.fortran_vec (), perm,
-               Complex (droptol), Complex (thresh), milu, udiag);
-
-      if (! error_state)
-        {
-          param_list.append (octave_value (L.cols ()));
-          SparseComplexMatrix eye =
-            feval ("speye", param_list)(0).sparse_complex_matrix_value ();
-          if (milu == "row")
-            {
-              if (nargout == 3)
-                {
-                  retval(2) = eye.index (idx_vector::colon, perm);
-                  retval(1) = U.index (idx_vector::colon, perm);
-                }
-              else if (nargout == 2)
-                retval(1) = U;
-              retval(0) = L + eye;
-            }
-          else
-            {
-              if (nargout == 3)
-                {
-                  retval(2) = eye.index (perm, idx_vector::colon);
-                  retval(1) = U;
-                  retval(0) = L.index (perm, idx_vector::colon) + eye;
-                }
-              else
-                {
-                  retval(1) = U;
-                  retval(0) = L + eye.index (perm, idx_vector::colon);
-                }
-            }
-        }
-    }
-
-  return retval;
-}
-
-/*
-## No test needed for internal helper function.
-%!assert (1)
-*/
-

--- a/libinterp/dldfcn/module-files	Sat Feb 28 07:42:26 2015 -0800
+++ b/libinterp/dldfcn/module-files	Fri Feb 27 19:44:28 2015 -0500
@@ -1,11 +1,8 @@
 # FILE|CPPFLAGS|LDFLAGS|LIBRARIES
 __delaunayn__.cc|$(QHULL_CPPFLAGS)|$(QHULL_LDFLAGS)|$(QHULL_LIBS)
-__dsearchn__.cc
 __eigs__.cc|$(ARPACK_CPPFLAGS) $(SPARSE_XCPPFLAGS)|$(ARPACK_LDFLAGS) $(SPARSE_XLDFLAGS)|$(ARPACK_LIBS) $(SPARSE_XLIBS) $(LAPACK_LIBS) $(BLAS_LIBS)
 __fltk_uigetfile__.cc|$(FLTK_CPPFLAGS) $(FT2_CPPFLAGS)|$(FLTK_LDFLAGS) $(FT2_LDFLAGS)|$(FLTK_LIBS) $(FT2_LIBS)
 __glpk__.cc|$(GLPK_CPPFLAGS)|$(GLPK_LDFLAGS)|$(GLPK_LIBS)
-__ichol__.cc
-__ilu__.cc
 __init_fltk__.cc|$(FLTK_CPPFLAGS) $(FT2_CPPFLAGS) $(FONTCONFIG_CPPFLAGS)|$(FLTK_LDFLAGS) $(FT2_LDFLAGS)|$(FLTK_LIBS) $(FT2_LIBS) $(OPENGL_LIBS)
 __init_gnuplot__.cc|$(FT2_CPPFLAGS) $(FONTCONFIG_CPPFLAGS)||
 __magick_read__.cc|$(MAGICK_CPPFLAGS)|$(MAGICK_LDFLAGS)|$(MAGICK_LIBS)
@@ -21,6 +18,5 @@
 qr.cc|$(QRUPDATE_CPPFLAGS) $(SPARSE_XCPPFLAGS)|$(QRUPDATE_LDFLAGS) $(SPARSE_XLDFLAGS)|$(QRUPDATE_LIBS) $(SPARSE_XLIBS)
 symbfact.cc|$(SPARSE_XCPPFLAGS)|$(SPARSE_XLDFLAGS)|$(SPARSE_XLIBS)
 symrcm.cc|$(SPARSE_XCPPFLAGS)|$(SPARSE_XLDFLAGS)|$(SPARSE_XLIBS)
-tsearch.cc
 audioread.cc|$(SNDFILE_CPPFLAGS)|$(SNDFILE_LDFLAGS)|$(SNDFILE_LIBS)
 audiodevinfo.cc|$(PORTAUDIO_CPPFLAGS)|$(PORTAUDIO_LDFLAGS)|$(PORTAUDIO_LIBS)

--- a/libinterp/dldfcn/tsearch.cc	Sat Feb 28 07:42:26 2015 -0800
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,188 +0,0 @@
-/*
-
-Copyright (C) 2002-2015 Andreas Stahel
-
-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/>.
-
-*/
-
-// Author: Andreas Stahel <Andreas.Stahel@hta-bi.bfh.ch>
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <iostream>
-#include <fstream>
-#include <string>
-
-#include "lo-ieee.h"
-#include "lo-math.h"
-
-#include "defun-dld.h"
-#include "error.h"
-#include "oct-obj.h"
-#include "parse.h"
-
-inline double max (double a, double b, double c)
-{
-  if (a < b)
-    return (b < c ? c : b);
-  else
-    return (a < c ? c : a);
-}
-
-inline double min (double a, double b, double c)
-{
-  if (a > b)
-    return (b > c ? c : b);
-  else
-    return (a > c ? c : a);
-}
-
-#define REF(x,k,i) x(static_cast<octave_idx_type>(elem((k), (i))) - 1)
-
-// for large data set the algorithm is very slow
-// one should presort (how?) either the elements of the points of evaluation
-// to cut down the time needed to decide which triangle contains the
-// given point
-
-// e.g., build up a neighbouring triangle structure and use a simplex-like
-// method to traverse it
-
-DEFUN_DLD (tsearch, args, ,
-           "-*- texinfo -*-\n\
-@deftypefn {Loadable Function} {@var{idx} =} tsearch (@var{x}, @var{y}, @var{t}, @var{xi}, @var{yi})\n\
-Search for the enclosing Delaunay convex hull.  For @code{@var{t} =\n\
-delaunay (@var{x}, @var{y})}, finds the index in @var{t} containing the\n\
-points @code{(@var{xi}, @var{yi})}.  For points outside the convex hull,\n\
-@var{idx} is NaN.\n\
-@seealso{delaunay, delaunayn}\n\
-@end deftypefn")
-{
-  const double eps=1.0e-12;
-
-  octave_value_list retval;
-  const int nargin = args.length ();
-  if (nargin != 5)
-    {
-      print_usage ();
-      return retval;
-    }
-
-  const ColumnVector x (args(0).vector_value ());
-  const ColumnVector y (args(1).vector_value ());
-  const Matrix elem (args(2).matrix_value ());
-  const ColumnVector xi (args(3).vector_value ());
-  const ColumnVector yi (args(4).vector_value ());
-
-  if (error_state)
-    return retval;
-
-  const octave_idx_type nelem = elem.rows ();
-
-  ColumnVector minx (nelem);
-  ColumnVector maxx (nelem);
-  ColumnVector miny (nelem);
-  ColumnVector maxy (nelem);
-  for (octave_idx_type k = 0; k < nelem; k++)
-    {
-      minx(k) = min (REF (x, k, 0), REF (x, k, 1), REF (x, k, 2)) - eps;
-      maxx(k) = max (REF (x, k, 0), REF (x, k, 1), REF (x, k, 2)) + eps;
-      miny(k) = min (REF (y, k, 0), REF (y, k, 1), REF (y, k, 2)) - eps;
-      maxy(k) = max (REF (y, k, 0), REF (y, k, 1), REF (y, k, 2)) + eps;
-    }
-
-  const octave_idx_type np = xi.length ();
-  ColumnVector values (np);
-
-  double x0, y0, a11, a12, a21, a22, det;
-  x0 = y0 = 0.0;
-  a11 = a12 = a21 = a22 = 0.0;
-  det = 0.0;
-
-  octave_idx_type k = nelem; // k is a counter of elements
-  for (octave_idx_type kp = 0; kp < np; kp++)
-    {
-      const double xt = xi(kp);
-      const double yt = yi(kp);
-
-      // check if last triangle contains the next point
-      if (k < nelem)
-        {
-          const double dx1 = xt - x0;
-          const double dx2 = yt - y0;
-          const double c1 = (a22 * dx1 - a21 * dx2) / det;
-          const double c2 = (-a12 * dx1 + a11 * dx2) / det;
-          if (c1 >= -eps && c2 >= -eps && (c1 + c2) <= (1 + eps))
-            {
-              values(kp) = double(k+1);
-              continue;
-            }
-        }
-
-      // it doesn't, so go through all elements
-      for (k = 0; k < nelem; k++)
-        {
-          OCTAVE_QUIT;
-          if (xt >= minx(k) && xt <= maxx(k) && yt >= miny(k) && yt <= maxy(k))
-            {
-              // element inside the minimum rectangle: examine it closely
-              x0  = REF (x, k, 0);
-              y0  = REF (y, k, 0);
-              a11 = REF (x, k, 1) - x0;
-              a12 = REF (y, k, 1) - y0;
-              a21 = REF (x, k, 2) - x0;
-              a22 = REF (y, k, 2) - y0;
-              det = a11 * a22 - a21 * a12;
-
-              // solve the system
-              const double dx1 = xt - x0;
-              const double dx2 = yt - y0;
-              const double c1 = (a22 * dx1 - a21 * dx2) / det;
-              const double c2 = (-a12 * dx1 + a11 * dx2) / det;
-              if ((c1 >= -eps) && (c2 >= -eps) && ((c1 + c2) <= (1 + eps)))
-                {
-                  values(kp) = double(k+1);
-                  break;
-                }
-            } //endif # examine this element closely
-        } //endfor # each element
-
-      if (k == nelem)
-        values(kp) = lo_ieee_nan_value ();
-
-    } //endfor # kp
-
-  retval(0) = values;
-
-  return retval;
-}
-
-/*
-%!shared x, y, tri
-%! x = [-1;-1;1];
-%! y = [-1;1;-1];
-%! tri = [1, 2, 3];
-%!assert (tsearch (x,y,tri,-1,-1), 1)
-%!assert (tsearch (x,y,tri, 1,-1), 1)
-%!assert (tsearch (x,y,tri,-1, 1), 1)
-%!assert (tsearch (x,y,tri,-1/3, -1/3), 1)
-%!assert (tsearch (x,y,tri, 1, 1), NaN)
-
-%!error tsearch ()
-*/