"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "Basic/MatrixOps/matrixops.pd" between
PDL-2.081.tar.gz and PDL-2.082.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.081)	:	matrixops.pd  (PDL-2.082)
skipping to change at line 325		skipping to change at line 325
provide the key but leave the value undefined, then the LU decomposition		provide the key but leave the value undefined, then the LU decomposition
goes in here; if you put an LU decomposition here, it will be used and		goes in here; if you put an LU decomposition here, it will be used and
the matrix will not be decomposed again.		the matrix will not be decomposed again.
=back		=back
=cut		=cut
*PDL::det = \&det;		*PDL::det = \&det;
sub det {		sub det {
my($x) = shift;		my ($x, $opt) = @_;
my($opt) = shift;
$opt = {} unless defined($opt);		$opt = {} unless defined($opt);
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 : 0 );
}		}
EOD		EOD
######################################################################		######################################################################
pp_add_exported('','determinant');		pp_add_exported('','determinant');
pp_addpm(<<'EOD');		pp_addpm(<<'EOD');
=head2 determinant		=head2 determinant
skipping to change at line 709		skipping to change at line 706
DEVEL NOTES:		DEVEL NOTES:
For now, there is no distinction between a complex eigenvalue and an		For now, there is no distinction between a complex eigenvalue and an
invalid eigenvalue, although the underlying code generates complex		invalid eigenvalue, although the underlying code generates complex
numbers. It might be useful to be able to return complex eigenvalues.		numbers. It might be useful to be able to return complex eigenvalues.
=for usage		=for usage
($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
=cut
',);		',);
######################################################################		######################################################################
### svd		### svd
pp_def(		pp_def(
"svd",		"svd",
HandleBad => 0,		HandleBad => 0,
Pars => 'a(n,m); [o]u(n,m); [o,phys]z(n); [o]v(n,n);',		Pars => 'a(n,m); [t]w(wsize); [o]u(n,m); [o,phys]z(n); [o]v(n,n);',
GenericTypes => ['D'],		GenericTypes => ['D'],
		RedoDimsCode => '
		if ($SIZE(m)<$SIZE(n))
		$CROAK("svd requires input ndarrays to have m >= n; you have m=%td a
		nd n=%td. Try inputting the transpose. See the docs for svd.",$SIZE(m),$SIZE(n)
		);
		$SIZE(wsize) = $SIZE(n) * ($SIZE(m) + $SIZE(n));
		',
Code => '		Code => '
extern void SVD( double *W, double *Z, int nRow, int nCol );		extern void SVD( double *W, double *Z, int nRow, int nCol );
int sm = $SIZE (m), sn = $SIZE (n), i;		double *t = $P(w);
if (sm<sn) {
$CROAK("svd requires input ndarrays to have m >= n; you have m=
%d and n=%d. Try inputting the transpose. See the docs for svd.",sm,sn);
}
double *w, *t, zv;
t = w = (double *) malloc(sn*(sm+sn)*sizeof(double));
loop (m) %{		loop (m) %{
loop(n) %{		loop(n) %{
*t++ = $a ();		*t++ = $a();
%}		%}
%}		%}
SVD(w, $P (z), sm, sn);		SVD($P(w), $P(z), $SIZE(m), $SIZE(n));
t = w;
loop (n) %{		loop (n) %{
zv = sqrt($z ());		$z() = sqrt($z());
$z () = zv;
%}		%}
		t = $P(w);
loop (m) %{		loop (m) %{
loop (n) %{		loop (n) %{
$u () = *t++/$z ();		$u() = *t++/$z();
%}		%}
%}		%}
loop (n) %{		loop (n1) %{
for (i=0;i<sn;i++) {		loop (n0) %{
$v (n0=>i, n1=>n) = *t++;		$v() = *t++;
}		%}
%}		%}
free(w);
',		',
, Doc => q{		, Doc => q{
=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.
skipping to change at line 800		skipping to change at line 792
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 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
}		}
=cut
},);		},);
######################################################################		######################################################################
pp_add_exported('','lu_decomp');		pp_add_exported('','lu_decomp');
pp_addpm(<<'EOD');		pp_addpm(<<'EOD');
=head2 lu_decomp		=head2 lu_decomp
=for sig		=for sig
End of changes. 16 change blocks.
31 lines changed or deleted		21 lines changed or added

---