--- a/NEWS	Thu Jun 30 18:20:37 2016 +1000
+++ b/NEWS	Mon Jul 25 14:50:29 2016 +0100
@@ -1,6 +1,8 @@
 Summary of important user-visible changes for version 4.2:
 ---------------------------------------------------------
 
+ ** condest now works with a normest1 compatible syntax.
+
  ** The parser has been extended to accept, but ignore, underscore
     characters in numbers.  This facilitates writing more legible code
     by using '_' as a thousands separator or to group nibbles into bytes
@@ -110,6 +112,7 @@
       hash
       im2double
       localfunctions
+      normest1
       ode45
       odeget
       odeset
@@ -130,6 +133,7 @@
       mahalanobis          | mahal in Octave-Forge statistics pkg
       md5sum               | hash
       octve_config_info    | __octave_config_info__
+      onenormest           | normest1
       sleep                | pause
       usleep               | pause
       wavread              | audioread

--- a/doc/interpreter/sparse.txi	Thu Jun 30 18:20:37 2016 +1000
+++ b/doc/interpreter/sparse.txi	Mon Jul 25 14:50:29 2016 +0100
@@ -486,7 +486,7 @@
 
 @item Linear algebra:
   @dfn{condest}, @dfn{eigs}, @dfn{matrix_type},
-  @dfn{normest}, @dfn{sprank}, @dfn{spaugment}, @dfn{svds}
+  @dfn{normest}, @dfn{normest1}, @dfn{sprank}, @dfn{spaugment}, @dfn{svds}
 
 @item Iterative techniques:
   @dfn{ichol}, @dfn{ilu}, @dfn{pcg}, @dfn{pcr}
@@ -830,7 +830,7 @@
 
 @DOCSTRING(normest)
 
-@DOCSTRING(onenormest)
+@DOCSTRING(normest1)
 
 @DOCSTRING(condest)
 

--- a/scripts/deprecated/module.mk	Thu Jun 30 18:20:37 2016 +1000
+++ b/scripts/deprecated/module.mk	Mon Jul 25 14:50:29 2016 +0100
@@ -19,6 +19,7 @@
   scripts/deprecated/nfields.m \
   scripts/deprecated/octave_config_info.m \
   scripts/deprecated/octave_tmp_file_name.m \
