LinearAlgebra

Usage

use LinearAlgebra;

or

import LinearAlgebra;

Submodules

Sparse

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 for LinearAlgebra 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 the k’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 and B can be a matrix, vector, or scalar.

Note

When A is a vector and B is a matrix, this function implicitly computes dot(transpose(A), B), which may not be as efficient as passing A and B in the reverse order.

Note

Dense matrix-matrix and matrix-vector multiplication will utilize the BLAS module for improved performance, if available. Compile with --set blasImpl=off to opt out of the BLAS implementation.

proc _array.dot(A: []) where isDenseArr(this) && isDenseArr(A): Compute the dot-product.

Note

Dense matrix-matrix and matrix-vector multiplication will utilize the BLAS module for improved performance, if available. Compile with --set blasImpl=off to opt out of the BLAS implementation.

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 if lapackImpl is off.

proc matPow(A: [], b) where isNumeric(b): Return the matrix A to the b’th power, where b is a positive integral type.

Note

matPow will utilize the BLAS module for improved performance, if available. Compile with --set blasImpl=off to opt out of the BLAS implementation.

proc cross(A: [?Adom] ?eltType, B: [?Bdom] eltType): Return cross-product of 3-element vectors A and B with domain of A.

proc diag(A: [?Adom] ?eltType, k = 0): Return a Vector containing the diagonal elements of A if the argument A is of rank 2. Return a diagonal Matrix whose diagonal contains elements of A if argument A is of rank 1.

proc tril(A: [?D] ?eltType, k = 0)

Return lower triangular part of matrix, above the diagonal + k, where k = 0 includes the diagonal, and k = -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, where k = 0 includes the diagonal, and k = 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, where k = 0 does not include the diagonal, and k = 1 includes the diagonal.

proc isTriu(A: [?D] ?eltType, k = 0) : bool: Return true if matrix is upper triangular above the diagonal + k, where k = 0 does not include the diagonal, and k = -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 default: Default - depends on array dimensions. See norm for details.

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 system L * x = b where L 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 system U * x = b where U 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 system A * 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 vector x 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. If cond < 0.0 (defaults to -1.0), the threshold will be set to max((...A.shape)) * epsilon, where epsilon is the machine precision for A.eltType.

Returns a tuple of (x, residues, rank, s), where:

x is the least-squares solution with shape of b

residues is:

the square of the 2-norm for each column in b - a x if M > N and A.rank == n.

a scalar if b is 1-D

rank is the effective rank of A

s is the singular values of a

Throws:: LinearAlgebraError – If the number of rows in A does not match the size of b, or if either is empty.

Note

This procedure depends on the LAPACK module, and will generate a compiler error if lapackImpl is none.

proc cholesky(A: [] ?t, lower = true) where A.rank == 2 && isLAPACKType(t) && usingLAPACK: Perform a Cholesky factorization on matrix A. A must be square. Argument lower indicates whether to return the lower or upper triangular factor. Matrix A is not modified. Returns an array with the same shape as argument A with the lower or upper triangular Cholesky factorization of A.

Note

This procedure depends on the LAPACK module, and will generate a compiler error if lapackImpl is off.

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 is true, or upper-triangular part of the matrix only. If overwrite is true, on exiting, this part of the matrix, including the diagonal is overwritten.

Throws:: LinearAlgebraError – If A is not square, or if the eigenvalue computation does not converge.

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 if lapackImpl is off.

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 is true, 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 the i’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.

Throws:: LinearAlgebraError – If A is not square, or if the eigenvalue computation does not converge.

Note

This procedure depends on the LAPACK module, and will generate a compiler error if lapackImpl is off.

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 if lapackImpl is off.

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 is true then the “left” eigenvectors are computed. The return value is a tuple with two elements: (eigenvalues, leftEigenvectors)
If right is true then the “right” eigenvectors are computed. The return value is a tuple with two elements: (eigenvalues, rightEigenvectors)
If left and right are both true then both eigenvectors are computed. The return value is a tuple with three elements: (eigenvalues, leftEigenvectors, rightEigenvectors)
If left and right are both false 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 if lapackImpl is off.

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), where s 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 of A, while Vh inherits the column indexing of A. The columns of U, rows of Vh, and s all share the same 0-based indexing.

Note

This procedure depends on the LAPACK module, and will generate a compiler error if lapackImpl is off.

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 when maxiter is reached or error is smaller than tol, whichever comes first. Return the number of iterations performed.

Note

X is passed as a reference, meaning the initial solution guess can be stored in X 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 matrix B. If the size of A is x * y and of B is a * 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.