MatrixOps.pm (PDL-2.082) | : | MatrixOps.pm (PDL-2.083) | ||
---|---|---|---|---|
skipping to change at line 17 | skipping to change at line 17 | |||
our %EXPORT_TAGS = (Func=>\@EXPORT_OK); | our %EXPORT_TAGS = (Func=>\@EXPORT_OK); | |||
use PDL::Core; | use PDL::Core; | |||
use PDL::Exporter; | use PDL::Exporter; | |||
use DynaLoader; | use DynaLoader; | |||
our @ISA = ( 'PDL::Exporter','DynaLoader' ); | our @ISA = ( 'PDL::Exporter','DynaLoader' ); | |||
push @PDL::Core::PP, __PACKAGE__; | push @PDL::Core::PP, __PACKAGE__; | |||
bootstrap PDL::MatrixOps ; | bootstrap PDL::MatrixOps ; | |||
#line 8 "matrixops.pd" | #line 9 "matrixops.pd" | |||
use strict; | use strict; | |||
use warnings; | use warnings; | |||
=head1 NAME | =head1 NAME | |||
PDL::MatrixOps -- Some Useful Matrix Operations | PDL::MatrixOps -- Some Useful Matrix Operations | |||
=head1 SYNOPSIS | =head1 SYNOPSIS | |||
skipping to change at line 131 | skipping to change at line 131 | |||
=cut | =cut | |||
use Carp; | use Carp; | |||
use strict; | use strict; | |||
#line 134 "MatrixOps.pm" | #line 134 "MatrixOps.pm" | |||
=head1 FUNCTIONS | =head1 FUNCTIONS | |||
=cut | =cut | |||
#line 123 "matrixops.pd" | #line 124 "matrixops.pd" | |||
=head2 identity | =head2 identity | |||
=for sig | =for sig | |||
Signature: (n; [o]a(n,n)) | Signature: (n; [o]a(n,n)) | |||
=for ref | =for ref | |||
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 | |||
skipping to change at line 161 | skipping to change at line 161 | |||
!(my $was_pdl = UNIVERSAL::isa($n,'PDL')) ? zeroes($n,$n) : | !(my $was_pdl = UNIVERSAL::isa($n,'PDL')) ? zeroes($n,$n) : | |||
$n->getndims == 0 ? zeroes($n->type, $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($n->type, @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" | |||
$was_pdl ? bless $out, ref($n) : $out; | $was_pdl ? bless $out, ref($n) : $out; | |||
} | } | |||
#line 178 "MatrixOps.pm" | ||||
#line 156 "matrixops.pd" | ||||
#line 157 "matrixops.pd" | ||||
=head2 stretcher | =head2 stretcher | |||
=for sig | =for sig | |||
Signature: (a(n); [o]b(n,n)) | Signature: (a(n); [o]b(n,n)) | |||
=for usage | =for usage | |||
$mat = stretcher($eigenvalues); | $mat = stretcher($eigenvalues); | |||
skipping to change at line 188 | skipping to change at line 186 | |||
=cut | =cut | |||
sub stretcher { | sub stretcher { | |||
my $in = shift; | my $in = shift; | |||
my $out = zeroes($in->dim(0),$in->dims); | my $out = zeroes($in->dim(0),$in->dims); | |||
my $tmp; # work around for perl -d "feature" | my $tmp; # work around for perl -d "feature" | |||
($tmp = $out->diagonal(0,1)) += $in; | ($tmp = $out->diagonal(0,1)) += $in; | |||
$out; | $out; | |||
} | } | |||
#line 208 "MatrixOps.pm" | ||||
#line 187 "matrixops.pd" | ||||
#line 188 "matrixops.pd" | ||||
=head2 inv | =head2 inv | |||
=for sig | =for sig | |||
Signature: (a(m,m); sv opt ) | Signature: (a(m,m); sv opt ) | |||
=for usage | =for usage | |||
$a1 = inv($a, {$opt}); | $a1 = inv($a, {$opt}); | |||
skipping to change at line 283 | skipping to change at line 279 | |||
} | } | |||
my $out = lu_backsub($lu,$perm,$par,identity($x))->transpose->sever; | my $out = lu_backsub($lu,$perm,$par,identity($x))->transpose->sever; | |||
return $out | return $out | |||
unless($x->is_inplace); | unless($x->is_inplace); | |||
$x .= $out; | $x .= $out; | |||
$x; | $x; | |||
} | } | |||
#line 307 "MatrixOps.pm" | ||||
#line 287 "matrixops.pd" | ||||
#line 288 "matrixops.pd" | ||||
=head2 det | =head2 det | |||
=for sig | =for sig | |||
Signature: (a(m,m); sv opt) | Signature: (a(m,m); sv opt) | |||
=for usage | =for usage | |||
$det = det($a,{opt}); | $det = det($a,{opt}); | |||
skipping to change at line 338 | skipping to change at line 332 | |||
my($lu,$perm,$par); | my($lu,$perm,$par); | |||
if(exists ($opt->{lu}) and (ref $opt->{lu} eq 'ARRAY')) { | if(exists ($opt->{lu}) and (ref $opt->{lu} eq 'ARRAY')) { | |||
($lu,$perm,$par) = @{$opt->{lu}}; | ($lu,$perm,$par) = @{$opt->{lu}}; | |||
} else { | } else { | |||
($lu,$perm,$par) = lu_decomp($x); | ($lu,$perm,$par) = lu_decomp($x); | |||
$opt->{lu} = [$lu,$perm,$par] | $opt->{lu} = [$lu,$perm,$par] | |||
if(exists($opt->{lu})); | if(exists($opt->{lu})); | |||
} | } | |||
defined $lu ? $lu->diagonal(0,1)->prodover * $par : PDL->zeroes(sbyte,1); | defined $lu ? $lu->diagonal(0,1)->prodover * $par : PDL->zeroes(sbyte,1); | |||
} | } | |||
#line 365 "MatrixOps.pm" | ||||
#line 346 "matrixops.pd" | ||||
#line 347 "matrixops.pd" | ||||
=head2 determinant | =head2 determinant | |||
=for sig | =for sig | |||
Signature: (a(m,m)) | Signature: (a(m,m)) | |||
=for usage | =for usage | |||
$det = determinant($x); | $det = determinant($x); | |||
skipping to change at line 422 | skipping to change at line 414 | |||
determinant($x->slice("0:".($i-1).",1:-1")-> | determinant($x->slice("0:".($i-1).",1:-1")-> | |||
append($x->slice(($i+1).":-1,1:-1"))); | append($x->slice(($i+1).":-1,1:-1"))); | |||
} | } | |||
# Do beginning and end submatrices | # Do beginning and end submatrices | |||
$sum += $x->slice("(0),(0)") * determinant($x->slice('1:-1,1:-1')); | $sum += $x->slice("(0),(0)") * determinant($x->slice('1:-1,1:-1')); | |||
$sum -= $x->slice("(-1),(0)") * determinant($x->slice('0:-2,1:-1')) * (1 - 2*( $n % 2)); | $sum -= $x->slice("(-1),(0)") * determinant($x->slice('0:-2,1:-1')) * (1 - 2*( $n % 2)); | |||
return $sum; | return $sum; | |||
} | } | |||
#line 452 "MatrixOps.pm" | #line 427 "MatrixOps.pm" | |||
#line 958 "../../blib/lib/PDL/PP.pm" | ||||
=head2 eigens_sym | =head2 eigens_sym | |||
=for sig | =for sig | |||
Signature: ([phys]a(m); [o,phys]ev(n,n); [o,phys]e(n)) | Signature: ([phys]a(m); [o,phys]ev(n,n); [o,phys]e(n)) | |||
=for ref | =for ref | |||
Eigenvalues and -vectors of a symmetric square matrix. If passed | Eigenvalues and -vectors of a symmetric square matrix. If passed | |||
skipping to change at line 468 | skipping to change at line 458 | |||
($ev, $e) = eigens_sym($x); # e-vects & e-values | ($ev, $e) = eigens_sym($x); # e-vects & e-values | |||
$e = eigens_sym($x); # just eigenvalues | $e = eigens_sym($x); # just eigenvalues | |||
=for bad | =for bad | |||
eigens_sym ignores the bad-value flag of the input ndarrays. | eigens_sym ignores the bad-value flag of the input ndarrays. | |||
It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. | It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. | |||
=cut | =cut | |||
#line 506 "MatrixOps.pm" | ||||
#line 959 "../../blib/lib/PDL/PP.pm" | ||||
sub PDL::eigens_sym { | sub PDL::eigens_sym { | |||
my ($x) = @_; | my ($x) = @_; | |||
my (@d) = $x->dims; | my (@d) = $x->dims; | |||
barf "Need real square matrix for eigens_sym" | barf "Need real square matrix for eigens_sym" | |||
if $#d < 1 or $d[0] != $d[1]; | if $#d < 1 or $d[0] != $d[1]; | |||
my ($n) = $d[0]; | my ($n) = $d[0]; | |||
my ($sym) = 0.5*($x + $x->transpose); | my ($sym) = 0.5*($x + $x->transpose); | |||
my ($err) = PDL::max(abs($sym)); | my ($err) = PDL::max(abs($sym)); | |||
barf "Need symmetric component non-zero for eigens_sym" | barf "Need symmetric component non-zero for eigens_sym" | |||
skipping to change at line 504 | skipping to change at line 491 | |||
)->copy; | )->copy; | |||
my $ev = PDL->zeroes($sym->dims); | my $ev = PDL->zeroes($sym->dims); | |||
my $e = PDL->zeroes($sym->index(0)->dims); | my $e = PDL->zeroes($sym->index(0)->dims); | |||
&PDL::_eigens_sym_int($lt, $ev, $e); | &PDL::_eigens_sym_int($lt, $ev, $e); | |||
return $ev->transpose, $e | return $ev->transpose, $e | |||
if(wantarray); | if(wantarray); | |||
$e; #just eigenvalues | $e; #just eigenvalues | |||
} | } | |||
#line 545 "MatrixOps.pm" | ||||
#line 960 "../../blib/lib/PDL/PP.pm" | ||||
*eigens_sym = \&PDL::eigens_sym; | *eigens_sym = \&PDL::eigens_sym; | |||
#line 552 "MatrixOps.pm" | ||||
#line 958 "../../blib/lib/PDL/PP.pm" | ||||
=head2 eigens | =head2 eigens | |||
=for sig | =for sig | |||
Signature: ([phys]a(m); [o,phys]ev(l,n,n); [o,phys]e(l,n)) | Signature: ([phys]a(m); [o,phys]ev(l,n,n); [o,phys]e(l,n)) | |||
=for ref | =for ref | |||
Real eigenvalues and -vectors of a real square matrix. | Real eigenvalues and -vectors of a real square matrix. | |||
skipping to change at line 581 | skipping to change at line 562 | |||
($ev, $e) = eigens($x); # e'vects & e'vals | ($ev, $e) = eigens($x); # e'vects & e'vals | |||
$e = eigens($x); # just eigenvalues | $e = eigens($x); # just eigenvalues | |||
=for bad | =for bad | |||
eigens ignores the bad-value flag of the input ndarrays. | eigens ignores the bad-value flag of the input ndarrays. | |||
It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. | It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. | |||
=cut | =cut | |||
#line 631 "MatrixOps.pm" | ||||
#line 959 "../../blib/lib/PDL/PP.pm" | ||||
sub PDL::eigens { | sub PDL::eigens { | |||
my ($x) = @_; | my ($x) = @_; | |||
my (@d) = $x->dims; | my (@d) = $x->dims; | |||
my $n = $d[0]; | my $n = $d[0]; | |||
barf "Need real square matrix for eigens" | barf "Need real square matrix for eigens" | |||
if $#d < 1 or $d[0] != $d[1]; | if $#d < 1 or $d[0] != $d[1]; | |||
my $deviation = PDL::max(abs($x - $x->transpose))/PDL::max(abs($x)); | my $deviation = PDL::max(abs($x - $x->transpose))/PDL::max(abs($x)); | |||
if ( $deviation <= 1e-5 ) { | if ( $deviation <= 1e-5 ) { | |||
#taken from eigens_sym code | #taken from eigens_sym code | |||
skipping to change at line 630 | skipping to change at line 608 | |||
my $ev = PDL->zeroes(2, $x->dims); | my $ev = PDL->zeroes(2, $x->dims); | |||
my $e = PDL->zeroes(2, $x->index(0)->dims); | my $e = PDL->zeroes(2, $x->index(0)->dims); | |||
&PDL::_eigens_int($x->clump(0,1), $ev, $e); | &PDL::_eigens_int($x->clump(0,1), $ev, $e); | |||
return $ev->index(0)->transpose->sever, $e->index(0)->sever | return $ev->index(0)->transpose->sever, $e->index(0)->sever | |||
if(wantarray); | if(wantarray); | |||
return $e->index(0)->sever; #just eigenvalues | return $e->index(0)->sever; #just eigenvalues | |||
} | } | |||
} | } | |||
#line 683 "MatrixOps.pm" | ||||
#line 960 "../../blib/lib/PDL/PP.pm" | ||||
*eigens = \&PDL::eigens; | *eigens = \&PDL::eigens; | |||
#line 690 "MatrixOps.pm" | ||||
#line 958 "../../blib/lib/PDL/PP.pm" | ||||
=head2 svd | =head2 svd | |||
=for sig | =for sig | |||
Signature: (a(n,m); [t]w(wsize); [o]u(n,m); [o,phys]z(n); [o]v(n,n)) | Signature: (a(n,m); [t]w(wsize); [o]u(n,m); [o,phys]z(n); [o]v(n,n)) | |||
=for usage | =for usage | |||
($u, $s, $v) = svd($x); | ($u, $s, $v) = svd($x); | |||
skipping to change at line 699 | skipping to change at line 671 | |||
$r2 *= $s; # implicit broadcasting for cheap mult. | $r2 *= $s; # implicit broadcasting for cheap mult. | |||
$x .= $r2 x $r1; # a gets r2 x ess x r1 | $x .= $r2 x $r1; # a gets r2 x ess x r1 | |||
} | } | |||
=for bad | =for bad | |||
svd ignores the bad-value flag of the input ndarrays. | svd ignores the bad-value flag of the input ndarrays. | |||
It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. | It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. | |||
=cut | =cut | |||
#line 761 "MatrixOps.pm" | ||||
#line 960 "../../blib/lib/PDL/PP.pm" | ||||
*svd = \&PDL::svd; | *svd = \&PDL::svd; | |||
#line 768 "MatrixOps.pm" | ||||
#line 800 "matrixops.pd" | #line 801 "matrixops.pd" | |||
=head2 lu_decomp | =head2 lu_decomp | |||
=for sig | =for sig | |||
Signature: (a(m,m); [o]lu(m,m); [o]perm(m); [o]parity) | Signature: (a(m,m); [o]lu(m,m); [o]perm(m); [o]parity) | |||
=for ref | =for ref | |||
LU decompose a matrix, with row permutation | LU decompose a matrix, with row permutation | |||
skipping to change at line 880 | skipping to change at line 848 | |||
return if !$notbig->isempty; | 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 | |||
wantarray ? ($out,$permute,$parity) : $out; | wantarray ? ($out,$permute,$parity) : $out; | |||
} | } | |||
#line 947 "MatrixOps.pm" | ||||
#line 979 "matrixops.pd" | ||||
#line 980 "matrixops.pd" | ||||
=head2 lu_decomp2 | =head2 lu_decomp2 | |||
=for sig | =for sig | |||
Signature: (a(m,m); [o]lu(m,m)) | Signature: (a(m,m); [o]lu(m,m)) | |||
=for ref | =for ref | |||
LU decompose a matrix, with no row permutation | LU decompose a matrix, with no row permutation | |||
skipping to change at line 998 | skipping to change at line 964 | |||
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 | |||
wantarray ? ($out,$perm,$par) : $out; | wantarray ? ($out,$perm,$par) : $out; | |||
} | } | |||
#line 1069 "MatrixOps.pm" | ||||
#line 1101 "matrixops.pd" | ||||
#line 1102 "matrixops.pd" | ||||
=head2 lu_backsub | =head2 lu_backsub | |||
=for sig | =for sig | |||
Signature: (lu(m,m); perm(m); b(m)) | Signature: (lu(m,m); perm(m); b(m)) | |||
=for ref | =for ref | |||
Solve A x = B for matrix A, by back substitution into A's LU decomposition. | Solve A x = B for matrix A, by back substitution into A's LU decomposition. | |||
skipping to change at line 1213 | skipping to change at line 1177 | |||
)->sumover; | )->sumover; | |||
($tmp = $out->index($r1)) /= $ludiag->index($r1)->dummy(0); # TODO: check broadcast dims | ($tmp = $out->index($r1)) /= $ludiag->index($r1)->dummy(0); # TODO: check broadcast dims | |||
} | } | |||
if ($y->is_inplace) { | if ($y->is_inplace) { | |||
$y->setdims([$out->dims]) if !PDL::all($y->shape == $out->shape); # assgn n eeds same shape | $y->setdims([$out->dims]) if !PDL::all($y->shape == $out->shape); # assgn n eeds same shape | |||
$y .= $out; | $y .= $out; | |||
} | } | |||
$out; | $out; | |||
} | } | |||
#line 1287 "MatrixOps.pm" | #line 1219 "MatrixOps.pm" | |||
#line 958 "../../blib/lib/PDL/PP.pm" | ||||
=head2 simq | =head2 simq | |||
=for sig | =for sig | |||
Signature: ([phys]a(n,n); [phys]b(n); [o,phys]x(n); int [o,phys]ips(n); int fl ag) | Signature: ([phys]a(n,n); [phys]b(n); [o,phys]x(n); int [o,phys]ips(n); int fl ag) | |||
=for ref | =for ref | |||
Solution of simultaneous linear equations, C<a x = b>. | Solution of simultaneous linear equations, C<a x = b>. | |||
skipping to change at line 1251 | skipping to change at line 1213 | |||
See also L</lu_backsub>, which does the same thing with a slightly | See also L</lu_backsub>, which does the same thing with a slightly | |||
less opaque interface. | less opaque interface. | |||
=for bad | =for bad | |||
simq ignores the bad-value flag of the input ndarrays. | simq ignores the bad-value flag of the input ndarrays. | |||
It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. | It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. | |||
=cut | =cut | |||
#line 1333 "MatrixOps.pm" | ||||
#line 960 "../../blib/lib/PDL/PP.pm" | ||||
*simq = \&PDL::simq; | *simq = \&PDL::simq; | |||
#line 1340 "MatrixOps.pm" | ||||
#line 958 "../../blib/lib/PDL/PP.pm" | ||||
=head2 squaretotri | =head2 squaretotri | |||
=for sig | =for sig | |||
Signature: (a(n,n); b(m)) | Signature: (a(n,n); [o]b(m)) | |||
=for ref | =for ref | |||
Convert a symmetric square matrix to triangular vector storage. | Convert a lower-triangular square matrix to triangular vector storage. | |||
Ignores upper half of input. | ||||
=for bad | =for bad | |||
squaretotri does not process bad values. | squaretotri does not process bad values. | |||
It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. | It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. | |||
=cut | =cut | |||
#line 1368 "MatrixOps.pm" | ||||
#line 960 "../../blib/lib/PDL/PP.pm" | ||||
*squaretotri = \&PDL::squaretotri; | *squaretotri = \&PDL::squaretotri; | |||
#line 1375 "MatrixOps.pm" | ||||
#line 1396 "matrixops.pd" | #line 1390 "matrixops.pd" | |||
=head1 AUTHOR | =head1 AUTHOR | |||
Copyright (C) 2002 Craig DeForest (deforest@boulder.swri.edu), | Copyright (C) 2002 Craig DeForest (deforest@boulder.swri.edu), | |||
R.J.R. Williams (rjrw@ast.leeds.ac.uk), Karl Glazebrook | R.J.R. Williams (rjrw@ast.leeds.ac.uk), Karl Glazebrook | |||
(kgb@aaoepp.aao.gov.au). There is no warranty. You are allowed to | (kgb@aaoepp.aao.gov.au). There is no warranty. You are allowed to | |||
redistribute and/or modify this work under the same conditions as PDL | redistribute and/or modify this work under the same conditions as PDL | |||
itself. If this file is separated from the PDL distribution, then the | itself. If this file is separated from the PDL distribution, then the | |||
PDL copyright notice should be included in this file. | PDL copyright notice should be included in this file. | |||
=cut | =cut | |||
#line 1394 "MatrixOps.pm" | #line 1307 "MatrixOps.pm" | |||
# Exit with OK status | # Exit with OK status | |||
1; | 1; | |||
End of changes. 33 change blocks. | ||||
63 lines changed or deleted | 16 lines changed or added |