Provided by: libpdl-linearalgebra-perl_0.12-2build2_amd64 bug

NAME

       PDL::LinearAlgebra - Linear Algebra utils for PDL

SYNOPSIS

        use PDL::LinearAlgebra;

        $a = random (100,100);
        ($U, $s, $V) = mdsvd($a);

DESCRIPTION

       This module provides a convenient interface to PDL::LinearAlgebra::Real and
       PDL::LinearAlgebra::Complex. Its primary purpose is educational.  You have to know that
       routines defined here are not optimized, particularly in term of memory. Since Blas and
       Lapack use a column major ordering scheme some routines here need to transpose matrices
       before calling fortran routines and transpose back (see the documentation of each
       routine). If you need optimized code use directly  PDL::LinearAlgebra::Real and
       PDL::LinearAlgebra::Complex. It's planned to "port" this module to PDL::Matrix such that
       transpositions will not be necessary, the major problem is that two new modules need to be
       created PDL::Matrix::Real and PDL::Matrix::Complex.

FUNCTIONS

   setlaerror
       Sets action type when an error is encountered, returns previous type. Available values are
       NO, WARN and BARF (predefined constants).  If, for example, in computation of the inverse,
       singularity is detected, the routine can silently return values from computation (see
       manuals), warn about singularity or barf. BARF is the default value.

        # h : x -> g(f(x))

        $a = sequence(5,5);
        $err = setlaerror(NO);
        ($b, $info)= f($a);
        setlaerror($err);
        $info ? barf "can't compute h" : return g($b);

   getlaerror
       Gets action type when an error is encountered.

               0 => NO,
               1 => WARN,
               2 => BARF

   t
        PDL = t(PDL, SCALAR(conj))
        conj : Conjugate Transpose = 1 | Transpose = 0, default = 1;

       Convenient function for transposing real or complex 2D array(s).  For PDL::Complex, if
       conj is true returns conjugate transposed array(s) and doesn't support dataflow.  Supports
       threading.

   issym
        PDL = issym(PDL, SCALAR|PDL(tol),SCALAR(hermitian))
        tol : tolerance value, default: 1e-8 for double else 1e-5
        hermitian : Hermitian = 1 | Symmetric = 0, default = 1;

       Checks symmetricity/Hermitianicity of matrix.  Supports threading.

   diag
       Returns i-th diagonal if matrix in entry or matrix with i-th diagonal with entry. I-th
       diagonal returned flows data back&forth.  Can be used as lvalue subs if your perl supports
       it.  Supports threading.

        PDL = diag(PDL, SCALAR(i), SCALAR(vector)))
        i      : i-th diagonal, default = 0
        vector : create diagonal matrices by threading over row vectors, default = 0

        my $a = random(5,5);
        my $diag  = diag($a,2);
        # If your perl support lvaluable subroutines.
        $a->diag(-2) .= pdl(1,2,3);
        # Construct a (5,5,5) PDL (5 matrices) with
        # diagonals from row vectors of $a
        $a->diag(0,1)

   tritosym
       Returns symmetric or Hermitian matrix from lower or upper triangular matrix.  Supports
       inplace and threading.  Uses tricpy or ctricpy from Lapack.

        PDL = tritosym(PDL, SCALAR(uplo), SCALAR(conj))
        uplo : UPPER = 0 | LOWER = 1, default = 0
        conj : Hermitian = 1 | Symmetric = 0, default = 1;

        # Assume $a is symmetric triangular
        my $a = random(10,10);
        my $b = tritosym($a);

   positivise
       Returns entry pdl with changed sign by row so that average of positive sign > 0.  In other
       words threads among dimension 1 and row  =  -row if sum(sign(row)) < 0.  Works inplace.

        my $a = random(10,10);
        $a -= 0.5;
        $a->xchg(0,1)->inplace->positivise;

   mcrossprod
       Computes the cross-product of two matrix: A' x  B.  If only one matrix is given, takes B
       to be the same as A.  Supports threading.  Uses crossprod or ccrossprod.

        PDL = mcrossprod(PDL(A), (PDL(B))

        my $a = random(10,10);
        my $crossproduct = mcrossprod($a);

   mrank
       Computes the rank of a matrix, using a singular value decomposition.  from Lapack.

        SCALAR = mrank(PDL, SCALAR(TOL))
        TOL:   tolerance value, default : mnorm(dims(PDL),'inf') * mnorm(PDL) * EPS

        my $a = random(10,10);
        my $b = mrank($a, 1e-5);

   mnorm
       Computes norm of real or complex matrix Supports threading.

        PDL(norm) = mnorm(PDL, SCALAR(ord));
        ord :
               0|'inf' : Infinity norm
               1|'one' : One norm
               2|'two' : norm 2 (default)
               3|'fro' : frobenius norm

        my $a = random(10,10);
        my $norm = mnorm($a);

   mdet
       Computes determinant of a general square matrix using LU factorization.  Supports
       threading.  Uses getrf or cgetrf from Lapack.

        PDL(determinant) = mdet(PDL);

        my $a = random(10,10);
        my $det = mdet($a);

   mposdet
       Compute determinant of a symmetric or Hermitian positive definite square matrix using
       Cholesky factorization.  Supports threading.  Uses potrf or cpotrf from Lapack.

        (PDL, PDL) = mposdet(PDL, SCALAR)
        SCALAR : UPPER = 0 | LOWER = 1, default = 0

        my $a = random(10,10);
        my $det = mposdet($a);

   mcond
       Computes the condition number (two-norm) of a general matrix.

       The condition number in two-n is defined:

               norm (a) * norm (inv (a)).

       Uses a singular value decomposition.  Supports threading.

        PDL = mcond(PDL)

        my $a = random(10,10);
        my $cond = mcond($a);

   mrcond
       Estimates the reciprocal condition number of a general square matrix using LU
       factorization in either the 1-norm or the infinity-norm.

       The reciprocal condition number is defined:

               1/(norm (a) * norm (inv (a)))

       Supports threading.  Works on transposed array(s)

        PDL = mrcond(PDL, SCALAR(ord))
        ord :
               0 : Infinity norm (default)
               1 : One norm

        my $a = random(10,10);
        my $rcond = mrcond($a,1);

   morth
       Returns an orthonormal basis of the range space of matrix A.

        PDL = morth(PDL(A), SCALAR(tol))
        tol : tolerance for determining rank, default: 1e-8 for double else 1e-5

        my $a = sequence(10,10);
        my $ortho = morth($a, 1e-8);

   mnull
       Returns an orthonormal basis of the null space of matrix A.  Works on transposed array.

        PDL = mnull(PDL(A), SCALAR(tol))
        tol : tolerance for determining rank, default: 1e-8 for double else 1e-5

        my $a = sequence(10,10);
        my $null = mnull($a, 1e-8);

   minv
       Computes inverse of a general square matrix using LU factorization. Supports inplace and
       threading.  Uses getrf and getri or cgetrf and cgetri from Lapack and returns "inverse,
       info" in array context.

        PDL(inv)  = minv(PDL)

        my $a = random(10,10);
        my $inv = minv($a);

   mtriinv
       Computes inverse of a triangular matrix. Supports inplace and threading.  Uses trtri or
       ctrtri from Lapack.  Returns "inverse, info" in array context.

        (PDL, PDL(info))) = mtriinv(PDL, SCALAR(uplo), SCALAR|PDL(diag))
        uplo : UPPER = 0 | LOWER = 1, default = 0
        diag : UNITARY DIAGONAL = 1, default = 0

        # Assume $a is upper triangular
        my $a = random(10,10);
        my $inv = mtriinv($a);

   msyminv
       Computes inverse of a symmetric square matrix using the Bunch-Kaufman diagonal pivoting
       method.  Supports inplace and threading.  Uses sytrf and sytri or csytrf and csytri from
       Lapack and returns "inverse, info" in array context.

        (PDL, (PDL(info))) = msyminv(PDL, SCALAR|PDL(uplo))
        uplo : UPPER = 0 | LOWER = 1, default = 0

        # Assume $a is symmetric
        my $a = random(10,10);
        my $inv = msyminv($a);

   mposinv
       Computes inverse of a symmetric positive definite square matrix using Cholesky
       factorization.  Supports inplace and threading.  Uses potrf and potri or cpotrf and cpotri
       from Lapack and returns "inverse, info" in array context.

        (PDL, (PDL(info))) = mposinv(PDL, SCALAR|PDL(uplo))
        uplo : UPPER = 0 | LOWER = 1, default = 0

        # Assume $a is symmetric positive definite
        my $a = random(10,10);
        $a = $a->crossprod($a);
        my $inv = mposinv($a);

   mpinv
       Computes pseudo-inverse (Moore-Penrose) of a general matrix.  Works on transposed array.

        PDL(pseudo-inv)  = mpinv(PDL, SCALAR(tol))
        TOL:   tolerance value, default : mnorm(dims(PDL),'inf') * mnorm(PDL) * EPS

        my $a = random(5,10);
        my $inv = mpinv($a);

   mlu
       Computes LU factorization.  Uses getrf or cgetrf from Lapack and returns L, U, pivot and
       info.  Works on transposed array.

        (PDL(l), PDL(u), PDL(pivot), PDL(info)) = mlu(PDL)

        my $a = random(10,10);
        ($l, $u, $pivot, $info) = mlu($a);

   mchol
       Computes Cholesky decomposition of a symmetric matrix also knows as symmetric square root.
       If inplace flag is set, overwrite  the leading upper or lower triangular part of A else
       returns triangular matrix. Returns "cholesky, info" in array context.  Supports threading.
       Uses potrf or cpotrf from Lapack.

        PDL(Cholesky) = mchol(PDL, SCALAR)
        SCALAR : UPPER = 0 | LOWER = 1, default = 0

        my $a = random(10,10);
        $a = crossprod($a, $a);
        my $u  = mchol($a);

   mhessen
       Reduces a square matrix to Hessenberg form H and orthogonal matrix Q.

       It reduces a general matrix A to upper Hessenberg form H by an orthogonal similarity
       transformation:

               Q' x A x Q = H

       or

               A = Q x H x Q'

       Uses gehrd and orghr or cgehrd and cunghr from Lapack and returns "H" in scalar context
       else "H" and "Q".  Works on transposed array.

        (PDL(h), (PDL(q))) = mhessen(PDL)

        my $a = random(10,10);
        ($h, $q) = mhessen($a);

   mschur
       Computes Schur form, works inplace.

               A = Z x T x Z'

       Supports threading for unordered eigenvalues.  Uses gees or cgees from Lapack and returns
       schur(T) in scalar context.  Works on tranposed array(s).

        ( PDL(schur), (PDL(eigenvalues), (PDL(left schur vectors), PDL(right schur vectors), $sdim), $info) ) = mschur(PDL(A), SCALAR(schur vector),SCALAR(left eigenvector), SCALAR(right eigenvector),SCALAR(select_func), SCALAR(backtransform), SCALAR(norm))
        schur vector        : Schur vectors returned, none = 0 | all = 1 | selected = 2, default = 0
        left eigenvector    : Left eigenvectors returned, none = 0 | all = 1 | selected = 2, default = 0
        right eigenvector   : Right eigenvectors returned, none = 0 | all = 1 | selected = 2, default = 0
        select_func         : Select_func is used to select eigenvalues to sort
                              to the top left of the Schur form.
                              An eigenvalue is selected if PerlInt select_func(PDL::Complex(w)) is true;
                              Note that a selected complex eigenvalue may no longer
                              satisfy select_func(PDL::Complex(w)) = 1 after ordering, since
                              ordering may change the value of complex eigenvalues
                              (especially if the eigenvalue is ill-conditioned).
                              All eigenvalues/vectors are selected if select_func is undefined.
        backtransform       : Whether or not backtransforms eigenvectors to those of A.
                              Only supported if schur vectors are computed, default = 1.
        norm                : Whether or not computed eigenvectors are normalized to have Euclidean norm equal to
                              1 and largest component real, default = 1

        Returned values     :
                              Schur form T (SCALAR CONTEXT),
                              eigenvalues,
                              Schur vectors (Z) if requested,
                              left eigenvectors if requested
                              right eigenvectors if requested
                              sdim: Number of eigenvalues selected if select_func is defined.
                              info: Info output from gees/cgees.

        my $a = random(10,10);
        my $schur  = mschur($a);
        sub select{
               my $m = shift;
               # select "discrete time" eigenspace
               return $m->Cabs < 1 ? 1 : 0;
        }
        my ($schur,$eigen, $svectors,$evectors)  = mschur($a,1,1,0,\&select);

   mschurx
       Computes Schur form, works inplace.  Uses geesx or cgeesx from Lapack and returns schur(T)
       in scalar context.  Works on transposed array.

        ( PDL(schur) (,PDL(eigenvalues))  (, PDL(schur vectors), HASH(result)) ) = mschurx(PDL, SCALAR(schur vector), SCALAR(left eigenvector), SCALAR(right eigenvector),SCALAR(select_func), SCALAR(sense), SCALAR(backtransform), SCALAR(norm))
        schur vector        : Schur vectors returned, none = 0 | all = 1 | selected = 2, default = 0
        left eigenvector    : Left eigenvectors returned, none = 0 | all = 1 | selected = 2, default = 0
        right eigenvector   : Right eigenvectors returned, none = 0 | all = 1 | selected = 2, default = 0
        select_func         : Select_func is used to select eigenvalues to sort
                              to the top left of the Schur form.
                              An eigenvalue is selected if PerlInt select_func(PDL::Complex(w)) is true;
                              Note that a selected complex eigenvalue may no longer
                              satisfy select_func(PDL::Complex(w)) = 1 after ordering, since
                              ordering may change the value of complex eigenvalues
                              (especially if the eigenvalue is ill-conditioned).
                              All  eigenvalues/vectors are selected if select_func is undefined.
        sense               : Determines which reciprocal condition numbers will be computed.
                               0: None are computed
                               1: Computed for average of selected eigenvalues only
                               2: Computed for selected right invariant subspace only
                               3: Computed for both
                               If select_func is undefined, sense is not used.
        backtransform       : Whether or not backtransforms eigenvectors to those of A.
                              Only supported if schur vector are computed, default = 1
        norm                : Whether or not computed eigenvectors are normalized to have Euclidean norm equal to
                              1 and largest component real, default = 1

        Returned values     :
                              Schur form T (SCALAR CONTEXT),
                              eigenvalues,
                              Schur vectors if requested,
                              HASH{VL}: left eigenvectors if requested
                              HASH{VR}: right eigenvectors if requested
                              HASH{info}: info output from gees/cgees.
                              if select_func is defined:
                               HASH{n}: number of eigenvalues selected,
                               HASH{rconde}: reciprocal condition numbers for the average of
                               the selected eigenvalues if requested,
                               HASH{rcondv}: reciprocal condition numbers for the selected
                               right invariant subspace if requested.

        my $a = random(10,10);
        my $schur  = mschurx($a);
        sub select{
               my $m = shift;
               # select "discrete time" eigenspace
               return $m->Cabs < 1 ? 1 : 0;
        }
        my ($schur,$eigen, $vectors,%ret)  = mschurx($a,1,0,0,\&select);

   mgschur
       Computes generalized Schur decomposition of the pair (A,B).

               A = Q x S x Z'
               B = Q x T x Z'

       Uses gges or cgges from Lapack.  Works on transposed array.

        ( PDL(schur S), PDL(schur T), PDL(alpha), PDL(beta), HASH{result}) = mgschur(PDL(A), PDL(B), SCALAR(left schur vector),SCALAR(right schur vector),SCALAR(left eigenvector), SCALAR(right eigenvector), SCALAR(select_func), SCALAR(backtransform), SCALAR(scale))
        left schur vector   : Left Schur vectors returned, none = 0 | all = 1 | selected = 2, default = 0
        right schur vector  : Right Schur vectors returned, none = 0 | all = 1 | selected = 2, default = 0
        left eigenvector    : Left eigenvectors returned, none = 0 | all = 1 | selected = 2, default = 0
        right eigenvector   : Right eigenvectors returned, none = 0 | all = 1 | selected = 2, default = 0
        select_func         : Select_func is used to select eigenvalues to sort.
                              to the top left of the Schur form.
                              An eigenvalue w = wr(j)+sqrt(-1)*wi(j) is selected if
                              PerlInt select_func(PDL::Complex(alpha),PDL | PDL::Complex (beta)) is true;
                              Note that a selected complex eigenvalue may no longer
                              satisfy select_func = 1 after ordering, since
                              ordering may change the value of complex eigenvalues
                              (especially if the eigenvalue is ill-conditioned).
                              All eigenvalues/vectors are selected if select_func is undefined.
        backtransform       : Whether or not backtransforms eigenvectors to those of (A,B).
                              Only supported if right and/or left schur vector are computed,
        scale               : Whether or not computed eigenvectors are scaled so the largest component
                              will have abs(real part) + abs(imag. part) = 1, default = 1

        Returned values     :
                              Schur form S,
                              Schur form T,
                              alpha,
                              beta (eigenvalues = alpha/beta),
                              HASH{info}: info output from gges/cgges.
                              HASH{SL}: left Schur vectors if requested
                              HASH{SR}: right Schur vectors if requested
                              HASH{VL}: left eigenvectors if requested
                              HASH{VR}: right eigenvectors if requested
                              HASH{n} : Number of eigenvalues selected if select_func is defined.

        my $a = random(10,10);
        my $b = random(10,10);
        my ($S,$T) = mgschur($a,$b);
        sub select{
               my ($alpha,$beta) = @_;
               return $alpha->Cabs < abs($beta) ? 1 : 0;
        }
        my ($S, $T, $alpha, $beta, %res)  = mgschur( $a, $b, 1, 1, 1, 1,\&select);

   mgschurx
       Computes generalized Schur decomposition of the pair (A,B).

               A = Q x S x Z'
               B = Q x T x Z'

       Uses ggesx or cggesx from Lapack. Works on transposed array.

        ( PDL(schur S), PDL(schur T), PDL(alpha), PDL(beta), HASH{result}) = mgschurx(PDL(A), PDL(B), SCALAR(left schur vector),SCALAR(right schur vector),SCALAR(left eigenvector), SCALAR(right eigenvector), SCALAR(select_func), SCALAR(sense), SCALAR(backtransform), SCALAR(scale))
        left schur vector   : Left Schur vectors returned, none = 0 | all = 1 | selected = 2, default = 0
        right schur vector  : Right Schur vectors returned, none = 0 | all = 1 | selected = 2, default = 0
        left eigenvector    : Left eigenvectors returned, none = 0 | all = 1 | selected = 2, default = 0
        right eigenvector   : Right eigenvectors returned, none = 0 | all = 1 | selected = 2, default = 0
        select_func         : Select_func is used to select eigenvalues to sort.
                              to the top left of the Schur form.
                              An eigenvalue w = wr(j)+sqrt(-1)*wi(j) is selected if
                              PerlInt select_func(PDL::Complex(alpha),PDL | PDL::Complex (beta)) is true;
                              Note that a selected complex eigenvalue may no longer
                              satisfy select_func = 1 after ordering, since
                              ordering may change the value of complex eigenvalues
                              (especially if the eigenvalue is ill-conditioned).
                              All eigenvalues/vectors are selected if select_func is undefined.
        sense               : Determines which reciprocal condition numbers will be computed.
                               0: None are computed
                               1: Computed for average of selected eigenvalues only
                               2: Computed for selected deflating subspaces only
                               3: Computed for both
                               If select_func is undefined, sense is not used.

        backtransform       : Whether or not backtransforms eigenvectors to those of (A,B).
                              Only supported if right and/or left schur vector are computed, default = 1
        scale               : Whether or not computed eigenvectors are scaled so the largest component
                              will have abs(real part) + abs(imag. part) = 1, default = 1

        Returned values     :
                              Schur form S,
                              Schur form T,
                              alpha,
                              beta (eigenvalues = alpha/beta),
                              HASH{info}: info output from gges/cgges.
                              HASH{SL}: left Schur vectors if requested
                              HASH{SR}: right Schur vectors if requested
                              HASH{VL}: left eigenvectors if requested
                              HASH{VR}: right eigenvectors if requested
                              HASH{rconde}: reciprocal condition numbers for average of selected eigenvalues if requested
                              HASH{rcondv}: reciprocal condition numbers for selected deflating subspaces if requested
                              HASH{n} : Number of eigenvalues selected if select_func is defined.

        my $a = random(10,10);
        my $b = random(10,10);
        my ($S,$T) = mgschurx($a,$b);
        sub select{
               my ($alpha,$beta) = @_;
               return $alpha->Cabs < abs($beta) ? 1 : 0;
        }
        my ($S, $T, $alpha, $beta, %res)  = mgschurx( $a, $b, 1, 1, 1, 1,\&select,3);

   mqr
       Computes QR decomposition.  For complex number needs object of type PDL::Complex.  Uses
       geqrf and orgqr or cgeqrf and cungqr from Lapack and returns "Q" in scalar context. Works
       on transposed array.

        (PDL(Q), PDL(R), PDL(info)) = mqr(PDL, SCALAR)
        SCALAR : ECONOMIC = 0 | FULL = 1, default = 0

        my $a = random(10,10);
        my ( $q, $r )  = mqr($a);
        # Can compute full decomposition if nrow > ncol
        $a = random(5,7);
        ( $q, $r )  = $a->mqr(1);

   mrq
       Computes RQ decomposition.  For complex number needs object of type PDL::Complex.  Uses
       gerqf and orgrq or cgerqf and cungrq from Lapack and returns "Q" in scalar context. Works
       on transposed array.

        (PDL(R), PDL(Q), PDL(info)) = mrq(PDL, SCALAR)
        SCALAR : ECONOMIC = 0 | FULL = 1, default = 0

        my $a = random(10,10);
        my ( $r, $q )  = mrq($a);
        # Can compute full decomposition if nrow < ncol
        $a = random(5,7);
        ( $r, $q )  = $a->mrq(1);

   mql
       Computes QL decomposition.  For complex number needs object of type PDL::Complex.  Uses
       geqlf and orgql or cgeqlf and cungql from Lapack and returns "Q" in scalar context. Works
       on transposed array.

        (PDL(Q), PDL(L), PDL(info)) = mql(PDL, SCALAR)
        SCALAR : ECONOMIC = 0 | FULL = 1, default = 0

        my $a = random(10,10);
        my ( $q, $l )  = mql($a);
        # Can compute full decomposition if nrow > ncol
        $a = random(5,7);
        ( $q, $l )  = $a->mql(1);

   mlq
       Computes LQ decomposition.  For complex number needs object of type PDL::Complex.  Uses
       gelqf and orglq or cgelqf and cunglq from Lapack and returns "Q" in scalar context. Works
       on transposed array.

        ( PDL(L), PDL(Q), PDL(info) ) = mlq(PDL, SCALAR)
        SCALAR : ECONOMIC = 0 | FULL = 1, default = 0

        my $a = random(10,10);
        my ( $l, $q )  = mlq($a);
        # Can compute full decomposition if nrow < ncol
        $a = random(5,7);
        ( $l, $q )  = $a->mlq(1);

   msolve
       Solves linear system of equations using LU decomposition.

               A * X = B

       Returns X in scalar context else X, LU, pivot vector and info.  B is overwritten by X if
       its inplace flag is set.  Supports threading.  Uses gesv or cgesv from Lapack.  Works on
       transposed arrays.

        (PDL(X), (PDL(LU), PDL(pivot), PDL(info))) = msolve(PDL(A), PDL(B) )

        my $a = random(5,5);
        my $b = random(10,5);
        my $X = msolve($a, $b);

   msolvex
       Solves linear system of equations using LU decomposition.

               A * X = B

       Can optionnally equilibrate the matrix.  Uses gesvx or cgesvx from Lapack.  Works on
       transposed arrays.

        (PDL, (HASH(result))) = msolvex(PDL(A), PDL(B), HASH(options))
        where options are:
        transpose:     solves A' * X = B
                       0: false
                       1: true
        equilibrate:   equilibrates A if necessary.
                       form equilibration is returned in HASH{'equilibration'}:
                               0: no equilibration
                               1: row equilibration
                               2: column equilibration
                       row scale factors are returned in HASH{'row'}
                       column scale factors are returned in HASH{'column'}
                       0: false
                       1: true
        LU:            returns lu decomposition in HASH{LU}
                       0: false
                       1: true
        A:             returns scaled A if equilibration was done in HASH{A}
                       0: false
                       1: true
        B:             returns scaled B if equilibration was done in HASH{B}
                       0: false
                       1: true
        Returned values:
                       X (SCALAR CONTEXT),
                       HASH{'pivot'}:
                        Pivot indice from LU factorization
                       HASH{'rcondition'}:
                        Reciprocal condition of the matrix
                       HASH{'ferror'}:
                        Forward error bound
                       HASH{'berror'}:
                        Componentwise relative backward error
                       HASH{'rpvgrw'}:
                        Reciprocal pivot growth factor
                       HASH{'info'}:
                        Info: output from gesvx

        my $a = random(10,10);
        my $b = random(5,10);
        my %options = (
                       LU=>1,
                       equilibrate => 1,
                       );
        my( $X, %result) = msolvex($a,$b,%options);

   mtrisolve
       Solves linear system of equations with triangular matrix A.

               A * X = B  or A' * X = B

       B is overwritten by X if its inplace flag is set.  Supports threading.  Uses trtrs or
       ctrtrs from Lapack.  Work on transposed array(s).

        (PDL(X), (PDL(info)) = mtrisolve(PDL(A), SCALAR(uplo), PDL(B), SCALAR(trans), SCALAR(diag))
        uplo   : UPPER  = 0 | LOWER = 1
        trans  : NOTRANSPOSE  = 0 | TRANSPOSE = 1, default = 0
        uplo   : UNITARY DIAGONAL = 1, default = 0

        # Assume $a is upper triagonal
        my $a = random(5,5);
        my $b = random(5,10);
        my $X = mtrisolve($a, 0, $b);

   msymsolve
       Solves linear system of equations using diagonal pivoting method with symmetric matrix A.

               A * X = B

       Returns X in scalar context else X, block diagonal matrix D (and the multipliers), pivot
       vector an info. B is overwritten by X if its inplace flag is set.  Supports threading.
       Uses sysv or csysv from Lapack.  Works on transposed array(s).

        (PDL(X), ( PDL(D), PDL(pivot), PDL(info) ) ) = msymsolve(PDL(A), SCALAR(uplo), PDL(B) )
        uplo : UPPER  = 0 | LOWER = 1, default = 0

        # Assume $a is symmetric
        my $a = random(5,5);
        my $b = random(5,10);
        my $X = msymsolve($a, 0, $b);

   msymsolvex
       Solves linear system of equations using diagonal pivoting method with symmetric matrix A.

               A * X = B

       Uses sysvx or csysvx from Lapack. Works on transposed array.

        (PDL, (HASH(result))) = msymsolvex(PDL(A), SCALAR (uplo), PDL(B), SCALAR(d))
        uplo : UPPER  = 0 | LOWER = 1, default = 0
        d    : whether return diagonal matrix d and pivot vector
               FALSE  = 0 | TRUE = 1, default = 0
        Returned values:
                       X (SCALAR CONTEXT),
                       HASH{'D'}:
                        Block diagonal matrix D (and the multipliers) (if requested)
                       HASH{'pivot'}:
                        Pivot indice from LU factorization (if requested)
                       HASH{'rcondition'}:
                        Reciprocal condition of the matrix
                       HASH{'ferror'}:
                        Forward error bound
                       HASH{'berror'}:
                        Componentwise relative backward error
                       HASH{'info'}:
                        Info: output from sysvx

        # Assume $a is symmetric
        my $a = random(10,10);
        my $b = random(5,10);
        my ($X, %result) = msolvex($a, 0, $b);

   mpossolve
       Solves linear system of equations using Cholesky decomposition with symmetric positive
       definite matrix A.

               A * X = B

       Returns X in scalar context else X, U or L and info.  B is overwritten by X if its inplace
       flag is set.  Supports threading.  Uses posv or cposv from Lapack.  Works on transposed
       array(s).

        (PDL, (PDL, PDL, PDL)) = mpossolve(PDL(A), SCALAR(uplo), PDL(B) )
        uplo : UPPER  = 0 | LOWER = 1, default = 0

        # asume $a is symmetric positive definite
        my $a = random(5,5);
        my $b = random(5,10);
        my $X = mpossolve($a, 0, $b);

   mpossolvex
       Solves linear system of equations using Cholesky decomposition with symmetric positive
       definite matrix A

               A * X = B

       Can optionnally equilibrate the matrix.  Uses posvx or cposvx from Lapack.  Works on
       transposed array(s).

        (PDL, (HASH(result))) = mpossolvex(PDL(A), SCARA(uplo), PDL(B), HASH(options))
        uplo : UPPER  = 0 | LOWER = 1, default = 0
        where options are:
        equilibrate:   equilibrates A if necessary.
                       form equilibration is returned in HASH{'equilibration'}:
                               0: no equilibration
                               1: equilibration
                       scale factors are returned in HASH{'scale'}
                       0: false
                       1: true
        U|L:           returns Cholesky factorization in HASH{U} or HASH{L}
                       0: false
                       1: true
        A:             returns scaled A if equilibration was done in HASH{A}
                       0: false
                       1: true
        B:             returns scaled B if equilibration was done in HASH{B}
                       0: false
                       1: true
        Returned values:
                       X (SCALAR CONTEXT),
                       HASH{'rcondition'}:
                        Reciprocal condition of the matrix
                       HASH{'ferror'}:
                        Forward error bound
                       HASH{'berror'}:
                        Componentwise relative backward error
                       HASH{'info'}:
                        Info: output from gesvx

        # Assume $a is symmetric positive definite
        my $a = random(10,10);
        my $b = random(5,10);
        my %options = (U=>1,
                       equilibrate => 1,
                       );
        my ($X, %result) = msolvex($a, 0, $b,%opt);

   mlls
       Solves overdetermined or underdetermined real linear systems using QR or LQ factorization.

       If M > N in the M-by-N matrix A, returns the residual sum of squares too.  Uses gels or
       cgels from Lapack.  Works on transposed arrays.

        PDL(X) = mlls(PDL(A), PDL(B), SCALAR(trans))
        trans : NOTRANSPOSE  = 0 | TRANSPOSE/CONJUGATE = 1, default = 0

        $a = random(4,5);
        $b = random(3,5);
        ($x, $res) = mlls($a, $b);

   mllsy
       Computes the minimum-norm solution to a real linear least squares problem using a complete
       orthogonal factorization.

       Uses gelsy or cgelsy from Lapack. Works on tranposed arrays.

        ( PDL(X), ( HASH(result) ) ) = mllsy(PDL(A), PDL(B))
        Returned values:
                       X (SCALAR CONTEXT),
                       HASH{'A'}:
                        complete orthogonal factorization of A
                       HASH{'jpvt'}:
                        details of columns interchanges
                       HASH{'rank'}:
                        effective rank of A

        my $a = random(10,10);
        my $b = random(10,10);
        $X = mllsy($a, $b);

   mllss
       Computes the minimum-norm solution to a real linear least squares problem using a singular
       value decomposition.

       Uses gelss or gelsd from Lapack.  Works on transposed arrays.

        ( PDL(X), ( HASH(result) ) )= mllss(PDL(A), PDL(B), SCALAR(method))
        method: specifie which method to use (see Lapack for further details)
               '(c)gelss' or '(c)gelsd', default = '(c)gelsd'
        Returned values:
                       X (SCALAR CONTEXT),
                       HASH{'V'}:
                        if method = (c)gelss, the right singular vectors, stored columnwise
                       HASH{'s'}:
                        singular values from SVD
                       HASH{'res'}:
                        if A has full rank the residual sum-of-squares for the solution
                       HASH{'rank'}:
                        effective rank of A
                       HASH{'info'}:
                        info output from method

        my $a = random(10,10);
        my $b = random(10,10);
        $X = mllss($a, $b);

   mglm
       Solves a general Gauss-Markov Linear Model (GLM) problem.  Supports threading.  Uses ggglm
       or cggglm from Lapack. Works on transposed arrays.

        (PDL(x), PDL(y)) = mglm(PDL(a), PDL(b), PDL(d))
        where d is the left hand side of the GLM equation

        my $a = random(8,10);
        my $b = random(7,10);
        my $d = random(10);
        my ($x, $y) = mglm($a, $b, $d);

   mlse
       Solves a linear equality-constrained least squares (LSE) problem.  Uses gglse or cgglse
       from Lapack. Works on transposed arrays.

        (PDL(x), PDL(res2)) = mlse(PDL(a), PDL(b), PDL(c), PDL(d))
        where
        c      : The right hand side vector for the
                 least squares part of the LSE problem.
        d      : The right hand side vector for the
                 constrained equation.
        x      : The solution of the LSE problem.
        res2   : The residual sum of squares for the solution
                 (returned only in array context)

        my $a = random(5,4);
        my $b = random(5,3);
        my $c = random(4);
        my $d = random(3);
        my ($x, $res2) = mlse($a, $b, $c, $d);

   meigen
       Computes eigenvalues and, optionally, the left and/or right eigenvectors of a general
       square matrix (spectral decomposition).  Eigenvectors are normalized (Euclidean norm = 1)
       and largest component real.  The eigenvalues and eigenvectors returned are object of type
       PDL::Complex.  If only eigenvalues are requested, info is returned in array context.
       Supports threading.  Uses geev or cgeev from Lapack.  Works on transposed arrays.

        (PDL(values), (PDL(LV),  (PDL(RV)), (PDL(info))) = meigen(PDL, SCALAR(left vector), SCALAR(right vector))
        left vector  : FALSE = 0 | TRUE = 1, default = 0
        right vector : FALSE = 0 | TRUE = 1, default = 0

        my $a = random(10,10);
        my ( $eigenvalues, $left_eigenvectors, $right_eigenvectors )  = meigen($a,1,1);

   meigenx
       Computes eigenvalues, one-norm and, optionally, the left and/or right eigenvectors of a
       general square matrix (spectral decomposition).  Eigenvectors are normalized (Euclidean
       norm = 1) and largest component real.  The eigenvalues and eigenvectors returned are
       object of type PDL::Complex.  Uses geevx or cgeevx from Lapack.  Works on transposed
       arrays.

        (PDL(value), (PDL(lv),  (PDL(rv)), HASH(result)), HASH(result)) = meigenx(PDL, HASH(options))
        where options are:
        vector:     eigenvectors to compute
                       'left':  computes left eigenvectors
                       'right': computes right eigenvectors
                       'all':   computes left and right eigenvectors
                        0:     doesn't compute (default)
        rcondition: reciprocal condition numbers to compute (returned in HASH{'rconde'} for eigenvalues and HASH{'rcondv'} for eigenvectors)
                       'value':  computes reciprocal condition numbers for eigenvalues
                       'vector': computes reciprocal condition numbers for eigenvectors
                       'all':    computes reciprocal condition numbers for eigenvalues and eigenvectors
                        0:      doesn't compute (default)
        error:      specifie whether or not it computes the error bounds (returned in HASH{'eerror'} and HASH{'verror'})
                    error bound = EPS * One-norm / rcond(e|v)
                    (reciprocal condition numbers for eigenvalues or eigenvectors must be computed).
                       1: returns error bounds
                       0: not computed
        scale:      specifie whether or not it diagonaly scales the entry matrix
                    (scale details returned in HASH : 'scale')
                       1: scales
                       0: Doesn't scale (default)
        permute:    specifie whether or not it permutes row and columns
                    (permute details returned in HASH{'balance'})
                       1: permutes
                       0: Doesn't permute (default)
        schur:      specifie whether or not it returns the Schur form (returned in HASH{'schur'})
                       1: returns Schur form
                       0: not returned
        Returned values:
                   eigenvalues (SCALAR CONTEXT),
                   left eigenvectors if requested,
                   right eigenvectors if requested,
                   HASH{'norm'}:
                       One-norm of the matrix
                   HASH{'info'}:
                       Info: if > 0, the QR algorithm failed to compute all the eigenvalues
                       (see syevx for further details)

        my $a = random(10,10);
        my %options = ( rcondition => 'all',
                    vector => 'all',
                    error => 1,
                    scale => 1,
                    permute=>1,
                    shur => 1
                    );
        my ( $eigenvalues, $left_eigenvectors, $right_eigenvectors, %result)  = meigenx($a,%options);
        print "Error bounds for eigenvalues:\n $eigenvalues\n are:\n". transpose($result{'eerror'}) unless $info;

   mgeigen
       Computes generalized eigenvalues and, optionally, the left and/or right generalized
       eigenvectors for a pair of N-by-N real nonsymmetric matrices (A,B) .  The alpha from ratio
       alpha/beta is object of type PDL::Complex.  Supports threading. Uses ggev or cggev from
       Lapack.  Works on transposed arrays.

        ( PDL(alpha), PDL(beta), ( PDL(LV),  (PDL(RV) ), PDL(info)) = mgeigen(PDL(A),PDL(B) SCALAR(left vector), SCALAR(right vector))
        left vector  : FALSE = 0 | TRUE = 1, default = 0
        right vector : FALSE = 0 | TRUE = 1, default = 0

        my $a = random(10,10);
        my $b = random(10,10);
        my ( $alpha, $beta, $left_eigenvectors, $right_eigenvectors )  = mgeigen($a, $b,1, 1);

   mgeigenx
       Computes generalized eigenvalues, one-norms and, optionally, the left and/or right
       generalized eigenvectors for a pair of N-by-N real nonsymmetric matrices (A,B).  The alpha
       from ratio alpha/beta is object of type PDL::Complex.  Uses ggevx or cggevx from Lapack.
       Works on transposed arrays.

        (PDL(alpha), PDL(beta), PDL(lv),  PDL(rv), HASH(result) ) = mgeigenx(PDL(a), PDL(b), HASH(options))
        where options are:
        vector:     eigenvectors to compute
                       'left':  computes left eigenvectors
                       'right': computes right eigenvectors
                       'all':   computes left and right eigenvectors
                        0:     doesn't compute (default)
        rcondition: reciprocal condition numbers to compute (returned in HASH{'rconde'} for eigenvalues and HASH{'rcondv'} for eigenvectors)
                       'value':  computes reciprocal condition numbers for eigenvalues
                       'vector': computes reciprocal condition numbers for eigenvectors
                       'all':    computes reciprocal condition numbers for eigenvalues and eigenvectors
                        0:      doesn't compute (default)
        error:      specifie whether or not it computes the error bounds (returned in HASH{'eerror'} and HASH{'verror'})
                    error bound = EPS * sqrt(one-norm(a)**2 + one-norm(b)**2) / rcond(e|v)
                    (reciprocal condition numbers for eigenvalues or eigenvectors must be computed).
                       1: returns error bounds
                       0: not computed
        scale:      specifie whether or not it diagonaly scales the entry matrix
                    (scale details returned in HASH : 'lscale' and 'rscale')
                       1: scales
                       0: doesn't scale (default)
        permute:    specifie whether or not it permutes row and columns
                    (permute details returned in HASH{'balance'})
                       1: permutes
                       0: Doesn't permute (default)
        schur:      specifie whether or not it returns the Schur forms (returned in HASH{'aschur'} and HASH{'bschur'})
                    (right or left eigenvectors must be computed).
                       1: returns Schur forms
                       0: not returned
        Returned values:
                   alpha,
                   beta,
                   left eigenvectors if requested,
                   right eigenvectors if requested,
                   HASH{'anorm'}, HASH{'bnorm'}:
                       One-norm of the matrix A and B
                   HASH{'info'}:
                       Info: if > 0, the QR algorithm failed to compute all the eigenvalues
                       (see syevx for further details)

        $a = random(10,10);
        $b = random(10,10);
        %options = (rcondition => 'all',
                    vector => 'all',
                    error => 1,
                    scale => 1,
                    permute=>1,
                    shur => 1
                    );
        ($alpha, $beta, $left_eigenvectors, $right_eigenvectors, %result)  = mgeigenx($a, $b,%options);
        print "Error bounds for eigenvalues:\n $eigenvalues\n are:\n". transpose($result{'eerror'}) unless $info;

   msymeigen
       Computes eigenvalues and, optionally eigenvectors of a real symmetric square or complex
       Hermitian matrix (spectral decomposition).  The eigenvalues are computed from lower or
       upper triangular matrix.  If only eigenvalues are requested, info is returned in array
       context.  Supports threading and works inplace if eigenvectors are requested.  From
       Lapack, uses syev or syevd for real and cheev or cheevd for complex.  Works on transposed
       array(s).

        (PDL(values), (PDL(VECTORS)), PDL(info)) = msymeigen(PDL, SCALAR(uplo), SCALAR(vector), SCALAR(method))
        uplo : UPPER  = 0 | LOWER = 1, default = 0
        vector : FALSE = 0 | TRUE = 1, default = 0
        method : 'syev' | 'syevd' | 'cheev' | 'cheevd', default = 'syevd'|'cheevd'

        # Assume $a is symmetric
        my $a = random(10,10);
        my ( $eigenvalues, $eigenvectors )  = msymeigen($a,0,1, 'syev');

   msymeigenx
       Computes eigenvalues and, optionally eigenvectors of a symmetric square matrix (spectral
       decomposition).  The eigenvalues are computed from lower or upper triangular matrix and
       can be selected by specifying a range. From Lapack, uses syevx or syevr for real and
       cheevx or cheevr for complex. Works on transposed arrays.

        (PDL(value), (PDL(vector)), PDL(n), PDL(info), (PDL(support)) ) = msymeigenx(PDL, SCALAR(uplo), SCALAR(vector), HASH(options))
        uplo : UPPER  = 0 | LOWER = 1, default = 0
        vector : FALSE = 0 | TRUE = 1, default = 0
        where options are:
        range_type:    method for selecting eigenvalues
                       indice:  range of indices
                       interval: range of values
                       0: find all eigenvalues and optionally all vectors
        range:         PDL(2), lower and upper bounds interval or smallest and largest indices
                       1<=range<=N for indice
        abstol:        specifie error tolerance for eigenvalues
        method:        specifie which method to use (see Lapack for further details)
                       'syevx' (default)
                       'syevr'
                       'cheevx' (default)
                       'cheevr'
        Returned values:
                       eigenvalues (SCALAR CONTEXT),
                       eigenvectors if requested,
                       total number of eigenvalues found (n),
                       info
                       issupz or ifail (support) according to method used and returned info,
                       for (sy|che)evx returns support only if info != 0

        # Assume $a is symmetric
        my $a = random(10,10);
        my $overflow = lamch(9);
        my $range = cat pdl(0),$overflow;
        my $abstol = pdl(1.e-5);
        my %options = (range_type=>'interval',
                       range => $range,
                       abstol => $abstol,
                       method=>'syevd');
        my ( $eigenvalues, $eigenvectors, $n, $isuppz )  = msymeigenx($a,0,1, %options);

   msymgeigen
       Computes eigenvalues and, optionally eigenvectors of a real generalized symmetric-definite
       or Hermitian-definite eigenproblem.  The eigenvalues are computed from lower or upper
       triangular matrix If only eigenvalues are requested, info is returned in array context.
       Supports threading. From Lapack, uses sygv or sygvd for real or chegv or chegvd for
       complex.  Works on transposed array(s).

        (PDL(values), (PDL(vectors)), PDL(info)) = msymgeigen(PDL(a), PDL(b),SCALAR(uplo), SCALAR(vector), SCALAR(type), SCALAR(method))
        uplo : UPPER  = 0 | LOWER = 1, default = 0
        vector : FALSE = 0 | TRUE = 1, default = 0
        type :
               1: A * x = (lambda) * B * x
               2: A * B * x = (lambda) * x
               3: B * A * x = (lambda) * x
               default = 1
        method : 'sygv' | 'sygvd' for real or  ,'chegv' | 'chegvd' for complex,  default = 'sygvd' | 'chegvd'

        # Assume $a is symmetric
        my $a = random(10,10);
        my $b = random(10,10);
        $b = $b->crossprod($b);
        my ( $eigenvalues, $eigenvectors )  = msymgeigen($a, $b, 0, 1, 1, 'sygv');

   msymgeigenx
       Computes eigenvalues and, optionally eigenvectors of a real generalized symmetric-definite
       or Hermitian eigenproblem.  The eigenvalues are computed from lower or upper triangular
       matrix and can be selected by specifying a range. Uses sygvx or cheevx from Lapack. Works
       on transposed arrays.

        (PDL(value), (PDL(vector)), PDL(info), PDL(n), (PDL(support)) ) = msymeigenx(PDL(a), PDL(b), SCALAR(uplo), SCALAR(vector), HASH(options))
        uplo : UPPER  = 0 | LOWER = 1, default = 0
        vector : FALSE = 0 | TRUE = 1, default = 0
        where options are:
        type :         Specifies the problem type to be solved
                       1: A * x = (lambda) * B * x
                       2: A * B * x = (lambda) * x
                       3: B * A * x = (lambda) * x
                       default = 1
        range_type:    method for selecting eigenvalues
                       indice:  range of indices
                       interval: range of values
                       0: find all eigenvalues and optionally all vectors
        range:         PDL(2), lower and upper bounds interval or smallest and largest indices
                       1<=range<=N for indice
        abstol:        specifie error tolerance for eigenvalues
        Returned values:
                       eigenvalues (SCALAR CONTEXT),
                       eigenvectors if requested,
                       total number of eigenvalues found (n),
                       info
                       ifail according to returned info (support).

        # Assume $a is symmetric
        my $a = random(10,10);
        my $b = random(10,10);
        $b = $b->crossprod($b);
        my $overflow = lamch(9);
        my $range = cat pdl(0),$overflow;
        my $abstol = pdl(1.e-5);
        my %options = (range_type=>'interval',
                       range => $range,
                       abstol => $abstol,
                       type => 1);
        my ( $eigenvalues, $eigenvectors, $n, $isuppz )  = msymgeigenx($a, $b, 0,1, %options);

   mdsvd
       Computes SVD using Coppen's divide and conquer algorithm.  Return singular values in
       scalar context else left (U), singular values, right (V' (hermitian for complex)) singular
       vectors and info.  Supports threading.  If only singulars values are requested, info is
       only returned in array context.  Uses gesdd or cgesdd from Lapack.

        (PDL(U), (PDL(s), PDL(V)), PDL(info)) = mdsvd(PDL, SCALAR(job))
        job :  0 = computes only singular values
               1 = computes full SVD (square U and V)
               2 = computes SVD (singular values, right and left singular vectors)
               default = 1

        my $a = random(5,10);
        my ($u, $s, $v) = mdsvd($a);

   msvd
       Computes SVD.  Can compute singular values, either U or V or neither.  Return singulars
       values in scalar context else left (U), singular values, right (V' (hermitian for complex)
       singulars vector and info.  Supports threading.  If only singulars values are requested,
       info is returned in array context.  Uses gesvd or cgesvd from Lapack.

        ( (PDL(U)), PDL(s), (PDL(V), PDL(info)) = msvd(PDL, SCALAR(jobu), SCALAR(jobv))
        jobu : 0 = Doesn't compute U
               1 = computes full SVD (square U)
               2 = computes right singular vectors
               default = 1
        jobv : 0 = Doesn't compute V
               1 = computes full SVD (square V)
               2 = computes left singular vectors
               default = 1

        my $a = random(10,10);
        my ($u, $s, $v) = msvd($a);

   mgsvd
       Computes generalized (or quotient) singular value decomposition.  If the effective rank of
       (A',B')' is 0 return only unitary V, U, Q.  For complex number, needs object of type
       PDL::Complex.  Uses ggsvd or cggsvd from Lapack. Works on transposed arrays.

        (PDL(sa), PDL(sb), %ret) = mgsvd(PDL(a), PDL(b), %HASH(options))
        where options are:
        V:    whether or not computes V (boolean, returned in HASH{'V'})
        U:    whether or not computes U (boolean, returned in HASH{'U'})
        Q:    whether or not computes Q (boolean, returned in HASH{'Q'})
        D1:   whether or not computes D1 (boolean, returned in HASH{'D1'})
        D2:   whether or not computes D2 (boolean, returned in HASH{'D2'})
        0R:   whether or not computes 0R (boolean, returned in HASH{'0R'})
        R:    whether or not computes R (boolean, returned in HASH{'R'})
        X:    whether or not computes X (boolean, returned in HASH{'X'})
        all:  whether or not computes all the above.
        Returned value:
                sa,sb          : singular value pairs of A and B (generalized singular values = sa/sb)
                $ret{'rank'}   : effective numerical rank of (A',B')'
                $ret{'info'}   : info from (c)ggsvd

        my $a = random(5,5);
        my $b = random(5,7);
        my ($c, $s, %ret) = mgsvd($a, $b, X => 1);

AUTHOR

       Copyright (C) Grégory Vanuxem 2005-2007.

       This library is free software; you can redistribute it and/or modify it under the terms of
       the artistic license as specified in the Artistic file.