PDL::LinearAlgebra - Linear Algebra utils for PDL
use PDL::LinearAlgebra; $a = random (100,100); ($U, $s, $V) = mdsvd($a);
This module provides a convenient interface to PDL::LinearAlgebra::Real and PDL::LinearAlgebra::Complex.
Set action type when 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.
$a = sequence(5,5); $err = setlaerror(NO); ($inv, $info)= minv($a); if ($info){ # Change the diagonal (the inverse doesn't exist but it's an example) $a->diagonal(0,1)+=1e-8; ($inv, $info)= minv($a); } if ($info){ print "Can't compute the inverse\n"; } else{ print "Inverse of \$a is $inv"; } setlaerror($err);
Get error type.
0 => NO, 1 => WARN, 2 => BARF
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 transpose array(s) and doesn't support dataflow. Supports threading.
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;
Check symmetricity/Hermitianicity of matrix. Supports threading.
Return 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)
Return 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);
Return entry pdl with changed sign by row so that average of positive sign > 0. In other words thread 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;
Compute the cross-product of two matrix: A' x B. If only one matrix is given, take 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);
Compute 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);
Compute 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 = mnrom($a);
Compute 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);
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);
Compute the condition number (two-norm) of a general matrix.
The condition number (two-norm) 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);
Estimate 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.
PDL = mrcond(PDL, SCALAR(ord)) ord : 0 : Infinity norm (default) 1 : One norm
my $a = random(10,10); my $rcond = mrcond($a,1);
Return 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 = random(10,10); my $ortho = morth($a, 1e-8);
Return an orthonormal basis of the null space of matrix A.
PDL = mnull(PDL(A), SCALAR(tol)) tol : tolerance for determining rank, default: 1e-8 for double else 1e-5
my $a = random(10,10); my $null = mnull($a, 1e-8);
Compute inverse of a general square matrix using LU factorization. Supports inplace and threading. Uses getrf and getri or cgetrf and cgetri from Lapack and return inverse, info in array context.
inverse, info
PDL(inv) = minv(PDL)
my $a = random(10,10); my $inv = minv($a);
Compute 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);
Compute 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);
Compute 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);
Compute pseudo-inverse (Moore-Penrose) of a general matrix.
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);
Compute LU factorization. Uses getrf or cgetrf from Lapack and return L, U, pivot and info.
(PDL(l), PDL(u), PDL(pivot), PDL(info)) = mlu(PDL)
my $a = random(10,10); ($l, $u, $pivot, $info) = mlu($a);
Compute 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.
cholesky, info
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);
Reduce 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.
H
Q
(PDL(h), (PDL(q))) = mhessen(PDL)
my $a = random(10,10); ($h, $q) = mhessen($a);
Compute 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.
( 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);
Compute Schur form, works inplace. Uses gees or cgees from Lapack and returns schur(T) in scalar context.
( 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);
Compute 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.
( 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);
Uses ggesx or cggesx from Lapack.
( 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);
Compute 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.
(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 row > col $a = random(5,7); ( $q, $r ) = $a->mqr(1);
Compute 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.
(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 row < col $a = random(5,7); ( $r, $q ) = $a->mrq(1);
Compute 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.
(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 row > col $a = random(5,7); ( $q, $l ) = $a->mql(1);
Compute 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.
( 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 row < col $a = random(5,7); ( $l, $q ) = $a->mlq(1);
Solve linear sytem of equations using LU decomposition.
A * X = B
Returns X in scalar context else X, LU, pivot vector and info. B is overwrited by X if its inplace flag is set. Supports threading. Uses gesv or cgesv from Lapack.
(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);
Can optionnally equilibrate the matrix. Uses gesvx or cgesvx from Lapack.
(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);
Solve linear sytem of equations with triangular matrix A.
A * X = B or A' * X = B
B is overwrited by X if its inplace flag is set. Supports threading. Uses trtrs or ctrtrs from Lapack.
(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);
Solve linear sytem of equations using diagonal pivoting method with symmetric matrix A.
Returns X in scalar context else X, block diagonal matrix D (and the multipliers), pivot vector an info. B is overwrited by X if its inplace flag is set. Supports threading. Uses sysv or csysv from Lapack.
(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);
Uses sysvx or csysvx from Lapack.
(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); mu( $X, %result) = msolvex($a, 0, $b);
Solve linear sytem of equations using Cholesky decomposition with symmetric positive definite matrix A.
Returns X in scalar context else X, U or L and info. B is overwrited by X if its inplace flag is set. Supports threading. Uses posv or cposv from Lapack.
(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);
Solve linear sytem of equations using Cholesky decomposition with symmetric positive definite matrix A
Can optionnally equilibrate the matrix. Uses posvx or cposvx from Lapack.
(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);
Solve 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.
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);
Compute the minimum-norm solution to a real linear least squares problem using a complete orthogonal factorization.
Uses gelsy or cgelsy from Lapack.
( 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);
Compute the minimum-norm solution to a real linear least squares problem using a singular value decomposition.
Uses gelss or gelsd from Lapack.
( 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);
Solve a general Gauss-Markov Linear Model (GLM) problem. Supports threading. Uses ggglm or cggglm from Lapack.
(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);
Solve a linear equality-constrained least squares (LSE) problem. Uses gglse or cgglse from Lapack.
(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);
Compute 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.
(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);
Compute 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.
(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;
Compute 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 are object of type PDL::Complex. Supports threading. Uses ggev or cggev from Lapack.
( 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);
Compute 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 are object of type PDL::Complex. Uses ggevx or cggevx from Lapack.
(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;
Compute 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 work inplace if eigenvectors are requested. From Lapack, uses syev or syevd for real and cheev or cheevd for complex.
(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');
Compute 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.
(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);
Compute 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.
(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');
Compute 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.
(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);
Compute 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);
Compute SVD. Can compute singular value 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);
Compute 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.
(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);
Others things
rectangular diag usage is_inplace and function which modify entry matrix avoid xchg threading support automatically create PDL inplace operation and memory check s after he/she/it and matrix(s) PDL type, verify float/double eig_det qr_det (g)schur(x): if conjugate pair non generalized pb: $seldim ?? (cf: generalized) return conjugate pair if only selected? port to PDL::Matrix
Copyright (C) Grégory Vanuxem 2005.
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.
1 POD Error
The following errors were encountered while parsing the POD:
Non-ASCII character seen before =encoding in 'Grégory'. Assuming CP1252
To install PDL::LinearAlgebra, copy and paste the appropriate command in to your terminal.
cpanm
cpanm PDL::LinearAlgebra
CPAN shell
perl -MCPAN -e shell install PDL::LinearAlgebra
For more information on module installation, please visit the detailed CPAN module installation guide.