/ 
Meschach-0.03
River stage zero No dependents
/ Version

NAME

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]

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.



* 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.

USAGE

To load Meschach in a script : either :

 use PDL::Meschach;           # Only "Friendly" functions.

or   

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

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, 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.
qrsolve
 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.
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.