I found a strange behavior in lu_backsub from PDL::MatrixOps:
pdl> my $M=pdl([1,2],[3,4]);
pdl> my $M1=pdl([3,4],[1,2]); # interchange two rows
pdl> my $y=pdl(1,1);
pdl> my $x=lu_backsub($M->copy->lu_decomp, $y);
pdl> my $x1=lu_backsub($M1->copy->lu_decomp, $y);
pdl> p $x->info;
PDL: Double D [2,1]
pdl> p $x1->info;
PDL: Double D [2]
In both cases I'm solving the same 2x2 system of equations using a 1D
vector as the right hand side, but in the second case I interchanged
rows, so I expected the same solution. In the first case I get back a 'row
vector', i.e., a matrix with one row. However, in the second case I
got a real 1D vector, not a matrix.
My guess is that the problem is related to the code around line 1298
of matrixops.pd, which contains some warnings:
if($nontrivial) {
if($y->is_inplace) {
$y .= $y->dummy(1,$y->dim(0))->index($perm->dummy(1,1))->sever; #
TODO: check threading
$out = $y;
} else {
$out = $y->dummy(1,$y->dim(0))->index($perm->dummy(1,1))->sever; #
TODO: check threading
}
} else {
# should check for more matrix dims to thread over
# but ignore the issue for now
$out = ($y->is_inplace ? $y : $y->copy);
}
print STDERR "lu_backsub: starting with \$out" . pdl($out->dims) . "\n"
if $PDL::debug;
Best regards,
Luis
--
o
W. Luis Mochán, | tel:(52)(777)329-1734 /<(*)
Instituto de Ciencias Físicas, UNAM | fax:(52)(777)317-5388 `>/ /\
Av. Universidad s/n CP 62210 | (*)/\/ \
Cuernavaca, Morelos, México | moc...@fis.unam.mx /\_/\__/
GPG: 791EB9EB, C949 3F81 6D9B 1191 9A16 C2DF 5F0A C52B 791E B9EB
_______________________________________________
pdl-general mailing list
pdl-general@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/pdl-general
