/ 
PDL-2.003
River stage three • 116 direct dependents • 124 total dependents
/ Basic/Math/math.pd

NAME

PDL::Math - extended mathematical operations and special functions

SYNOPSIS

use PDL::Math;

use PDL::Graphics::TriD;
imag3d [SURF2D,bessj0(rvals(zeroes(50,50))/2)];

DESCRIPTION

This module extends PDL with more advanced mathematical functions than provided by standard Perl.

All the functions have one input pdl, and one output, unless otherwise stated.

The functions are usually available from the system maths library, however if they are not (determined when PDL is compiled) a version from the Cephes math library is used.

The standard Bessel Functions

log gamma function

This returns 2 piddles -- the first set gives the log(gamma) values, while the second set, of integer values, gives the sign of the gamma function. This is useful for determining factorials, amongst other things.

Singular value decomposition of array.

($u, $s, $v) = svd($a);

',);

# XXX The next two really need driver routines...

pp_def("eigens", Pars => '[phys]a(m); [o,phys]ev(n,n); [o,phys]e(n)', GenericTypes => [D], Code => ' extern void eigens( double *A, double *RR, double *E, int N ); register int sn = $SIZE(n); if($SIZE(m) != (sn * (sn + 1))/2) { barf("Wrong sized args for eigens"); } eigens($P(a), $P(ev), $P(e), sn); ', PMCode =>' sub PDL::eigens { my ($a) = @_; my (@d) = $a->dims; barf "Need real square matrix for eigens" if $#d != 1 or $d[0] != $d[1]; my ($n) = $d[0]; my ($sym) = 0.5*($a + $a->mv(0,1)); my ($err) = PDL::max(abs($sym)); barf "Need symmetric component non-zero for eigens" if $err == 0; $err = PDL::max(abs($a-$sym))/$err; warn "Using symmetrized version of the matrix in eigens" if $err > 1e-5; # Get lower triangular form -- my ($lt) = PDL::where($sym,$sym->xvals <= $sym->yvals)->copy; my ($ev) = PDL->zeroes($n,$n); my ($e) = PDL->zeroes($n); &PDL::_eigens_int($lt, $ev, $e); return $ev, $e; } ' , Doc => ' =for ref

Eigenvalues and -vectors of a symmetric square matrix. If passed an asymmetric matrix, the routine will warn and symmetrize it.

($e, $ev) = eigens($a);

',);

# XXX Destroys a!!! # To use the new a again, must store both a and ips. pp_def("simq", Pars => '[phys]a(n,n); [phys]b(n); [o,phys]x(n); int [o,phys]ips(n)', OtherPars => 'int flag;', GenericTypes => [D], Code => ' extern int simq( double *A, double *B, double *X, int n, int flag, int *IPS ); simq($P(a),$P(b),$P(x),$SIZE(n),$COMP(flag),$P(ips)); ', Doc => ' =for ref

Solution of simultaneous linear equations, a x = b.

a is an n x n matrix (i.e., a vector of length n*n), stored row-wise: that is, a(i,j) = a[ij], where ij = i*n + j. While this is the transpose of the normal column-wise storage, this corresponds to normal PDL usage. The contents of matrix a may be altered (but may be required for subsequent calls with flag = -1).

b, x, ips are vectors of length n.

Set flag=0 to solve. Set flag=-1 to do a new back substitution for different b vector using the same a matrix previously reduced when flag=0 (the ips vector generated in the previous solution is also required).

');

pp_def("squaretotri", Pars => 'a(n,n); b(m)', Code => ' register int mna=0, nb=0, ns = $SIZE(n); #if (PERL_VERSION >= 5) && (PERL_SUBVERSION >= 57) dXSARGS; #endif

            threadloop %{
		if($SIZE(m) != (ns * (ns+1))/2) {
			barf("Wrong sized args for squaretotri");
		}
		loop(m) %{
                       $b() = $a(n0 => mna, n1 => nb);
                      mna++; if(mna > nb) {mna = 0; nb ++;}
		%}
             %}
	',
	Doc => '
=for ref

Convert a symmetric square matrix to triangular vector storage

BUGS

Hasn't been tested on all platforms to ensure Cephes versions are picked up automatically and used correctly.

AUTHOR

Copyright (C) R.J.R. Williams 1997 (rjrw@ast.leeds.ac.uk), Karl Glazebrook (kgb@aaoepp.aao.gov.au) and Tuomas J. Lukka (Tuomas.Lukka@helsinki.fi).

All rights reserved. There is no warranty. You are allowed to redistribute this software / documentation under certain conditions. For details, see the file COPYING in the PDL distribution. If this file is separated from the PDL distribution, the copyright notice should be included in the file.

2 POD Errors

The following errors were encountered while parsing the POD:

Around line 63:

=pod directives shouldn't be over one line long! Ignoring all 3 lines of content

Around line 70:

=pod directives shouldn't be over one line long! Ignoring all 3 lines of content