Re: rcond, condest, and a block 1-norm estimator

I believe the rcond function should take rcond from the LAPACK
factorization, at least I believe is what matlab does. Matlab calls the
normest1 function, that is basically your block_onenorm_est, and then
multiples the estimated inverse norm by norm(A,1) to get the condition
number estimate.

condest is one of the annoying missing sparse functions in Octave and I
appreciate that your code treats sparse matrices as a special case. I'd
hesitate to include it though with the rcond function with its current
name...

D.




Jason Riedy wrote:
> Appended are the rcond(), condest(), and block_onenorm_est()
> routines I've been using. I don't think they exactly match the
> MATLAB(tm) routines, but I also don't care. ;) The files also
> are living at
>   http://www.cs.berkeley.edu/~ejr/dev/linalg/
> for now.
> =

> They have been moderately beaten around, but there still may be
> bugs lurking. I'm not 100% sure I'm handling complex matrices
> correctly, but I also haven't run into any noticable problems.
> =

> Improvements and bug fixes gleefully accepted.
> =

> Jason
> =

> =

> =

> ------------------------------------------------------------------------
> =

> ## Copyright (C) 2007, Regents of the University of California
> ## -*- mode: octave; -*-
> ##
> ## This program is free software distributed under the "modified" or
> ## 3-clause BSD license appended to this file.
> =

> function [est, v, w] =3D rcond (varargin)
>   ## -*- texinfo -*-
>   ## @deftypefn {Function File} {[@var{est}, @var{v}] =3D } rcond (@var{A=
}, @var{t} =3D min (size (A, 1), 5)) =

>   ## @deftypefnx {Function File} {[@var{est}, @var{v}] =3D } rcond (@var{=
A}, @var{SOLVE}, @var{SOVLE_T}, @var{t} =3D min (n, 5))
>   ## @deftypefnx {Function File} {[@var{est}, @var{v}] =3D } rcond (@var{=
APPLY}, @var{APPLY_T}, @var{SOLVE}, @var{SOVLE_T}, @var{n}, @var{t} =3D min=
 (n, 5))
>   ##
>   ## Estimate the inverse of the 1-norm condition number of a matrix
>   ## matrix @var{A} using @var{t} test vectors pusing the randomized
>   ## 1-norm estimator in block_onenorm_est.  If the matrix is not
>   ## explicit, e.g. when estimating the condition number of @var{A}
>   ## given an LU factorization, rcond 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}.
>   ## @item SOLVE
>   ## @code{A \ b} for a matrix @code{b} of size @var{n} by @var{t}.
>   ## @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{rcond} uses a randomized algorithm to approximate
>   ## the 1-norms.
>   ##
>   ## @code{rcond} returns the inverse of the 1-norm condition estimate
>   ## @var{est} and a vector @var{v} satisfying @code{norm
>   ## (@var{A}*@var{v}, 1) =3D=3D norm (@var{A}, 1) * norm (@var{v}, 1) *
>   ## @var{est}}. When @var{est} is large, @var{v} is an approximate null
>   ## vector.
>   ##
>   ## References: =

>   ## @itemize
>   ## @item Nicholas J. Higham and Fran=E7oise Tisseur, "A Block Algorithm
>   ## for Matrix 1-Norm Estimation, with an Application to 1-Norm
>   ## Pseudospectra." SIMAX vol 21, no 4, pp 1185-1201.
>   ## http://dx.doi.org/10.1137/S0895479899356080
>   ## @item Nicholas J. Higham and Fran=E7oise Tisseur, "A Block Algorithm
>   ## for Matrix 1-Norm Estimation, with an Application to 1-Norm
>   ## Pseudospectra." http://citeseer.ist.psu.edu/223007.html
>   ## @end itemize
>   ##
>   ## @seealso{block_onenorm_est, condest, norm, cond}
>   ##
>   ## @end deftypefn
> =

>   ## Author: Jason Riedy <[email protected]>
>   ## Keywords: linear-algebra norm estimation
>   ## Version: 0.2
> =

