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.
Meschach should function even if it is PDL_F; but it hasn't been tested.
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.
