Version 0.03 alpha
PDL::Meschach Etienne Grossmann, 11/11/1996 [etienne@isr.ist.utl.pt] meschach 1.2 David E. Stewart [david.stewart@anu.edu.au] Zbigniew Leyk [zbigniew.leyk@anu.edu.au]
PDL::Meschach links PDL to a few matrix functions from meschach 1.2 : Diagonal,upper,lower triangle extraction ... Matrix exponentiation ... LU , Cholesky , QR Factorisation and associated Linear equation solvers. Symmetric matrix eigenvector/eigenvalue extraction. Singular value decomposition.
The functionalities are available either through "friendly" functions, e.g. : # Solve A.x == b using LU decomposition $x = lusolve($A,$b); or through "raw" functions : $LU = $A + 0 ; # Copy $A lufac_($LU,$Perm) # LU decomposition overwrites $LU lusolve_($x,$b,$LU,$Perm) That may be more efficient in terms of memory allocation. All "raw" function names have a trailing underscore.
* Dimensions are expressed ALWAYS as COLUMN, ROW instead of the usual matrix row, column. * pdl members "incs" and "offs" are not yet handled correctly, and neither can they be detected reliably. The pdls passed to Meschach routines are assumed to hold data consecutively arranged in memory.
use PDL::Meschach; # Only "Friendly" functions. or use PDL::Meschach qw( :Raw ) # "Raw" functions too.
eval "use Meschach qw( :All )"; if($@ ne ""){ print "Meschach NOT AVAILABLE : \n$@\n"; } else { print "Meschach found\n";} somewhere in it, e.g. before the line : eval "use Term::ReadLine";
$T = ut( $A ); Puts the upper triangle of $A in $T. $T = ut( $Col, $Row ); Sets to 1 the upper-triangle of T, to 0 its strictly lower triangle. Raw function : ut_($T,$A); ut_($T,$Col,$Row); The output is the same type as th input. For lower triangle, use lt_, lt.
$Vec = diag( $Mat ); diag_( $Vec, $Mat); Both put into $Vec the diagonal of $Mat. output is PDL_D (!). $Mat = diag ( $Vec [,$cols, $rows] ); diag_($Vec,$Mat); Make $Mat a diagonal matrix with $Vec as diagonal. $Mat may be of arbitrary size if $cols, $rows are used : $Mat = diag ( $Vec,$c ); # $c x $c $Mat = diag ( $Vec,$c, $r ); # $c x $m output is PDL_F (!)
$Mat = ident($Cols [,$Rows = $Cols ] ) Returns a PDL_D Identity Matrix.
$Out = mpow( $In, $Pow ); mpow_( $Out, $In, $Pow [,$Coerce] ); Integer (negative or positive) powers of a matrix, If $Coerce is true (default), the result is Real typed. Otherwise, $$Out{Datatype} is unchanged.
$Out = inv( $In ); inv_( $Out, $In [,$Coerce] ); Inverse of a matrix.
Solve $A x $x == $b , by LU Factorization, with pivoting. ($LU,$Perm,$x) = lusolve( $b, $A ); $LU,$Perm describe the factorization. They may be re-used (which spares some computation). $LU is a matrix the same size as $A. $Perm is an integer vector describing the pivoting used in the factorization. $x = lusolve($b, $LU, $Perm ); The third argument is what decides lusolve not to factor. In all cases ($LU,$Perm,$x) is returned. An estimate of the conditioning (the ratio of the greatest and smallest eigenvalues) is returned by : lucond($LU,$Perm); Raw method : $LU = $A + 0 ; # Copy $A. lufac_( $LU, $Perm ); # Factorize. lusolve_( $x, $b, $LU, $Perm ); # Solve.
Solve $A x $x == $b , by QR Factorization (A = Q.R where Q is orthogonal, and R upper triangular). The "QR" functions are almost equivalent to the "LU" ones. ($QR,$V,$x) = qrsolve( $b, $A ); $x = qrsolve($b, $QR, $V ); $QR,$V describe the factorization. The "R" matrix is represented by the strictly higher triangle in $QR, and its diagonal is $V. $Q is represented by the lower triangle. The third argument is what decides qrsolve not to factor. In all cases ($QR,$V,$x) is returned. Raw method : $QR = $A + 0 ; # Copy $A. qrfac_( $QR, $V ); # Factorize. qrsolve_( $x, $b, $QR, $V ); # Solve. An estimate of the conditioning (the ratio of the greatest and smallest eigenvalues) is returned by : qrcond($QR); # only QR is passed.
If A is positive definite, solving $A x $x == $b by Cholesky Factorization may be more efficient : ($CH,$x) = chsolve( $b, $A ); $CH may be re-used : $x = chsolve ( $b , $CH , 1 ); If the third argument is true, chsolve considers that it is already in factored form. Raw method : $CH = $A + 0 ; # Copy chfac_($CH); # Factor chsolve_($x,$b,$CH); # Solve
Symmetric Matrix Eigenvalues/vectors ( $Mat, $Vec ) = symmeig( $A ); $Mat contains the eigenvectors of $A, $Vec contains the eigenvalues of $A, Raw method : symmeig_( $Mat, $Vec, $A ); symmeig_ returns true upon success.
Singular Value Decomposition of a ncol,nrow matrix : ($U,$V,$l) = svd( $A ) ; $U Contains the left "singular vectors" (nrow,nrow matrix). $V Contains the right "singular vectors" (ncol,ncol matrix). $l Contains the "singular values" (min(ncol,nrow) vector). Raw method : svd_( $U, $V, $l, $A ); or svd_( $l, $A ); # Only the "singular values". svd_ returns true upon success.
"Real" is the floating point datatype used by Meschach. It probably always corresponds to PDL_D.
Meschach should function even if it is PDL_F; but it hasn't been tested.
Similar conversion is done for meschach's "u_int" and a PDL integer type (should be PDL_L).
The equivalence between types is done in the BOOT: part of Meschach.xs
The results of most functions are "Real"-typed pdls.
Sometimes (e.g. the $Perm argument in lufac_ ) it is an integer-typed pdls.
"Raw" functions may conserve the type of their return-argument when they accept a $coerce argument.
To install PDL::Meschach, copy and paste the appropriate command in to your terminal.
cpanm
cpanm PDL::Meschach
CPAN shell
perl -MCPAN -e shell install PDL::Meschach
For more information on module installation, please visit the detailed CPAN module installation guide.