new file mode 100644
--- /dev/null
+++ b/scripts/linear-algebra/onenormest.m
@@ -0,0 +1,265 @@
+## Copyright (C) 2007, 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 {Function File} {[@var{est}, @var{v}, @var{w}, @var{iter}] =} onenormest (@var{A}, @var{t}) 
+## @deftypefnx {Function File} {[@var{est}, @var{v}, @var{w}, @var{iter}] =} onenormest (@var{apply}, @var{apply_t}, @var{n}, @var{t})
+##
+## Apply 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 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 Nicholas J. Higham and Françoise Tisseur, "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 Nicholas J. Higham and Françoise Tisseur, "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.
+##
+## Relicensed to GPL for inclusion in Octave.
+
+## Author: Jason Riedy <ejr@cs.berkeley.edu>
+## Keywords: linear-algebra norm estimation
+## Version: 0.2
+
+function [est, v, w, iter] = onenormest (varargin)
+
+  if (size (varargin, 2) < 1 || size (varargin, 2) > 4)
+    print_usage ();
+  endif
+
+  default_t = 5;
+  itmax = 10;
+
+  if (ismatrix (varargin{1}))
+    n = size (varargin{1}, 1);
+    if n != size (varargin{1}, 2),
+      error ("onenormest: matrix must be square.");
+    endif
+    apply = @(x) varargin{1} * x;
+    apply_t = @(x) varargin{1}' * x;
+    if (size (varargin) > 1)
+      t = varargin{2};
+    else
+      t = min (n, default_t);
+    endif
+  else
+    if (size (varargin, 2) < 3)
+      print_usage();
+    endif
+    n = varargin{3};
+    apply = varargin{1};
+    apply_t = varargin{2};
+    if (size (varargin) > 3)
+      t = varargin{4};
+    else
+      t = default_t;
+    endif
+  endif
+
+  ## Initial test vectors X.
+  X = rand (n, t);
+  X = X ./ (ones (n,1) * sum (abs (X), 1));
+
+  been_there = zeros (n, 1); # Track if a vertex has been visited.
+  est_old = 0; # To check if the estimate has increased.
+  S = zeros (n, t); # Normalized vector of signs.  The normalization is 
+
+  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
+%!  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)