The following programs fails to solve a simple 2x2 system of equations; I believe something is wrong with either lu_decomp and/or lu_backsub from PDL::MatrixOps. (or am I making a something wrong?). Maybe it is related with amd64.
Program:
#!/usr/bin/perl -w
use strict;
use PDL;
use PDL::NiceSlice;
my $m=pdl([[2,1],[1,2]]); #make a simple matrix
my ($lu,$perm,$par)=lu_decomp($m); #decompose
my $x=pdl([[1],[0]]); #make a simple column vector
#back substitute directly and transpossing
my $y=lu_backsub($lu,$perm,$par,$x);
my $y1=lu_backsub($lu,$perm,$par,$x->transpose)->transpose;
print "Solve \$m x \$y= \$x\n";
print "Original matrix times result: ", $m x $y;
print "Original matrix times alternative result: ", $m x $y1;
print "One of the results should have yielded the original \$x: $x";
Output:
Solve $m x $y= $x
Original matrix times result:
[
[ 1]
[0.5]
]
Original matrix times alternative result:
[
[ 0.38888889]
[-0.13888889]
]
One of the results should have yielded the original $x:
[
[1]
[0]
]
The output from perldl -V is attached
W_Luis
2008-07-21
output from perldl -V
Chris Marshall
2008-07-27
Logged In: YES
user_id=44920
Originator: NO
The signature of the lu_backsub() routine is:
(lu(m,m); perm(m); b(m))
Which means that the calling args for the two tests
should be
(lu(2,2); perm(2); x(1,2))
and
(lu(2,2); perm(2); x(2,1))
neither of which matches the required signature.
I tried using $x = pdl([1,0]) instead with call
of
(lu(2,2); perm(2); x(2))
which does match and got for the result
perldl> $y = lu_backsub($lu,$perm,$d,$x)
perldl> p $y
[
[ 0.66666667 -0.33333333]
[ 0.66666667 -0.33333333]
]
perldl> p $y x $m
[
[1 0]
[1 0]
]
which does give the correct solution
although it looks like there is a problem
with the threading.
I tried padding the $lu and $perm
arrays with extra singleton dimension
which should trivially thread and got
the incorrect answer in case #2 above.
When I used
$lu = $lu->dummy(-1,2)-sever
$perm = $perm->dummy(-1,2)->sever
$x = $x->dummy(-1,2)->sever
then I got the repeated "correct"
answer as above (i.e. it has two
identical solutions but the solutions
are 2x2 rather than 2x1).
SUMMARY: the lu_decomp() appears to work correctly.
The lu_backsub() gives correct (but redundant) solns
if the arguments have the correct dimensionality.
There may be a problem with trivial threading or
the threading in general. Needs more looking at...
Chris Marshall
2008-07-27
lu_backsub() pp implementation and test
Chris Marshall
2008-07-27
Logged In: YES
user_id=44920
Originator: NO
I implemented a PDL::PP version of the NR in C back substitution code (see attached Inline::Pdlpp code) which appeared to work for the 2x2 problem of this bug report. When I tried it on a problem with a nontrivial permutation, it does not work. No more time to work on this but maybe this can help someone finish this off. The PP version of the code, lu_backsub2(), seems to handle the loop/thread dimensions but something is breaking when a permutation is obtained in the lu_decomp().
File Added: lu_back.pl
Chris Marshall
2008-07-29
Logged In: YES
user_id=44920
Originator: NO
I took another look at the NR C code I translated into a PP lu_backup2() routine. The LU decomp still looks solid. The problem remaining appears to be the back substitution. I think it may be tied into row/col-major issues and dimension ordering---something that can be tricky to keep straight at times.
Chris Marshall
2008-09-18
Another alternative here would be to implement the GSL versions of these functions. That might be a better alternative in the long run.
W_Luis
2008-09-18
I actually solved my problem using the Slatec routines (geco and gesl), but the bug in MatrixOps still worries me.
Chris Marshall
2008-10-06
This needs to be resolved for PDL-2.4.4 if possible.
Chris Marshall
2008-10-06
Chris Marshall
2008-10-13
Bug fixed in CVS.
Thanks for reporting the problem!