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 |