Mercurial > octave

changeset 23099:d44537a50f4b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
__ilu__ and __ichol__: use liboctave functions instead of feval (bug #50105) * __ichol__.cc(__ichol0__, __icholt__): Remove parameter checking, as it is done in ichol.m, all input arguments must be given. Replace feval calls to "tril" by "Ftril". * __ichol__.cc(__icholt__): Replace feval "norm" calls for column-wise l1-Norm by liboctave's "xcolnorms". * __ilu__.cc(__ilu0__, __iluc__, __ilutp__): Remove parameter checking, as it is done in ilu.m, all input arguments must be given. Replace feval calls to "tril" and "triu" by "Ftril" and "Ftriu". Replace feval calls to "speye" by using SparseMatrix constructor calls. Replace feval "norm" calls for column-wise l2-Norm by liboctave's "xcolnorms" and for row-wise l2-Norm by "xrownorms". * ilu.m: Avoid undefined values for "ilutp" and one return value. Refactor tests for "ilutp". Added xtest, where failure was just commented out.
author Kai T. Ohlhus <k.ohlhus@gmail.com>
date Fri, 27 Jan 2017 16:39:15 +0100
parents 03f817ed37c5
children a40e586b66ba
files libinterp/corefcn/__ichol__.cc libinterp/corefcn/__ilu__.cc scripts/sparse/ilu.m
diffstat 3 files changed, 124 insertions(+), 211 deletions(-) [+]
line wrap: on
line diff
--- a/libinterp/corefcn/__ichol__.cc	Thu Jan 26 23:46:47 2017 -0500
+++ b/libinterp/corefcn/__ichol__.cc	Fri Jan 27 16:39:15 2017 +0100
@@ -1,7 +1,7 @@
 /*
 
-Copyright (C) 2014-2016 Eduardo Ramos Fernández <eduradical951@gmail.com>
-Copyright (C) 2013-2016 Kai T. Ohlhus <k.ohlhus@gmail.com>
+Copyright (C) 2014-2017 Eduardo Ramos Fernández <eduradical951@gmail.com>
+Copyright (C) 2013-2017 Kai T. Ohlhus <k.ohlhus@gmail.com>
 
 This file is part of Octave.
 
@@ -26,10 +26,12 @@
 #endif
 
 #include "oct-locbuf.h"
+#include "oct-norm.h"
 
 #include "defun.h"
 #include "error.h"
-#include "parse.h"
+
+#include "builtin-defun-decls.h"
 
 // Secondary functions for complex and real case used in ichol algorithms.
 Complex ichol_mult_complex (Complex a, Complex b)
@@ -181,37 +183,31 @@
 
 DEFUN (__ichol0__, args, ,
        doc: /* -*- texinfo -*-
-@deftypefn  {} {@var{L} =} __ichol0__ (@var{A})
-@deftypefnx {} {@var{L} =} __ichol0__ (@var{A}, @var{michol})
+@deftypefn {} {@var{L} =} __ichol0__ (@var{A}, @var{michol})
 Undocumented internal function.
 @end deftypefn */)
 {
-  std::string michol = "off";
-  if (args.length ())
-    michol = args(1).string_value ();
+  if (args.length () != 2)
+    print_usage ();
+
+  std::string michol = args(1).string_value ();
 
   // In ICHOL0 algorithm the zero-pattern of the input matrix is preserved
   // so its 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 arg_list;
   if (! args(0).is_complex_type ())
     {
-      SparseMatrix sm = args(0).sparse_matrix_value ();
-      arg_list.append (sm);
-      sm = octave::feval ("tril", arg_list)(0).sparse_matrix_value ();
+      SparseMatrix sm = Ftril (args(0))(0).sparse_matrix_value ();
       ichol_0 <SparseMatrix, double, ichol_mult_real,
                ichol_checkpivot_real> (sm, michol);
-
       return ovl (sm);
     }
   else
     {
-      SparseComplexMatrix sm = args(0).sparse_complex_matrix_value ();
-      arg_list.append (sm);
-      sm = octave::feval ("tril", arg_list)(0).sparse_complex_matrix_value ();
+      SparseComplexMatrix sm =
+        Ftril (args(0))(0).sparse_complex_matrix_value ();
       ichol_0 <SparseComplexMatrix, Complex, ichol_mult_complex,
                ichol_checkpivot_complex> (sm, michol);
-
       return ovl (sm);
     }
 }