> %!demo
> %!  N =3D 100;
> %!  A =3D randn (N) + eye (N);
> %!  rcond (A)
> %!  [L,U,P] =3D lu (A);
> %!  rcond (A, @(x) U\ (L\ (P*x)), @(x) P'*(L'\ (U'\x)))
> %!  rcond (@(x) A*x, @(x) A'*x, @(x) U\ (L\ (P*x)), @(x) P'*(L'\ (U'\x)),=
 N)
> %!  1 / (norm (inv (A), 1) * norm (A, 1))
> =

>   if size (varargin, 2) < 1 || size (varargin, 2) > 5,
>     usage("rcond: Incorrect arguments.");
>   endif
> =

>   DEFAULT_T =3D 5;
>   ITMAX =3D 10;
> =

>   if ismatrix (varargin{1}),
>     n =3D size (varargin{1}, 1);
>     if n !=3D size (varargin{1}, 2),
>       usage("Matrix must be square.");
>     endif
>     A =3D varargin{1};
> =

>     if size (varargin, 2) > 1,
>       if isscalar (varargin{2}),
> 	t =3D varargin{2};
>       else
> 	if size (varargin, 2) < 3,
> 	  usage("Must supply both SOLVE and SOLVE_T.");
> 	else
> 	  SOLVE =3D varargin{2};
> 	  SOLVE_T =3D varargin{3};
> 	  if size (varargin, 2) > 3,
> 	    t =3D varargin{4};
> 	  endif
> 	endif
>       endif
>     endif
>   else
>     if size (varargin, 2) < 5,
>       usage("Implicit form of rcond requires at least 5 arguments.");
>     endif
>     APPLY =3D varargin{1};
>     APPLY_T =3D varargin{2};
>     SOLVE =3D varargin{3};
>     SOLVE_T =3D varargin{4};
>     n =3D varargin{5};
>     if !isscalar (n),
>       usage("Dimension argument of implicit form must be scalar.");
>     endif
>     if size (varargin, 2) > 5,
>       t =3D varargin{6};
>     endif
>   endif
> =

>   if !exist ("t", "var"),
>     t =3D min (n, DEFAULT_T);
>   endif
> =

>   if !exist ("SOLVE", "var"),
>     if issparse (A),
>       colord =3D colamd(A);
>       Pc =3D speye(n);
>       Pc(:,colord) =3D Pc;
>       [L,U,P] =3D lu(A(:,Pc));
>       SOLVE =3D @(x) Pc' * (U\ (L\ (P*x)));
>       SOLVE_T =3D @(x) P'*(L'\ (U'\ (Pc*x)));
>     else
>       [L,U,P] =3D lu(A);
>       SOLVE =3D @(x) U\ (L\ (P*x));
>       SOLVE_T =3D @(x) P' * (L'\ (U'\x));
>     endif
>   endif
> =

>   if exist ("A", "var"),
>     Anorm =3D norm(A, 1);
>   else
>     Anorm =3D block_onenorm_est(APPLY, APPLY_T, n, t);
>   endif
> =

>   [Ainv_norm, v, w] =3D block_onenorm_est(SOLVE, SOLVE_T, n, t);
> =

>   est =3D 1/Anorm;
>   est /=3D Ainv_norm;
>   v =3D w / norm (w, 1);
> =

> endfunction
> =

> ## Yes, these test bounds are really loose.  There's
> ## enough randomization to trigger odd cases with hilb().
> =

> %!test
> %!  N =3D 6;
> %!  A =3D hilb (N);
> %!  cA =3D 1 / rcond (A);
> %!  cA_test =3D norm (inv (A), 1) * norm (A, 1);
> %!  assert (cA, cA_test, 2**-12);
> =

> %!test
> %!  N =3D 6;
> %!  A =3D hilb (N);
> %!  SOLVE =3D @(x) A\x; SOLVE_T =3D @(x) A'\x;
> %!  cA =3D 1 / rcond (A, SOLVE, SOLVE_T);
> %!  cA_test =3D norm (inv (A), 1) * norm (A, 1);
> %!  assert (cA, cA_test, 2**-12);
> =

> %!test
> %!  N =3D 6;
> %!  A =3D hilb (N);
> %!  APPLY =3D @(x) A*x; APPLY_T =3D @(x) A'*x;
> %!  SOLVE =3D @(x) A\x; SOLVE_T =3D @(x) A'\x;
> %!  cA =3D 1 / rcond (APPLY, APPLY_T, SOLVE, SOLVE_T, N);
> %!  cA_test =3D norm (inv (A), 1) * norm (A, 1);
> %!  assert (cA, cA_test, 2**-6);
> =

> %!test
> %!  N =3D 12;
> %!  A =3D hilb (N);
> %!  [rcondA, v] =3D rcond (A);
> %!  x =3D A*v;
> %!  assert (norm(x, inf), 0, eps);
> =

> ## 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 distribut=
ion.
> ##     * Neither the name of the University of California, Berkeley nor t=
he
> ##       names of its contributors may be used to endorse or promote prod=
ucts
> ##       derived from this software without specific prior written permis=
sion.
> ##
> ## THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AN=
D ANY
> ## EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPL=
IED
> ## WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
> ## DISCLAIMED. IN NO EVENT SHALL THE REGENTS AND CONTRIBUTORS BE LIABLE F=
OR ANY
> ## DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAM=
AGES
> ## (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SER=
VICES;
> ## LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSE=
D 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.
> =

> =

> ------------------------------------------------------------------------
> =

> ## Copyright (C) 2007, Regents of the University of California
> ## -*- mode: octave; -*-
> ##
> ## This program is free software distributed under the "modified" or
> ## 3-clause BSD license appended to this file.
> =

> function [est, v, w, iter] =3D block_onenorm_est (varargin)
>   ## -*- texinfo -*-
>   ## @deftypefn {Function File} {[@var{est}, @var{v}, @var{w}, @var{iter}=
] =3D } block_onenorm_est (@var{A}, @var{t} =3D min (size (A, 1), 5)) =

>   ## @deftypefnx {Function File} {[@var{est}, @var{v}, @var{w}, @var{iter=
}] =3D } block_onenorm_est (@var{APPLY}, @var{APPLY_T}, @var{n}, @var{t} =
=3D min (n, 5))
>   ##
>   ## Apply Higham and Tisseur's randomized block 1-norm estimator to
>   ## matrix @var{A} using @var{t} test vectors.  If the matrix is not
>   ## explicit, e.g. when estimating the norm of @code{inv(@var{A})} given=
 an
