Hey Benoit,
The class looks good! I have only very few comments. First, I am wondering why useThreshold is returning a reference to the LU object and that actually only according to the signature since the current implementation returns nothing. I would simply make the return type void. Second, for the threshold retrieval, I woul make only 'RealScalar threshold() const' public and move the defaultThreshold() method to the private section or even consider removing it. When the default threshold is set, the thresold() function does return it - in case it is not set, I don't see why a 'public' user should be interested in it.
Finally, I am wondering why all the variables are 'protected'. Do we expect somebody to inherit from LU? If that were the case it should have a virtual destructor. There is not harm in declaring the variables as protected - I was simply wondering.
Cheers,
Hauke
On Sun, Oct 18, 2009 at 7:14 AM, Benoit Jacob
<jacob.benoit.1@xxxxxxxxx> wrote:
Hi,
This is mostly for Gael and Hauke (if you're busy, skip to 5. below)
who I bugged recently about that, and a followup to the goog old
thread "LU precision tuning" on this list.
It's been very very tedious but i think we're now on the good side of
the API refactoring in LU (and once I've finished LU the rest should
be the same).
See my fork http://bitbucket.org/bjacob/refactor-solve-api/
1) What Eigen 2.0 des.
In Eigen 2.0, we determine the rank in the LU decomposition itself. As
soon as some pivot is small enough to be considered zero (we have a
threshold for that, see "LU precision tuning" on this list), we stop
the LU decomposition. In LU::solve(), we also use that threshold to
determine whether a solution exists. If no solution exists according
to that test, then we don't try to return a solution.
2) What LAPACK does.
In dgesvx, there are 2 cases. If a pivot is exactly zero, so the
matrix is exactly singular, dgesvx bails out and doesn't give a
solution. This means that dgesvx can't be used on simple 2x2 examples
like the matrix (1&0\\0&0) and the right-hand-side (1\\0). If no pivot
is exactly zero, dgesvx just computes a solution without worrying
about accuracy if a pivot was near zero.
3) What the new Eigen (in the fork) *cough* 3.0 *cough* does.
As previously discussed, i realized that we should try as hard as
possible to get a solution, even approximate; the only case in which
we need to bail out is if continuing would generate inf/nan values
(that would be worthless to the user anyway); and the precision check
is something that the user can always do by himself afterwards. So
that means that we get much closer to LAPACK. On the other hand, I
don't see the need to restrict to nonsingular matrices like LAPACK
does. If some pivot is zero, we have an equation of the form 0*X=Y, we
can just ignore Y and set X=0, which can give a useful approximate
solution if Y was itself meant to be zero.
4) precision of the rank determination (problem solved, was a PEBKAC)
As you may know, in the past we've found it very difficult to have
good accuracy for rank() methods in e.g. LU. We always had a small
amount of failures in the LU test.
I just found the reason. Our LU with full pivoting algo consists in
taking the biggest element in the remaining corner and moving it to
become the next pivot. From there we assumed that the pivots were
automatically sorted in decreasing order. That was approximately true,
but not exactly true. In practice it was possible to get cases where
the (n+1)-th pivot was several times bigger than the n-th pivot. So
now we no longer make this assumption, we manually sort the pivots as
needed (not possible to do in the decomposition itself), and bingo! I
can now run "test_lu r1000" and it succeeds!
5) The threshold to determine when a pivot is nonzero is, by default,
the same as before (Higham formula). The user can override that by
calling usePrecision(value) and I introduced a special type Default_t
so one can come back to the default by doing useThreshold(Default). Ah
yes and nevermind what i said about virtual methods, I had been
smoking some good raccoon fur blend (a Canadian Thanksgiving
tradition).
Benoit