@@ -416,55 +412,32 @@
 
 DEFUN (__icholt__, args, ,
        doc: /* -*- texinfo -*-
-@deftypefn  {} {@var{L} =} __icholt__ (@var{A})
-@deftypefnx {} {@var{L} =} __icholt__ (@var{A}, @var{droptol})
-@deftypefnx {} {@var{L} =} __icholt__ (@var{A}, @var{droptol}, @var{michol})
+@deftypefn {} {@var{L} =} __icholt__ (@var{A}, @var{droptol}, @var{michol})
 Undocumented internal function.
 @end deftypefn */)
 {
-  int nargin = args.length ();
-  // Default values of parameters
-  std::string michol = "off";
-  double droptol = 0;
+  if (args.length () != 3)
+    print_usage ();
 
-  // Don't repeat input validation of arguments done in ichol.m
-  if (nargin >= 2)
-    droptol = args(1).double_value ();
+  double droptol = args(1).double_value ();
+  std::string michol = args(2).string_value ();
 
-  if (nargin == 3)
-    michol = args(2).string_value ();
-
-  octave_value_list arg_list;
   if (! args(0).is_complex_type ())
     {
-      Array <double> cols_norm;
       SparseMatrix L;
-      arg_list.append (args(0).sparse_matrix_value ());
-      SparseMatrix sm_l =
-        octave::feval ("tril", arg_list)(0).sparse_matrix_value ();
-      arg_list(0) = sm_l;
-      arg_list(1) = 1;
-      arg_list(2) = "cols";
-      cols_norm = octave::feval ("norm", arg_list)(0).vector_value ();
-      arg_list.clear ();
+      SparseMatrix sm_l = Ftril (args(0))(0).sparse_matrix_value ();
       ichol_t <SparseMatrix,
                double, ichol_mult_real, ichol_checkpivot_real>
-               (sm_l, L, cols_norm.fortran_vec (), droptol, michol);
+               (sm_l, L, xcolnorms (sm_l, 1).fortran_vec (), droptol, michol);
 
       return ovl (L);
     }
   else
     {
-      Array <Complex> cols_norm;
       SparseComplexMatrix L;
-      arg_list.append (args(0).sparse_complex_matrix_value ());
       SparseComplexMatrix sm_l =
-        octave::feval ("tril", arg_list)(0).sparse_complex_matrix_value ();
-      arg_list(0) = sm_l;
-      arg_list(1) = 1;
-      arg_list(2) = "cols";
-      cols_norm = octave::feval ("norm", arg_list)(0).complex_vector_value ();
-      arg_list.clear ();
+        Ftril (args(0))(0).sparse_complex_matrix_value ();
+      Array <Complex> cols_norm = xcolnorms (sm_l, 1);
       ichol_t <SparseComplexMatrix,
                Complex, ichol_mult_complex, ichol_checkpivot_complex>
                (sm_l, L, cols_norm.fortran_vec (),
--- a/libinterp/corefcn/__ilu__.cc	Thu Jan 26 23:46:47 2017 -0500
+++ b/libinterp/corefcn/__ilu__.cc	Fri Jan 27 16:39:15 2017 +0100
@@ -1,7 +1,7 @@
 /*
 
-Copyright (C) 2014-2016 Eduardo Ramos Fernández <eduradical951@gmail.com>
-Copyright (C) 2013-2016 Kai T. Ohlhus <k.ohlhus@gmail.com>
+Copyright (C) 2014-2017 Eduardo Ramos Fernández <eduradical951@gmail.com>
+Copyright (C) 2013-2017 Kai T. Ohlhus <k.ohlhus@gmail.com>
 
 This file is part of Octave.
 
@@ -26,10 +26,12 @@
 #endif
 
 #include "oct-locbuf.h"
+#include "oct-norm.h"
 
 #include "defun.h"
 #include "error.h"
-#include "parse.h"
+
+#include "builtin-defun-decls.h"
 
 // This function implements the IKJ and JKI variants of Gaussian elimination to
 // perform the ILUTP decomposition.  The behavior is controlled by milu
@@ -135,20 +137,17 @@
 
 DEFUN (__ilu0__, args, ,
        doc: /* -*- texinfo -*-
-@deftypefn  {} {[@var{L}, @var{U}] =} __ilu0__ (@var{A})
-@deftypefnx {} {[@var{L}, @var{U}] =} __ilu0__ (@var{A}, @var{milu})
-@deftypefnx {} {[@var{L}, @var{U}, @var{P}] =} __ilu0__ (@var{A}, @dots{})
+@deftypefn  {} {[@var{L}, @var{U}] =} __ilu0__ (@var{A}, @var{milu})
+@deftypefnx {} {[@var{L}, @var{U}, @var{P}] =} __ilu0__ (@var{A}, @var{milu})
 Undocumented internal function.
 @end deftypefn */)
 {
-  int nargin = args.length ();
-
-  if (nargin < 1 || nargin > 2)
+  if (args.length () != 2)
     print_usage ();
 
   octave_value_list retval (2);
 
-  std::string milu;
+  std::string milu = args(1).string_value ();
 
   // In ILU0 algorithm the zero-pattern of the input matrix is preserved so
   // its structure does not change during the algorithm.  The same input
@@ -157,28 +156,23 @@
   if (! args(0).is_complex_type ())
     {
       SparseMatrix sm = args(0).sparse_matrix_value ();
+      SparseMatrix speye (DiagMatrix (sm.cols (), sm.cols (), 1.0));
+
       ilu_0 <SparseMatrix, double> (sm, milu);
 
-      arg_list.append (sm);
-      retval(1) = octave::feval ("triu", arg_list)(0).sparse_matrix_value ();
-      SparseMatrix eye =
-        octave::feval ("speye", ovl (sm.cols ()))(0).sparse_matrix_value ();
-      arg_list.append (-1);
-      retval(0) = eye +
-        octave::feval ("tril", arg_list)(0).sparse_matrix_value ();
+      retval(0) = speye + Ftril (ovl (sm, -1))(0).sparse_matrix_value ();
+      retval(1) = Ftriu (ovl (sm))(0).sparse_matrix_value ();
     }
   else
     {
       SparseComplexMatrix sm = args(0).sparse_complex_matrix_value ();
+      SparseMatrix speye (DiagMatrix (sm.cols (), sm.cols (), 1.0));
+
       ilu_0 <SparseComplexMatrix, Complex> (sm, milu);
 
-      arg_list.append (sm);
-      retval(1) = octave::feval ("triu", arg_list)(0).sparse_complex_matrix_value ();
-      SparseComplexMatrix eye =
-        octave::feval ("speye", ovl (sm.cols ()))(0).sparse_complex_matrix_value ();
-      arg_list.append (-1);
-      retval(0) = eye +
-        octave::feval ("tril", arg_list)(0).sparse_complex_matrix_value ();
+      retval(0) = speye +
+        Ftril (ovl (sm, -1))(0).sparse_complex_matrix_value ();
+      retval(1) = Ftriu (ovl (sm))(0).sparse_complex_matrix_value ();
     }
 
   return retval;
@@ -468,72 +462,49 @@
 
 DEFUN (__iluc__, args, ,
        doc: /* -*- texinfo -*-
-@deftypefn  {} {[@var{L}, @var{U}] =} __iluc__ (@var{A})
-@deftypefnx {} {[@var{L}, @var{U}] =} __iluc__ (@var{A}, @var{droptol})
-@deftypefnx {} {[@var{L}, @var{U}] =} __iluc__ (@var{A}, @var{droptol}, @var{milu})
+@deftypefn {} {[@var{L}, @var{U}] =} __iluc__ (@var{A}, @var{droptol}, @var{milu})
 Undocumented internal function.
 @end deftypefn */)
 {
-  int nargin = args.length ();
-
-  if (nargin < 1 || nargin > 3)
+  if (args.length () != 3)
     print_usage ();
 
-  std::string milu = "off";
-  double droptol = 0;
+  double droptol = args(1).double_value ();
+  std::string milu = args(2).string_value ();
 
-  // 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 arg_list;
   if (! args(0).is_complex_type ())
     {
-      Array<double> cols_norm, rows_norm;
-      arg_list.append (args(0).sparse_matrix_value ());
-      SparseMatrix sm_u = octave::feval ("triu", arg_list)(0).sparse_matrix_value ();
-      arg_list.append (-1);
-      SparseMatrix sm_l = octave::feval ("tril", arg_list)(0).sparse_matrix_value ();
-      arg_list(1) = "rows";
-      rows_norm = octave::feval ("norm", arg_list)(0).vector_value ();
-      arg_list(1) = "cols";
-      cols_norm = octave::feval ("norm", arg_list)(0).vector_value ();
-      arg_list.clear ();
+      SparseMatrix sm = args(0).sparse_matrix_value ();
+      SparseMatrix sm_u = Ftriu (ovl (sm))(0).sparse_matrix_value ();
+      SparseMatrix sm_l = Ftril (ovl (sm, -1))(0).sparse_matrix_value ();
       SparseMatrix U, L;
+
       ilu_crout <SparseMatrix, double> (sm_l, sm_u, L, U,
-                                        cols_norm.fortran_vec (),
-                                        rows_norm.fortran_vec (),
+                                        xcolnorms (sm).fortran_vec (),
+                                        xrownorms (sm).fortran_vec (),
                                         droptol, milu);
-      arg_list.append (octave_value (L.cols ()));
-      SparseMatrix eye = octave::feval ("speye", arg_list)(0).sparse_matrix_value ();
-      return ovl (L + eye, U);
+
+      SparseMatrix speye (DiagMatrix (L.cols (), L.cols (), 1.0));
+
+      return ovl (L + speye, U);
     }
   else
     {
-      Array<Complex> cols_norm, rows_norm;
-      arg_list.append (args(0).sparse_complex_matrix_value ());
-      SparseComplexMatrix sm_u =
-        octave::feval ("triu", arg_list)(0).sparse_complex_matrix_value ();
-      arg_list.append (-1);
-      SparseComplexMatrix sm_l =
-        octave::feval ("tril", arg_list)(0).sparse_complex_matrix_value ();
-      arg_list(1) = "rows";
-      rows_norm = octave::feval ("norm", arg_list)(0).complex_vector_value ();
-      arg_list(1) = "cols";
-      cols_norm = octave::feval ("norm", arg_list)(0).complex_vector_value ();
-      arg_list.clear ();
+      SparseComplexMatrix sm = args(0).sparse_complex_matrix_value ();
+      SparseComplexMatrix sm_u = Ftriu (ovl (sm))(0).sparse_complex_matrix_value ();
+      SparseComplexMatrix sm_l = Ftril (ovl (sm, -1))(0).sparse_complex_matrix_value ();
       SparseComplexMatrix U, L;
-      ilu_crout <SparseComplexMatrix, Complex>
-                (sm_l, sm_u, L, U, cols_norm.fortran_vec (),
-                 rows_norm.fortran_vec (), Complex (droptol), milu);
+      Array<Complex> cols_norm = xcolnorms (sm);
+      Array<Complex> rows_norm = xrownorms (sm);
 
-      arg_list.append (octave_value (L.cols ()));
-      SparseComplexMatrix eye =
-        octave::feval ("speye", arg_list)(0).sparse_complex_matrix_value ();
-      return ovl (L + eye, U);
+      ilu_crout <SparseComplexMatrix, Complex> (sm_l, sm_u, L, U,
+                                                cols_norm.fortran_vec (),
+                                                rows_norm.fortran_vec (),
+                                                Complex (droptol), milu);
+
+      SparseMatrix speye (DiagMatrix (L.cols (), L.cols (), 1.0));
+
+      return ovl (L + speye, U);
     }
 }
 
@@ -926,138 +897,102 @@
 
 DEFUN (__ilutp__, args, nargout,
        doc: /* -*- texinfo -*-
-@deftypefn  {} {[@var{L}, @var{U}] =} __ilutp__ (@var{A})
-@deftypefnx {} {[@var{L}, @var{U}] =} __ilutp__ (@var{A}, @var{droptol})
-@deftypefnx {} {[@var{L}, @var{U}] =} __ilutp__ (@var{A}, @var{droptol}, @var{thresh})
-@deftypefnx {} {[@var{L}, @var{U}] =} __ilutp__ (@var{A}, @var{droptol}, @var{thresh}, @var{milu})
-@deftypefnx {} {[@var{L}, @var{U}] =} __ilutp__ (@var{A}, @var{droptol}, @var{thresh}, @var{milu}, @var{udiag})
-@deftypefnx {} {[@var{L}, @var{U}, @var{P}] =} __ilutp__ (@var{A}, @dots{})
+@deftypefn  {} {[@var{L}, @var{U}] =} __ilutp__ (@var{A}, @var{droptol}, @var{thresh}, @var{milu}, @var{udiag})
+@deftypefnx {} {[@var{L}, @var{U}, @var{P}] =} __ilutp__ (@dots{})
 Undocumented internal function.
 @end deftypefn */)
 {
-  int nargin = args.length ();
-
-  if (nargin < 1 || nargin > 5)
+  if (args.length () != 5)
     print_usage ();
 
   octave_value_list retval;
-  std::string milu = "";
-  double droptol = 0;
-  double thresh = 1;
-  double udiag = 0;
-
-  // 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 ();
+  double droptol = args(1).double_value ();
+  double thresh = args(2).double_value ();
+  std::string milu = args(3).string_value ();
+  double udiag = args(4).double_value ();
 
   octave_value_list arg_list;
   octave_idx_type nnz_u, nnz_l;
   if (! args(0).is_complex_type ())
     {
+      SparseMatrix sm = args(0).sparse_matrix_value ();
+      SparseMatrix U, L;
+      nnz_u = (Ftriu (ovl (sm))(0).sparse_matrix_value ()).nnz ();
+      nnz_l = (Ftril (ovl (sm, -1))(0).sparse_matrix_value ()).nnz ();
       Array <double> rc_norm;
-      SparseMatrix sm = args(0).sparse_matrix_value ();
-      arg_list.append (sm);
-      nnz_u = (octave::feval ("triu", arg_list)(0).sparse_matrix_value ()).nnz ();
-      arg_list.append (-1);
-      nnz_l = (octave::feval ("tril", arg_list)(0).sparse_matrix_value ()).nnz ();
       if (milu == "row")
-        arg_list (1) = "rows";
+        rc_norm = xrownorms (sm);
       else
-        arg_list (1) = "cols";
-      rc_norm = octave::feval ("norm", arg_list)(0).vector_value ();
-      arg_list.clear ();
+        rc_norm = xcolnorms (sm);
       Array <octave_idx_type> perm (dim_vector (sm.cols (), 1));
-      SparseMatrix U, L;
+
       ilu_tp <SparseMatrix, double> (sm, L, U, nnz_u, nnz_l,
                                      rc_norm.fortran_vec (),
                                      perm, droptol, thresh, milu, udiag);
-      arg_list.append (octave_value (L.cols ()));
-      SparseMatrix eye =
-        octave::feval ("speye", arg_list)(0).sparse_matrix_value ();
+
+      SparseMatrix speye (DiagMatrix (L.cols (), L.cols (), 1.0));
       if (milu == "row")
         {
+          retval(0) = L + speye;
           if (nargout == 3)
             {
-              retval(2) = eye.index (idx_vector::colon, perm);
               retval(1) = U.index (idx_vector::colon, perm);
+              retval(2) = speye.index (idx_vector::colon, perm);
             }
-          else if (nargout == 2)
+          else
             retval(1) = U;
-          retval(0) = L + eye;
         }
       else
         {
+          retval(1) = U;
           if (nargout == 3)
             {
-              retval(2) = eye.index (perm, idx_vector::colon);
-              retval(1) = U;
-              retval(0) = L.index (perm, idx_vector::colon) + eye;
+              retval(0) = L.index (perm, idx_vector::colon) + speye;
+              retval(2) = speye.index (perm, idx_vector::colon);
             }
           else
-            {
-              retval(1) = U;
-              retval(0) = L + eye.index (perm, idx_vector::colon);
-            }
+            retval(0) = L + speye.index (perm, idx_vector::colon);
         }
     }
   else
     {
+      SparseComplexMatrix sm = args(0).sparse_complex_matrix_value ();
+      SparseComplexMatrix U, L;
+      nnz_u = (Ftriu (ovl (sm))(0).sparse_complex_matrix_value ()).nnz ();
+      nnz_l = (Ftril (ovl (sm, -1))(0).sparse_complex_matrix_value ()).nnz ();
       Array <Complex> rc_norm;
-      SparseComplexMatrix sm = args(0).sparse_complex_matrix_value ();
-      arg_list.append (sm);
-      nnz_u =
-        octave::feval ("triu", arg_list)(0).sparse_complex_matrix_value ().nnz ();
-      arg_list.append (-1);
-      nnz_l =
-        octave::feval ("tril", arg_list)(0).sparse_complex_matrix_value ().nnz ();
       if (milu == "row")
-        arg_list(1) = "rows";
+        rc_norm = xrownorms (sm);
       else
-        arg_list(1) = "cols";
-      rc_norm = octave::feval ("norm", arg_list)(0).complex_vector_value ();
-      arg_list.clear ();
+        rc_norm = xcolnorms (sm);
       Array <octave_idx_type> perm (dim_vector (sm.cols (), 1));
-      SparseComplexMatrix U, L;
+
       ilu_tp <SparseComplexMatrix, Complex>
               (sm, L, U, nnz_u, nnz_l, rc_norm.fortran_vec (), perm,
                Complex (droptol), Complex (thresh), milu, udiag);
 
-      arg_list.append (octave_value (L.cols ()));
-      SparseComplexMatrix eye =
-        octave::feval ("speye", arg_list)(0).sparse_complex_matrix_value ();
+      SparseMatrix speye (DiagMatrix (L.cols (), L.cols (), 1.0));
       if (milu == "row")
         {
+          retval(0) = L + speye;
           if (nargout == 3)
             {
-              retval(2) = eye.index (idx_vector::colon, perm);
               retval(1) = U.index (idx_vector::colon, perm);
+              retval(2) = speye.index (idx_vector::colon, perm);
             }
           else if (nargout == 2)
             retval(1) = U;
-          retval(0) = L + eye;
         }
       else
         {
+          retval(1) = U;
           if (nargout == 3)
             {
-              retval(2) = eye.index (perm, idx_vector::colon);
-              retval(1) = U;
-              retval(0) = L.index (perm, idx_vector::colon) + eye;
+              retval(0) = L.index (perm, idx_vector::colon) + speye;
+              retval(2) = speye.index (perm, idx_vector::colon);
             }
           else
-            {
-              retval(1) = U;
-              retval(0) = L + eye.index (perm, idx_vector::colon);
-            }
+            retval(0) = L + speye.index (perm, idx_vector::colon);
         }
     }
 
--- a/scripts/sparse/ilu.m	Thu Jan 26 23:46:47 2017 -0500
+++ b/scripts/sparse/ilu.m	Fri Jan 27 16:39:15 2017 +0100
@@ -1,5 +1,5 @@
-## Copyright (C) 2014-2016 Eduardo Ramos Fernández <eduradical951@gmail.com>
-## Copyright (C) 2013-2016 Kai T. Ohlhus <k.ohlhus@gmail.com>
+## Copyright (C) 2014-2017 Eduardo Ramos Fernández <eduradical951@gmail.com>
+## Copyright (C) 2013-2017 Kai T. Ohlhus <k.ohlhus@gmail.com>
 ##
 ## This file is part of Octave.
 ##
@@ -242,12 +242,12 @@
           P = speye (length (A));
         endif
     case "ilutp"
-        if (nargout == 2)
-          [L, U]  = __ilutp__ (A, opts.droptol, opts.thresh,
-                                  opts.milu, opts.udiag);
-        elseif (nargout == 3)
-          [L, U, P]  = __ilutp__ (A, opts.droptol, opts.thresh,
-                                     opts.milu, opts.udiag);
+        if (nargout == 3)
+          [L, U, P] = __ilutp__ (A, opts.droptol, opts.thresh,
+                                    opts.milu, opts.udiag);
+        else
+          [L, U] = __ilutp__ (A, opts.droptol, opts.thresh,
+                                 opts.milu, opts.udiag);
         endif
   endswitch
 
@@ -459,25 +459,30 @@
 
 ## Specific tests for ilutp
 
-%!shared a1, a2
-%! a1 = sparse ([0 0 4 3 1; 5 1 2.3 2 4.5; 0 0 0 2 1;0 0 8 0 2.2; 0 0 9 9 1 ]);
-%! a2 = sparse ([3 1 0 0 4; 3 1 0 0 -2;0 0 8 0 0; 0 4 0 4 -4.5; 0 -1 0 0 1]);
+%!shared A
+%! A = sparse ([0 0 4 3 1; 5 1 2.3 2 4.5; 0 0 0 2 1;0 0 8 0 2.2; 0 0 9 9 1 ]);
+%! 
 %!test
 %! opts.udiag = 1;
 %! opts.type = "ilutp";
 %! opts.droptol = 0.2;
-%! [L, U, P] = ilu (a1, opts);
+%! [L, U, P] = ilu (A, opts);
 %! assert (norm (U, "fro"), 17.4577, 1e-4);
 %! assert (norm (L, "fro"), 2.4192, 1e-4);
-%! opts.udiag = 0;
-%! #fail ("ilu (a1, opts)");
 %!
 %!test
 %! opts.type = "ilutp";
+%! opts.udiag = 0;
+%! opts.droptol = 0.2;
+%! fail ("ilu (A, opts)", "ilu: encountered a pivot equal to 0");
+
+%!xtest
+%! A = sparse ([3 1 0 0 4; 3 1 0 0 -2;0 0 8 0 0; 0 4 0 4 -4.5; 0 -1 0 0 1]);
+%! opts.type = "ilutp";
 %! opts.droptol = 0;
 %! opts.thresh = 0;
 %! opts.milu = "row";
-%! #fail ("ilu (a2, opts)");
+%! [L, U, P] = ilu (A, opts);  % Matlab R2016b passes, no pivoting necessary
 
 ## Tests for input validation
 %!shared A_tiny