>   ## LU factorization, block_onenorm_est 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) =3D @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=E7oise Tisseur, "A Block Algorithm
>   ## for Matrix 1-Norm Estimation, with an Application to 1-Norm
>   ## Pseudospectra." SIMAX vol 21, no 4, pp 1185-1201.
>   ## http://dx.doi.org/10.1137/S0895479899356080
>   ## @item Nicholas J. Higham and Fran=E7oise Tisseur, "A Block Algorithm
>   ## for Matrix 1-Norm Estimation, with an Application to 1-Norm
>   ## Pseudospectra." http://citeseer.ist.psu.edu/223007.html
>   ## @end itemize
>   ##
>   ## @seealso{condest, norm, cond}
>   ##
>   ## @end deftypefn
> =

>   ## Author: Jason Riedy <[email protected]>
>   ## Keywords: linear-algebra norm estimation
>   ## Version: 0.2
> =

> %!demo
> %!  N =3D 100;
> %!  A =3D randn(N) + eye(N);
> %!  [L,U,P] =3D lu(A);
> %!  nm1inv =3D block_onenorm_est(@(x) U\(L\(P*x)), @(x) P'*(L'\(U'\x)), \
> %!                             N, 30)
> %!  norm(inv(A), 1)
> =

>   if size (varargin, 2) < 1 || size (varargin, 2) > 4,
>     print_usage();
>   endif
> =

>   DEFAULT_T =3D 5;
>   ITMAX =3D 10;
> =

>   if ismatrix (varargin{1}),
>     n =3D size (varargin{1}, 1);
>     if n !=3D size (varargin{1}, 2),
>       error("Matrix must be square.");
>     endif
>     APPLY =3D @(x) varargin{1} * x;
>     APPLY_T =3D @(x) varargin{1}' * x;
>     if size (varargin) > 1,
>       t =3D varargin{2};
>     else
>       t =3D min (n, DEFAULT_T);
>     endif
>   else
>     if size (varargin, 2) < 3,
>       print_usage();
>     endif
>     n =3D varargin{3};
>     APPLY =3D varargin{1};
>     APPLY_T =3D varargin{2};
>     if size (varargin) > 3,
>       t =3D varargin{4};
>     else
>       t =3D DEFAULT_T;
>     endif
>   endif
> =

>   ## Initial test vectors X.
>   X =3D rand (n, t);
>   X =3D X ./ (ones (n,1) * sum (abs (X), 1));
> =

>   been_there =3D zeros (n, 1); # Track if a vertex has been visited.
>   est_old =3D 0; # To check if the estimate has increased.
>   S =3D zeros (n, t); # Normalized vector of signs.  The normalization is =

> =

>   for iter=3D1:ITMAX+1,
>     Y =3D feval (APPLY, X);
> =

>     ## Find the initial estimate as the largest A*x.
>     [est, ind_best] =3D max (sum (abs (Y), 1));
>     if (est > est_old || iter =3D=3D 2),
>       w =3D Y(:,ind_best);
>     endif
>     if (iter >=3D 2 && est < est_old),
>       ## No improvement, so stop.
>       est =3D est_old;
>       break;
>     endif
> =

>     est_old =3D est;
>     S_old =3D S;
>     if (iter > ITMAX),
>       ## Gone too far.  Stop.
>       break;
>     endif
> =

>     S =3D 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 =3D 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 =3D sum (partest);
>       replacements =3D 2*(rand (n,numpar) < 0.5) - 1;
>       S(:,partest) =3D replacements;
>     endif
>     ## Now test for parallel vectors within S.
>     partest =3D any ( (S' * S - eye (t)) =3D=3D n );
>     if any (partest),
>       numpar =3D sum (partest);
>       replacements =3D 2*(rand (n,numpar) < 0.5) - 1;
>       S(:,partest) =3D replacements;
>     endif
>     =

>     Z =3D feval (APPLY_T, S);
> =

>     ## Now find the largest non-previously-visted index per
>     ## vector.
>     h =3D max (abs (Z),2);
>     [mh, mhi] =3D max (h);
>     if iter >=3D 2 && mhi =3D=3D ind_best,
>       ## Hit a cycle, stop.
>       break;
>     endif
>     [h, ind] =3D sort (h, 'descend');
>     if t > 1,
>       firstind =3D ind(1:t);
>       if all (been_there(firstind)),
> 	## Visited all these before, so stop.
> 	break;
>       endif
>       ind=3Dind(!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 =3D zeros (n, t);
>     for zz =3D 1:t,
>       X(ind(zz),zz) =3D 1;
>     endfor
>     been_there(ind(1:t)) =3D 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 =3D zeros (n, 1);
>   v(ind_best) =3D 1;
> endfunction
> =

> %!test
> %!  N =3D 10;
> %!  A =3D ones (N);
> %!  [nm1, v1, w1] =3D block_onenorm_est (A);
> %!  [nminf, vinf, winf] =3D block_onenorm_est (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 =3D 10;
> %!  A =3D ones (N);
> %!  [nm1, v1, w1] =3D block_onenorm_est (@(x) A*x, @(x) A'*x, N, 3);
> %!  [nminf, vinf, winf] =3D block_onenorm_est (@(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 =3D 5;
> %!  A =3D hilb (N);
> %!  [nm1, v1, w1] =3D block_onenorm_est (A);
> %!  [nminf, vinf, winf] =3D block_onenorm_est (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 =3D 100;
> %!  A =3D rand (N);
> %!  [nm1, v1, w1] =3D block_onenorm_est (A);
> %!  [nminf, vinf, winf] =3D block_onenorm_est (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)
> =

> ## 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 distribut=
ion.
> ##     * Neither the name of the University of California, Berkeley nor t=
he
> ##       names of its contributors may be used to endorse or promote prod=
ucts
> ##       derived from this software without specific prior written permis=
sion.
> ##
> ## THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AN=
D ANY
> ## EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPL=
IED
> ## WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
> ## DISCLAIMED. IN NO EVENT SHALL THE REGENTS AND CONTRIBUTORS BE LIABLE F=
OR ANY
> ## DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAM=
AGES
> ## (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SER=
VICES;
> ## LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSE=
D 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.
> =

> =

> ------------------------------------------------------------------------
> =

> ## Copyright (C) 2007, Regents of the University of California
> ## -*- mode: octave; -*-
> ##
> ## This program is free software distributed under the "modified" or
> ## 3-clause BSD license appended to this file.
> =

> function [est, v] =3D condest (varargin)
>   ## -*- texinfo -*-
>   ## @deftypefn {Function File} {[@var{est}, @var{v}] =3D } condest (@var=
{A}, @var{t} =3D min (size (A, 1), 5)) =

>   ## @deftypefnx {Function File} {[@var{est}, @var{v}] =3D } condest (@va=
r{A}, @var{SOLVE}, @var{SOVLE_T}, @var{t} =3D min (n, 5))
>   ## @deftypefnx {Function File} {[@var{est}, @var{v}] =3D } condest (@va=
r{APPLY}, @var{APPLY_T}, @var{SOLVE}, @var{SOVLE_T}, @var{n}, @var{t} =3D m=
in (n, 5))
>   ##
>   ## Estimate the 1-norm condition number of a matrix matrix @var{A}
>   ## using @var{t} test vectors pusing the randomized 1-norm estimator in
>   ## block_onenorm_est.  If the matrix is not explicit, e.g. when
>   ## estimating the condition number of @var{A} given an LU
>   ## factorization, 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}.
>   ## @item SOLVE
>   ## @code{A \ b} for a matrix @code{b} of size @var{n} by @var{t}.
>   ## @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} returns the 1-norm condition estimate @var{est} and
>   ## a vector @var{v} satisfying @code{norm (@var{A}*@var{v}, 1) =3D=3D n=
orm
>   ## (@var{A}, 1) * norm (@var{v}, 1) / @var{est}}. When @var{est} is
>   ## large, @var{v} is an approximate null vector.
>   ##
>   ## References: =

>   ## @itemize
>   ## @item Nicholas J. Higham and Fran=E7oise Tisseur, "A Block Algorithm
>   ## for Matrix 1-Norm Estimation, with an Application to 1-Norm
>   ## Pseudospectra." SIMAX vol 21, no 4, pp 1185-1201.
>   ## http://dx.doi.org/10.1137/S0895479899356080
>   ## @item Nicholas J. Higham and Fran=E7oise Tisseur, "A Block Algorithm
>   ## for Matrix 1-Norm Estimation, with an Application to 1-Norm
>   ## Pseudospectra." http://citeseer.ist.psu.edu/223007.html
>   ## @end itemize
>   ##
>   ## @seealso{rcond, block_onenorm_est, norm, cond}
>   ##
>   ## @end deftypefn
> =

>   ## Author: Jason Riedy <[email protected]>
>   ## Keywords: linear-algebra norm estimation
>   ## Version: 0.2
> =

> %!demo
> %!  N =3D 100;
> %!  A =3D randn (N) + eye (N);
> %!  condest (A)
> %!  [L,U,P] =3D 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)
> =

>   [est, v] =3D rcond(varargin{:});
>   est =3D 1/est;
> endfunction
> =

> ## These tests are the same as in rcond().
> =

> %!test
> %!  N =3D 6;
> %!  A =3D hilb (N);
> %!  cA =3D condest (A);
> %!  rcA =3D 1 / rcond (A);
> %!  assert (cA, rcA, 2**-12);
> =

> ## 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 distribut=
ion.
> ##     * Neither the name of the University of California, Berkeley nor t=
he
> ##       names of its contributors may be used to endorse or promote prod=
ucts
> ##       derived from this software without specific prior written permis=
sion.
> ##
> ## THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AN=
D ANY
> ## EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPL=
IED
> ## WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
> ## DISCLAIMED. IN NO EVENT SHALL THE REGENTS AND CONTRIBUTORS BE LIABLE F=
OR ANY
> ## DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAM=
AGES
> ## (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SER=
VICES;
> ## LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSE=
D 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.
> =

> =

> ------------------------------------------------------------------------
> =

> _______________________________________________
> Octave-sources mailing list
> [email protected]
> https://www.cae.wisc.edu/mailman/listinfo/octave-sources