"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "Basic/MatrixOps/matrixops.pd" between
PDL-2.080.tar.gz and PDL-2.081.tar.gz

About: PDL (Perl Data Language) aims to turn perl into an efficient numerical language for scientific computing (similar to IDL and MatLab).

matrixops.pd  (PDL-2.080):matrixops.pd  (PDL-2.081)
skipping to change at line 140 skipping to change at line 140
Return an identity matrix of the specified size. If you hand in a Return an identity matrix of the specified size. If you hand in a
scalar, its value is the size of the identity matrix; if you hand in a scalar, its value is the size of the identity matrix; if you hand in a
dimensioned PDL, the 0th dimension is the first two dimensions of the dimensioned PDL, the 0th dimension is the first two dimensions of the
matrix, with higher dimensions preserved. matrix, with higher dimensions preserved.
=cut =cut
sub identity { sub identity {
my $n = shift; my $n = shift;
my $out = my $out =
!UNIVERSAL::isa($n,'PDL') ? zeroes($n,$n) : !(my $was_pdl = UNIVERSAL::isa($n,'PDL')) ? zeroes($n,$n) :
$n->getndims == 0 ? zeroes($n->at(0),$n->at(0)) : $n->getndims == 0 ? zeroes($n->type, $n->at(0),$n->at(0)) :
undef; undef;
if (!defined $out) { if (!defined $out) {
my @dims = $n->dims; my @dims = $n->dims;
$out = zeroes(@dims[0, 0, 2..$#dims]); $out = zeroes($n->type, @dims[0, 0, 2..$#dims]);
} }
(my $tmp = $out->diagonal(0,1))++; # work around perl -d "feature" (my $tmp = $out->diagonal(0,1))++; # work around perl -d "feature"
$out; $was_pdl ? bless $out, ref($n) : $out;
} }
EOD EOD
###################################################################### ######################################################################
pp_add_exported('','stretcher'); pp_add_exported('','stretcher');
pp_addpm(<<'EOD'); pp_addpm(<<'EOD');
=head2 stretcher =head2 stretcher
=for sig =for sig
skipping to change at line 926 skipping to change at line 926
$parity = $in->slice('(0),(0)')->ones; $parity = $in->slice('(0),(0)')->ones;
} }
my($scales) = $in->copy->abs->maximum; # elementwise by rows my($scales) = $in->copy->abs->maximum; # elementwise by rows
if(($scales==0)->sum) { if(($scales==0)->sum) {
return undef; return undef;
} }
# Some holding tanks # Some holding tanks
my($tmprow) = $out->slice('(0)')->double->zeroes; my($tmprow) = $out->slice('(0)')->zeroes;
$tmprow = $tmprow->double if $tmprow->type < double;
my($tmpval) = $tmprow->slice('(0)')->sever; my($tmpval) = $tmprow->slice('(0)')->sever;
my($col,$row); my($col,$row);
for $col(0..$n1) { for $col(0..$n1) {
for $row(1..$n1) { for $row(1..$n1) {
my($klim) = $row<$col ? $row : $col; my($klim) = $row<$col ? $row : $col;
if($klim > 0) { if($klim > 0) {
$klim--; $klim--;
my($el) = $out->index2d($col,$row); my($el) = $out->index2d($col,$row);
$el -= ( $out->slice("($col),0:$klim") * $el -= ( $out->slice("($col),0:$klim") *
skipping to change at line 968 skipping to change at line 969
$sl1 = $permute->index($whc); $sl1 = $permute->index($whc);
$sl2 = $permute->index($col); $sl2 = $permute->index($col);
$tmpval .= $sl1; $sl1 .= $sl2; $sl2 .= $tmpval; $tmpval .= $sl1; $sl1 .= $sl2; $sl2 .= $tmpval;
{ my $tmp; { my $tmp;
($tmp = $parity->where($wh>0)) *= -1.0; ($tmp = $parity->where($wh>0)) *= -1.0;
} }
} }
# Sidestep near-singularity (NR does this; not sure if it is helpful) # LAPACK cgetrf does not try fix singularity so nor do we, even though
NR does
my $notbig = $big->where(abs($big) < $TINY); my $notbig = $big->where(abs($big) < $TINY);
$notbig .= $TINY * (1.0 - 2.0*($notbig < 0)); return if !$notbig->isempty;
# Divide by the diagonal element (which is now the largest element) # Divide by the diagonal element (which is now the largest element)
my $tout; my $tout;
($tout = $out->slice("($col),".($col+1).":$n1")) /= $big->slice('*1'); ($tout = $out->slice("($col),".($col+1).":$n1")) /= $big->slice('*1');
} # end of pivoting part } # end of pivoting part
} # end of column loop } # end of column loop
if(wantarray) { wantarray ? ($out,$permute,$parity) : $out;
return ($out,$permute,$parity);
}
$out;
} }
EOD EOD
###################################################################### ######################################################################
pp_add_exported('','lu_decomp2'); pp_add_exported('','lu_decomp2');
pp_addpm(<<'EOD'); pp_addpm(<<'EOD');
=head2 lu_decomp2 =head2 lu_decomp2
skipping to change at line 1103 skipping to change at line 1100
# Figure a_ij, with no pivoting # Figure a_ij, with no pivoting
if($col < $n1) { if($col < $n1) {
# Divide the rest of the column by the diagonal element # Divide the rest of the column by the diagonal element
my $tmp; # work around for perl -d "feature" my $tmp; # work around for perl -d "feature"
($tmp = $out->slice("($col),".($col+1).":$n1")) /= $diagonal->index($col)- >dummy(0,$n1-$col); ($tmp = $out->slice("($col),".($col+1).":$n1")) /= $diagonal->index($col)- >dummy(0,$n1-$col);
} }
} # end of column loop } # end of column loop
if(wantarray) { wantarray ? ($out,$perm,$par) : $out;
return ($out,$perm,$par);
}
$out;
} }
EOD EOD
###################################################################### ######################################################################
pp_add_exported('','lu_backsub'); pp_add_exported('','lu_backsub');
pp_addpm(<<'EOD'); pp_addpm(<<'EOD');
=head2 lu_backsub =head2 lu_backsub
 End of changes. 8 change blocks. 
16 lines changed or deleted 11 lines changed or added

Home  |  About  |  Features  |  All  |  Newest  |  Dox  |  Diffs  |  RSS Feeds  |  Screenshots  |  Comments  |  Imprint  |  Privacy  |  HTTP(S)