LinearAlgebra¶
Usage
use LinearAlgebra;
or
import LinearAlgebra;
Submodules
A high-level interface to linear algebra operations and procedures.
Compiling with Linear Algebra¶
Programs using the LinearAlgebra
module can be built with no additional
dependencies if they do not use any procedures that rely on BLAS
or
LAPACK
. Procedure dependencies are specified in procedure documentation
below.
If a program calls a procedure that depends on BLAS
or LAPACK
, the
headers and library will need to be available during compilation, typically
through compiler flags and/or environment variables.
Some procedures have implementations both with and without dependencies. By
default, the implementation with dependencies will be selected. Users can
explicitly opt out of using the BLAS
and LAPACK
dependent
implementations by setting the blasImpl
and lapackImpl
flags to off
.
Building programs with no dependencies
// example1.chpl
var A = Matrix([0.0, 1.0, 1.0],
[1.0, 0.0, 1.0],
[1.0, 1.0, 0.0]);
var I = eye(3,3);
var B = A + I;
The program above has no dependencies and can therefore be compiled without
the BLAS
or LAPACK
headers and libraries available:
chpl example1.chpl
If this program had used a procedure with a dependency such as
cholesky
(depends on LAPACK
), compilation without LAPACK
headers and libraries available would result in a compilation error.
Building programs with dependencies
// example2.chpl
var A = Matrix([2.0, 1.0],
[1.0, 2.0]);
var (eigenvalues, eigenvectors) = eig(A, right=true);
The program above uses eig
, which depends on LAPACK
.
Compilation without LAPACK
headers and libraries available would result in
a compilation error.
Following the instructions from the LAPACK
module
documentation, the above program could be compiled if LAPACK
is available
on the system and specified with the following compilation flags:
# Building with LAPACK dependency
chpl -I/path/to/lapack/header/file \
-L/path/to/gfortran/library/file -lgfortran \
-L/path/to/lapack/library/files -llapacke -llapack -lrefblas \
example2.chpl
Building programs with optional dependencies
// example3.chpl
var A = Matrix(3,5);
A = 2;
var AA = A.dot(A.T);
The program above uses dot
, which has two available implementations:
one that depends on BLAS
and one that is written in Chapel. The
program will default to using the more performant BLAS
implementation
of matrix-matrix multiplication.
Following the instructions from the
BLAS
module documentation, the above program could be compiled if
BLAS
is available on the system and specified with the following
compilation flags:
# Building with BLAS dependency
chpl -I$PATH_TO_CBLAS_DIR \
-L$PATH_TO_BLAS_LIBS -lblas \
example3.chpl
Note
Users can set environment variables like LDFLAGS
for -L
arguments and
CFLAGS
for -I
arguments, to avoid throwing these flags every time.
Additionally, the required linker flags (-l
) may vary depending on the
BLAS
and LAPACK
implementations being used.
To opt out of using the BLAS
implementation, users can add the --set
blasImpl=off
flag, so that BLAS
is no longer a dependency:
# Building with BLAS dependency explicitly disabled
chpl --set blasImpl=off example3.chpl
Similarly, users can opt out of of LAPACK
implementations with the --set
lapackImpl=off
flag. Setting both flags to off
will always choose the
Chapel implementation when available, and will emit a compiler error
when no native implementation is available:
# Building with all dependencies explicitly disabled
chpl --set lapackImpl=off --set blasImpl=off example3.chpl
See the documentation of BLAS
or LAPACK
for
more details on these flags.
Linear Algebra Interface¶
Matrix and Vector representation
Matrices and vectors are represented as Chapel arrays, which means the convenience constructors provided in this module are not necessary to use this module’s functions. This also means that matrices (2D arrays) and vectors (1D arrays) will work with any other Chapel library that works with arrays.
Note
This documentation uses the terms matrix to refer to 2D arrays, and vector to refer to 1D arrays.
Domain offsets
All functions that return arrays will inherit their domains from the input array if possible. Otherwise they will return arrays with 0-based indices.
Row vs Column vectors
Row and column vectors are both represented as 1D arrays and are
indistinguishable in Chapel. In the dot
function, matrix-vector
multiplication assumes a column vector, vector-matrix multiplication assumes a
row vector, vector-vector multiplication is always treated as an inner-product,
as the function name implies.
An outer product can be computed with the outer
function.
Domain maps
All of the functions in this module only support
DefaultRectangular
arrays (the default domain map), unless explicitly
specified in the function’s documentation. Other domain maps
are supported through submodules, such LinearAlgebra.Sparse
for the
CS
layout.
- class LinearAlgebraError : Error¶
Base
Error
type forLinearAlgebra
errors.- var info : string¶
Stores message to be emitted upon uncaught throw.
- proc Vector(length, type eltType = real)¶
Return a vector (1D array) over domain
{0..<length}
.
- proc Vector(space: range, type eltType = real)
Return a vector (1D array) over domain
{space}
.
- proc Vector(Dom: domain(1), type eltType = real)
Return a vector (1D array) over domain
Dom
.
- proc Vector(A: [?Dom] ?Atype, type eltType = Atype)
Return a vector (1D array) with domain and values of
A
.
- proc Vector(x: ?t, Scalars ...?n, type eltType) where isNumericType(t)
Return a vector (1D array), given 2 or more numeric values.
If type is omitted, it will be inferred from the first argument.
- proc Matrix(rows, type eltType = real) where isIntegral(rows)¶
Return a square matrix (2D array) over domain
{0..<rows, 0..<rows}
.
- proc Matrix(rows, cols, type eltType = real) where isIntegral(rows) && isIntegral(cols)
Return a matrix (2D array) over domain
{0..<rows, 0..<cols}
.
- proc Matrix(space: range, type eltType = real)
Return a square matrix (2D array) over domain
{space, space}
.
- proc Matrix(rowSpace: range, colSpace: range, type eltType = real)
Return a matrix (2D array) over domain
{rowSpace, colSpace}
.
- proc Matrix(Dom: domain, type eltType = real) where Dom.rank == 2
Return a matrix (2D array) over domain
Dom
.
- proc Matrix(A: [?Dom] ?Atype, type eltType = Atype) where isDenseMatrix(A)
Return a matrix (2D array) with domain and values of
A
.A
can be sparse (CS) or dense.
- proc Matrix(const Arrays ...?n, type eltType)
Return a matrix (2D array), given 2 or more vectors, such that the vectors form the rows of the matrix. In other words, the vectors are concatenated such that the
ith
vector corresponds to the matrix slice:A[i, ..]
.If type is omitted, it will be inferred from the first array.
For example:
var A = Matrix([1, 2, 3], [4, 5, 6], [7, 8, 9]); /* Produces the 3x3 matrix of integers: 1 2 3 4 5 6 7 8 9 */
- proc eye(m: integral, type eltType = real)¶
Return a square identity matrix over domain
{0..<m, 0..<m}
.
- proc eye(m: integral, n: integral, type eltType = real)
Return an identity matrix over domain
{0..<m, 0..<n}
.
- proc eye(Dom: domain(2), type eltType = real)
Return an identity matrix over domain
Dom
.
- proc setDiag(ref X: [?D] ?eltType, in k: int = 0, val: eltType = 0) where isDenseMatrix(X)¶
Sets the value of a diagonal in a matrix in-place. If the matrix is sparse, indices on the diagonal will be added to its domain.
k > 0
, represents an upper diagonal starting from thek
’th column,k == 0
represents the main diagonal,k < 0
represents a lower diagonal starting from the-k
’th row.k
is 0-indexed.
- proc transpose(A: [?Dom] ?eltType) where isDenseMatrix(A)¶
Transpose vector, matrix, or domain.
Note
Since row vectors and columns vectors are indistinguishable, passing a vector to this function will return that vector unchanged.
- proc _array.T where isDenseMatrix(this)¶
Transpose vector or matrix.
- proc _array.plus(A: [?Adom] ?eltType) where isDenseArr(A) && isDenseArr(this)¶
Element-wise addition. Same as
A + B
.
- proc _array.minus(A: [?Adom] ?eltType) where isDenseArr(A) && isDenseArr(this)¶
Element-wise subtraction. Same as
A - B
.
- proc _array.times(A: [?Adom]) where isDenseArr(A) && isDenseArr(this)¶
Element-wise multiplication. Same as
A * B
.
- proc _array.elementDiv(A: [?Adom]) where isDenseArr(A) && isDenseArr(this)¶
Element-wise division. Same as
A / B
.
- proc dot(A: [?Adom] ?eltType, B: [?Bdom] eltType) where isDenseArr(A) && isDenseArr(B)¶
Generic matrix multiplication,
A
andB
can be a matrix, vector, or scalar.Note
When
A
is a vector andB
is a matrix, this function implicitly computesdot(transpose(A), B)
, which may not be as efficient as passingA
andB
in the reverse order.
- proc _array.dot(A: []) where isDenseArr(this) && isDenseArr(A)¶
Compute the dot-product.
- proc inner(const ref A: [?Adom] ?eltType, const ref B: [?Bdom])¶
Inner product of 2 vectors.
- proc outer(A: [?Adom] ?eltType, B: [?Bdom] eltType)¶
Outer product of 2 vectors.
- proc inv(ref A: [?Adom] ?eltType, overwrite = false) where usingLAPACK¶
Returns the inverse of the square matrix
A
.Note
This procedure depends on the
LAPACK
module, and will generate a compiler error iflapackImpl
isoff
.
- proc matPow(A: [], b) where isNumeric(b)¶
Return the matrix
A
to theb
’th power, whereb
is a positive integral type.
- proc cross(A: [?Adom] ?eltType, B: [?Bdom] eltType)¶
Return cross-product of 3-element vectors
A
andB
with domain ofA
.
- proc diag(A: [?Adom] ?eltType, k = 0)¶
Return a Vector containing the diagonal elements of
A
if the argumentA
is of rank 2. Return a diagonal Matrix whose diagonal contains elements ofA
if argumentA
is of rank 1.
- proc tril(A: [?D] ?eltType, k = 0)¶
Return lower triangular part of matrix, above the diagonal +
k
, wherek = 0
includes the diagonal, andk = -1
does not include the diagonal. For example:var A = Matrix(4, 4, eltType=int); A = 1; tril(A); /* Returns: 1 0 0 0 1 1 0 0 1 1 1 0 1 1 1 1 */ tril(A, 1); /* Returns: 1 1 0 0 1 1 1 0 1 1 1 1 1 1 1 1 */ tril(A, -1); /* Returns: 0 0 0 0 1 0 0 0 1 1 0 0 1 1 1 0 */
- proc triu(A: [?D] ?eltType, k = 0)¶
Return upper triangular part of matrix, above the diagonal +
k
, wherek = 0
includes the diagonal, andk = 1
does not include the diagonal. For example:var A = Matrix(4, 4, eltType=int); A = 1; triu(A); /* Returns: 1 1 1 1 0 1 1 1 0 0 1 1 0 0 0 1 */ triu(A, 1); /* Returns: 0 1 1 1 0 0 1 1 0 0 0 1 0 0 0 0 */ triu(A, -1); /* Returns: 1 1 1 1 1 1 1 1 0 1 1 1 0 0 1 1 */
- proc isDiag(A: [?D] ?eltType) where isDenseMatrix(A)¶
Return true if matrix is diagonal.
- proc isZero(A: [?D] ?eltType) where isDenseMatrix(A)¶
Return true if matrix is the additive identity (zero matrix).
- proc isEye(A: [?D] ?eltType) where isDenseMatrix(A)¶
Return true if matrix is the multiplicative identity (identity matrix).
- proc isHermitian(A: [?D]) where isDenseMatrix(A)¶
Return true if matrix is Hermitian.
- proc isSymmetric(A: [?D]) where isDenseMatrix(A)¶
Return true if matrix is symmetric.
- proc isTril(A: [?D] ?eltType, k = 0) : bool¶
Return true if matrix is lower triangular below the diagonal +
k
, wherek = 0
does not include the diagonal, andk = 1
includes the diagonal.
- proc isTriu(A: [?D] ?eltType, k = 0) : bool¶
Return true if matrix is upper triangular above the diagonal +
k
, wherek = 0
does not include the diagonal, andk = -1
includes the diagonal.
- proc isSquare(A: [?D])¶
Return true if matrix is square.
- proc trace(A: [?D] ?eltType)¶
Return the trace (sum of diagonal elements) of
A
.
- proc lu(A: [?Adom] ?eltType)¶
Compute an LU factorization of square matrix A using partial pivoting, such that A = P * L * U where P is a permutation matrix. Return a tuple of size 2 (LU, ipiv).
L and U are stored in the same matrix LU where the unit diagonal elements of L are not stored.
ipiv contains the pivot indices such that row i of A was interchanged with row ipiv(i).
Note
Arrays with any offset are supported, and LU and ipiv inherit the indexing.
- proc det(A: [?Adom] ?eltType)¶
Return the determinant of a square matrix.
Note
This procedure performs LU factorization to compute the determinant. In certain cases, e.g. having a lower/upper triangular matrix, it is more desirable to compute the determinant manually.
- proc norm(x: [], param p = normType.default)¶
Compute the default norm on x.
For a 1D array this is the 2-norm, for a 2D array, this is the Frobenius norm.
- Return type:
x.eltType
- enum normType { default, norm1, norm2, normInf, normFrob }¶
Indicates the different types of norms supported by
norm
.- enum constant norm1¶
1-norm
- enum constant norm2¶
2-norm
- enum constant normInf¶
Infinity norm
- enum constant normFrob¶
Frobenius norm
- proc _norm(x: [?D], param p: normType) where x.rank == 2¶
Compute the norm indicated by p on the 2D array x.
p cannot be normType.norm2.
- Return type:
x.eltType
- proc solve_tril(const ref L: [?Ldom] ?eltType, const ref b: [?bdom] eltType, unit_diag = true)¶
Return the solution
x
to the linear systemL * x = b
whereL
is a lower triangular matrix. Setting unit_diag to true will assume the diagonal elements as 1 and will not be referenced within this procedure.Note
Arrays with any offset are supported, and
x
inherits the indexing.
- proc solve_triu(const ref U: [?Udom] ?eltType, const ref b: [?bdom] eltType)¶
Return the solution
x
to the linear systemU * x = b
whereU
is an upper triangular matrix.Note
Arrays with any offset are supported, and
x
inherits the indexing.
- proc solve(const ref A: [?Adom] ?eltType, const ref b: [?bdom] eltType)¶
Return the solution
x
to the linear systemA * x = b
.Note
Arrays with any offset are supported, and
x
inherits the indexing.
- proc leastSquares(A: [] ?t, b: [] t, cond = -1.0) throws where usingLAPACK && isLAPACKType(t)¶
Compute least-squares solution to
A * x = b
. Compute a vectorx
such that the 2-norm|b - A x|
is minimized.cond
is the cut-off threshold such that singular values will be considered 0.0. Ifcond < 0.0
(defaults to-1.0
), the threshold will be set tomax((...A.shape)) * epsilon
, whereepsilon
is the machine precision forA.eltType
.Returns a tuple of
(x, residues, rank, s)
, where:x
is the least-squares solution with shape ofb
residues
is:the square of the 2-norm for each column in
b - a x
ifM > N
andA.rank == n
.a scalar if
b
is 1-D
rank
is the effective rank ofA
s
is the singular values ofa
Note
This procedure depends on the
LAPACK
module, and will generate a compiler error iflapackImpl
isnone
.
- proc cholesky(A: [] ?t, lower = true) where A.rank == 2 && isLAPACKType(t) && usingLAPACK¶
Perform a Cholesky factorization on matrix
A
.A
must be square. Argumentlower
indicates whether to return the lower or upper triangular factor. MatrixA
is not modified. Returns an array with the same shape as argumentA
with the lower or upper triangular Cholesky factorization ofA
.Note
This procedure depends on the
LAPACK
module, and will generate a compiler error iflapackImpl
isoff
.
- proc eigvalsh(ref A: [] ?t, lower = true, param overwrite = false) throws where A.domain.rank == 2 && usingLAPACK¶
Find the eigenvalues of a real-symmetric/complex-hermitian matrix
A
.A
must be square.The algorithms uses either the lower-triangular (if
lower
istrue
, or upper-triangular part of the matrix only. Ifoverwrite
is true, on exiting, this part of the matrix, including the diagonal is overwritten.Note
This procedure currently just returns all eigenvalues. To selectively return certain eigenvalues, the user should call the LAPACK routine directly.
Note
This procedure depends on the
LAPACK
module, and will generate a compiler error iflapackImpl
isoff
.
- proc eigh(ref A: [] ?t, lower = true, param eigvalsOnly = false, param overwrite = false) throws where A.domain.rank == 2 && usingLAPACK¶
Find the eigenvalues and eigenvectors of a real-symmetric/complex-hermitian matrix
A
.A
must be square.The algorithms uses either the lower-triangular (if
lower
istrue
, or upper-triangular part of the matrix only.If
overwrite
is true, the matrix is overwritten with the eigenvectors and only the eigenvalues are returned, otherwise the original matrix is preserved.The eigenvectors are stored in the columns of the returned matrix i.e.
A[..,i]
is thei
’th eigenvector.Note
This procedure currently returns all eigenvalues and eigenvectors. To selectively return certain eigenvalues/eigenvectors, the user should call the LAPACK routine directly.
Note
This procedure depends on the
LAPACK
module, and will generate a compiler error iflapackImpl
isoff
.
- proc eigvals(A: [] ?t) where A.domain.rank == 2 && usingLAPACK¶
Find the eigenvalues of matrix
A
.A
must be square.Note
This procedure depends on the
LAPACK
module, and will generate a compiler error iflapackImpl
isoff
.
- proc eig(A: [] ?t, param left = false, param right = false) where A.domain.rank == 2 && usingLAPACK¶
Find the eigenvalues and eigenvectors of matrix
A
.A
must be square.If
left
istrue
then the “left” eigenvectors are computed. The return value is a tuple with two elements:(eigenvalues, leftEigenvectors)
If
right
istrue
then the “right” eigenvectors are computed. The return value is a tuple with two elements:(eigenvalues, rightEigenvectors)
If
left
andright
are bothtrue
then both eigenvectors are computed. The return value is a tuple with three elements:(eigenvalues, leftEigenvectors, rightEigenvectors)
If
left
andright
are bothfalse
only the eigenvalues are computed, and returned as a single array.
Note
This procedure depends on the
LAPACK
module, and will generate a compiler error iflapackImpl
isoff
.
- proc svd(A: [?Adom] ?t) throws where isLAPACKType(t) && usingLAPACK && Adom.rank == 2 && Adom.strides == strideKind.one¶
Singular Value Decomposition.
Factorizes the m x n matrix
A
such that:\[\mathbf{A} = \textbf{U} \cdot \Sigma \cdot \mathbf{V^H}\]where
\(\mathbf{U}\) is an m x m unitary matrix,
\(\Sigma\) is a diagonal m x n matrix,
\(\mathbf{V}\) is an n x n unitary matrix, and \(\mathbf{V^H}\) is the Hermitian transpose.
This procedure returns a tuple of
(U, s, Vh)
, wheres
is a vector containing the diagonal elements of \(\Sigma\), known as the singular values.For example:
var A = Matrix([3, 2, 2], [2, 3, -2], eltType=real); var (U, s, Vh) = svd(A);
- Throws:
LinearAlgebraError – if the SVD computation does not converge or an illegal argument, such as a matrix containing a
NAN
value, is given.
Note
A temporary copy of
A
will be created within this computation.Note
Arrays with strided domains are not supported.
Note
Arrays whose domains have nonzero offsets are supported.
U
inherits the row indexing ofA
, whileVh
inherits the column indexing ofA
. The columns ofU
, rows ofVh
, ands
all share the same 0-based indexing.Note
This procedure depends on the
LAPACK
module, and will generate a compiler error iflapackImpl
isoff
.
- proc jacobi(A: [?Adom] ?eltType, ref X: [?Xdom] eltType, b: [Xdom] eltType, tol = 0.0001, maxiter = 1000)¶
Compute the approximate solution to
A * x = b
using the Jacobi method. iteration will stop whenmaxiter
is reached or error is smaller thantol
, whichever comes first. Return the number of iterations performed.Note
X
is passed as a reference, meaning the initial solution guess can be stored inX
before calling the procedure, and the approximate solution will be stored in the same array.Dense and CSR arrays are supported.
- proc kron(A: [?ADom] ?eltType, B: [?BDom] eltType)¶
Return the Kronecker Product of matrix
A
and matrixB
. If the size of A isx * y
and of B isa * b
then size of the resulting matrix will be(x * a) * (y * b)
.
- proc expm(A: [], param useExactOneNorm = true) throws¶
Matrix exponential using Pade approximation. This method returns a square matrix which is the matrix exponential of
A
.- Arguments:
A : A – Expects a square matrix.
useExactOneNorm :
bool
– boolean value specifying if the onenorm has to be exact. Defaults to true.
- Throws:
LinearAlgebraError – If the input matrix is not a square matrix.
- Returns:
Matrix exponential of the given matrix.
- Return type:
A
Note
[1] Awad H. Al-Mohy and Nicholas J. Higham (2009) “A New Scaling and Squaring Algorithm for the Matrix Exponential.” SIAM Journal on Matrix Analysis and Applications. 31 (3). pp. 970-989. ISSN 1095-7162
- proc sincos(A: []) throws¶
This method returns both sine and cosine of the matrix
A
.- Arguments:
A : A – Expects a square matrix.
- Returns:
Matrix a tuple of sin and cosine of the given matrix.
- Return type:
(A, A)
- Throws:
LinearAlgebraError – If the input matrix is not a square matrix.
- proc sinm(A: []) throws¶
This method returns the sine of the matrix
A
.- Arguments:
A : A – Expects a square matrix.
- Returns:
Matrix returns the sine of the given matrix.
- Return type:
A
- Throws:
LinearAlgebraError – If the input matrix is not a square matrix.
- proc cosm(A: []) throws¶
This method returns the cosine of the matrix
A
.- Arguments:
A : A – Expects a square matrix.
- Returns:
Matrix returns the cosine of the given matrix.
- Return type:
A
- Throws:
LinearAlgebraError – If the input matrix is not a square matrix.