The London Perl and Raku Workshop takes place on 26th Oct 2024. If your company depends on Perl, please consider sponsoring and/or attending.
 / 
Meschach-0.01
River stage zero No dependents
/ Version

NAME

Version 0.01 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]

DESCRIPTION

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.

INTRODUCTION

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.

CAVEAT

Dimensions are expressed 

       ALWAYS as COLUMN, ROW 

instead of the usual matrix row, column.

USAGE

Load Meschach in a script using either :

 use Meschach;           # Only "Friendly" functions.

or   

 use Meschach qw( :Raw ) # "Raw" functions too.


 Note : "use PDL::Meschach;"  fails. Bug in Makefile.PL ? 

To use Meschach in perldl, insert


       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";

AVAILABLE FUNCTIONS

ut, lt
   $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.
diag
   $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 (!)
ident
$Mat = ident($Cols [,$Rows = $Cols ] )

Returns a PDL_D Identity Matrix.
mpow
   $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.
inv
   $Out = inv( $In );

   inv_( $Out, $In [,$Coerce] );

Inverse of a matrix.
lusolve
Solve $A x $x == $b , by LU Factorization  

   ($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. 

Raw method :

   $LU = $A + 0 ;                   # Copy $A 
   lufac_( $LU, $Perm );            # Factorize. 
   lusolve_( $x, $b, $LU, $Perm );  # Solve
chsolve
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
symmeig
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.
svd
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.

TYPE CONVERSION

  • "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.