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