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

And David Bateman writes:
> I believe the rcond function should take rcond from the LAPACK
> factorization, at least I believe is what matlab does. [...]

For various reasons, I'm not referring to MATLAB(tm) itself or
any docs other than what is on-line and trivially findable.  I
have no idea if MATLAB's rcond applies to sparse matrices.  I
know LAPACK's does not.

That said, I have been losing the argument to use the block
estimator in LAPACK...  We may weaken the routine slightly
instead.  For "small" problems, computing RCOND is very, very
expensive right now.  The current norm estimator uses a very safe
but very slow triangular solve.  We likely will replace that with
the fast solve, followed by a check for overflow and re-running
with the slow solve if necessary.  That will catch and re-scale
overflows, but it won't catch possibly avoidable underflows
during estimation.

So it's up to y'all what rcond() should be.  The routine I
included computes the inverse of the 1-norm condition number of A
(iirc), and LAPACK's RCOND returned by xyySVX is the inverse of
the inf-norm condition number of the matrix that actually is
factored.  That matrix may have been equilibrated, so determining
how LAPACK's RCOND relates to your original matrix requires
examining a bunch of parameters MATLAB's rcond() doesn't return.
Without experimenting, I can't tell what MATLAB's rcond() really
does.

The only use for LAPACK's RCOND directly is to see if the
factorization algorithm came close to solving with a nearly a
singular matrix when forming a Schur complement.  It tells you
very little about the actual problem you are trying to solve or
even about any errors when using the factorization; RCOND may be
measured with an entirely different norm (scaled by
equilibration) than your problem.

And if you pick up our extra-precise refinement drivers (once we
actually ship the damned things, see LAWN 165), you will have
more and more relevant condition numbers from which to chose.

> 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.

>From MathWorks' pages, it seems normest has a very different
interface.  Hence the cumbersome block_onenorm_est name.

But your description is not the RCOND returned from LAPACK's
linear system drivers.  ;)

> condest is one of the annoying missing sparse functions in
> Octave...

Amen.  Limits Octave's applicability for many colleagues.

> I'd hesitate to include it though with the rcond function with
> its current name...

I have little ego about these things.  Use, change, and rename as
you see fit.  Just thought I should send along my good-enough
versions as a starting point.

I'm not in a position to get into the copyright assignment
argument here right now, so feel free to make sufficient,
copyrightable changes that are assigned and licensed
appropriately.  ;) These are just the routines I'm using right
now.

Jason