Following the spirited discussion of dummy() last week and the report in sf.net bug #3021567 where lu_backsub is returning an answer with too many dimensions:
http://sourceforge.net/tracker/?func=detail&aid=3021567&group_id=612&atid=100612 I would like to suggest this as a good exercise for a practical analysis of the use of dimensions, dummy and threading in one of the PDL Core routines. I would like to challenge any interested students of the PDL way to see if they can work through the code of the lu_backsub method to discover where the problem is that introduces the undesired extra dimensions in the output. HINT: working through a call to the lu_backsub method a line at a time using the perldl or pdl2 shell and determining where things act differently than you would expect could be one way to approach the problem. Happy PDL-ing! Chris P.S. I've attached the lu_backsub source from MatrixOps.pm to the above bug ticket and I'm including it below for those who can't wait that long... sub lu_backsub { my ($lu, $perm, $b, $par); if(@_==3) { ($lu, $perm, $b) = @_; } elsif(@_==4) { ($lu, $perm, $par, $b) = @_; } barf("lu_backsub: LU decomposition is undef -- probably from a singular matrix.\n") unless defined($lu); barf("Usage: \$x = lu_backsub(\$lu,\$perm,\$b); all must be PDLs\n") unless(UNIVERSAL::isa($lu,'PDL') && UNIVERSAL::isa($perm,'PDL') && UNIVERSAL::isa($b,'PDL')); my $n = $b->dim(0); my $n1 = $n; $n1--; # Permute the vector and make a copy if necessary. my $out; my $nontrivial = !(($perm==(PDL->xvals($perm->dims)))->all); if($nontrivial) { if($b->is_inplace) { $b .= $b->dummy(1,$b->dim(0))->index($perm)->sever; $out = $b; } else { $out = $b->dummy(1,$b->dim(0))->index($perm)->sever; } } else { $out = ($b->is_inplace ? $b : $b->copy); } # Make sure threading over lu happens OK... if($out->ndims < $lu->ndims) { do { $out = $out->dummy(-1,$lu->dim($out->ndims)); } while($out->ndims < $lu->ndims); $out = $out->sever; } ## Do forward substitution into L my $row; my $r1; for $row(1..$n1) { $r1 = $row-1; my $tmp; # work around perl -d "feature ($tmp = $out->index($row)) -= ($lu->(0:$r1,$row) * $out->(0:$r1) )->sumover; } ## Do backward substitution into U, and normalize by the diagonal my $ludiag = $lu->diagonal(0,1)->dummy(1,$n); { my $tmp; # work around for perl -d "feature" ($tmp = $out->index($n1)) /= $ludiag->index($n1); } for ($row=$n1; $row>0; $row--) { $r1 = $row-1; my $tmp; # work around for perl -d "feature" ($tmp = $out->index($r1)) -= ($lu->($row:$n1,$r1) * $out->($row:$n1) )->sumover; ($tmp = $out->index($r1)) /= $ludiag->index($r1); } $out; } _______________________________________________ Perldl mailing list [email protected] http://mailman.jach.hawaii.edu/mailman/listinfo/perldl
