.. default-domain:: chpl .. module:: LinearAlgebra :synopsis: A high-level interface to linear algebra operations and procedures. LinearAlgebra ============= **Usage** .. code-block:: chapel use LinearAlgebra; or .. code-block:: chapel import LinearAlgebra; **Submodules** .. toctree:: :maxdepth: 1 :glob: LinearAlgebra/* A high-level interface to linear algebra operations and procedures. Compiling with Linear Algebra ----------------------------- Programs using the :mod:`LinearAlgebra` module can be built with no additional dependencies `if` they do not use any procedures that rely on :mod:`BLAS` or :mod:`LAPACK`. Procedure dependencies are specified in procedure documentation below. If a program calls a procedure that depends on :mod:`BLAS` or :mod:`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 :mod:`BLAS` and :mod:`LAPACK` dependent implementations by setting the ``blasImpl`` and ``lapackImpl`` flags to ``off``. **Building programs with no dependencies** .. code-block:: chpl // 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: .. code-block:: bash chpl example1.chpl If this program had used a procedure with a dependency such as :proc:`cholesky` (depends on ``LAPACK``), compilation without ``LAPACK`` headers and libraries available would result in a compilation error. **Building programs with dependencies** .. code-block:: chpl // example2.chpl var A = Matrix([2.0, 1.0], [1.0, 2.0]); var (eigenvalues, eigenvectors) = eig(A, right=true); The program above uses :proc:`eig`, which depends on :mod:`LAPACK`. Compilation without ``LAPACK`` headers and libraries available would result in a compilation error. Following the instructions from the :mod:`LAPACK` module documentation, the above program could be compiled if ``LAPACK`` is available on the system and specified with the following compilation flags: .. code-block:: bash # 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** .. code-block:: chpl // example3.chpl var A = Matrix(3,5); A = 2; var AA = A.dot(A.T); The program above uses :proc:`dot`, which has two available implementations: one that depends on :mod:`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 :mod:`BLAS` module documentation, the above program could be compiled if ``BLAS`` is available on the system and specified with the following compilation flags: .. code-block:: bash # 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: .. code-block:: bash # 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: .. code-block:: bash # Building with all dependencies explicitly disabled chpl --set lapackImpl=off --set blasImpl=off example3.chpl See the documentation of :mod:`BLAS` or :mod:`LAPACK` for more details on these flags. .. _LinearAlgebraInterface: 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 :proc:`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 :proc:`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. .. attribute:: var info: string Stores message to be emitted upon uncaught throw. .. function:: proc Vector(length, type eltType = real) Return a vector (1D array) over domain ``{0.. 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. .. function:: 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. .. method:: proc _array.T where isDenseMatrix(this) Transpose vector or matrix. .. method:: proc _array.plus(A: [?Adom] ?eltType) where isDenseArr(A) && isDenseArr(this) Element-wise addition. Same as ``A + B``. .. method:: proc _array.minus(A: [?Adom] ?eltType) where isDenseArr(A) && isDenseArr(this) Element-wise subtraction. Same as ``A - B``. .. method:: proc _array.times(A: [?Adom]) where isDenseArr(A) && isDenseArr(this) Element-wise multiplication. Same as ``A * B``. .. method:: proc _array.elementDiv(A: [?Adom]) where isDenseArr(A) && isDenseArr(this) Element-wise division. Same as ``A / B``. .. function:: 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 :mod:`BLAS` module for improved performance, if available. Compile with ``--set blasImpl=off`` to opt out of the :mod:`BLAS` implementation. .. method:: proc _array.dot(A: []) where isDenseArr(this) && isDenseArr(A) Compute the dot-product. .. note:: Dense matrix-matrix and matrix-vector multiplication will utilize the :mod:`BLAS` module for improved performance, if available. Compile with ``--set blasImpl=off`` to opt out of the :mod:`BLAS` implementation. .. function:: proc inner(const ref A: [?Adom] ?eltType, const ref B: [?Bdom]) Inner product of 2 vectors. .. function:: proc outer(A: [?Adom] ?eltType, B: [?Bdom] eltType) Outer product of 2 vectors. .. function:: proc inv(ref A: [?Adom] ?eltType, overwrite = false) where usingLAPACK Returns the inverse of the square matrix ``A``. .. note:: This procedure depends on the :mod:`LAPACK` module, and will generate a compiler error if ``lapackImpl`` is ``off``. .. function:: 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 :mod:`BLAS` module for improved performance, if available. Compile with ``--set blasImpl=off`` to opt out of the :mod:`BLAS` implementation. .. function:: proc cross(A: [?Adom] ?eltType, B: [?Bdom] eltType) Return cross-product of 3-element vectors ``A`` and ``B`` with domain of ``A``. .. function:: 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. .. function:: 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: .. code-block:: chapel 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 */ .. function:: 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: .. code-block:: chapel 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 */ .. function:: proc isDiag(A: [?D] ?eltType) where isDenseMatrix(A) Return `true` if matrix is diagonal. .. function:: proc isZero(A: [?D] ?eltType) where isDenseMatrix(A) Return `true` if matrix is the additive identity (zero matrix). .. function:: proc isEye(A: [?D] ?eltType) where isDenseMatrix(A) Return `true` if matrix is the multiplicative identity (identity matrix). .. function:: proc isHermitian(A: [?D]) where isDenseMatrix(A) Return `true` if matrix is Hermitian. .. function:: proc isSymmetric(A: [?D]) where isDenseMatrix(A) Return `true` if matrix is symmetric. .. function:: 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. .. function:: 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. .. function:: proc isSquare(A: [?D]) Return `true` if matrix is square. .. function:: proc trace(A: [?D] ?eltType) Return the trace (sum of diagonal elements) of ``A``. .. function:: 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. .. function:: 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. .. function:: 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. :rtype: x.eltType .. enum:: enum normType { default, norm1, norm2, normInf, normFrob } Indicates the different types of norms supported by :proc:`norm`. .. enumconstant:: enum constant default Default - depends on array dimensions. See :proc:`norm` for details. .. enumconstant:: enum constant norm1 1-norm .. enumconstant:: enum constant norm2 2-norm .. enumconstant:: enum constant normInf Infinity norm .. enumconstant:: enum constant normFrob Frobenius norm .. function:: 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`. :rtype: x.eltType .. function:: 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. .. function:: 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. .. function:: 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. .. function:: 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`` .. note:: This procedure depends on the :mod:`LAPACK` module, and will generate a compiler error if ``lapackImpl`` is ``none``. .. function:: 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 :mod:`LAPACK` module, and will generate a compiler error if ``lapackImpl`` is ``off``. .. function:: 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. .. 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 :mod:`LAPACK` module, and will generate a compiler error if ``lapackImpl`` is ``off``. .. function:: 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. .. note:: This procedure depends on the :mod:`LAPACK` module, and will generate a compiler error if ``lapackImpl`` is ``off``. .. function:: 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 :mod:`LAPACK` module, and will generate a compiler error if ``lapackImpl`` is ``off``. .. function:: 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 :mod:`LAPACK` module, and will generate a compiler error if ``lapackImpl`` is ``off``. .. function:: 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: .. math:: \mathbf{A} = \textbf{U} \cdot \Sigma \cdot \mathbf{V^H} where - :math:`\mathbf{U}` is an `m x m` unitary matrix, - :math:`\Sigma` is a diagonal `m x n` matrix, - :math:`\mathbf{V}` is an `n x n` unitary matrix, and :math:`\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 :math:`\Sigma`, known as the singular values. For example: .. code-block:: chapel 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 :mod:`LAPACK` module, and will generate a compiler error if ``lapackImpl`` is ``off``. .. function:: 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. .. function:: 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)``. .. function:: 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``. :arg A: Expects a square matrix. :type A: `A` :arg useExactOneNorm: boolean value specifying if the onenorm has to be exact. Defaults to `true`. :type useExactOneNorm: bool :throws LinearAlgebraError: If the input matrix is not a square matrix. :returns: Matrix exponential of the given matrix. :rtype: `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 .. function:: proc sincos(A: []) throws This method returns both sine and cosine of the matrix ``A``. :arg A: Expects a square matrix. :type A: `A` :returns: Matrix a tuple of sin and cosine of the given matrix. :rtype: (`A`, `A`) :throws LinearAlgebraError: If the input matrix is not a square matrix. .. function:: proc sinm(A: []) throws This method returns the sine of the matrix ``A``. :arg A: Expects a square matrix. :type A: `A` :returns: Matrix returns the sine of the given matrix. :rtype: `A` :throws LinearAlgebraError: If the input matrix is not a square matrix. .. function:: proc cosm(A: []) throws This method returns the cosine of the matrix ``A``. :arg A: Expects a square matrix. :type A: `A` :returns: Matrix returns the cosine of the given matrix. :rtype: `A` :throws LinearAlgebraError: If the input matrix is not a square matrix.