+  scripts/deprecated/onenormest.m \
   scripts/deprecated/playaudio.m \
   scripts/deprecated/saveaudio.m \
   scripts/deprecated/setaudio.m \

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/deprecated/onenormest.m	Mon Jul 25 14:50:29 2016 +0100
@@ -0,0 +1,307 @@
+## Copyright (C) 2007-2015 Regents of the University of California
+##
+## 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/>.
+
+## -*- texinfo -*-
+## @deftypefn  {} {[@var{est}, @var{v}, @var{w}, @var{iter}] =} onenormest (@var{A}, @var{t})
+## @deftypefnx {} {[@var{est}, @var{v}, @var{w}, @var{iter}] =} onenormest (@var{apply}, @var{apply_t}, @var{n}, @var{t})
+##
+## @code{onenormest} is deprecated and will be removed in Octave version 4.4.
+## Use @code{normest1} for the equivalent functionality.
+##
+## Apply @nospell{Higham and Tisseur's} randomized block 1-norm estimator to
+## matrix @var{A} using @var{t} test vectors.
+##
+## If @var{t} exceeds 5, then only 5 test vectors are used.
+##
+## If the matrix is not explicit, e.g., when estimating the norm of
+## @code{inv (@var{A})} given an LU@tie{}factorization, @code{onenormest}
+## applies @var{A} and its conjugate transpose through a pair of functions
+## @var{apply} and @var{apply_t}, respectively, to a dense matrix of size
+## @var{n} by @var{t}.  The implicit version requires an explicit dimension
+## @var{n}.
+##
+## Returns the norm estimate @var{est}, two vectors @var{v} and @var{w} related
+## by norm @code{(@var{w}, 1) = @var{est} * norm (@var{v}, 1)}, and the number
+## of iterations @var{iter}.  The number of iterations is limited to 10 and is
+## at least 2.
+##
+## References:
+##
+## @itemize
+## @item
+## @nospell{N.J. Higham and F. Tisseur}, @cite{A Block Algorithm
+## for Matrix 1-Norm Estimation, with an Application to 1-Norm
+## Pseudospectra}. SIMAX vol 21, no 4, pp 1185-1201.
+## @url{http://dx.doi.org/10.1137/S0895479899356080}
+##
+## @item
+## @nospell{N.J. Higham and F. Tisseur}, @cite{A Block Algorithm
+## for Matrix 1-Norm Estimation, with an Application to 1-Norm
+## Pseudospectra}. @url{http://citeseer.ist.psu.edu/223007.html}
+## @end itemize
+##
+## @seealso{condest, norm, cond}
+## @end deftypefn
+
+## Code originally licensed under:
+##
+## Copyright (c) 2007, Regents of the University of California
+## All rights reserved.
+##
+## Redistribution and use in source and binary forms, with or without
+## modification, are permitted provided that the following conditions
+## are met:
+##
+##    * Redistributions of source code must retain the above copyright
+##      notice, this list of conditions and the following disclaimer.
+##
+##    * Redistributions in binary form must reproduce the above
+##      copyright notice, this list of conditions and the following
+##      disclaimer in the documentation and/or other materials provided
+##      with the distribution.
+##
+##    * Neither the name of the University of California, Berkeley nor
+##      the names of its contributors may be used to endorse or promote
+##      products derived from this software without specific prior
+##      written permission.
+##
+## THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS''
+## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+## TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+## PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS AND
+## CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+## SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+## LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
+## USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+## ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
+## OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+## SUCH DAMAGE.
+
+## Author: Jason Riedy <ejr@cs.berkeley.edu>
+## Keywords: linear-algebra norm estimation
+## Version: 0.2
+
+function [est, v, w, iter] = onenormest (varargin)
+
+  persistent warned = false;
+  if (! warned)
+    warned = true;
+    warning ("Octave:deprecated-function",
+             "onenormest is obsolete and will be removed from a future version of Octave, please use normest1 instead");
+  endif
+
+
+  if (nargin < 1 || nargin > 4)
+    print_usage ();
+  endif
+
+  default_t = 5;
+  itmax = 10;
+
+  if (isnumeric (varargin{1}))
+    [n, nc] = size (varargin{1});
+    if (n != nc)
+      error ("onenormest: matrix must be square");
+    endif
+    apply = @(x) varargin{1} * x;
+    apply_t = @(x) varargin{1}' * x;
+    if (nargin > 1)
+      t = varargin{2};
+    else
+      t = min (n, default_t);
+    endif
+    issing = isa (varargin{1}, "single");
+  else
+    if (nargin < 3)
+      print_usage ();
+    endif
+    apply = varargin{1};
+    apply_t = varargin{2};
+    n = varargin{3};
+    if (nargin > 3)
+      t = varargin{4};
+    else
+      t = default_t;
+    endif
+    issing = isa (n, "single");
+  endif
+
+  ## Initial test vectors X.
+  X = rand (n, t);
+  X ./= ones (n,1) * sum (abs (X), 1);
+
+  ## Track if a vertex has been visited.
+  been_there = zeros (n, 1);
+
+  ## To check if the estimate has increased.
+  est_old = 0;
+
+  ## Normalized vector of signs.
+  S = zeros (n, t);
+
+  if (issing)
+    myeps = eps ("single");
+    X = single (X);
+  else
+    myeps = eps;
+  endif
+
+  for iter = 1 : itmax + 1
+    Y = feval (apply, X);
+
+    ## Find the initial estimate as the largest A*x.
+    [est, ind_best] = max (sum (abs (Y), 1));
+    if (est > est_old || iter == 2)
+      w = Y(:,ind_best);
+    endif
+    if (iter >= 2 && est < est_old)
+      ## No improvement, so stop.
+      est = est_old;
+      break;
+    endif
+
+    est_old = est;
+    S_old = S;
+    if (iter > itmax),
+      ## Gone too far.  Stop.
+      break;
+    endif
+
+    S = sign (Y);
+
+    ## Test if any of S are approximately parallel to previous S
+    ## vectors or current S vectors.  If everything is parallel,
+    ## stop.  Otherwise, replace any parallel vectors with
+    ## rand{-1,+1}.
+    partest = any (abs (S_old' * S - n) < 4*eps*n);
+    if (all (partest))
+      ## All the current vectors are parallel to old vectors.
+      ## We've hit a cycle, so stop.
+      break;
+    endif
+    if (any (partest))
+      ## Some vectors are parallel to old ones and are cycling,
+      ## but not all of them.  Replace the parallel vectors with
+      ## rand{-1,+1}.
+      numpar = sum (partest);
+      replacements = 2*(rand (n,numpar) < 0.5) - 1;
+      S(:,partest) = replacements;
+    endif
+    ## Now test for parallel vectors within S.
+    partest = any ((S' * S - eye (t)) == n);
+    if (any (partest))
+      numpar = sum (partest);
+      replacements = 2*(rand (n,numpar) < 0.5) - 1;
+      S(:,partest) = replacements;
+    endif
+
+    Z = feval (apply_t, S);
+
+    ## Now find the largest non-previously-visted index per vector.
+    h = max (abs (Z),2);
+    [mh, mhi] = max (h);
+    if (iter >= 2 && mhi == ind_best)
+      ## Hit a cycle, stop.
+      break;
+    endif
+    [h, ind] = sort (h, 'descend');
+    if (t > 1)
+      firstind = ind(1:t);
+      if (all (been_there(firstind)))
+        ## Visited all these before, so stop.
+        break;
+      endif
+      ind = ind(! been_there(ind));
+      if (length (ind) < t)
+        ## There aren't enough new vectors, so we're practically
+        ## in a cycle.  Stop.
+        break;
+      endif
+    endif
+
+    ## Visit the new indices.
+    X = zeros (n, t);
+    for zz = 1 : t
+      X(ind(zz),zz) = 1;
+    endfor
+    been_there(ind(1 : t)) = 1;
+  endfor
+
+  ## The estimate est and vector w are set in the loop above.
+  ## The vector v selects the ind_best column of A.
+  v = zeros (n, 1);
+  v(ind_best) = 1;
+
+endfunction
+
+
+%!demo
+%! N = 100;
+%! A = randn (N) + eye (N);
+%! [L,U,P] = lu (A);
+%! nm1inv = onenormest (@(x) U\(L\(P*x)), @(x) P'*(L'\(U'\x)), N, 30)
+%! norm (inv (A), 1)
+
+%!test
+%! warning ("off", "Octave:deprecated-function", "local");
+%! N = 10;
+%! A = ones (N);
+%! [nm1, v1, w1] = onenormest (A);
+%! [nminf, vinf, winf] = onenormest (A', 6);
+%! assert (nm1, N, -2*eps);
+%! assert (nminf, N, -2*eps);
+%! assert (norm (w1, 1), nm1 * norm (v1, 1), -2*eps);
+%! assert (norm (winf, 1), nminf * norm (vinf, 1), -2*eps);
+
+%!test
+%! warning ("off", "Octave:deprecated-function", "local");
+%! N = 10;
+%! A = ones (N);
+%! [nm1, v1, w1] = onenormest (@(x) A*x, @(x) A'*x, N, 3);
+%! [nminf, vinf, winf] = onenormest (@(x) A'*x, @(x) A*x, N, 3);
+%! assert (nm1, N, -2*eps);
+%! assert (nminf, N, -2*eps);
+%! assert (norm (w1, 1), nm1 * norm (v1, 1), -2*eps);
+%! assert (norm (winf, 1), nminf * norm (vinf, 1), -2*eps);
+
+%!test
+%! warning ("off", "Octave:deprecated-function", "local");
+%! N = 5;
+%! A = hilb (N);
+%! [nm1, v1, w1] = onenormest (A);
+%! [nminf, vinf, winf] = onenormest (A', 6);
+%! assert (nm1, norm (A, 1), -2*eps);
+%! assert (nminf, norm (A, inf), -2*eps);
+%! assert (norm (w1, 1), nm1 * norm (v1, 1), -2*eps);
+%! assert (norm (winf, 1), nminf * norm (vinf, 1), -2*eps);
+
+## Only likely to be within a factor of 10.
+%!test
+%! warning ("off", "Octave:deprecated-function", "local");
+%! old_state = rand ("state");
+%! restore_state = onCleanup (@() rand ("state", old_state));
+%! rand ("state", 42);  # Initialize to guarantee reproducible results
+%! N = 100;
+%! A = rand (N);
+%! [nm1, v1, w1] = onenormest (A);
+%! [nminf, vinf, winf] = onenormest (A', 6);
+%! assert (nm1, norm (A, 1), -.1);
+%! assert (nminf, norm (A, inf), -.1);
+%! assert (norm (w1, 1), nm1 * norm (v1, 1), -2*eps);
+%! assert (norm (winf, 1), nminf * norm (vinf, 1), -2*eps);

--- a/scripts/linear-algebra/condest.m	Thu Jun 30 18:20:37 2016 +1000
+++ b/scripts/linear-algebra/condest.m	Mon Jul 25 14:50:29 2016 +0100
@@ -1,4 +1,5 @@
 ## Copyright (C) 2007-2015 Regents of the University of California
+## Copyright (C) 2016 Marco Caliari
 ##
 ## This file is part of Octave.
 ##
@@ -20,8 +21,8 @@
 ## @deftypefn  {} {} condest (@var{A})
 ## @deftypefnx {} {} condest (@var{A}, @var{t})
 ## @deftypefnx {} {[@var{est}, @var{v}] =} condest (@dots{})
-## @deftypefnx {} {[@var{est}, @var{v}] =} condest (@var{A}, @var{solve}, @var{solve_t}, @var{t})
-## @deftypefnx {} {[@var{est}, @var{v}] =} condest (@var{apply}, @var{apply_t}, @var{solve}, @var{solve_t}, @var{n}, @var{t})
+## @deftypefnx {} {[@var{est}, @var{v}] =} condest (@var{A}, @var{solvefun}, @var{t}, @var{p1}, @var{p2}, @dots{})
+## @deftypefnx {} {[@var{est}, @var{v}] =} condest (@var{afun}, @var{solvefun}, @var{t}, @var{p1}, @var{p2}, @dots{})
 ##
 ## Estimate the 1-norm condition number of a matrix @var{A} using @var{t} test
 ## vectors using a randomized 1-norm estimator.
@@ -32,23 +33,39 @@
 ## number of @var{A} given an LU@tie{}factorization, @code{condest} uses the
 ## following functions:
 ##
-## @table @var
-## @item apply
-## @code{A*x} for a matrix @code{x} of size @var{n} by @var{t}.
-##
-## @item apply_t
-## @code{A'*x} for a matrix @code{x} of size @var{n} by @var{t}.
+## @itemize @minus
+## @item @var{afun} which should returns
+## @itemize @bullet
+## @item
+## the dimension @var{n} of @var{a}, if @var{flag} is "dim"
+## @item
+## true if @var{a} is a real operator, if @var{flag} is "real"
+## @item
+## the result @code{@var{a} * @var{x}}, if @var{flag} is "notransp"
+## @item
+## the result @code{@var{a}' * @var{x}}, if @var{flag} is "transp"
+## @end itemize
+## @item @var{solvefun} which should returns
+## @itemize @bullet
+## @item
+## the dimension @var{n} of @var{a}, if @var{flag} is "dim"
+## @item
+## true if @var{a} is a real operator, if @var{flag} is "real"
+## @item
+## the result @code{@var{a} \ @var{x}}, if @var{flag} is "notransp"
+## @item
+## the result @code{@var{a}' \ @var{x}}, if @var{flag} is "transp"
+## @end itemize
+## @end itemize
 ##
-## @item solve
-## @code{A \ b} for a matrix @code{b} of size @var{n} by @var{t}.
+## The parameters @var{p1}, @var{p2}, @dots{} are arguments of
+## @code{@var{afun} (@var{flag}, @var{x}, @var{p1}, @var{p2}, @dots{})}
+## and @code{@var{solvefun} (@var{flag}, @var{x}, @var{p1}, @var{p2},
+## @dots{})}.
 ##
-## @item solve_t
-## @code{A' \ b} for a matrix @code{b} of size @var{n} by @var{t}.
-## @end table
-##
-## The implicit version requires an explicit dimension @var{n}.
-##
-## @code{condest} uses a randomized algorithm to approximate the 1-norms.
+## @code{condest} uses a randomized algorithm to approximate the
+## 1-norms. Therefore, if consistent results are required, the "state" of the
+## random generator should be fixed before invoking @code{normest1}.
 ##
 ## @code{condest} returns the 1-norm condition estimate @var{est} and a vector
 ## @var{v} satisfying @code{norm (A*v, 1) == norm (A, 1) * norm
@@ -70,7 +87,7 @@
 ## Pseudospectra}. @url{http://citeseer.ist.psu.edu/223007.html}
 ## @end itemize
 ##
-## @seealso{cond, norm, onenormest}
+## @seealso{cond, norm, normest1}
 ## @end deftypefn
 
 ## Code originally licensed under:
@@ -120,10 +137,129 @@
 
   default_t = 5;
 
+  if ((nargin == 3 && is_function_handle (varargin{3}))
+      || (nargin == 4 && is_function_handle (varargin{3})
+          && isnumeric (varargin{4})))
+    ## onenormest syntax, deprecated in 4.2
+    [est, v] = condest_legacy (varargin{:});
+    return
+  elseif ((nargin >= 5) && is_function_handle (varargin{4}))
+    ## onenormest syntax, deprecated in 4.2
+    [est, v] = condest_legacy (varargin{:});
+    return
+  endif
+
+  have_A = false;
+  have_t = false;
+  have_apply_normest1 = false;
+  have_solve_normest1 = false;
+
+  if (isnumeric (varargin{1}))
+    A = varargin{1};
+    if (! issquare (A))
+      error ("condest: matrix must be square");
+    endif
+    n = rows (A);
+    have_A = true;
+    if (nargin > 1)
+      if (is_function_handle (varargin{2}))
+        solve = varargin{2};
+        have_solve_normest1 = true;
+        if (nargin > 2)
+          t = varargin{3};
+          have_t = true;
+        endif
+      else
+        t = varargin{2};
+        have_t = true;
+        real_op = isreal (A);
+      endif
+    else
+      real_op = isreal (A);
+    endif
+  else # varargin{1} is function handle
+    apply = varargin{1};
+    if (nargin > 1)
+      solve = varargin{2};
+      have_apply_normest1 = true;
+      have_solve_normest1 = true;
+      n = apply ("dim", [], varargin{4:end});
+      if (nargin > 2)
+        t = varargin{3};
+        have_t = true;
+      endif
+    else
+      error("condest: wrong number of input parameters");
+    endif
+  endif
+
+  if (! have_t)
+    t = min (n, default_t);
+  endif
+
+  if (! have_solve_normest1)
+     ## prepare solve in normest1 form
+    if (issparse (A))
+      [L, U, P, Pc] = lu (A);
+      solve = @(flag, x) solve_sparse (flag, x, n, real_op, L, U, P, Pc);
+    else
+      [L, U, P] = lu (A);
+      solve = @(flag, x) solve_not_sparse (flag, x, n, real_op, L, U, P);
+    endif
+  endif
+
+  if (have_A)
+    Anorm = norm (A, 1);
+  else
+    Anorm = normest1 (apply, t, [], varargin{4:end});
+  endif
+  [Ainv_norm, v, w] = normest1 (solve, t, [], varargin{4:end});
+
+  est = Anorm * Ainv_norm;
+  v = w / norm (w, 1);
+
+endfunction
+
+function value = solve_sparse (flag, x, n, real_op, L , U , P , Pc)
+  switch flag
+    case "dim"
+      value = n;
+    case "real"
+      value = real_op;
+    case "notransp"
+      value = P' * (L' \ (U' \ (Pc * x)));
+    case "transp"
+      value = Pc' * (U \ (L \ (P * x)));
+  endswitch
+endfunction
+
+function value = solve_not_sparse (flag, x, n, real_op, L, U, P)
+  switch flag
+    case "dim"
+      value = n;
+    case "real"
+      value = real_op;
+    case "notransp"
+      value = P' * (L' \ (U' \ x));
+    case "transp"
+      value = U \ (L \ (P * x));
+  endswitch
+endfunction
+
+function [est, v] = condest_legacy (varargin) # to be removed after 4.2
+
+  persistent warned = false;
+  if (! warned)
+    warned = true;
+    warning ("Octave:deprecated-function",
+             "condest: this syntax is deprecated, call condest (A, SOLVEFUN, T, P1, P2, ...) instead.");
+  endif
+
+  default_t = 5;
+
   have_A = false;
   have_t = false;
   have_solve = false;
-
   if (isnumeric (varargin{1}))
     A = varargin{1};
     if (! issquare (A))
@@ -182,6 +318,10 @@
     endif
   endif
 
+  ## We already warned about this usage being deprecated.  Don't
+  ## warn again about onenormest.
+  warning ("off", "Octave:deprecated-function", "local");
+
   if (have_A)
     Anorm = norm (A, 1);
   else
@@ -195,19 +335,40 @@
 
 endfunction
 
-
-%!demo
-%! N = 100;
-%! A = randn (N) + eye (N);
-%! condest (A)
-%! [L,U,P] = lu (A);
-%! condest (A, @(x) U \ (L \ (P*x)), @(x) P'*(L' \ (U'\x)))
-%! condest (@(x) A*x, @(x) A'*x, @(x) U \ (L \ (P*x)), @(x) P'*(L' \ (U'\x)), N)
-%! norm (inv (A), 1) * norm (A, 1)
-
 ## Yes, these test bounds are really loose.  There's
 ## enough randomization to trigger odd cases with hilb().
 
+%!function value = apply_fun (flag, x, A, m)
+%!  if (nargin == 3)
+%!    m = 1;
+%!  endif
+%!  switch flag
+%!    case "dim"
+%!      value = length (A);
+%!    case "real"
+%!      value = isreal (A);
+%!    case "notransp"
+%!      value = x; for i = 1:m, value = A * value;, endfor
+%!    case "transp"
+%!      value = x; for i = 1:m, value = A' * value;, endfor
+%!  endswitch
+%!endfunction
+%!function value = solve_fun (flag, x, A, m)
+%!  if (nargin == 3)
+%!    m = 1;
+%!  endif
+%!  switch flag
+%!    case "dim"
+%!      value = length (A);
+%!    case "real"
+%!      value = isreal (A);
+%!    case "notransp"
+%!      value = x; for i = 1:m, value = A \ value;, endfor;
+%!    case "transp"
+%!      value = x; for i = 1:m, value = A' \ value;, endfor;
+%!  endswitch
+%!endfunction
+
 %!test
 %! N = 6;
 %! A = hilb (N);
@@ -215,7 +376,8 @@
 %! cA_test = norm (inv (A), 1) * norm (A, 1);
 %! assert (cA, cA_test, -2^-8);
 
-%!test
+%!test # to be removed after 4.2
+%! warning ("off", "Octave:deprecated-function", "local");
 %! N = 6;
 %! A = hilb (N);
 %! solve = @(x) A\x; solve_t = @(x) A'\x;
@@ -223,7 +385,8 @@
 %! cA_test = norm (inv (A), 1) * norm (A, 1);
 %! assert (cA, cA_test, -2^-8);
 
-%!test
+%!test # to be removed after 4.2
+%! warning ("off", "Octave:deprecated-function", "local");
 %! N = 6;
 %! A = hilb (N);
 %! apply = @(x) A*x; apply_t = @(x) A'*x;
@@ -232,6 +395,16 @@
 %! cA_test = norm (inv (A), 1) * norm (A, 1);
 %! assert (cA, cA_test, -2^-6);
 
+%!test # to be removed after 4.2
+%! warning ("off", "Octave:deprecated-function", "local");
+%! N = 6;
+%! A = hilb (N);
+%! apply = @(x) A*x; apply_t = @(x) A'*x;
+%! solve = @(x) A\x; solve_t = @(x) A'\x;
+%! cA = condest (apply, apply_t, solve, solve_t, N, 2);
+%! cA_test = norm (inv (A), 1) * norm (A, 1);
+%! assert (cA, cA_test, -2^-6);
+
 %!test
 %! warning ("off", "Octave:nearly-singular-matrix", "local");
 %! N = 12;
@@ -240,3 +413,26 @@
 %! x = A*v;
 %! assert (norm (x, inf), 0, eps);
 
+%!test
+%! N = 6;
+%! A = hilb (N);
+%! solve = @(flag, x) solve_fun (flag, x, A);
+%! cA = condest (A, solve);
+%! cA_test = norm (inv (A), 1) * norm (A, 1);
+%! assert (cA, cA_test, -2^-6);
+%!test
+%! N = 6;
+%! A = hilb (N);
+%! apply = @(flag, x) apply_fun (flag, x, A);
+%! solve = @(flag, x) solve_fun (flag, x, A);
+%! cA = condest (apply, solve);
+%! cA_test = norm (inv (A), 1) * norm (A, 1);
+%! assert (cA, cA_test, -2^-6);
+
+%!test # parameters for apply and solve functions
+%! N = 6;
+%! A = hilb (N);
+%! m = 2;
+%! cA = condest (@apply_fun, @solve_fun, [], A, m);
+%! cA_test = norm (inv (A^2), 1) * norm (A^2, 1);
+%! assert (cA, cA_test, -2^-6);

--- a/scripts/linear-algebra/module.mk	Thu Jun 30 18:20:37 2016 +1000
+++ b/scripts/linear-algebra/module.mk	Mon Jul 25 14:50:29 2016 +0100
@@ -21,8 +21,8 @@
   scripts/linear-algebra/linsolve.m \
   scripts/linear-algebra/logm.m \
   scripts/linear-algebra/normest.m \
+  scripts/linear-algebra/normest1.m \
   scripts/linear-algebra/null.m \
-  scripts/linear-algebra/onenormest.m \
   scripts/linear-algebra/orth.m \
   scripts/linear-algebra/planerot.m \
   scripts/linear-algebra/qzhess.m \

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/linear-algebra/normest1.m	Mon Jul 25 14:50:29 2016 +0100
@@ -0,0 +1,356 @@
+## Copyright (C) 2016 Marco Caliari
+##
+## 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/>.
+
+## -*- texinfo -*-
+## @deftypefn  {Function File} {@var{c} =} normest1 (@var{a})
+## @deftypefnx {Function File} {@var{c} =} normest1 (@var{a}, @var{t})
+## @deftypefnx {Function File} {@var{c} =} normest1 (@var{a}, @var{t}, @var{x0})
+## @deftypefnx {Function File} {@var{c} =} normest1 (@var{afun}, @var{t}, @var{x0}, @var{p1}, @var{p2}, @dots{})
+## @deftypefnx {Function File} {[@var{c}, @var{v}] =} normest1 (@var{a}, @dots{})
+## @deftypefnx {Function File} {[@var{c}, @var{v}, @var{w}] =} normest1 (@var{a}, @dots{})
+## @deftypefnx {Function File} {[@var{c}, @var{v}, @var{w}, @var{it}] =} normest1 (@var{a}, @dots{})
+## Estimate the 1-norm of the matrix @var{a} using a block algorithm.
+##
+## For a medium size matrix @var{a}, @code{norm (@var{a}, 1)} should be
+## used instead. For a large sparse matrix, when only an estimate of the norm
+## is needed, @code{normest1 (@var{a})} might be faster. Moreover, it can be
+## used for the estimate of the 1-norm of a linear
+## operator @var{a} when matrix-vector products @code{@var{a} * @var{x}} and
+## @code{@var{a}' * @var{x}} can be cheaply computed. In this case,
+## instead of the matrix @var{a}, a function
+## @code{@var{afun} (@var{flag}, @var{x})} can be used. It should return:
+##
+## @itemize @bullet
+## @item
+## the dimension @var{n} of @var{a}, if @var{flag} is @qcode{"dim"}
+## @item
+## true if @var{a} is a real operator, if @var{flag} is @qcode{"real"}
+## @item
+## the result @code{@var{a} * @var{x}}, if @var{flag} is @qcode{"notransp"}
+## @item
+## the result @code{@var{a}' * @var{x}}, if @var{flag} is @qcode{"transp"}
+## @end itemize
+##
+## A typical case is @var{a} defined by @code{@var{b} ^ @var{m}},
+## in which the result @code{@var{a} * @var{x}} can be computed without
+## even forming explicitely @code{@var{b} ^ @var{m}} by:
+##
+## @example
+## @var{y} = @var{x};
+## for @var{i} = 1:@var{m}
+##   @var{y} = @var{b} * @var{y};
+## endfor
+## @end example
+##
+## The parameters @var{p1}, @var{p2}, @dots{} are arguments of
+## @code{@var{afun} (@var{flag}, @var{x}, @var{p1}, @var{p2}, @dots{})}.
+##
+## The default value for @var{t} is 2. The algorithm requires
+## matrix-matrix products with sizes @var{n} x @var{n} and
+## @var{n} x @var{t}.
+##
+## The initial matrix @var{x0} should have columns of unit 1-norm.
+## The default initial matrix @var{x0} has the first column
+## @code{ones (@var{n}, 1) / @var{n}}
+## and, if @var{t} >  1, the remaining columns with random elements
+## @code{-1 / @var{n}}, @code{1 / @var{n}}, divided by @var{n}.
+## Therefore, if consistent results are required, the "state" of the
+## random generator should be fixed before invoking @code{normest1}.
+##
+## On output, @var{c} is the desired estimate, @var{v} and @var{w}
+## vectors such that @code{@var{w} = @var{a} * @var{v}}, with
+## @code{norm (@var{w}, 1)} = @code{@var{c} * norm (@var{v}, 1)}.
+## @var{it} contains in @code{@var{it}(1)} the number of iterations
+## (the maximum number is hardcoded to 5) and in  @code{@var{it}(2)}
+## the total number of products @code{@var{a} * @var{x}} or
+## @code{@var{a}' * @var{x}} performed by the algorithm.
+##
+## Reference: @nospell{N. J. Higham and F. Tisseur},
+## @cite{A block algorithm for matrix 1-norm estimation, with and
+## application to 1-norm pseudospectra}, SIAM J. Matrix Anal. Appl.,
+## pp. 1185--1201, Vol 21, No. 4, 2000.
+##
+## @seealso{normest, rand}
+## @end deftypefn
+
+## Ideally, we would set t and X to their default values but Matlab
+## compatibility would require we set the default even when they are
+## empty.
+function [est, v, w, k] = normest1 (A, t = [], X = [], varargin)
+
+  if (nargin < 1)
+    print_usage ();
+  endif
+
+  if (isempty (t))
+    t = 2;
+  endif
+
+  ## FIXME: t < 0 should print trace information
+  if (isnumeric (A) && issquare (A))
+    Aisnum = true;
+    n = rows (A);
+    if ((n <= 4) || (t == n))
+      ## compute directly
+      [est, idx] = max (sum (abs (A), 1), [] ,2);
+      v = zeros (n, 1);
+      v(idx) = 1;
+      w = A(:, idx);
+      k = [0, 1];
+      return
+    else
+      realm = isreal (A);
+    endif
+  elseif (is_function_handle (A))
+    Aisnum = false;
+    n = A ("dim", [], varargin{:});
+    realm = A ("real", [], varargin{:});
+    Afun = @(x) A ("notransp", x, varargin{:});
+    A1fun = @(x) A ("transp", x, varargin{:});
+  else
+    error ("normest1: A must be a square matrix or a function handle");
+  endif
+
+  t = min (t, n);
+
+  if (isempty (X))
+    X = [ones(n, 1), sign(2 * rand(n, t - 1) - 1)];
+    i = 2;
+    imax = min (t, 2 ^ (n - 1));
+    ## There are at most 2^(n-1) unparallel columns, see later.
+    while (i <= imax)
+      if (any (abs (X(:,i)' * X(:,1:i-1)) == n))
+        ## column i is parallel to a colum 1:i-1. Change it.
+        X(:,i) = sign (2 * rand (n, 1) - 1);
+      else
+        i++;
+      endif
+    endwhile
+    X /= n;
+  elseif (columns (X) < t)
+    error ("normest1: X should have %d columns", t);
+  endif
+
+  itmax = 5;
+  ind_hist = zeros (n, 1);
+  est_old = 0;
+  ind = zeros (n, 1);
+  S = zeros (n, t);
+  k = [0; 0];
+  conv = false;
+  while ((! conv) && (k(1) < itmax))
+    k(1)++;
+    if (Aisnum)
+      Y = A * X;
+    else
+      Y = Afun (X);
+    endif
+    k(2)++;
+    [est, j] = max (sum (abs (Y), 1), [], 2);
+    if ((est > est_old) || (k(1) == 2))
+      ind_best = ind(j);
+      w = Y(:, j); # there is an error in Algorithm 2.4
+    endif
+    if ((est <= est_old) && (k(1) >= 2)) # (1) of Algorithm 2.4
+      est = est_old;
+      break; # while
+    endif
+    est_old = est;
+    Sold = S;
+    S = sign (Y);
+    S(S==0) = 1;
+    possible_break = false;
+    if (realm)
+      ## test parallel (only real case)
+      if (all (any (abs (Sold' * S) == n))) # (2) of Algorithm 2.4
+        ## all columns of S parallel to a column of Sold, exit
+        possible_break = true;
+        conv = true;
+      else
+        if (t > 1)
+          ## at least two columns of S are not parallel
+          i = 1;
+          imax = min (t, 2 ^ (n - 1));
+          while (i <= imax)
+            ## The maximum number of parallel columns of length n with entries
+            ## {-1,1} is 2^(n-1). Therefore, if the number of columns of S is
+            ## greater than 2^(n-1), for sure some of them are parallel to some
+            ## columns of Sold. Don't even try to change them (i <= 2^(n-1)).
+            ## Now, check if S(:,i) is parallel to any previous column of S
+            p = (any (abs (S(:,i)' * S(:,1:i-1)) == n));
+            if (p || (any (abs (S(:,i)' * Sold) == n)))
+              ## i-th column of S parallel to a previous or to a
+              ## column of Sold: change it.
+              S(:,i) = sign (2*rand (n, 1)-1);
+            else
+              i++;
+            endif
+          endwhile
+        endif
+      endif
+    endif
+    if (! possible_break)
+      if (Aisnum)
+        Z = A' * S;
+      else
+        Z = A1fun (S); # (3) of Algorithm 2.4
+      endif
+      k(2)++;
+      h = max (abs (Z), [], 2);
+      ind = (1:n)';
+      if ((k(1) >= 2) && (max (h, [], 1) == h(ind_best))) # (4) of Alg. 2.4
+        break; # while
+      endif
+      [h, ind] = sort (h, "descend");
+      if (t > 1)
+       if (all (ind_hist(ind(1:t)))) # (5) of Algorithm 2.4
+          break; # while
+        endif
+        ind = ind(! ind_hist(ind));
+        ## length(ind) could be less than t, especially if t is not << n.
+        ## This is not considered in point (5) of Algorithm 2.4.
+        tmax = min (numel (ind), t);
+        ind = ind(1:tmax);
+      else
+        tmax = 1;
+      endif
+      X = zeros (n, tmax);
+      X(sub2ind (size (X), ind(1:tmax), (1:tmax)')) = 1;
+      ind_hist(ind(1:tmax)) = 1;
+    endif
+  endwhile
+  v = zeros (n, 1);
+  v(ind_best) = 1;
+endfunction
+
+%!function z = afun_A (flag, x, A, n)
+%! switch flag
+%! case {"dim"}
+%!   z = n;
+%! case {"real"}
+%!   z = isreal (A);
+%! case {"transp"}
+%!   z = A' * x;
+%! case {"notransp"}
+%!   z = A * x;
+%! endswitch
+%!endfunction
+%!function z = afun_A_P (flag, x, A, m)
+%! switch flag
+%! case "dim"
+%!   z = length (A);
+%! case "real"
+%!   z = isreal (A);
+%! case "transp"
+%!   z = x; for i = 1:m, z = A' * z;, endfor
+%! case "notransp"
+%!   z = x; for i = 1:m, z = A * z;, endfor
+%! endswitch
+%!endfunction
+
+%!test
+%! A = reshape ((1:16)-8, 4, 4);
+%! assert (normest1 (A), norm (A, 1), eps)
+
+## test t=1
+%!test
+%! A = rand (4); # for positive matrices always work
+%! assert (normest1 (A, 1), norm (A,1), 2 * eps)
+
+## test t=3
+%!test
+%! A = [-0.21825   0.16598   0.19388   0.75297
+%!      -1.47732   0.78443  -1.04254   0.42240
+%!       1.39857  -0.34046   2.28617   0.68089
+%!       0.31205   1.50529  -0.75804  -1.22476];
+%! X = [1,1,-1;1,1,1;1,1,-1;1,-1,-1]/3;
+%! assert (normest1 (A, 3, X), norm (A, 1), 2 * eps)
+
+## test afun
+%!test
+%! A = rand (10);
+%! n = length (A);
+%! afun = @(flag, x) afun_A (flag, x, A, n);
+%! assert (normest1 (afun), norm (A, 1), 2 * eps)
+
+## test afun with parameters
+%!test
+%! A = rand (10);
+%! assert (normest1 (@afun_A_P, [], [], A, 3), norm (A ^ 3, 1), 1000 * eps)
+
+## test output
+%!test
+%! A = reshape (1:16,4,4);
+%! [c, v, w, it] = normest1 (A);
+%! assert (norm (w, 1), c * norm (v, 1), eps)
+
+## test output
+%!test
+%! A = rand (100);
+%! A(A <= 1/3) = -1;
+%! A(A > 1/3 & A <= 2/3) = 0;
+%! A(A > 2/3) = 1;
+%! [c, v, w, it] = normest1 (A, 10);
+%! assert (w, A * v, eps)
+
+%!test
+%! A = rand (5);
+%! c = normest1 (A, 6);
+%! assert (c, norm (A,1), eps)
+
+%!test
+%! A = rand (5);
+%! c = normest1 (A, 2, ones (5, 2) / 5);
+%! assert (c, norm (A,1), eps)
+
+%!test
+%! N = 10;
+%! A = ones (N);
+%! [nm1, v1, w1] = normest1 (A);
+%! [nminf, vinf, winf] = normest1 (A', 6);
+%! assert (nm1, N, -2*eps)
+%! assert (nminf, N, -2*eps)
+%! assert (norm (w1, 1), nm1 * norm (v1, 1), -2*eps)
+%! assert (norm (winf, 1), nminf * norm (vinf, 1), -2*eps)
+
+%!test
+%! N = 5;
+%! A = hilb (N);
+%! [nm1, v1, w1] = normest1 (A);
+%! [nminf, vinf, winf] = normest1 (A', 6);
+%! assert (nm1, norm (A, 1), -2*eps)
+%! assert (nminf, norm (A, inf), -2*eps)
+%! assert (norm (w1, 1), nm1 * norm (v1, 1), -2*eps)
+%! assert (norm (winf, 1), nminf * norm (vinf, 1), -2*eps)
+
+## Only likely to be within a factor of 10.
+%!test
+%! old_state = rand ("state");
+%! unwind_protect
+%!   rand ("state", 42);  # Initialize to guarantee reproducible results
+%!   N = 100;
+%!   A = rand (N);
+%!   [nm1, v1, w1] = normest1 (A);
+%!   [nminf, vinf, winf] = normest1 (A', 6);
+%!   assert (nm1, norm (A, 1), -.1)
+%!   assert (nminf, norm (A, inf), -.1)
+%!   assert (norm (w1, 1), nm1 * norm (v1, 1), -2*eps)
+%!   assert (norm (winf, 1), nminf * norm (vinf, 1), -2*eps)
+%! unwind_protect_cleanup
+%!   rand ("state", old_state);
+%! end_unwind_protect

--- a/scripts/linear-algebra/onenormest.m	Thu Jun 30 18:20:37 2016 +1000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,293 +0,0 @@
-## Copyright (C) 2007-2015 Regents of the University of California
-##
-## 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/>.
-
-## -*- texinfo -*-
-## @deftypefn  {} {[@var{est}, @var{v}, @var{w}, @var{iter}] =} onenormest (@var{A}, @var{t})
-## @deftypefnx {} {[@var{est}, @var{v}, @var{w}, @var{iter}] =} onenormest (@var{apply}, @var{apply_t}, @var{n}, @var{t})
-##
-## Apply @nospell{Higham and Tisseur's} randomized block 1-norm estimator to
-## matrix @var{A} using @var{t} test vectors.
-##
-## If @var{t} exceeds 5, then only 5 test vectors are used.
-##
-## If the matrix is not explicit, e.g., when estimating the norm of
-## @code{inv (@var{A})} given an LU@tie{}factorization, @code{onenormest}
-## applies @var{A} and its conjugate transpose through a pair of functions
-## @var{apply} and @var{apply_t}, respectively, to a dense matrix of size
-## @var{n} by @var{t}.  The implicit version requires an explicit dimension
-## @var{n}.
-##
-## Returns the norm estimate @var{est}, two vectors @var{v} and @var{w} related
-## by norm @code{(@var{w}, 1) = @var{est} * norm (@var{v}, 1)}, and the number
-## of iterations @var{iter}.  The number of iterations is limited to 10 and is
-## at least 2.
-##
-## References:
-##
-## @itemize
-## @item
-## @nospell{N.J. Higham and F. Tisseur}, @cite{A Block Algorithm
-## for Matrix 1-Norm Estimation, with an Application to 1-Norm
-## Pseudospectra}. SIMAX vol 21, no 4, pp 1185-1201.
-## @url{http://dx.doi.org/10.1137/S0895479899356080}
-##
-## @item
-## @nospell{N.J. Higham and F. Tisseur}, @cite{A Block Algorithm
-## for Matrix 1-Norm Estimation, with an Application to 1-Norm
-## Pseudospectra}. @url{http://citeseer.ist.psu.edu/223007.html}
-## @end itemize
-##
-## @seealso{condest, norm, cond}
-## @end deftypefn
-
-## Code originally licensed under:
-##
-## Copyright (c) 2007, Regents of the University of California
-## All rights reserved.
-##
-## Redistribution and use in source and binary forms, with or without
-## modification, are permitted provided that the following conditions
-## are met:
-##
-##    * Redistributions of source code must retain the above copyright
-##      notice, this list of conditions and the following disclaimer.
-##
-##    * Redistributions in binary form must reproduce the above
-##      copyright notice, this list of conditions and the following
-##      disclaimer in the documentation and/or other materials provided
-##      with the distribution.
-##
-##    * Neither the name of the University of California, Berkeley nor
-##      the names of its contributors may be used to endorse or promote
-##      products derived from this software without specific prior
-##      written permission.
-##
-## THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS''
-## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
-## TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
-## PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS AND
-## CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
-## SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
-## LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
-## USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
-## ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
-## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
-## OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
-## SUCH DAMAGE.
-
-## Author: Jason Riedy <ejr@cs.berkeley.edu>
-## Keywords: linear-algebra norm estimation
-## Version: 0.2
-
-function [est, v, w, iter] = onenormest (varargin)
-
-  if (nargin < 1 || nargin > 4)
-    print_usage ();
-  endif
-
-  default_t = 5;
-  itmax = 10;
-
-  if (isnumeric (varargin{1}))
-    [n, nc] = size (varargin{1});
-    if (n != nc)
-      error ("onenormest: matrix must be square");
-    endif
-    apply = @(x) varargin{1} * x;
-    apply_t = @(x) varargin{1}' * x;
-    if (nargin > 1)
-      t = varargin{2};
-    else
-      t = min (n, default_t);
-    endif
-    issing = isa (varargin{1}, "single");
-  else
-    if (nargin < 3)
-      print_usage ();
-    endif
-    apply = varargin{1};
-    apply_t = varargin{2};
-    n = varargin{3};
-    if (nargin > 3)
-      t = varargin{4};
-    else
-      t = default_t;
-    endif
-    issing = isa (n, "single");
-  endif
-
-  ## Initial test vectors X.
-  X = rand (n, t);
-  X ./= ones (n,1) * sum (abs (X), 1);
-
-  ## Track if a vertex has been visited.
-  been_there = zeros (n, 1);
-
-  ## To check if the estimate has increased.
-  est_old = 0;
-
-  ## Normalized vector of signs.
-  S = zeros (n, t);
-
-  if (issing)
-    myeps = eps ("single");
-    X = single (X);
-  else
-    myeps = eps;
-  endif
-
-  for iter = 1 : itmax + 1
-    Y = feval (apply, X);
-
-    ## Find the initial estimate as the largest A*x.
-    [est, ind_best] = max (sum (abs (Y), 1));
-    if (est > est_old || iter == 2)
-      w = Y(:,ind_best);
-    endif
-    if (iter >= 2 && est < est_old)
-      ## No improvement, so stop.
-      est = est_old;
-      break;
-    endif
-
-    est_old = est;
-    S_old = S;
-    if (iter > itmax),
-      ## Gone too far.  Stop.
-      break;
-    endif
-
-    S = sign (Y);
-
-    ## Test if any of S are approximately parallel to previous S
-    ## vectors or current S vectors.  If everything is parallel,
-    ## stop.  Otherwise, replace any parallel vectors with
-    ## rand{-1,+1}.
-    partest = any (abs (S_old' * S - n) < 4*eps*n);
-    if (all (partest))
-      ## All the current vectors are parallel to old vectors.
-      ## We've hit a cycle, so stop.
-      break;
-    endif
-    if (any (partest))
-      ## Some vectors are parallel to old ones and are cycling,
-      ## but not all of them.  Replace the parallel vectors with
-      ## rand{-1,+1}.
-      numpar = sum (partest);
-      replacements = 2*(rand (n,numpar) < 0.5) - 1;
-      S(:,partest) = replacements;
-    endif
-    ## Now test for parallel vectors within S.
-    partest = any ((S' * S - eye (t)) == n);
-    if (any (partest))
-      numpar = sum (partest);
-      replacements = 2*(rand (n,numpar) < 0.5) - 1;
-      S(:,partest) = replacements;
-    endif
-
-    Z = feval (apply_t, S);
-
-    ## Now find the largest non-previously-visted index per vector.
-    h = max (abs (Z),2);
-    [mh, mhi] = max (h);
-    if (iter >= 2 && mhi == ind_best)
-      ## Hit a cycle, stop.
-      break;
-    endif
-    [h, ind] = sort (h, 'descend');
-    if (t > 1)
-      firstind = ind(1:t);
-      if (all (been_there(firstind)))
-        ## Visited all these before, so stop.
-        break;
-      endif
-      ind = ind(! been_there(ind));
-      if (length (ind) < t)
-        ## There aren't enough new vectors, so we're practically
-        ## in a cycle.  Stop.
-        break;
-      endif
-    endif
-
-    ## Visit the new indices.
-    X = zeros (n, t);
-    for zz = 1 : t
-      X(ind(zz),zz) = 1;
-    endfor
-    been_there(ind(1 : t)) = 1;
-  endfor
-
-  ## The estimate est and vector w are set in the loop above.
-  ## The vector v selects the ind_best column of A.
-  v = zeros (n, 1);
-  v(ind_best) = 1;
-
-endfunction
-
-
-%!demo
-%! N = 100;
-%! A = randn (N) + eye (N);
-%! [L,U,P] = lu (A);
-%! nm1inv = onenormest (@(x) U\(L\(P*x)), @(x) P'*(L'\(U'\x)), N, 30)
-%! norm (inv (A), 1)
-
-%!test
-%! N = 10;
-%! A = ones (N);
-%! [nm1, v1, w1] = onenormest (A);
-%! [nminf, vinf, winf] = onenormest (A', 6);
-%! assert (nm1, N, -2*eps);
-%! assert (nminf, N, -2*eps);
-%! assert (norm (w1, 1), nm1 * norm (v1, 1), -2*eps);
-%! assert (norm (winf, 1), nminf * norm (vinf, 1), -2*eps);
-
-%!test
-%! N = 10;
-%! A = ones (N);
-%! [nm1, v1, w1] = onenormest (@(x) A*x, @(x) A'*x, N, 3);
-%! [nminf, vinf, winf] = onenormest (@(x) A'*x, @(x) A*x, N, 3);
-%! assert (nm1, N, -2*eps);
-%! assert (nminf, N, -2*eps);
-%! assert (norm (w1, 1), nm1 * norm (v1, 1), -2*eps);
-%! assert (norm (winf, 1), nminf * norm (vinf, 1), -2*eps);
-
-%!test
-%! N = 5;
-%! A = hilb (N);
-%! [nm1, v1, w1] = onenormest (A);
-%! [nminf, vinf, winf] = onenormest (A', 6);
-%! assert (nm1, norm (A, 1), -2*eps);
-%! assert (nminf, norm (A, inf), -2*eps);
-%! assert (norm (w1, 1), nm1 * norm (v1, 1), -2*eps);
-%! assert (norm (winf, 1), nminf * norm (vinf, 1), -2*eps);
-
-## Only likely to be within a factor of 10.
-%!test
-%! old_state = rand ("state");
-%! restore_state = onCleanup (@() rand ("state", old_state));
-%! rand ("state", 42);  # Initialize to guarantee reproducible results
-%! N = 100;
-%! A = rand (N);
-%! [nm1, v1, w1] = onenormest (A);
-%! [nminf, vinf, winf] = onenormest (A', 6);
-%! assert (nm1, norm (A, 1), -.1);
-%! assert (nminf, norm (A, inf), -.1);
-%! assert (norm (w1, 1), nm1 * norm (v1, 1), -2*eps);
-%! assert (norm (winf, 1), nminf * norm (vinf, 1), -2*eps);
-