LinearAlgebra¶
Usage
use 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 none
.
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_LAPACKE_INCLUDE_DIR \
-L$PATH_TO_LIBGFORTRAN -lgfortran \
-L$PATH_TO_LAPACK_BINARIES -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=none
flag, so that BLAS
is no longer a dependency:
# Building with BLAS dependency explicitly disabled
chpl --set blasImpl=none example3.chpl
Similarly, users can opt out of of LAPACK
implementations with the --set
lapackImpl=none
flag. Setting both flags to none
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=none --set blasImpl=none 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 1-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
-
var
-
proc
Vector
(length, type eltType = real)¶ Return a vector (1D array) over domain
{1..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) 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)¶ Return a square matrix (2D array) over domain
{1..rows, 1..rows}
-
proc
Matrix
(rows, cols, type eltType = real) Return a matrix (2D array) over domain
{1..rows, 1..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) Return a matrix (2D array) over domain
Dom
-
proc
Matrix
(A: [?Dom] ?Atype, type eltType = Atype) Return a matrix (2D array) with domain and values of
A
.A
can be sparse (CS) or dense.
-
proc
Matrix
(const Arrays: ?t ...?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, type eltType = real)¶ Return a square identity matrix over domain
{1..m, 1..m}
-
proc
eye
(m, n, type eltType = real) Return an identity matrix over domain
{1..m, 1..n}
-
proc
eye
(Dom: domain(2), type eltType = real) Return an identity matrix over domain
Dom
-
proc
transpose
(A: [?Dom] ?eltType)¶ 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
¶ Transpose vector or matrix
-
proc
_array.
plus
(A: [?Adom] ?eltType)¶ Element-wise addition. Same as
A + B
.
-
proc
_array.
minus
(A: [?Adom] ?eltType)¶ Element-wise subtraction. Same as
A - B
.
-
proc
_array.
times
(A: [?Adom])¶ Element-wise multiplication. Same as
A * B
.
-
proc
_array.
elementDiv
(A: [?Adom])¶ Element-wise division. Same as
A / B
.
-
proc
dot
(A: [?Adom] ?eltType, B: [?Bdom] eltType)¶ 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: [])¶ Compute the dot-product
-
proc
inner
(A: [?Adom], B: [?Bdom])¶ Inner product of 2 vectors.
-
proc
outer
(A: [?Adom] ?eltType, B: [?Bdom] eltType)¶ Outer product of 2 vectors.
-
proc
matPow
(A: [], b)¶ Return the matrix
A
to thebth
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)¶ Return true if matrix is diagonal
-
proc
isHermitian
(A: [?D])¶ Return true if matrix is Hermitian
-
proc
isSymmetric
(A: [?D]): bool¶ 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
cholesky
(A: [] ?t, lower = true)¶ 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
isnone
.
-
proc
eigvals
(A: [] ?t)¶ 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
isnone
.
-
proc
eig
(A: [] ?t, param left = false, param right = false)¶ 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
isnone
.- If
-
proc
svd
(A: [?Adom] ?t) throws¶ 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);
LinearAlgebraError
will be thrown if the SVD computation does not converge or an illegal argument, such as a matrix containing aNAN
value, is given.Note
A temporary copy of
A
will be created within this computation.Note
This procedure depends on the
LAPACK
module, and will generate a compiler error iflapackImpl
isnone
.
-
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)