That is what I was aiming to achieve (and have captured in hopefully-thorough tests). My thinking is that the linear algebra stuff is founded on matrices, with at least 2 dimensions. Therefore, getting something back that has at least 2 dimensions (and always at least 2) seems to me the least surprising outcome (especially given your previous note was apparently prompted by surprise).
What do you (and others) think? Best regards, Ed From: Luis Mochan<mailto:[email protected]> Sent: 01 June 2021 01:46 To: perldl<mailto:[email protected]> Subject: Re: [Pdl-general] [Pdl-devel] lu_backsub 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
_______________________________________________ pdl-general mailing list [email protected] https://lists.sourceforge.net/lists/listinfo/pdl-general
