"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "GENERATED/PDL/MatrixOps.pm" between
PDL-2.074.tar.gz and PDL-2.075.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.pm  (PDL-2.074):MatrixOps.pm  (PDL-2.075)
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 9 "matrixops.pd" #line 8 "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
$inv = $x->inv; $inv = $x->inv;
skipping to change at line 71 skipping to change at line 72
matrix itself, counting rightwards and downwards from the upper left matrix itself, counting rightwards and downwards from the upper left
corner. This means that if you print a PDL that contains a matrix, corner. This means that if you print a PDL that contains a matrix,
the matrix appears correctly on the screen, but if you index a matrix the matrix appears correctly on the screen, but if you index a matrix
element, you use the indices in the reverse order that you would in a element, you use the indices in the reverse order that you would in a
math textbook. If you prefer your matrices indexed in (row, column) math textbook. If you prefer your matrices indexed in (row, column)
order, you can try using the L<PDL::Matrix> object, which order, you can try using the L<PDL::Matrix> object, which
includes an implicit exchange of the first two dimensions but should includes an implicit exchange of the first two dimensions but should
be compatible with most of these matrix operations. TIMTOWDTI.) be compatible with most of these matrix operations. TIMTOWDTI.)
Matrices, row vectors, and column vectors can be multiplied with the 'x' Matrices, row vectors, and column vectors can be multiplied with the 'x'
operator (which is, of course, threadable): operator (which is, of course, broadcastable):
$m3 = $m1 x $m2; $m3 = $m1 x $m2;
$col_vec2 = $m1 x $col_vec1; $col_vec2 = $m1 x $col_vec1;
$row_vec2 = $row_vec1 x $m1; $row_vec2 = $row_vec1 x $m1;
$scalar = $row_vec x $col_vec; $scalar = $row_vec x $col_vec;
Because of the (column,row) addressing order, 1-D PDLs are treated as Because of the (column,row) addressing order, 1-D PDLs are treated as
_row_ vectors; if you want a _column_ vector you must add a dummy dimension: _row_ vectors; if you want a _column_ vector you must add a dummy dimension:
$rowvec = pdl(1,2); # row vector $rowvec = pdl(1,2); # row vector
$colvec = $rowvec->slice('*1'); # 1x2 column vector $colvec = $rowvec->slice('*1'); # 1x2 column vector
$matrix = pdl([[3,4],[6,2]]); # 2x2 matrix $matrix = pdl([[3,4],[6,2]]); # 2x2 matrix
$rowvec2 = $rowvec x $matrix; # right-multiplication by matrix $rowvec2 = $rowvec x $matrix; # right-multiplication by matrix
$colvec = $matrix x $colvec; # left-multiplication by matrix $colvec = $matrix x $colvec; # left-multiplication by matrix
$m2 = $matrix x $rowvec; # Throws an error $m2 = $matrix x $rowvec; # Throws an error
Implicit threading works correctly with most matrix operations, but Implicit broadcasting works correctly with most matrix operations, but
you must be extra careful that you understand the dimensionality. In you must be extra careful that you understand the dimensionality. In
particular, matrix multiplication and other matrix ops need nx1 PDLs particular, matrix multiplication and other matrix ops need nx1 PDLs
as row vectors and 1xn PDLs as column vectors. In most cases you must as row vectors and 1xn PDLs as column vectors. In most cases you must
explicitly include the trailing 'x1' dimension in order to get the expected explicitly include the trailing 'x1' dimension in order to get the expected
results when you thread over multiple row vectors. results when you broadcast over multiple row vectors.
When threading over matrices, it's very easy to get confused about When broadcasting over matrices, it's very easy to get confused about
which dimension goes where. It is useful to include comments with which dimension goes where. It is useful to include comments with
every expression, explaining what you think each dimension means: every expression, explaining what you think each dimension means:
$x = xvals(360)*3.14159/180; # (angle) $x = xvals(360)*3.14159/180; # (angle)
$rot = cat(cat(cos($x),sin($x)), # rotmat: (col,row,angle) $rot = cat(cat(cos($x),sin($x)), # rotmat: (col,row,angle)
cat(-sin($x),cos($x))); cat(-sin($x),cos($x)));
=head1 ACKNOWLEDGEMENTS =head1 ACKNOWLEDGEMENTS
MatrixOps includes algorithms and pre-existing code from several MatrixOps includes algorithms and pre-existing code from several
skipping to change at line 124 skipping to change at line 125
small-to-mid sized matrices. The algorithms may not scale well to small-to-mid sized matrices. The algorithms may not scale well to
large matrices (hundreds by hundreds) or to near singular matrices. large matrices (hundreds by hundreds) or to near singular matrices.
If there is something you want that is not here, please add and If there is something you want that is not here, please add and
document it! document it!
=cut =cut
use Carp; use Carp;
use strict; use strict;
#line 133 "MatrixOps.pm" #line 134 "MatrixOps.pm"
=head1 FUNCTIONS =head1 FUNCTIONS
=cut =cut
#line 124 "matrixops.pd" #line 123 "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
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
skipping to change at line 159 skipping to change at line 161
!UNIVERSAL::isa($n,'PDL') ? zeroes($n,$n) : !UNIVERSAL::isa($n,'PDL') ? zeroes($n,$n) :
$n->getndims == 0 ? zeroes($n->at(0),$n->at(0)) : $n->getndims == 0 ? zeroes($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(@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; $out;
} }
#line 176 "MatrixOps.pm" #line 178 "MatrixOps.pm"
#line 157 "matrixops.pd" #line 156 "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 186 skipping to change at line 188
=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 205 "MatrixOps.pm" #line 208 "MatrixOps.pm"
#line 188 "matrixops.pd" #line 187 "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 281 skipping to change at line 283
} }
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 304 "MatrixOps.pm" #line 308 "MatrixOps.pm"
#line 289 "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 339 skipping to change at line 341
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 : 0 ); ( (defined $lu) ? $lu->diagonal(0,1)->prodover * $par : 0 );
} }
#line 364 "MatrixOps.pm" #line 369 "MatrixOps.pm"
#line 351 "matrixops.pd" #line 350 "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);
=for ref =for ref
Determinant of a square matrix, using recursive descent (threadable). Determinant of a square matrix, using recursive descent (broadcastable).
This is the traditional, robust recursive determinant method taught in This is the traditional, robust recursive determinant method taught in
most linear algebra courses. It scales like C<O(n!)> (and hence is most linear algebra courses. It scales like C<O(n!)> (and hence is
pitifully slow for large matrices) but is very robust because no pitifully slow for large matrices) but is very robust because no
division is involved (hence no division-by-zero errors for singular division is involved (hence no division-by-zero errors for singular
matrices). It's also threadable, so you can find the determinants of matrices). It's also broadcastable, so you can find the determinants of
a large collection of matrices all at once if you want. a large collection of matrices all at once if you want.
Matrices up to 3x3 are handled by direct multiplication; larger matrices Matrices up to 3x3 are handled by direct multiplication; larger matrices
are handled by recursive descent to the 3x3 case. are handled by recursive descent to the 3x3 case.
The LU-decomposition method L</det> is faster in isolation for The LU-decomposition method L</det> is faster in isolation for
single matrices larger than about 4x4, and is much faster if you end up single matrices larger than about 4x4, and is much faster if you end up
reusing the LU decomposition of C<$a> (NOTE: check performance and reusing the LU decomposition of C<$a> (NOTE: check performance and
threading benchmarks with new code). broadcasting benchmarks with new code).
=cut =cut
*PDL::determinant = \&determinant; *PDL::determinant = \&determinant;
sub determinant { sub determinant {
my($x) = shift; my($x) = shift;
my($n); my($n);
return undef unless( return undef unless(
UNIVERSAL::isa($x,'PDL') && UNIVERSAL::isa($x,'PDL') &&
$x->getndims >= 2 && $x->getndims >= 2 &&
($n = $x->dim(0)) == $x->dim(1) ($n = $x->dim(0)) == $x->dim(1)
); );
return $x->clump(2) if($n==1); return $x->clump(2) if($n==1);
if($n==2) { if($n==2) {
my($y) = $x->clump(2); my($y) = $x->clump(2);
return $y->index(0)*$y->index(3) - $y->index(1)*$y->index(2); return $y->index(0)*$y->index(3) - $y->index(1)*$y->index(2);
} }
if($n==3) { if($n==3) {
my($y) = $x->clump(2); my($y) = $x->clump(2);
my $y3 = $y->index(3); my $y3 = $y->index(3);
my $y4 = $y->index(4); my $y4 = $y->index(4);
my $y5 = $y->index(5); my $y5 = $y->index(5);
my $y6 = $y->index(6); my $y6 = $y->index(6);
my $y7 = $y->index(7); my $y7 = $y->index(7);
my $y8 = $y->index(8); my $y8 = $y->index(8);
return ( return (
$y->index(0) * ( $y4 * $y8 - $y5 * $y7 ) $y->index(0) * ( $y4 * $y8 - $y5 * $y7 )
+ $y->index(1) * ( $y5 * $y6 - $y3 * $y8 ) + $y->index(1) * ( $y5 * $y6 - $y3 * $y8 )
skipping to change at line 424 skipping to change at line 425
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 451 "MatrixOps.pm" #line 456 "MatrixOps.pm"
#line 1059 "../../blib/lib/PDL/PP.pm" #line 1058 "../../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
an asymmetric matrix, the routine will warn and symmetrize it, by taking an asymmetric matrix, the routine will warn and symmetrize it, by taking
the average value. That is, it will solve for 0.5*($a+$a->transpose). the average value. That is, it will solve for 0.5*($a+$a->transpose).
It's threadable, so if C<$a> is 3x3x100, it's treated as 100 separate 3x3 It's broadcastable, so if C<$a> is 3x3x100, it's treated as 100 separate 3x3
matrices, and both C<$ev> and C<$e> get extra dimensions accordingly. matrices, and both C<$ev> and C<$e> get extra dimensions accordingly.
If called in scalar context it hands back only the eigenvalues. Ultimately, If called in scalar context it hands back only the eigenvalues. Ultimately,
it should switch to a faster algorithm in this case (as discarding the it should switch to a faster algorithm in this case (as discarding the
eigenvectors is wasteful). eigenvectors is wasteful).
The algorithm used is due to J. vonNeumann, which was a rediscovery of The algorithm used is due to J. vonNeumann, which was a rediscovery of
L<Jacobi's Method|http://en.wikipedia.org/wiki/Jacobi_eigenvalue_algorithm> . L<Jacobi's Method|http://en.wikipedia.org/wiki/Jacobi_eigenvalue_algorithm> .
The eigenvectors are returned in COLUMNS of the returned PDL. That The eigenvectors are returned in COLUMNS of the returned PDL. That
skipping to change at line 470 skipping to change at line 471
($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 504 "MatrixOps.pm" #line 510 "MatrixOps.pm"
#line 1060 "../../blib/lib/PDL/PP.pm" #line 1059 "../../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"
if $err == 0; if $err == 0;
$err = PDL::max(abs($x-$sym))/$err; $err = PDL::max(abs($x-$sym))/$err;
warn "Using symmetrized version of the matrix in eigens_sym" warn "Using symmetrized version of the matrix in eigens_sym"
if $err > 1e-5 && $PDL::debug; if $err > 1e-5 && $PDL::debug;
## Get lower diagonal form ## Get lower diagonal form
## Use whichND/indexND because whereND doesn't exist (yet?) and ## Use whichND/indexND because whereND doesn't exist (yet?) and
## the combo is threadable (unlike where). Note that for historical ## the combo is broadcastable (unlike where). Note that for historical
## reasons whichND needs a scalar() around it to give back a ## reasons whichND needs a scalar() around it to give back a
## nice 2xn PDL index. ## nice 2xn PDL index.
my $lt = PDL::indexND($sym, my $lt = PDL::indexND($sym,
scalar(PDL::whichND(PDL->xvals($n,$n) <= scalar(PDL::whichND(PDL->xvals($n,$n) <=
PDL->yvals($n,$n))) PDL->yvals($n,$n)))
)->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 542 "MatrixOps.pm" #line 549 "MatrixOps.pm"
#line 1060 "../../blib/lib/PDL/PP.pm"
#line 1061 "../../blib/lib/PDL/PP.pm"
*eigens_sym = \&PDL::eigens_sym; *eigens_sym = \&PDL::eigens_sym;
#line 548 "MatrixOps.pm" #line 556 "MatrixOps.pm"
#line 1059 "../../blib/lib/PDL/PP.pm" #line 1058 "../../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 550 skipping to change at line 552
experimental! For asymmetric matrices, nearly all observed matrices experimental! For asymmetric matrices, nearly all observed matrices
with real eigenvalues produce incorrect results, due to errors of the with real eigenvalues produce incorrect results, due to errors of the
sslib algorithm. If your assymmetric matrix returns all NaNs, do not sslib algorithm. If your assymmetric matrix returns all NaNs, do not
assume that the values are complex. Also, problems with memory access assume that the values are complex. Also, problems with memory access
is known in this library. is known in this library.
Not all square matrices are diagonalizable. If you feed in a Not all square matrices are diagonalizable. If you feed in a
non-diagonalizable matrix, then one or more of the eigenvectors will non-diagonalizable matrix, then one or more of the eigenvectors will
be set to NaN, along with the corresponding eigenvalues. be set to NaN, along with the corresponding eigenvalues.
C<eigens> is threadable, so you can solve 100 eigenproblems by C<eigens> is broadcastable, so you can solve 100 eigenproblems by
feeding in a 3x3x100 array. Both C<$ev> and C<$e> get extra dimensions according ly. feeding in a 3x3x100 array. Both C<$ev> and C<$e> get extra dimensions according ly.
If called in scalar context C<eigens> hands back only the eigenvalues. This If called in scalar context C<eigens> hands back only the eigenvalues. This
is somewhat wasteful, as it calculates the eigenvectors anyway. is somewhat wasteful, as it calculates the eigenvectors anyway.
The eigenvectors are returned in COLUMNS of the returned PDL (ie the The eigenvectors are returned in COLUMNS of the returned PDL (ie the
the 0 dimension). That makes it slightly easier to access individual the 0 dimension). That makes it slightly easier to access individual
eigenvectors, since the 0th dim of the output PDL runs across the eigenvectors, since the 0th dim of the output PDL runs across the
eigenvectors and the 1st dim runs across their components. eigenvectors and the 1st dim runs across their components.
skipping to change at line 582 skipping to change at line 584
($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 627 "MatrixOps.pm" #line 636 "MatrixOps.pm"
#line 1060 "../../blib/lib/PDL/PP.pm" #line 1059 "../../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 631 skipping to change at line 633
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 678 "MatrixOps.pm" #line 688 "MatrixOps.pm"
#line 1060 "../../blib/lib/PDL/PP.pm"
#line 1061 "../../blib/lib/PDL/PP.pm"
*eigens = \&PDL::eigens; *eigens = \&PDL::eigens;
#line 684 "MatrixOps.pm" #line 695 "MatrixOps.pm"
#line 1059 "../../blib/lib/PDL/PP.pm" #line 1058 "../../blib/lib/PDL/PP.pm"
=head2 svd =head2 svd
=for sig =for sig
Signature: (a(n,m); [o]u(n,m); [o,phys]z(n); [o]v(n,n)) Signature: (a(n,m); [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);
=for ref =for ref
Singular value decomposition of a matrix. Singular value decomposition of a matrix.
C<svd> is threadable. C<svd> is broadcastable.
Given an m x n matrix C<$a> that has m rows and n columns (m >= n), Given an m x n matrix C<$a> that has m rows and n columns (m >= n),
C<svd> computes matrices C<$u> and C<$v>, and a vector of the singular C<svd> computes matrices C<$u> and C<$v>, and a vector of the singular
values C<$s>. Like most implementations, C<svd> computes what is values C<$s>. Like most implementations, C<svd> computes what is
commonly referred to as the "thin SVD" of C<$a>, such that C<$u> is m commonly referred to as the "thin SVD" of C<$a>, such that C<$u> is m
x n, C<$v> is n x n, and there are <=n singular values. As long as m x n, C<$v> is n x n, and there are <=n singular values. As long as m
>= n, the original matrix can be reconstructed as follows: >= n, the original matrix can be reconstructed as follows:
($u,$s,$v) = svd($x); ($u,$s,$v) = svd($x);
$ess = zeroes($x->dim(0),$x->dim(0)); $ess = zeroes($x->dim(0),$x->dim(0));
skipping to change at line 689 skipping to change at line 692
EXAMPLE EXAMPLE
The computing literature has loads of examples of how to use SVD. The computing literature has loads of examples of how to use SVD.
Here's a trivial example (used in L<PDL::Transform::map|PDL::Transform/map>) Here's a trivial example (used in L<PDL::Transform::map|PDL::Transform/map>)
of how to make a matrix less, er, singular, without changing the of how to make a matrix less, er, singular, without changing the
orientation of the ellipsoid of transformation: orientation of the ellipsoid of transformation:
{ my($r1,$s,$r2) = svd $x; { my($r1,$s,$r2) = svd $x;
$s++; # fatten all singular values $s++; # fatten all singular values
$r2 *= $s; # implicit threading 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 755 "MatrixOps.pm" #line 767 "MatrixOps.pm"
#line 1060 "../../blib/lib/PDL/PP.pm"
#line 1061 "../../blib/lib/PDL/PP.pm"
*svd = \&PDL::svd; *svd = \&PDL::svd;
#line 761 "MatrixOps.pm" #line 774 "MatrixOps.pm"
#line 817 "matrixops.pd" #line 815 "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 731 skipping to change at line 735
$lu = lu_decomp($x, $perm, $par); # $perm and $par are outputs! $lu = lu_decomp($x, $perm, $par); # $perm and $par are outputs!
lu_decomp($x->inplace,$perm,$par); # Everything in place. lu_decomp($x->inplace,$perm,$par); # Everything in place.
=for description =for description
C<lu_decomp> returns an LU decomposition of a square matrix, C<lu_decomp> returns an LU decomposition of a square matrix,
using Crout's method with partial pivoting. It's ported using Crout's method with partial pivoting. It's ported
from I<Numerical Recipes>. The partial pivoting keeps it from I<Numerical Recipes>. The partial pivoting keeps it
numerically stable but means a little more overhead from numerically stable but means a little more overhead from
threading. broadcasting.
C<lu_decomp> decomposes the input matrix into matrices L and C<lu_decomp> decomposes the input matrix into matrices L and
U such that LU = A, L is a subdiagonal matrix, and U is a U such that LU = A, L is a subdiagonal matrix, and U is a
superdiagonal matrix. By convention, the diagonal of L is superdiagonal matrix. By convention, the diagonal of L is
all 1's. all 1's.
The single output matrix contains all the variable elements The single output matrix contains all the variable elements
of both the L and U matrices, stacked together. Because the of both the L and U matrices, stacked together. Because the
method uses pivoting (rearranging the lower part of the method uses pivoting (rearranging the lower part of the
matrix for better numerical stability), you have to permute matrix for better numerical stability), you have to permute
skipping to change at line 769 skipping to change at line 773
those are handled OK but give a zero determinant (and hence those are handled OK but give a zero determinant (and hence
can't be inverted). can't be inverted).
C<lu_decomp> uses pivoting, which rearranges the values in the C<lu_decomp> uses pivoting, which rearranges the values in the
matrix for more numerical stability. This makes it really matrix for more numerical stability. This makes it really
good for large and even near-singular matrices. There is good for large and even near-singular matrices. There is
a non-pivoting version C<lu_decomp2> available which is a non-pivoting version C<lu_decomp2> available which is
from 5 to 60 percent faster for typical problems at from 5 to 60 percent faster for typical problems at
the expense of failing to compute a result in some cases. the expense of failing to compute a result in some cases.
Now that the C<lu_decomp> is threaded, it is the recommended Now that the C<lu_decomp> is broadcasted, it is the recommended
LU decomposition routine. It no longer falls back to C<lu_decomp2>. LU decomposition routine. It no longer falls back to C<lu_decomp2>.
C<lu_decomp> is ported from I<Numerical Recipes> to PDL. It C<lu_decomp> is ported from I<Numerical Recipes> to PDL. It
should probably be implemented in C. should probably be implemented in C.
=cut =cut
*PDL::lu_decomp = \&lu_decomp; *PDL::lu_decomp = \&lu_decomp;
sub lu_decomp { sub lu_decomp {
skipping to change at line 882 skipping to change at line 886
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) { if(wantarray) {
return ($out,$permute,$parity); return ($out,$permute,$parity);
} }
$out; $out;
} }
#line 943 "MatrixOps.pm" #line 957 "MatrixOps.pm"
#line 1000 "matrixops.pd" #line 998 "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 1003 skipping to change at line 1007
($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) { if(wantarray) {
return ($out,$perm,$par); return ($out,$perm,$par);
} }
$out; $out;
} }
#line 1067 "MatrixOps.pm" #line 1082 "MatrixOps.pm"
#line 1125 "matrixops.pd" #line 1123 "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 1057 skipping to change at line 1061
# or with LAPACK # or with LAPACK
use PDL::LinearAlgebra::Real; use PDL::LinearAlgebra::Real;
getrf($lu=$A->copy, $ipiv=null, $info=null); getrf($lu=$A->copy, $ipiv=null, $info=null);
getrs($lu, 1, $x=$B->transpose->copy, $ipiv, $info=null); # again, need transp ose getrs($lu, 1, $x=$B->transpose->copy, $ipiv, $info=null); # again, need transp ose
$x=$x->inplace->transpose; $x=$x->inplace->transpose;
# or with GSL # or with GSL
use PDL::GSL::LINALG; use PDL::GSL::LINALG;
LU_decomp(my $lu=$A->copy, my $p=null, my $signum=null); LU_decomp(my $lu=$A->copy, my $p=null, my $signum=null);
# $B and $x, first dim is because GSL treats as vector, higher dims thread # $B and $x, first dim is because GSL treats as vector, higher dims broadcast
# so we transpose in and back out # so we transpose in and back out
LU_solve($lu, $p, $B->transpose, my $x=null); LU_solve($lu, $p, $B->transpose, my $x=null);
$x=$x->inplace->transpose; $x=$x->inplace->transpose;
# proof of the pudding is in the eating: # proof of the pudding is in the eating:
print $A x $x; print $A x $x;
=for description =for description
Given the LU decomposition of a square matrix (from L</lu_decomp>), Given the LU decomposition of a square matrix (from L</lu_decomp>),
C<lu_backsub> does back substitution into the matrix to solve C<lu_backsub> does back substitution into the matrix to solve
C<a x = b> for given vector C<b>. It is separated from the C<a x = b> for given vector C<b>. It is separated from the
C<lu_decomp> method so that you can call the cheap C<lu_backsub> C<lu_decomp> method so that you can call the cheap C<lu_backsub>
multiple times and not have to do the expensive LU decomposition multiple times and not have to do the expensive LU decomposition
more than once. more than once.
C<lu_backsub> acts on single vectors and threads in the usual C<lu_backsub> acts on single vectors and broadcasts in the usual
way, which means that it treats C<$y> as the I<transpose> way, which means that it treats C<$y> as the I<transpose>
of the input. If you want to process a matrix, you must of the input. If you want to process a matrix, you must
hand in the I<transpose> of the matrix, and then transpose hand in the I<transpose> of the matrix, and then transpose
the output when you get it back. that is because pdls are the output when you get it back. that is because pdls are
indexed by (col,row), and matrices are (row,column) by indexed by (col,row), and matrices are (row,column) by
convention, so a 1-D pdl corresponds to a row vector, not a convention, so a 1-D pdl corresponds to a row vector, not a
column vector. column vector.
If C<$lu> is dense and you have more than a few points to If C<$lu> is dense and you have more than a few points to
solve for, it is probably cheaper to find C<a^-1> with solve for, it is probably cheaper to find C<a^-1> with
skipping to change at line 1116 skipping to change at line 1120
unless defined($lu); unless defined($lu);
barf("Usage: \$x = lu_backsub(\$lu,\$perm,\$y); all must be PDLs\n") barf("Usage: \$x = lu_backsub(\$lu,\$perm,\$y); all must be PDLs\n")
unless(UNIVERSAL::isa($lu,'PDL') && unless(UNIVERSAL::isa($lu,'PDL') &&
UNIVERSAL::isa($perm,'PDL') && UNIVERSAL::isa($perm,'PDL') &&
UNIVERSAL::isa($y,'PDL')); UNIVERSAL::isa($y,'PDL'));
my $n = $y->dim(0); my $n = $y->dim(0);
my $n1 = $n; $n1--; my $n1 = $n; $n1--;
# Make sure threading dimensions are compatible. # Make sure broadcasting dimensions are compatible.
# There are two possible sources of thread dims: # There are two possible sources of broadcast dims:
# #
# (1) over multiple LU (i.e., $lu,$perm) instances # (1) over multiple LU (i.e., $lu,$perm) instances
# (2) over multiple B (i.e., $y) column instances # (2) over multiple B (i.e., $y) column instances
# #
# The full dimensions of the function call looks like # The full dimensions of the function call looks like
# #
# lu_backsub( lu(m,m,X), perm(m,X), b(m,Y) ) # lu_backsub( lu(m,m,X), perm(m,X), b(m,Y) )
# #
# where X is the list of extra LU dims and Y is # where X is the list of extra LU dims and Y is
# the list of extra B dims. We have several possible # the list of extra B dims. We have several possible
skipping to change at line 1147 skipping to change at line 1151
my $m = $ludims->slice("(0)"); # this is the sig dimension my $m = $ludims->slice("(0)"); # this is the sig dimension
unless ( ($ludims->slice(0) == $m) and ($ludims->slice(1) == $m) and unless ( ($ludims->slice(0) == $m) and ($ludims->slice(1) == $m) and
($permdims->slice(0) == $m) and ($bdims->slice(0) == $m)) { ($permdims->slice(0) == $m) and ($bdims->slice(0) == $m)) {
barf "lu_backsub: mismatched sig dimensions"; barf "lu_backsub: mismatched sig dimensions";
} }
my $lunumthr = $ludims->dim(0)-2; my $lunumthr = $ludims->dim(0)-2;
my $permnumthr = $permdims->dim(0)-1; my $permnumthr = $permdims->dim(0)-1;
my $bnumthr = $bdims->dim(0)-1; my $bnumthr = $bdims->dim(0)-1;
unless ( ($lunumthr == $permnumthr) and ($ludims->slice("1:-1") == $permdims) ->all ) { unless ( ($lunumthr == $permnumthr) and ($ludims->slice("1:-1") == $permdims) ->all ) {
barf "lu_backsub: \$lu and \$perm thread dims not equal! \n"; barf "lu_backsub: \$lu and \$perm broadcast dims not equal! \n";
} }
# (2) If X == Y then default threading is ok # (2) If X == Y then default broadcasting is ok
if ( ($bnumthr==$permnumthr) and ($bdims==$permdims)->all) { if ( ($bnumthr==$permnumthr) and ($bdims==$permdims)->all) {
print STDERR "lu_backsub: have explicit thread dims, goto THREAD_OK\n" if print STDERR "lu_backsub: have explicit broadcast dims, goto BROADCAST_OK\
$PDL::debug; n" if $PDL::debug;
goto THREAD_OK; goto BROADCAST_OK;
} }
# (3) If X == (x,Y) then add x dummy to lu,perm # (3) If X == (x,Y) then add x dummy to lu,perm
# (4) If ndims(X) > ndims(Y) then must have #3 # (4) If ndims(X) > ndims(Y) then must have #3
# (5) If ndims(X) < ndims(Y) then foreach # (5) If ndims(X) < ndims(Y) then foreach
# non-trivial leading dim in X (x0,x1,..) # non-trivial leading dim in X (x0,x1,..)
# insert dummy (x0,x1) into lu and perm # insert dummy (x0,x1) into lu and perm
# This means that threading occurs over all # This means that broadcasting occurs over all
# leading non-trivial (not length 1) dims of # leading non-trivial (not length 1) dims of
# B unless all the thread dims are explicitly # B unless all the broadcast dims are explicitly
# matched to the LU dims. # matched to the LU dims.
THREAD_OK: BROADCAST_OK:
# Permute the vector and make a copy if necessary. # Permute the vector and make a copy if necessary.
my $out = $y->dummy(1,$y->dim(0))->index($perm->dummy(1)); my $out = $y->dummy(1,$y->dim(0))->index($perm->dummy(1));
$out = $out->sever if !$y->is_inplace; $out = $out->sever if !$y->is_inplace;
print STDERR "lu_backsub: starting with \$out" . pdl($out->dims) . "\n" if $P DL::debug; print STDERR "lu_backsub: starting with \$out" . pdl($out->dims) . "\n" if $P DL::debug;
# Make sure threading over lu happens OK... # Make sure broadcasting over lu happens OK...
if($out->ndims < $lu->ndims-1) { if($out->ndims < $lu->ndims-1) {
print STDERR "lu_backsub: adjusting dims for \$out" . pdl($out->dims) . "\ n" if $PDL::debug; print STDERR "lu_backsub: adjusting dims for \$out" . pdl($out->dims) . "\ n" if $PDL::debug;
do { do {
$out = $out->dummy(-1,$lu->dim($out->ndims+1)); $out = $out->dummy(-1,$lu->dim($out->ndims+1));
} while($out->ndims < $lu->ndims-1); } while($out->ndims < $lu->ndims-1);
} }
## Do forward substitution into L ## Do forward substitution into L
my $row; my $r1; my $row; my $r1;
skipping to change at line 1200 skipping to change at line 1204
my $tmp; # work around perl -d "feature my $tmp; # work around perl -d "feature
($tmp = $out->index($row)) -= ($lu->slice("0:$r1,$row") * ($tmp = $out->index($row)) -= ($lu->slice("0:$r1,$row") *
$out->slice("0:$r1") $out->slice("0:$r1")
)->sumover; )->sumover;
} }
## Do backward substitution into U, and normalize by the diagonal ## Do backward substitution into U, and normalize by the diagonal
my $ludiag = $lu->diagonal(0,1); my $ludiag = $lu->diagonal(0,1);
{ {
my $tmp; # work around for perl -d "feature" my $tmp; # work around for perl -d "feature"
($tmp = $out->index($n1)) /= $ludiag->index($n1)->dummy(0); # TODO: check threading ($tmp = $out->index($n1)) /= $ludiag->index($n1)->dummy(0); # TODO: check broadcasting
} }
for ($row=$n1; $row>0; $row--) { for ($row=$n1; $row>0; $row--) {
$r1 = $row-1; $r1 = $row-1;
my $tmp; # work around for perl -d "feature" my $tmp; # work around for perl -d "feature"
($tmp = $out->index($r1)) -= ($lu->slice("$row:$n1,$r1") * # TODO: check thread dims ($tmp = $out->index($r1)) -= ($lu->slice("$row:$n1,$r1") * # TODO: check broadcast dims
$out->slice("$row:$n1") $out->slice("$row:$n1")
)->sumover; )->sumover;
($tmp = $out->index($r1)) /= $ludiag->index($r1)->dummy(0); # TODO: check thread 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 1284 "MatrixOps.pm" #line 1300 "MatrixOps.pm"
#line 1059 "../../blib/lib/PDL/PP.pm" #line 1058 "../../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 1256 skipping to change at line 1260
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 1329 "MatrixOps.pm" #line 1346 "MatrixOps.pm"
#line 1060 "../../blib/lib/PDL/PP.pm"
#line 1061 "../../blib/lib/PDL/PP.pm"
*simq = \&PDL::simq; *simq = \&PDL::simq;
#line 1335 "MatrixOps.pm" #line 1353 "MatrixOps.pm"
#line 1059 "../../blib/lib/PDL/PP.pm" #line 1058 "../../blib/lib/PDL/PP.pm"
=head2 squaretotri =head2 squaretotri
=for sig =for sig
Signature: (a(n,n); b(m)) Signature: (a(n,n); b(m))
=for ref =for ref
Convert a symmetric square matrix to triangular vector storage. Convert a symmetric square matrix to triangular vector storage.
=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 1362 "MatrixOps.pm" #line 1381 "MatrixOps.pm"
#line 1060 "../../blib/lib/PDL/PP.pm"
#line 1061 "../../blib/lib/PDL/PP.pm"
*squaretotri = \&PDL::squaretotri; *squaretotri = \&PDL::squaretotri;
#line 1368 "MatrixOps.pm" #line 1388 "MatrixOps.pm"
#line 1420 "matrixops.pd" #line 1418 "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 1386 "MatrixOps.pm" #line 1407 "MatrixOps.pm"
# Exit with OK status # Exit with OK status
1; 1;
 End of changes. 72 change blocks. 
75 lines changed or deleted 81 lines changed or added

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