Cholseski update of rank one

[nice44]

Since they are still not in Lapack (or I'm blind?), and since I had a need for it, I have a made port of Choleski update/downdate from the original LINPACK dchud and dchdd routines.
http://www.netlib.org/linpack/dchud.f
http://www.netlib.org/linpack/dchdd.f
I followed the suggestions and guidelines of Matthias Seeger from:
Low Rank Updates for the Cholesky Decomposition - Matthias Seeger
http://infoscience.epfl.ch/record/16146 ... update.pdf
(I noted that he previouslyposted on the subject in 2007)
The main changes are :
- uses of BLAS drot and drotg for Given's rotations
- forces positive elements on the diagonal of R
- accepts upper or lower Cholesky factor R
I also changed sqrt(1-norm**2) into a sin(acos(norm)) or cos(asin(norm) so as to be more stable if norm is close to 1.
This might be questionable because there is no use of inverse trigonometric functions inthe whole Lapack, but if you have a better way...

I named the routines dchur1 and dchdr1 for choleski update/downdate of rank 1, but of course, this can be changed at will...
You can find the code attached, it is released under MIT license.
Tell me if it is of any interest.

Cheers

Nicolas Cellier

[admin]

Thank you Nicolas,
As it is, it is not possible to integrate this code directly into LAPACK.
- The routine needs a better description, name ,...
- We need the double complex version of the code (we have generator to get the single precision and the complex one)
- We need to have some TESTING that goes with it
Also, if we decide to integrate the routine, we may want to look for a code that use a "blocking" technic. But it would be fine to start with yours.

Nicolas, if you feel like finishing the polishing, we will be more than happy to count you as one of the several LAPACK contributors.
If you want more information about LAPACK Contribution, please take a look at http://www.netlib.org/lapack/#_strong_c ... ors_strong and especially the LAPACK Program Style document (Sorry, it has not been updated for a while).
Let us know if you have any questions
Julie

[sven]

Dear Nicolas,

Just a small point, but sqrt(1-norm**2) can be computed stably as sqrt((1+norm)*(1-norm)).

I agree with Julie; updating by rank k instead of rank one would allow the use of Householder transformations and block partitioning and hence the use of Level 3 BLAS. But of course we appreciate that this would be much more work. Daniel Kressner has software in this spirit for updating the QR factorization at

http://www.math.ethz.ch/~kressner/qrupdate.php

Best wishes,

Sven Hammarling.

[nice44]

Julie, yes the names would have to be changed.
I have some tests for the real/double cases but haven't any complex version available.
Also, I didn't start thinking about blocking.
I needed a rank-1 update in the real case, the version from Seeger was not available, so I decided to share this little piece.
My added value was minimal. I might have a deeper look if there is any interest, but I don't promise anything.

Sven, thanks for the sqrt(1-n**2) tips, that's dead simple!
And yes, of course, updates of rank n can be performed with QR (or its LQ transpose in case of lower triangular factors I think).
Do you know how to handle a rank n "downdate"?

Nicolas

[sven]

Dear Nicolas,

The downdating problem needs care. Some ideas, in the context of QR are given in:

[url]eprints.ma.man.ac.uk/1192/[/url]

The best reference for much of the background is probably

A Bjorck. Numerical Methods for Least Squares Problems. SIAM, 1996.

Best wishes,

Sven.

[nice44]

Sven, thanks, this is a much more advanced analysis!
What are you plans about developed code?