LinearAlgebra¶
Usage
use LinearAlgebra;
or
import LinearAlgebra;
Submodules
A highlevel 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_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 matrixmatrix 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 0based indices.
Row vs Column vectors
Row and column vectors are both represented as 1D arrays and are
indistinguishable in Chapel. In the dot
function, matrixvector
multiplication assumes a column vector, vectormatrix multiplication assumes a
row vector, vectorvector multiplication is always treated as an innerproduct,
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
{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) 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
{0..<rows, 0..<rows}

proc
Matrix
(rows, cols, type eltType = real) 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) 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: 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)¶ Sets the value of a diagonal in a matrix inplace. 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 thek``th row. ``k
is 0indexed.

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)¶ Elementwise addition. Same as
A + B
.

proc
_array.
minus
(A: [?Adom] ?eltType)¶ Elementwise subtraction. Same as
A  B
.

proc
_array.
times
(A: [?Adom])¶ Elementwise multiplication. Same as
A * B
.

proc
_array.
elementDiv
(A: [?Adom])¶ Elementwise 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 dotproduct

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)¶ Returns the inverse of
A
square matrix A.Note
This procedure depends on the
LAPACK
module, and will generate a compiler error iflapackImpl
isoff
.

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 crossproduct of 3element 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])¶ 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).

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 2norm, 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
:Default  depends on array dimensions. See
norm
for details.1norm
2norm
Infinity norm
Frobenius norm

proc
_norm
(x: [?D], param p: normType)¶ 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 system `` L * 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.

proc
solve_triu
(const ref U: [?Udom] ?eltType, const ref b: [?bdom] eltType)¶ Return the solution
x
to the linear system `` U * x = b `` whereU
is an upper triangular matrix.

proc
solve
(A: [?Adom] ?eltType, b: [?bdom] eltType)¶ Return the solution
x
to the linear systemA * x = b
.

proc
leastSquares
(A: [] ?t, b: [] t, cond = 1.0) throws¶ Compute leastsquares solution to
A * x = b
. Compute a vectorx
such that the 2normb  A x
is minimized.cond
is the cutoff threshold such that singular values will be considered 0.0. Ifcond < 0.0
(defaults to1.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 the leastsquares solution with shape ofb
residues
is:the square of the 2norm for each column in
b  a x
ifM > N
andA.rank == n
.a scalar if
b
is 1D
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)¶ 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
(A: [] ?t, lower = true, param overwrite = false) throws¶ Find the eigenvalues of a realsymmetric/complexhermitian matrix
A
.A
must be square.The algorithms uses either the lowertriangular (if
lower
istrue
, or uppertriangular 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
(A: [] ?t, lower = true, param eigvalsOnly = false, param overwrite = false) throws¶ Find the eigenvalues and eigenvectors of a realsymmetric/complexhermitian matrix
A
.A
must be square.The algorithms uses either the lowertriangular (if
lower
istrue
, or uppertriangular 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)¶ 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)¶ 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¶ 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
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)