--- a/scripts/sparse/ichol.m	Mon Aug 18 12:32:16 2014 +0100
+++ b/scripts/sparse/ichol.m	Tue Aug 26 15:27:21 2014 -0700
@@ -1,335 +1,273 @@
+## Copyright (C) 2014 Eduardo Ramos Fernández <eduradical951@gmail.com>
 ## Copyright (C) 2013 Kai T. Ohlhus <k.ohlhus@gmail.com>
-## Copyright (C) 2014 Eduardo Ramos Fernández <eduradical951@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/>.
+##
+## 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/>.
 
 ## -*- texinfo -*-
-## @deftypefn  {Function File} ichol (@var{A}, @var{opts})
+## @deftypefn  {Function File} {@var{L} =} ichol (@var{A})
 ## @deftypefnx {Function File} {@var{L} =} ichol (@var{A}, @var{opts})
 ##
-## @code{@var{L} = ichol (@var{A})} performs the incomplete Cholesky
-## factorization of A with zero-fill.
+## Compute the incomplete Cholesky factorization of the sparse square matrix
+## @var{A} with zero-fill.
 ##
-## @code{@var{L} = ichol (@var{A}, @var{opts})} performs the incomplete Cholesky
-## factorization of A with options specified by opts.
-##
-## By default, ichol references the lower triangle of A and produces lower
-## triangular factors.
+## By default, ichol references the lower triangle of @var{A} and produces
+## lower triangular factors.
 ##
 ## The factor given by this routine may be useful as a preconditioner for a
 ## system of linear equations being solved by iterative methods such as
-## PCG (Preconditioned conjugate gradient).
-##
-## ichol works only for sparse square matrices.
+## PCG (Preconditioned Conjugate Gradient).
 ##
-## The fields of @var{opts} must be named exactly as shown below. You can
-## include any number of these fields in the structure and define them in any
-## order. Any additional fields are ignored. Names and specifiers are case
-## sensitive.
+## The factorization may be modified by passing options in a structure
+## @var{opts}.  The option name is a field in the structure and the setting
+## is the value of field.  Names and specifiers are case sensitive.
 ##
 ## @table @asis
 ## @item type
 ## Type of factorization.
-## String indicating which flavor of incomplete Cholesky to perform. Valid
-## values of this field are @samp{nofill} and @samp{ict}. The
-## @samp{nofill} variant performs incomplete Cholesky with zero-fill [IC(0)].
-## The @samp{ict} variant performs incomplete Cholesky with threshold dropping
-## [ICT]. The default value is @samp{nofill}.
+## String indicating which flavor of incomplete Cholesky to perform.  Valid
+## values of this field are @qcode{"nofill"} and @qcode{"ict"}.  The
+## @qcode{"nofill"} variant performs incomplete Cholesky with zero-fill
+## [IC(0)].  The @qcode{"ict"} variant performs incomplete Cholesky with
+## threshold dropping [ICT].  The default value is @qcode{"nofill"}.
 ##
 ## @item droptol
-## Drop tolerance when type is @samp{ict}.
-## Nonnegative scalar used as a drop tolerance when performing ICT. Elements
+## Drop tolerance when type is @qcode{"ict"}.
+## Non-negative scalar used as a drop tolerance when performing ICT@.  Elements
 ## which are smaller in magnitude than a local drop tolerance are dropped from
 ## the resulting factor except for the diagonal element which is never dropped.
 ## The local drop tolerance at step j of the factorization is
-## @code{norm (@var{A}(j:end, j), 1) * droptol}. @samp{droptol} is ignored if
-## @samp{type} is @samp{nofill}. The default value is 0.
+## @code{norm (@var{A}(j:end, j), 1) * droptol}.  @code{droptol} is ignored if
+## @code{type} is @qcode{"nofill"}.  The default value is 0.
 ##
 ## @item michol
