"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "GENERATED/PDL/MatrixOps.pm" 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.pm  (PDL-2.080)	:	MatrixOps.pm  (PDL-2.081)
skipping to change at line 151		skipping to change at line 151
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;
}		}
#line 178 "MatrixOps.pm"		#line 178 "MatrixOps.pm"
#line 156 "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))
skipping to change at line 283		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 308 "MatrixOps.pm"		#line 307 "MatrixOps.pm"
#line 288 "matrixops.pd"		#line 287 "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 341		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 369 "MatrixOps.pm"		#line 368 "MatrixOps.pm"
#line 350 "matrixops.pd"		#line 349 "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 425		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 456 "MatrixOps.pm"		#line 455 "MatrixOps.pm"
#line 1058 "../../blib/lib/PDL/PP.pm"		#line 949 "../../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 471		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 510 "MatrixOps.pm"		#line 509 "MatrixOps.pm"
#line 1059 "../../blib/lib/PDL/PP.pm"		#line 950 "../../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 507		skipping to change at line 507
)->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 549 "MatrixOps.pm"		#line 548 "MatrixOps.pm"
#line 1060 "../../blib/lib/PDL/PP.pm"		#line 951 "../../blib/lib/PDL/PP.pm"
*eigens_sym = \&PDL::eigens_sym;		*eigens_sym = \&PDL::eigens_sym;
#line 556 "MatrixOps.pm"		#line 555 "MatrixOps.pm"
#line 1058 "../../blib/lib/PDL/PP.pm"		#line 949 "../../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 584		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 636 "MatrixOps.pm"		#line 635 "MatrixOps.pm"
#line 1059 "../../blib/lib/PDL/PP.pm"		#line 950 "../../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 633		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 688 "MatrixOps.pm"		#line 687 "MatrixOps.pm"
#line 1060 "../../blib/lib/PDL/PP.pm"		#line 951 "../../blib/lib/PDL/PP.pm"
*eigens = \&PDL::eigens;		*eigens = \&PDL::eigens;
#line 695 "MatrixOps.pm"		#line 694 "MatrixOps.pm"
#line 1058 "../../blib/lib/PDL/PP.pm"		#line 949 "../../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);
skipping to change at line 702		skipping to change at line 702
$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 767 "MatrixOps.pm"		#line 766 "MatrixOps.pm"
#line 1060 "../../blib/lib/PDL/PP.pm"		#line 951 "../../blib/lib/PDL/PP.pm"
*svd = \&PDL::svd;		*svd = \&PDL::svd;
#line 774 "MatrixOps.pm"		#line 773 "MatrixOps.pm"
#line 815 "matrixops.pd"		#line 814 "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 828		skipping to change at line 828
$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 870		skipping to change at line 871
$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;
}		}
#line 957 "MatrixOps.pm"		#line 952 "MatrixOps.pm"
#line 998 "matrixops.pd"		#line 993 "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 1002		skipping to change at line 999
# 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;
}		}
#line 1082 "MatrixOps.pm"		#line 1074 "MatrixOps.pm"
#line 1123 "matrixops.pd"		#line 1115 "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 1222		skipping to change at line 1216
)->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 1300 "MatrixOps.pm"		#line 1292 "MatrixOps.pm"
#line 1058 "../../blib/lib/PDL/PP.pm"		#line 949 "../../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 1260		skipping to change at line 1254
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 1346 "MatrixOps.pm"		#line 1338 "MatrixOps.pm"
#line 1060 "../../blib/lib/PDL/PP.pm"		#line 951 "../../blib/lib/PDL/PP.pm"
*simq = \&PDL::simq;		*simq = \&PDL::simq;
#line 1353 "MatrixOps.pm"		#line 1345 "MatrixOps.pm"
#line 1058 "../../blib/lib/PDL/PP.pm"		#line 949 "../../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 1381 "MatrixOps.pm"		#line 1373 "MatrixOps.pm"
#line 1060 "../../blib/lib/PDL/PP.pm"		#line 951 "../../blib/lib/PDL/PP.pm"
*squaretotri = \&PDL::squaretotri;		*squaretotri = \&PDL::squaretotri;
#line 1388 "MatrixOps.pm"		#line 1380 "MatrixOps.pm"
#line 1418 "matrixops.pd"		#line 1410 "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 1407 "MatrixOps.pm"		#line 1399 "MatrixOps.pm"
# Exit with OK status		# Exit with OK status
1;		1;
End of changes. 45 change blocks.
53 lines changed or deleted		48 lines changed or added

---