From: Gunter Winkler (guwi17_at_[hidden])
Date: 2006-10-24 14:25:56
Georg Baum schrieb:
> Am Montag, 23. Oktober 2006 20:40 schrieb Sorkin Dima:
>> 1) Is there any docs/tests for lu.hpp ?
> AFAIK none. If some exists they could be in the wiki, and IRC Gunter Winkler
> has an example on his page.
There is no documentation, because LU does exactly what its name
suggests: LU decomposition and LU solves.
>> 2) I saw that in lu.hpp the arguments
>> of lu_substitute can come in either
>> (src,dst) or (dst,src) order. Why two
>> implementations ?
> I did not know that two versions existed! This is a mistake IMO, it only
> creates confusion.
Yes. I think it is because all products and solves are twice implemented:
(x,y ... vectors, A .. matrix)
solve(A,x,y) -> Ax = y
solve(x,A,y) -> A^T x = y
unfortunately, this makes trouble when used with multiple right hand sides:
(X,Y ... matrices (usually column major), A ... matrix)
solve(A,X,Y) -> A X = Y
solve(X,A,Y) -> A^T X = Y
The compiler cannot distinguish both methods!
I suggest renaming the LU solves to
lu_substitute(A,X,Y), lu_substitute(A,x,y) -> solve L(Ux)=y
ul_substitute(A,X,Y), ul_substitute(A,x,y) -> U(Lx)=y
(thus lu_substitute(X,A,Y) would then be ul_substitute(trans(A),X,Y),
which is much clearer)