-## Indicates whether to perform modified incomplete Cholesky.
-## Indicates whether or not modified incomplete Cholesky [MIC] is performed.
-## The field may be @samp{on} or @samp{off}. When performing MIC, the diagonal
-## is compensated for dropped elements to enforce the relationship
+## Indicate whether modified incomplete Cholesky [MIC] is performed.
+## The field may be @qcode{"on"} or @qcode{"off"}.  When performing MIC, the
+## diagonal is compensated for dropped elements to enforce the relationship
 ## @code{@var{A} * @var{e} = @var{L} * @var{L}' * @var{e}} where
-## @code{@var{e} = ones (size (@var{A}, 2), 1))}. The default value is
-## @samp{off}.
+## @code{@var{e} = ones (columns (@var{A}), 1)}.  The default value is
+## @qcode{"off"}.
 ##
 ## @item diagcomp
 ## Perform compensated incomplete Cholesky with the specified coefficient.
-## Real nonnegative scalar used as a global diagonal shift @code{@var{alpha}}
-## in forming the incomplete Cholesky factor. That is, instead of performing
-## incomplete Cholesky on @code{@var{A}}, the factorization of
-## @code{@var{A} + @var{alpha} * diag (diag (@var{A}))} is formed. The default
-## value is 0.
+## The coefficient is a real non-negative scalar used as a global diagonal
+## shift @code{@var{alpha}} in forming the incomplete Cholesky factor.  That
+## is, instead of performing incomplete Cholesky on @code{@var{A}}, the
+## factorization of @code{@var{A} + @var{alpha} * diag (diag (@var{A}))} is
+## formed.  The default value is 0.
 ##
 ## @item shape
-## Determines which triangle is referenced and returned.
-## Valid values are @samp{upper} and @samp{lower}. If @samp{upper} is specified,
-## only the upper triangle of @code{@var{A}} is referenced and @code{@var{R}}
-## is constructed such that @code{@var{A}} is approximated by
-## @code{@var{R}' * @var{R}}. If @samp{lower} is specified, only the lower
-## triangle of @code{@var{A}} is referenced and @code{@var{L}} is constructed
-## such that @code{@var{A}} is approximated by @code{@var{L} * @var{L}'}. The
-## default value is @samp{lower}.
+## Determine which triangle is referenced and returned.
+## Valid values are @qcode{"upper"} and @qcode{"lower"}.  If @qcode{"upper"}
+## is specified, only the upper triangle of @code{@var{A}} is referenced and
+## @code{@var{R}} is constructed such that @code{@var{A}} is approximated by
+## @code{@var{R}' * @var{R}}.  If @qcode{"lower"} is specified, only the
+## lower triangle of @code{@var{A}} is referenced and @code{@var{L}} is
+## constructed such that @code{@var{A}} is approximated by @code{@var{L} *
+## @var{L}'}.  The default value is @qcode{"lower"}.
 ## @end table
 ##
-## EXAMPLES
+## Examples
 ##
 ## The following problem demonstrates how to factorize a sample symmetric
 ## positive definite matrix with the full Cholesky decomposition and with the
 ## incomplete one.
 ##
 ## @example
+## @group
 ## A = [ 0.37, -0.05,  -0.05,  -0.07;
 ##      -0.05,  0.116,  0.0,   -0.05;
 ##      -0.05,  0.0,    0.116, -0.05;
 ##      -0.07, -0.05,  -0.05,   0.202];
-## A = sparse(A);
-## nnz(tril (A))
+## A = sparse (A);
+## nnz (tril (A))
 ## ans =  9
-## L = chol(A, "lower");
+## L = chol (A, "lower");
 ## nnz (L)
 ## ans =  10
 ## norm (A - L * L', "fro") / norm (A, "fro")
 ## ans =  1.1993e-16
-## opts.type = 'nofill';
-## L = ichol(A,opts);
+## opts.type = "nofill";
+## L = ichol (A, opts);
 ## nnz (L)
 ## ans =  9
 ## norm (A - L * L', "fro") / norm (A, "fro")
 ## ans =  0.019736
+## @end group
 ## @end example
 ##
