I tested the new version and it seems to be consistent. So the current behavior is that if you feed lu_backsub with a 1D n-vector you consistently get out a 2D nX1 row matrix-vector, the same as if the input vector got an additional dummy dimension. Is this the expected behavior?
Best regards, Luis On Mon, May 31, 2021 at 09:55:30PM +0000, Ed . wrote: > Hi Luis, > > The difference you’re seeing is due to a sort of mirage; when you transpose > both outputs (suitable for matrix-multiplication) you’ll get the correct > results, which are a 1x2 column vector. > > However, the lu_backsub code was quite complicated and in ascertaining the > above, and also making it work “inplace” (did you know that the current > “copy” method when “inplace” just turns off “inplace” and returns the given > ndarray?) with tests, I had to make some fixes and changes. Now the dims of > the results will all be the same. > > As soon as I’m confident these changes haven’t broken any other stuff, I’ll > be releasing them. Thanks for the report! > > Best regards, > Ed > > From: Luis Mochan<mailto:[email protected]> > Sent: 29 May 2021 20:28 > To: perldl<mailto:[email protected]>; > perldl<mailto:[email protected]> > Subject: [Pdl-devel] lu_backsub > > 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 | [email protected] /\_/\__/ > GPG: 791EB9EB, C949 3F81 6D9B 1191 9A16 C2DF 5F0A C52B 791E B9EB > > > _______________________________________________ > pdl-devel mailing list > [email protected] > https://lists.sourceforge.net/lists/listinfo/pdl-devel > -- 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 | [email protected] /\_/\__/ GPG: 791EB9EB, C949 3F81 6D9B 1191 9A16 C2DF 5F0A C52B 791E B9EB _______________________________________________ pdl-general mailing list [email protected] https://lists.sourceforge.net/lists/listinfo/pdl-general