-## Another example for decomposition is finite difference matrix to solve a
-## boundary value problem on the unit square.
+## Another example for decomposition is a finite difference matrix used to
+## solve a boundary value problem on the unit square.
 ##
 ## @example
+## @group
 ## nx = 400; ny = 200;
 ## hx = 1 / (nx + 1); hy = 1 / (ny + 1);
-## Dxx = spdiags ([ones(nx, 1), -2 * ones(nx, 1), ones(nx, 1)], [-1 0 1 ], nx, nx) / (hx ^ 2);
-## Dyy = spdiags ([ones(ny, 1), -2 * ones(ny, 1), ones(ny, 1)], [-1 0 1 ], ny, ny) / (hy ^ 2);
+## Dxx = spdiags ([ones(nx, 1), -2*ones(nx, 1), ones(nx, 1)],
+##                [-1 0 1 ], nx, nx) / (hx ^ 2);
+## Dyy = spdiags ([ones(ny, 1), -2*ones(ny, 1), ones(ny, 1)],
+##                [-1 0 1 ], ny, ny) / (hy ^ 2);
 ## A = -kron (Dxx, speye (ny)) - kron (speye (nx), Dyy);
 ## nnz (tril (A))
 ## ans =  239400
-## opts.type = 'nofill';
+## opts.type = "nofill";
 ## L = ichol (A, opts);
 ## nnz (tril (A))
 ## ans =  239400
 ## norm (A - L * L', "fro") / norm (A, "fro")
 ## ans =  0.062327
+## @end group
 ## @end example
 ##
-## References for the implemented algorithms:
+## References for implemented algorithms:
 ##
-## [1] Saad, Yousef. "Preconditioning Techniques." Iterative Methods for Sparse Linear
-## Systems. PWS Publishing Company, 1996.
+## [1] @nospell{Y. Saad}. "Preconditioning Techniques." @cite{Iterative
+## Methods for Sparse Linear Systems}, PWS Publishing Company, 1996.
 ##
-## [2] Jones, Mark T. and Plassmann, Paul E.: An Improved Incomplete Cholesky
-## Factorization (1992).
+## [2] @nospell{M. Jones, P. Plassmann}: @cite{An Improved Incomplete
+## Cholesky Factorization}, 1992.
+## @seealso{chol, ilu, pcg}
 ## @end deftypefn
 
-function [L] = ichol (A, opts)
+function L = ichol (A, opts = struct ())
 
-  if ((nargin > 2) || (nargin < 1) || (nargout > 1))
+  if (nargin < 1 || nargin > 2 || nargout > 1)
     print_usage ();
   endif
 
-  % Check input matrix
-  if (~issparse(A) || ~issquare (A))
-    error ("ichol: Input A must be a non-empty sparse square matrix");
+  if (! (issparse (A) && issquare (A)))
+    error ("ichol: A must be a sparse square matrix");
   endif
 
-  % If A is empty and sparse then return empty L
-  % Compatibility with Matlab
+  if (! isstruct (opts))
+    error ("ichol: OPTS must be a structure.");
+  endif
+
+  ## If A is empty then return empty L for Matlab compatibility
   if (isempty (A)) 
     L = A;
     return;
   endif
 
-  % Check input structure, otherwise set default values
-  if (nargin == 2)
-    if (~isstruct (opts))
-      error ("ichol: Input \"opts\" must be a valid structure.");
-    endif
-  else
-    opts = struct ();
-  endif
-
-  if (~isfield (opts, "type"))
-    opts.type = "nofill"; % set default
+  ## Parse input options
+  if (! isfield (opts, "type"))
+    opts.type = "nofill";  # set default
   else
     type = tolower (getfield (opts, "type"));
-    if ((strcmp (type, "nofill") == 0)
-        && (strcmp (type, "ict") == 0))
-      error ("ichol: Invalid field \"type\" in input structure.");
-    else
-      opts.type = type;
+    if (! strcmp (type, "nofill") && ! strcmp (type, "ict"))
+      error ('ichol: TYPE must be "nofill" or "ict"');
     endif
+    opts.type = type;
   endif
 
-  if (~isfield (opts, "droptol"))
-    opts.droptol = 0; % set default
+  if (! isfield (opts, "droptol"))
+    opts.droptol = 0;      # set default
   else
-    if (~isscalar (opts.droptol) || (opts.droptol < 0))
-      error ("ichol: Invalid field \"droptol\" in input structure.");
+    if (! (isreal (opts.droptol) && isscalar (opts.droptol)
+           && opts.droptol >= 0))
+      error ("ichol: DROPTOL must be a non-negative real scalar");
     endif
   endif
 
   michol = "";
-  if (~isfield (opts, "michol"))
-    opts.michol = "off"; % set default
+  if (! isfield (opts, "michol"))
+    opts.michol = "off";   # set default
   else
     michol = tolower (getfield (opts, "michol"));
-    if ((strcmp (michol, "off") == 0) 
-        && (strcmp (michol, "on") == 0))
-      error ("ichol: Invalid field \"michol\" in input structure.");
-    else
-      opts.michol = michol;
+    if (! strcmp (michol, "off") && ! strcmp (michol, "on"))
+      error ('ichol: MICHOL must be "on" or "off"');
+    endif
+    opts.michol = michol;
+  endif
+
+  if (! isfield (opts, "diagcomp"))
+    opts.diagcomp = 0;     # set default
+  else
+    if (! (isreal (opts.diagcomp) && isscalar (opts.diagcomp)
+           && opts.diagcomp >= 0))
+      error ("ichol: DIAGCOMP must be a non-negative real scalar");
     endif
   endif
 
-  if (~isfield (opts, "diagcomp"))
-    opts.diagcomp = 0; % set default
+  if (! isfield (opts, "shape"))
+    opts.shape = "lower";  # set default
   else
-    if (~isscalar (opts.diagcomp) || (opts.diagcomp < 0))
-      error ("ichol: Invalid field \"diagcomp\" in input structure.");
+    shape = tolower (getfield (opts, "shape"));
+    if (! strcmp (shape, "lower") && ! strcmp (shape, "upper"))
+      error ('ichol: SHAPE must be "lower" or "upper"');
     endif
+    opts.shape = shape;
   endif
 
-  if (~isfield (opts, "shape"))
-    opts.shape = "lower"; % set default
-  else
-    shape = tolower (getfield (opts, "shape"));
-    if ((strcmp (shape, "lower") == 0) 
-        && (strcmp (shape, "upper") == 0))
-      error ("ichol: Invalid field \"shape\" in input structure.");
-    else
-      opts.shape = shape;
-    endif
-  endif
-
-  % Prepare input for specialized ICHOL
+  ## Prepare input for specialized ICHOL
   A_in = [];
   if (opts.diagcomp > 0)
     A += opts.diagcomp * diag (diag (A));
   endif
-  if (strcmp (opts.shape, "upper") == 1)
+  if (strcmp (opts.shape, "upper"))
     A_in = triu (A);
     A_in = A_in';
   else
     A_in = tril (A);
   endif
 
-  % Delegate to specialized ICHOL
+  ## Delegate to specialized ICHOL
   switch (opts.type)
     case "nofill"
       L  = __ichol0__ (A_in, opts.michol);
     case "ict"
       L = __icholt__ (A_in, opts.droptol, opts.michol);
-    otherwise
-      printf ("The input structure is invalid.\n");
   endswitch
 
-  if (strcmp (opts.shape, "upper") == 1)
+  if (strcmp (opts.shape, "upper"))
     L = L';
   endif
   
+endfunction
 
-endfunction
 
 %!shared A1, A2, A3, A4, A5, A6, A7
 %! A1 = [ 0.37, -0.05,  -0.05,  -0.07;
-%!      -0.05,  0.116,  0.0,   -0.05;
-%!      -0.05,  0.0,    0.116, -0.05;
-%!      -0.07, -0.05,  -0.05,   0.202];
-%! A1 = sparse(A1);
-%! A2 = gallery ('poisson', 30);
-%! A3 = gallery ('tridiag', 50);
+%!       -0.05,  0.116,  0.0,   -0.05;
+%!       -0.05,  0.0,    0.116, -0.05;
+%!       -0.07, -0.05,  -0.05,   0.202];
+%! A1 = sparse (A1);
+%! A2 = gallery ("poisson", 30);
+%! A3 = gallery ("tridiag", 50);
 %! nx = 400; ny = 200;
 %! hx = 1 / (nx + 1); hy = 1 / (ny + 1);
-%! Dxx = spdiags ([ones(nx, 1), -2 * ones(nx, 1), ones(nx, 1)], [-1 0 1 ], nx, nx) / (hx ^ 2);
-%! Dyy = spdiags ([ones(ny, 1), -2 * ones(ny, 1), ones(ny, 1)], [-1 0 1 ], ny, ny) / (hy ^ 2);
+%! Dxx = spdiags ([ones(nx, 1), -2*ones(nx, 1), ones(nx, 1)],
+%!                [-1 0 1 ], nx, nx) / (hx ^ 2);
+%! Dyy = spdiags ([ones(ny, 1), -2*ones(ny, 1), ones(ny, 1)],
+%!                [-1 0 1 ], ny, ny) / (hy ^ 2);
 %! A4 = -kron (Dxx, speye (ny)) - kron (speye (nx), Dyy);
-%! A5 = [ 0.37, -0.05,          -0.05,  -0.07;
-%!        -0.05,  0.116,          0.0,   -0.05 + 0.05i;
-%!        -0.05,  0.0,            0.116, -0.05;
-%!        -0.07, -0.05 - 0.05i,  -0.05,   0.202];
-%! A5 = sparse(A5);
-%! A6 = [ 0.37,    -0.05 - i, -0.05,  -0.07;
-%!        -0.05 + i, 0.116,     0.0,   -0.05;
-%!        -0.05,     0.0,       0.116, -0.05;
-%!        -0.07,    -0.05,     -0.05,   0.202];
-%! A6 = sparse(A6);
+%! A5 = [ 0.37, -0.05,         -0.05,  -0.07;
+%!       -0.05,  0.116,         0.0,   -0.05 + 0.05i;
+%!       -0.05,  0.0,           0.116, -0.05;
+%!       -0.07, -0.05 - 0.05i, -0.05,   0.202];
+%! A5 = sparse (A5);
+%! A6 = [ 0.37,     -0.05 - i, -0.05,  -0.07;
+%!       -0.05 + i,  0.116,     0.0,   -0.05;
+%!       -0.05,      0.0,       0.116, -0.05;
+%!       -0.07,     -0.05,     -0.05,   0.202];
+%! A6 = sparse (A6);
 %! A7 = A5;
 %! A7(1) = 2i;
-%!
 
-%!# Input validation tests
-
-%!test
-%!error ichol ([]);
-%!error ichol (0);
-%!error ichol (-0);
-%!error ichol (1);
-%!error ichol (-1);
-%!error ichol (i);
-%!error ichol (-i);
-%!error ichol (1 + 1i);
-%!error ichol (1 - 1i);
-%!error ichol (sparse (0));
-%!error ichol (sparse (-0));
-%!error ichol (sparse (-1));
-%!error ichol (sparse (-1));
-%!test
-%! opts.milu = 'foo';
-%!error L = ichol (A1, opts);
-%! opts.milu = 1;
-%!error L = ichol (A1, opts);
-%! opts.milu = [];
-%!error L = ichol (A1, opts);
-%!test
-%! opts.droptol = -1;
-%!error L = ichol (A1, opts);
-%! opts.droptol = 0.5i;
-%!error L = ichol (A1, opts);
-%! opts.droptol = [];
-%!error L = ichol (A1, opts);
-%!test
-%! opts.type = 'foo';
-%!error L = ichol (A1, opts);
-%! opts.type = 1;
-%!error L = ichol (A1, opts);
-%! opts.type = [];
-%!error L = ichol (A1, opts);
-%!test
-%! opts.shape = 'foo';
-%!error L = ichol (A1, opts);
-%! opts.shape = 1;
-%!error L = ichol (A1, opts);
-%! opts.shape = [];
-%!error L = ichol (A1, opts);
-%!test
-%! opts.diagcomp = -1;
-%!error L = ichol (A1, opts);
-%! opts.diagcomp = 0.5i;
-%!error L = ichol (A1, opts);
-%! opts.diagcomp = [];
-%!error L = ichol (A1, opts);
-
-%!# ICHOL0 tests
+## ICHOL0 tests
 
 %!test
 %! opts.type = "nofill";
@@ -401,13 +339,14 @@
 %! opts.michol = "on";
 %! L = ichol (A5, opts);
 %! assert (norm (A5 - L*L', "fro") / norm (A5, "fro"), 0.0276, 1e-4);
-%!test
-%% Negative pivot 
-%!error ichol (A6);
-%% Complex entry in the diagonal
-%!error ichol (A7);
 
-%%!ICHOLT tests
+## Negative pivot 
+%!error <negative pivot> ichol (A6)
+%!error ichol (A6)
+## Complex entry in the diagonal
+%!error <non-real pivot> ichol (A7)
+
+## ICHOLT tests
  
 %%!test
 %! opts.type = "ict";
@@ -475,9 +414,47 @@
 %! opts.michol = "on";
 %! L = ichol (A5, opts);
 %! assert (norm (A5 - L*L', "fro") / norm (A5, "fro"), 0.3231, 1e-4);
+
+%% Input validation tests
+
+%!error <A must be a sparse square matrix> ichol ([])
+%!error <A must be a sparse square matrix> ichol (0)
+%!error <pivot equal to 0> ichol (sparse (0))
+%!error <pivot equal to 0> ichol (sparse (-0))
+%!error <negative pivot> ichol (sparse (-1))
 %!test
-%% Negative pivot 
-%! opts.type = "ict";
-%!error ichol (A6, setup);
-%% Complex entry in the diagonal
-%!error ichol (A7, setup);
+%! opts.type = "foo";
+%! fail ("ichol (A1, opts)", 'TYPE must be "nofill"');
+%! opts.type = 1;
+%! fail ("ichol (A1, opts)", 'TYPE must be "nofill"');
+%! opts.type = [];
+%! fail ("ichol (A1, opts)", 'TYPE must be "nofill"');
+%!test
+%! opts.droptol = -1;
+%! fail ("ichol (A1, opts)", "DROPTOL must be a non-negative real scalar");
+%! opts.droptol = 0.5i;
+%! fail ("ichol (A1, opts)", "DROPTOL must be a non-negative real scalar");
+%! opts.droptol = [];
+%! fail ("ichol (A1, opts)", "DROPTOL must be a non-negative real scalar");
+%!test
+%! opts.michol = "foo";
+%! fail ("ichol (A1, opts)", 'MICHOL must be "on"');
+%! opts.michol = 1;
+%! fail ("ichol (A1, opts)", 'MICHOL must be "on"');
+%! opts.michol = [];
+%! fail ("ichol (A1, opts)", 'MICHOL must be "on"');
+%!test
+%! opts.diagcomp = -1;
+%! fail ("ichol (A1, opts)", "DIAGCOMP must be a non-negative real scalar");
+%! opts.diagcomp = 0.5i;
+%! fail ("ichol (A1, opts)", "DIAGCOMP must be a non-negative real scalar");
+%! opts.diagcomp = [];
+%! fail ("ichol (A1, opts)", "DIAGCOMP must be a non-negative real scalar");
+%!test
+%! opts.shape = "foo";
+%! fail ("ichol (A1, opts)", 'SHAPE must be "lower"');
+%! opts.shape = 1;
+%! fail ("ichol (A1, opts)", 'SHAPE must be "lower"');
+%! opts.shape = [];
+%! fail ("ichol (A1, opts)", 'SHAPE must be "lower"');
+