This package provides precision-agnostic, high-level linear algebra APIs for real and complex arguments in Modern Fortran. The APIs are similar to NumPy/SciPy operations, and leverage a Modern Fortran implementation of the Reference-LAPACK library.

A full and standardized implementation of the present library has been integrated into the Fortran Standard Library, and as such, most users should seek to access the functionality from stdlib. The present library is kept in place for those who seek a compact implementation of it.

Browse API

All procedures work with all types (real, complex) and kinds (32, 64, 128-bit floats).

chol - Cholesky factorization of a matrix (function).

Syntax

c = chol(a [, lower] [, other_zeroed])

Description

This function computes the Cholesky factorization of a real symmetric or complex Hermitian matrix $A$ :

$A = L L^T = U^T U$

where $L$ is a lower triangular matrix and $U$ is an upper triangular matrix. The function returns the factorized matrix as a new allocation, without modifying the input matrix.

Arguments

a: A real or complex matrix of size $[n,n]$ , representing the symmetric/Hermitian input matrix.
lower (optional): A logical flag indicating whether the lower ( $L$ ) or upper ( $U$ ) triangular factor should be computed. Defaults to lower = .true..
other_zeroed (optional): A logical flag determining whether the unused half of the returned matrix should be explicitly zeroed. Defaults to other_zeroed = .true..

Return value

c: A real or complex matrix of size $[n,n]$ , containing the Cholesky factors. The returned matrix is triangular (upper or lower, as selected).

Errors

Raises LINALG_VALUE_ERROR if the input matrix ihas invalid size.
Raises LINALG_ERROR if numerical instability prevents factorization.
If error handling is not provided, exceptions will trigger an error stop.

Notes

The function is based on LAPACK's POTRF routines.
This function allocates a new matrix to store the factorization. For an in-place version, use cholesky.

cholesky - Cholesky factorization of a matrix (subroutine).

Syntax

call cholesky(a [, c] [, lower] [, other_zeroed])

Description

This subroutine computes the Cholesky factorization of a real symmetric or complex Hermitian matrix $A$ :

$A = L L^T = U^T U$

where $L$ is a lower triangular matrix and $U$ is an upper triangular matrix. The factorization is performed in-place, modifying the input matrix a, or on a pre-allocated matrix c with the same type and kind as a.

Arguments

a: A real or complex matrix of size $[n,n]$ , representing the symmetric/Hermitian input matrix. if c is not provided, on return it contains the Cholesky factorization.
c (optional): A matrix of size $[n,n]$ , of the same type and kind as a, containing the Cholesky factorization. If provided, a is unchanged.
lower (optional): A logical flag indicating whether the lower ( $L$ ) or upper ( $U$ ) triangular factor should be computed. Defaults to lower = .true..
other_zeroed (optional): A logical flag determining whether the unused half of the matrix should be explicitly zeroed. Defaults to other_zeroed = .true..

Errors

Raises LINALG_VALUE_ERROR if the input matrix is not positive definite.
Raises LINALG_ERROR if numerical instability prevents factorization.
If error handling is not provided, exceptions will trigger an error stop.

Notes

The returned Cholesky factorization matrix is triangular (upper or lower, as selected).
The subroutine is based on LAPACK's POTRF routines.
This subroutine modifies the input matrix in-place. For a version that returns a newly allocated matrix, use chol.

eig - Eigendecomposition of a square matrix.

Syntax

call eig(a [, b] [, lambda] [, right] [, left] [, overwrite_a] [, overwrite_b] [, err])

Description

This interface provides methods for computing the eigenvalues and eigenvectors of a real or complex matrix. It supports both standard and generalized eigenvalue problems, allowing for the decomposition of a matrix A alone or a pair of matrices (A, B) in the generalized case.

Given a square matrix $A$ , this routine computes its eigenvalues $\lambda$ and, optionally, its right or left eigenvectors:

$A v = \lambda v$

where $v$ represents an eigenvector corresponding to eigenvalue $\lambda$ .

In the generalized eigenvalue problem case, the routine solves:

$A v = \lambda B v$

The computation supports both real and complex matrices. If requested, eigenvectors are returned as additional output arguments. The function provides options to allow in-place modification of A and B for performance optimization.

Note: The solution is based on LAPACK's GEEV and GGEV routines.

Arguments

a: A real or complex matrix of size $[n,n]$ , representing the input matrix to be decomposed.
b (optional): A real or complex matrix of size $[n,n]$ , same type and kind as a, representing the second matrix in the generalized eigenvalue problem.
lambda: A complex or real array of length $n$ , containing the computed eigenvalues.
right (optional): A complex matrix of size $[n,n]$ containing the right eigenvectors as columns.
left (optional): A complex matrix of size $[n,n]$ containing the left eigenvectors as columns.
overwrite_a (optional): A logical flag indicating whether A can be overwritten for performance optimization.
overwrite_b (optional): A logical flag indicating whether B can be overwritten (only in the generalized case).
err (optional): A type(la_state) variable to handle errors. If not provided, execution will stop on errors.

Errors

Raises LINALG_VALUE_ERROR if the matrices have invalid/incompatible sizes.
Raises LINALG_ERROR if the eigendecomposition fails.
If err is not provided, exceptions will trigger an error stop.

Notes

The computed eigenvectors are normalized.
If computing real eigenvalues, an error is returned if eigenvalues have nonzero imaginary parts.
This routine is based on LAPACK's GEEV and GGEV routines.
Overwriting A or B can improve performance but destroys the original matrix data.

eigh - Eigendecomposition of a real symmetric or complex Hermitian matrix.

Syntax

call eigh(a, lambda [, vectors] [, upper_a] [, overwrite_a] [, err])

Description

This interface provides methods for computing the eigenvalues and optionally the eigenvectors of a real symmetric or complex Hermitian matrix.

Given a real symmetric or complex Hermitian matrix $A$ , this routine computes its eigenvalues $\lambda$ and, optionally, its right eigenvectors:

$A v = \lambda v$

where $v$ represents an eigenvector corresponding to eigenvalue $\lambda$ .

The computation supports both real and complex matrices, and the eigenvectors, if requested, are returned as orthonormal vectors.

Note: The solution is based on LAPACK's SYEVD and HEEVD routines.

Arguments

a: A real or complex matrix of size $[n,n]$ , representing the input matrix to be decomposed. The matrix is overwritten with the eigenvalues on output.
lambda: A real array of length $n$ , with the same kind as a, containing the computed eigenvalues.
vectors (optional): A matrix of size $[n,n]$ , with the same type and kind as a, containing the right eigenvectors stored as columns.
overwrite_a (optional): A logical flag indicating whether the matrix A can be overwritten for performance optimization.
upper_a (optional): A logical flag indicating whether the upper half of matrix A should be used for computation (default is the lower half).
err (optional): A type(la_state) variable to handle errors. If not provided, execution will stop on errors.

Errors

Raises LINALG_VALUE_ERROR if the matrices have invalid/incompatible sizes.
Raises LINALG_ERROR if the eigendecomposition fails.
If err is not provided, exceptions will trigger an error stop.

Notes

The computed eigenvectors are orthonormal.
Eigenvalues are real for symmetric or Hermitian matrices.
This routine is based on LAPACK's SYEVD and HEEVD routines.
Overwriting the matrix A can improve performance but destroys the original matrix data.

eigvals - Eigenvalues of a square matrix (interface).

Syntax

lambda = eigvals(a [, b] [, err])

Description

This interface provides methods for computing the eigenvalues of a real or complex square matrix. It supports both standard and generalized eigenvalue problems.

In the standard eigenvalue problem, the function computes the eigenvalues $\lambda$ of the matrix $A$ such that:

$A v = \lambda v$

where $v$ is the eigenvector corresponding to eigenvalue $\lambda$ .

In the generalized eigenvalue problem, the function solves:

$A v = \lambda B v$

where $A$ and $B$ are the input matrices and $\lambda$ is the eigenvalue.

The function returns an array of eigenvalues computed for the input matrix $A$ , and optionally the matrix $B$ for the generalized case.

Note: The solution is based on LAPACK's GEEV and GGEV routines.

Arguments

a: A real or complex matrix of size $[n,n]$ , representing the input matrix to be decomposed.
b (optional, generalized case): A matrix of size $[n,n]$ and same type and kind as a, representing the second matrix in the generalized eigenvalue problem.
err (optional): A type(la_state) variable to handle errors. If not provided, execution will stop on errors.

Return value

lambda: A complex array of eigenvalues, computed from the input matrix $A$ (and $B$ if in the generalized case).

Errors

Raises LINALG_VALUE_ERROR if the matrices have invalid/incompatible sizes.
Raises LINALG_ERROR if the eigendecomposition fails.
If err is not provided, exceptions will trigger an error stop.

Notes

The eigenvalues are returned in a complex array, even for real matrices.
For the generalized eigenvalue problem, matrix $B$ must be provided and will be modified in-place.
This routine is based on LAPACK's GEEV and GGEV routines.

eigvalsh - Eigenvalues of a real symmetric or complex Hermitian matrix.

Syntax

lambda = eigvalsh(a [, upper_a] [, err])

Description

This interface provides methods for computing the eigenvalues of a real symmetric or complex Hermitian matrix. The function computes the eigenvalues of the matrix $A$ , and returns them in an array. The user can specify whether to use the upper or lower half of the matrix for computation.

The function solves the eigenvalue problem:

$A v = \lambda v$

where $v$ is the eigenvector corresponding to eigenvalue $\lambda$ .

The computation supports both real and complex matrices. Regardless, due to symmetry the eigenvalues are returned as an array of real values.

The user can specify whether to use the upper or lower half of the matrix $A$ for the computation (default: lower half).

Note: The solution is based on LAPACK's SYEV and HEEV routines.

Arguments

a: A real or complex matrix of size $[n,n]$ , representing the real symmetric or complex Hermitian matrix to be decomposed.
upper_a (optional): A logical flag indicating whether to use the upper half (.true.) or the lower half (.false.) of $A$ for the computation. The default is lower.
err (optional): A type(la_state) variable to handle errors. If not provided, execution will stop on errors.

Return value

lambda: A real array containing the computed eigenvalues of the matrix $A$ .

Errors

Raises LINALG_VALUE_ERROR if the matrix has invalid size.
Raises LINALG_ERROR if the eigendecomposition fails.
If err is not provided, execution will stop on errors.

Notes

The eigenvalues are returned in a real array, with the same kind as the input matrix a.
This routine is based on LAPACK's SYEV and HEEV routines.

solve - Solve a linear matrix equation or a linear system of equations.

Syntax

x = solve(a, b [, overwrite_a] [, err])

Description

Solve linear systems - one (b(:)) or many (b(:,:)).

Arguments

a: A real or complex coefficient matrix. If overwrite_a=.true., it is destroyed by the call.
b: A rank-1 (one system) or rank-2 (many systems) array of the same kind as a, containing the right-hand-side vector(s).
overwrite_a (optional, default = .false.): If .true., input matrix a will be used as temporary storage and overwritten, to avoid internal data allocation.
err (optional): A type(la_state) variable.

Return value

For a full-rank matrix, returns an array value that represents the solution to the linear system of equations.

Errors

Raises LINALG_ERROR if the matrix is singular to working precision.
Raises LINALG_VALUE_ERROR if the matrix and rhs vectors have invalid/incompatible sizes.
If err is not present, exceptions trigger an error stop.

lstsq - Compute a least squares solution to a system of linear equations.

Syntax

x = lstsq(a, b [, cond] [, overwrite_a] [, rank] [, err])

Description

Solves the least-squares problem for the system $A \cdot x = b$ , where $A$ is a square matrix of size $n \times n$ and $b$ is either a vector of size $n$ or a matrix of size $n \times nrhs$ . The function minimizes the 2-norm $\|b - A \cdot x\|$ by solving for $x$ .

The result $x$ is returned as an allocatable array, and it is either a vector (for a single right-hand side) or a matrix (for multiple right-hand sides).

Arguments

a: A real matrix of size $n \times n$ representing the coefficient matrix. If overwrite_a = .true., the contents of a may be modified during the computation.
b: A real vector or matrix representing the right-hand side. The size should be $n$ (for a single right-hand side) or $n \times nrhs$ (for multiple right-hand sides).
cond (optional): A cutoff for rank evaluation. Singular values $s(i)$ such that $s(i) \leq \text{cond} \cdot \max(s)$ are considered zero.
overwrite_a (optional, default = .false.): If .true., both a and b may be overwritten and destroyed during computation.
rank (optional): An integer variable that returns the rank of the matrix $A$ .
err (optional): A type(la_state) variable that returns the error state. If err is not provided, the function will stop execution on error.

Return value

Returns the solution array $x$ with size $n$ (for a single right-hand side) or $n \times nrhs$ (for multiple right-hand sides).

Errors

Raises LINALG_ERROR if the matrix $A$ is singular to working precision.
Raises LINALG_VALUE_ERROR if the matrix a and the right-hand side b have incompatible sizes.
If err is not provided, the function stops execution on error.

Notes

This function relies on LAPACK's least-squares solvers, such as GELSS.
If overwrite_a is enabled, the original contents of a and b may be lost.

det - Determinant of a scalar or rectangular matrix.

Syntax

d = det(a [, overwrite_a] [, err])

Description

This function computes the determinant of a square matrix $A$ . The matrix must be a real matrix of size $[m, n]$ , and the determinant is computed using an efficient factorization method (e.g., LU decomposition).

Arguments

a: A real matrix of size $[m, n]$ , representing the rectangular matrix for which the determinant is calculated. If overwrite_a, it is an inout argument and may be modified during computation.
overwrite_a (optional, default = .false.): A logical flag that determines whether the input matrix a can be overwritten. If .true., the matrix a may be destroyed and modified in place to save memory.
err (optional): A state return flag of type(la_state). If an error occurs and err is not provided, the function will stop execution.

Return value

The function returns a real scalar value representing the determinant of the input matrix $A$ , with the same kind as $A$ .

Errors

Raises LINALG_VALUE_ERROR if the matrix a is not square.
If err is not provided, the function will stop execution on errors.

Notes

The determinant of the matrix is computed using the LAPACK getrf backend.
If overwrite_a is enabled, the input matrix a will be destroyed during the computation process.

inv - Inverse of a square matrix.

Syntax

inv_a = inv(a [, err])

Description

This function computes the inverse $A^{-1}$ of a real or complex square matrix $A$ , provided that $A$ is non-singular. The inverse of a matrix is defined as:

$A A^{-1} = A^{-1} A = I$

where $I$ is the identity matrix of the same size as $A$ . The inverse exists only if $A$ is square and has full rank (i.e., all its singular values are nonzero).

The computation is performed using LU decomposition.

Arguments

a: A real or complex square matrix of size $[n,n]$ , representing the matrix to be inverted.
err (optional): A type(la_state) variable that returns the error state. If not provided, the function will stop execution on error.

Return value

inv_a: A real or complex square matrix of size $[n,n]$ , representing the inverse of a.

Errors

Raises LINALG_VALUE_ERROR if a is singular or has invalid size.
If err is not provided, exceptions will trigger an error stop.

Notes

This function computes the inverse using LAPACK's LU decomposition routine GETRF followed by GETRI.
The inverse should be used with caution in numerical computations. For solving linear systems, using solve is usually more stable and efficient than explicitly computing the inverse.

invert - In-place matrix inversion

Syntax

call invert(a [, err])

Description

This subroutine computes the inverse $A^{-1}$ of a real or complex square matrix $A$ in-place, modifying a directly. It uses the LU decomposition method via LAPACK's GETRF and GETRI routines.

Given a square matrix $A$ , the LU decomposition factorizes it as:

$A = P L U$

where:

$P$ is a permutation matrix,
$L$ is a lower triangular matrix with unit diagonal,
$U$ is an upper triangular matrix.

The inverse is then obtained by solving $A X = I$ using the LU factors.

Arguments

a: A real or complex square matrix of size $[n,n]$ . On output, it is replaced with its inverse $A^{-1}$ .
err (optional): A type(la_state) variable that returns the error state. If not provided, the function will stop execution on error.

Errors

Raises LINALG_ERROR if the matrix is singular.
Raises LINALG_VALUE_ERROR if a has invalid size.
If err is not provided, exceptions will trigger an error stop.

Notes

This subroutine modifies a in-place. If the original matrix needs to be preserved, use inv instead.
The determinant of a can be computed before inversion using det to check for singularity.
The computational complexity is $O(n^3)$ , making it expensive for large matrices.
It is recommended to use matrix factorizations (e.g., LU or QR) for solving linear systems instead of computing the inverse explicitly, as it is numerically more stable and efficient.

la_inverse::operator(.inv.) "operator(.inv.)" - Compute the inverse of a square matrix.

Syntax

inva = .inv. a

Description

This operator computes the inverse $A^{-1}$ of a square, non-singular real or complex matrix $A$ using an LU decomposition. The inversion satisfies:

$A A^{-1} = I$

where $I$ is the identity matrix of appropriate size.

This operator is functionally equivalent to inv but provides a more convenient syntax. It supports operator chaining, allowing multiple inversions within expressions:

Arguments

A: A real or complex square matrix of size $[n,n]$ , representing the input matrix to be inverted.

Return value

invA: A real or complex square matrix of size $[n,n]$ , and same kind as A representing its inverse.
If A is singular or the inversion fails, an empty matrix (size $[0,0]$ ) is returned instead of raising an error.

Errors

Unlike inv, this operator does not provide explicit error handling.
If A is singular or an error occurs during inversion, the function returns an empty matrix (size $[0,0]$ ) instead of raising an exception.
The caller should check the size of the returned matrix to determine if inversion was successful.

Notes

This operator internally calls LAPACK's LU decomposition routine GETRF followed by GETRI.
The chaining property allows for concise expressions but requires caution: if any intermediate inversion fails, subsequent operations may propagate errors due to empty matrix results.
If strict error handling is required, use inv instead.

pinv - Moore-Penrose pseudo-inverse of a matrix (function).

Syntax

pinva = pinv(a [, rtol] [, err])

Description

This function computes the Moore-Penrose pseudo-inverse $A^+$ of a real or complex matrix $A$ using Singular Value Decomposition (SVD). The pseudo-inverse provides a generalization of the inverse for non-square and singular matrices, making it useful for solving least-squares problems and underdetermined systems.

The computation is based on the singular value decomposition (SVD):

$A = U \Sigma V^T$

where $U$ and $V$ are orthogonal matrices, and $\Sigma$ is a diagonal matrix containing the singular values. The pseudo-inverse is computed as:

$A^+ = V \Sigma^+ U^T$

where $\Sigma^+$ is obtained by inverting the nonzero singular values.

Arguments

a: A real or complex matrix of size $[m,n]$ , representing the input matrix to be inverted.
rtol (optional): A real scalar specifying the relative tolerance for singular value truncation. Singular values smaller than rtol * max(singular_values(A)) are set to zero. If not provided, a default machine-precision-based tolerance is used.
err (optional): A type(la_state) variable that returns the error state. If not provided, the function will stop execution on error.

Return value

pinva: A real or complex matrix of size $[n,m]$ , representing the pseudo-inverse of a.

Errors

Raises LINALG_VALUE_ERROR if the SVD decomposition fails or the input matrix has invalid dimensions.
Raises LINALG_ERROR if numerical instability prevents inversion.
If err is not provided, exceptions will trigger an error stop.

Notes

This function computes the pseudo-inverse using LAPACK's SVD decomposition routine *GESVD.
The choice of rtol affects numerical stability and rank estimation: setting it too high may result in an inaccurate inverse, while setting it too low may amplify numerical noise.
This function returns a newly allocated matrix. For an in-place version, use pseudoinvert.

la_pseudoinverse::operator(.pinv.) "operator(.pinv.)" - Compute the Moore-Penrose pseudo-inverse of a matrix.

Syntax

pinva = .pinv. a

Description

This operator computes the Moore-Penrose pseudo-inverse $A^+$ of a real or complex matrix $A$ using Singular Value Decomposition (SVD). The pseudo-inverse is useful for solving least-squares problems and handling singular or underdetermined systems.

Given the singular value decomposition (SVD):

$A = U \Sigma V^T$

the pseudo-inverse is computed as:

$A^+ = V \Sigma^+ U^T$

where $\Sigma^+$ is the inverse of the nonzero singular values in $\Sigma$ .

Arguments

a: A real or complex matrix of size $[m,n]$ , representing the input matrix to be inverted.

Return value

pinva: A real or complex matrix of size $[n,m]$ , representing the pseudo-inverse of a.

Errors

Raises LINALG_VALUE_ERROR if the SVD decomposition fails or the input matrix has invalid dimensions.
Raises LINALG_ERROR if numerical instability prevents inversion.
If an error occurs, execution will stop.

Notes

This operator internally calls pinv and behaves identically.
The pseudo-inverse is computed using LAPACK's SVD decomposition routine GESVD.
This operator is a convenient shorthand for calling the functional interface pinv(a).

svd - Singular Value Decomposition (SVD) of a matrix

Syntax

call svd(a, s [, u] [, vt] [, overwrite_a] [, full_matrices] [, err])

Description

This subroutine computes the Singular Value Decomposition (SVD) of a matrix $A$ :

$A = U \cdot S \cdot V^T$

where:

$A$ is the input matrix of size $[m,n]$ .
$U$ is an orthogonal matrix of size $[m,m]$ (or $[m,k]$ for the reduced problem), containing the left singular vectors of $A$ .
$S$ is a diagonal matrix containing the singular values of size $[k,k]$ with $k = \min(m,n)$ .
$V^T$ is an orthogonal matrix of size $[n,n]$ (or $[k,n]$ for the reduced problem), containing the right singular vectors of $A^T$ .

The singular values are returned in the array $S$ , and optionally, the matrices $U$ and $V^T$ are computed and returned.

Arguments

a: A real matrix of size $[m,n]$ representing the input matrix $A$ . If overwrite_a = .true., this matrix may be modified during computation. This is an inout argument.
s: A real array of size $k = \min(m,n)$ , containing the singular values of $A$ . This is an output argument.
u: An optional real matrix of the same type and kind as a, representing the left singular vectors of $A$ . This has shape $[m,m]$ for the full problem or $[m,k]$ for the reduced problem. This is an output argument.
vt: An optional real matrix of the same type and kind as a, representing the right singular vectors of $A^T$ . This has shape $[n,n]$ for the full problem or $[k,n]$ for the reduced problem. This is an output argument.
overwrite_a: (Optional, default = .false.) A logical flag indicating whether the input matrix a may be overwritten during computation. If .true., a is overwritten to avoid additional memory allocation.
full_matrices: (Optional, default = .true.) A logical flag that determines whether to compute full-sized matrices $U$ and $V^T$ (shape $[m,m]$ and $[n,n]$ ). If .false., computes reduced matrices of shape $[m,k]$ and $[k,n]$ .
err: (Optional) A type(la_state) variable to capture the error state. If not provided, the function will stop execution on error.

Return value

The SVD of matrix $A$ is returned in the corresponding output arguments, with the singular values in s, and optionally, the matrices $U$ and $V^T$ .

Errors

Raises LINALG_VALUE_ERROR if the sizes of the matrices are incompatible with the full/reduced problem.
Raises LINALG_ERROR if there is insufficient storage space.
If err is not provided, exceptions will trigger an error stop.

Notes

This subroutine computes the Singular Value Decomposition using LAPACK's GESDD algorithm.
If overwrite_a is enabled, the input matrix a may be overwritten during computation.

svdvals - Singular Values Computation (function).

Syntax

s = svdvals(a [, err])

Description

This function computes the singular values of a real or complex matrix $A$ and returns them in a vector $s$ , where $s$ is an array of size $k = \min(m, n)$ . Singular values are non-negative values that provide important insights into the properties of the matrix, such as its rank and conditioning.

This function does not compute the full Singular Value Decomposition (SVD); instead, it directly calculates and returns only the singular values of matrix $A$ .

The Singular Value Decomposition of a matrix $A$ is expressed as:

$A = U \cdot S \cdot V^T$

where:

$A$ is the input matrix of size $[m,n]$ ,
$U$ and $V$ are orthogonal matrices,
$S$ is a diagonal matrix with the singular values of $A$ .

Arguments

a: A real or complex matrix of size $[m,n]$ , representing the input matrix whose singular values are to be computed.
err (optional): A type(la_state) variable that returns the error state. If not provided, the function will stop execution on error.

Return value

s: A real array containing the singular values of the matrix $A$ , with the same type and kind as the input matrix. The size of the array is $k = \min(m, n)$ .

Errors

Raises LINALG_VALUE_ERROR if the input matrix has invalid dimensions or if the SVD computation fails.
If err is not provided, exceptions will trigger an error stop.

Notes

This function only computes the singular values and does not compute the full SVD (i.e., matrices $U$ and $V$ are not computed).
The singular values are returned as a vector, sorted in decreasing order.
The function uses LAPACK's GESDD routine for singular value computation.

diag - Diagonal matrix.

Syntax

d = diag(n, source [, err]) for scalar input d = diag(source(:) [, err]) for array input

Description

This function generates a square diagonal matrix where the diagonal elements are populated either by a scalar value or an array of values. The size of the matrix is determined by the input parameter $n$ or the size of the input array. If a scalar is provided, the diagonal elements are all set to the same value. If an array is provided, its length determines the size of the matrix, and its elements are placed along the diagonal.

Arguments

n: The size of the square matrix (only used if a scalar is provided for the diagonal).
source:
- If a scalar, this value is used to populate all the diagonal elements of the matrix.
- If an array, the elements of the array are used to populate the diagonal of the matrix. The size of the array determines the matrix size.
err (optional): A state return flag of type(la_state). If an error occurs and err is not provided, the function will stop execution.

Return value

The function returns a matrix of size $n \times n$ , where the diagonal elements are either all equal to the scalar source or populated by the values from the input array.

Errors

Raises LINALG_VALUE_ERROR if the dimensions of the matrix are invalid or if the array size does not match the expected matrix size.
If err is not provided, the function will stop execution on errors.

Notes

The diagonal elements are set to the specified scalar or the array values in the order they appear in the input.
If the err parameter is provided, the error state of the function will be returned.

eye - Identity matrix.

Syntax

eye = eye(m [, n] [, mold] [, err])

Description

This function constructs an identity matrix of size $m \times n$ , where the diagonal elements are set to 1 and all off-diagonal elements are set to 0. If only the number of rows $m$ is provided, a square matrix of size $m \times m$ is returned. The matrix is populated with a real data type, by default real(real64), or a type specified by the user.

Arguments

m: The number of rows of the identity matrix.
n (optional): The number of columns of the identity matrix. If omitted, the matrix is square ( $m \times m$ ).
mold (optional): The data type to define the return type. Defaults to real(real64).
err (optional): A state return flag of type(la_state). If an error occurs and err is not provided, the function will stop execution.

Return value

The function returns a matrix of size $m \times n$ (or $m \times m$ if $n$ is omitted) with diagonal elements set to 1 and all other elements set to 0.

Errors

Raises LINALG_VALUE_ERROR if the dimensions of the matrix are invalid (e.g., negative values).
If err is not provided, the function will stop execution on errors.

Notes

The identity matrix is constructed with the specified data type, which defaults to real(real64) if no type is specified.
The mold scalar is used to provide a function return type.
If the err parameter is provided, the error state of the function will be returned.

qr - QR factorization of a matrix.

Syntax

call qr(a, q, r [, overwrite_a] [, storage] [, err])

Description

This subroutine computes the QR factorization of a real or complex matrix $A = Q \cdot R$ , where $Q$ is orthonormal and $R$ is upper-triangular. The matrix $A$ has size $[m,n]$ with $m \ge n$ . The result is returned in the output matrices $Q$ and $R$ , which have the same type and kind as $A$ .

Given $k = \min(m, n)$ , the matrix $A$ can be written as:

$A = \left( \begin{array}{cc} Q_1 & Q_2 \end{array} \right) \cdot \left( \begin{array}{cc} R_1 & 0 \end{array} \right)$

Because the lower rows of $R$ are zeros, a reduced problem $A = Q_1 R_1$ can be solved. The size of the input matrices determines which problem is solved:

For full matrices (shape(Q) == [m,m], shape(R) == [m,n]), the full problem is solved.
For reduced matrices (shape(Q) == [m,k], shape(R) == [k,n]), the reduced problem is solved.

Arguments

a: A real or complex matrix of size $[m,n]$ , representing the coefficient matrix. If overwrite_a = .false., this is an input argument. If overwrite_a = .true., it is an inout argument and is overwritten upon return.
q: A rank-2 array of the same type and kind as a, representing the orthonormal matrix $Q$ . This is an output argument with shape $[m,m]$ (for the full problem) or $[m,k]$ (for the reduced problem).
r: A rank-2 array of the same type and kind as a, representing the upper-triangular matrix $R$ . This is an output argument with shape $[m,n]$ (for the full problem) or $[k,n]$ (for the reduced problem).
storage (optional): A rank-1 array of the same type and kind as a, providing working storage for the solver. Its minimum size can be determined by a call to qr_space. This is an output argument.
overwrite_a (optional, default = .false.): A logical flag that determines whether the input matrix a can be overwritten. If .true., the matrix a is used as temporary storage and overwritten to avoid internal memory allocation. This is an input argument.
err (optional): A type(la_state) variable that returns the error state. If not provided, the function will stop execution on error.

Return value

The QR factorization matrices $Q$ and $R$ are returned in the corresponding arguments.

Errors

Raises LINALG_VALUE_ERROR if the sizes of the matrices are incompatible with the full/reduced problem.
Raises LINALG_ERROR if there is insufficient storage space.
If err is not provided, exceptions will trigger an error stop.

Notes

This subroutine computes the QR factorization using LAPACK's QR decomposition algorithm *GEQRF.
If overwrite_a is enabled, the input matrix a will be modified during computation.

qr_space - Workspace size for QR operations.

Syntax

call qr_space(a, lwork [, err])

Description

This subroutine computes the minimum workspace size required for performing QR factorization. The size of the workspace array needed for both QR factorization and solving the reduced problem is determined based on the input matrix $A$ .

The input matrix $A$ has size $[m,n]$ , and the output value $lwork$ represents the minimum size of the workspace array that should be allocated for QR operations.

Arguments

a: A real or complex matrix of size $[m,n]$ , representing the input matrix used to determine the required workspace size.
lwork: An integer variable that will return the minimum workspace size required for QR factorization.
err (optional): A type(la_state) variable that returns the error state. If not provided, the function will stop execution on error.

Return value

The workspace size $lwork$ that should be allocated before calling the QR factorization routine is returned.

Errors

Raises LINALG_ERROR if there is an issue determining the required workspace size.
If err is not provided, exceptions will trigger an error stop.

Notes

This subroutine is useful for preallocating memory for QR factorization in large systems.
It is important to ensure that the workspace size is correctly allocated before proceeding with QR factorization to avoid memory issues.

schur - Schur decomposition of a matrix.

Syntax

call schur(a, t, z [, eigvals] [, overwrite_a] [, storage] [, err])

Description

This subroutine computes the Schur decomposition of a real or complex matrix $A = Z T Z^H$ , where $Z$ is an orthonormal/unitary matrix, and $T$ is an upper-triangular or quasi-upper-triangular matrix. The matrix $A$ has size $[m,m]$ .

The decomposition produces:

$T$ , which is upper-triangular for complex matrices and quasi-upper-triangular for real matrices (with possible $2 \times 2$ blocks on the diagonal).
$Z$ , the transformation matrix, which is optional.
Optionally, the eigenvalues corresponding to the diagonal elements of $T$ .

If a pre-allocated workspace is provided, no internal memory allocations take place.

Arguments

a: A real or complex matrix of size $[m,m]$ . If overwrite_a = .false., this is an input argument. If overwrite_a = .true., it is an inout argument and is overwritten upon return.
t: A rank-2 array of the same type and kind as a, representing the Schur form of a. This is an output argument with shape $[m,m]$ .
z (optional): A rank-2 array of the same type and kind as a, representing the unitary/orthonormal transformation matrix $Z$ . This is an output argument with shape $[m,m]$ .
eigvals (optional): A complex array of size $[m]$ , representing the eigenvalues that appear on the diagonal of $T$ . This is an output argument.
storage (optional): A rank-1 array of the same type and kind as a, providing working storage for the solver. Its minimum size can be determined by a call to schur_space. This is an input argument.
overwrite_a (optional, default = .false.): A logical flag that determines whether the input matrix a can be overwritten. If .true., the matrix a is used as temporary storage and overwritten to avoid internal memory allocation. This is an input argument.
err (optional): A type(la_state) variable that returns the error state. If not provided, the function will stop execution on error.

Return value

The Schur decomposition matrices $T$ and optionally $Z$ are returned in the corresponding arguments.

Errors

Raises LINALG_VALUE_ERROR if the sizes of the matrices are incompatible.
Raises LINALG_ERROR if the algorithm did not converge.
If err is not provided, exceptions will trigger an error stop.

Notes

This subroutine computes the Schur decomposition using LAPACK's Schur decomposition routines (GEES).
Sorting options for eigenvalues can be requested, utilizing LAPACK's eigenvalue sorting mechanism.
If overwrite_a is enabled, the input matrix a will be modified during computation.

schur_space - Workspace size for Schur decomposition.

Syntax

call schur_space(a, lwork [, err])

Description

This subroutine computes the minimum workspace size required for performing Schur decomposition. The size of the workspace array needed is determined based on the input matrix $A$ .

The input matrix $A$ has size $[m,m]$ , and the output value $lwork$ represents the minimum size of the workspace array that should be allocated for Schur decomposition operations.

Arguments

a: A real or complex matrix of size $[m,m]$ , representing the input matrix used to determine the required workspace size.
lwork: An integer variable that will return the minimum workspace size required for Schur decomposition.
err (optional): A type(la_state) variable that returns the error state. If not provided, the function will stop execution on error.

Return value

The workspace size $lwork$ that should be allocated before calling the Schur decomposition routine is returned.

Errors

Raises LINALG_ERROR if there is an issue determining the required workspace size.
If err is not provided, exceptions will trigger an error stop.

Notes

This subroutine is useful for preallocating memory for Schur decomposition in large systems.
It is important to ensure that the workspace size is correctly allocated before proceeding with Schur decomposition to avoid memory issues.

BLAS, LAPACK

Modern Fortran modules with full explicit typing features are available as modules la_blas and la_lapack. The reference Fortran-77 library, forked from Release 3.10.1, was automatically processed and modernized. The following refactorings are applied:

All datatypes and accuracy constants standardized into a module (stdlib-compatible names)
Both libraries available for 32, 64 and 128-bit floats
Free format, lower-case style
implicit none(type, external) everywhere
all pure procedures where possible
intent added to all procedure arguments
Removed DO 10 .... 10 CONTINUE, replaced with do..end do loops or labelled loop_10: do ... cycle loop_10 ... end do loop_10 in case control statements are present
BLAS modularized into a single-file module
LAPACK modularized into a single-file module
All procedures prefixed (with stdlib_, currently).
F77-style parameters removed, and numeric constants moved to the top of each module.
Ambiguity in single vs. double precision constants (0.0, 0.d0, (1.0,0.0)) removed
preprocessor-based OpenMP directives retained.

The single-source module structure hopefully allows for cross-procedural inlining which is otherwise impossible without link-time optimization.

Building

An automated build is currently available via the Fortran Package Manager. To add fortran-lapack to your project, simply add it as a dependency:

[dependencies]

fortran-lapack = { git="https://github.com/perazz/fortran-lapack.git" }

fortran-lapack is compatible with the LAPACK API. If high-performance external BLAS/LAPACK libraries are available, it is sufficient to define macros

[dependencies]

fortran-lapack = { git="https://github.com/perazz/fortran-lapack.git", preprocess.cpp.macros=["LA_EXTERNAL_BLAS", "LA_EXTERNAL_LAPACK"] }

Extension to external BLAS/LAPACK libraries

Generic interfaces to most BLAS/LAPACK functions are exposed to modules la_blas and la_lapack. These interfaces drop the initial letter to wrap a precision-agnostic version. For example, axpy is a precision-agnostic interface to saxpy, daxpy, caxpy, zaxpy, qaxpy, waxpy. The naming convention is:

Type	32-bit	64-bit	128-bit
real	`s`	`d`	`q`
complex	`c`	`z`	`w`

All public interfaces in la_blas and la_lapack allow seamless linking against external libraries via a simple pre-processor flag. When an external library is available, just define macros LA_EXTERNAL_BLAS and LA_EXTERNAL_LAPACK. The kind-agnostic interface will just point to the external function. All such interfaces follow this template:

interface axpy
#ifdef LA_EXTERNAL_BLAS
    ! Use external library
    pure subroutine saxpy(n, a, x, incx, y, incy)
      import :: ik, sp
      integer, parameter :: wp = sp
      integer(ik), intent(in) :: n
      real(wp), intent(in) :: a
      real(wp), intent(in) :: x(*)
      integer(ik), intent(in) :: incx
      real(wp), intent(inout) :: y(*)
      integer(ik), intent(in) :: incy
    end subroutine saxpy
#else
    ! Use internal implementation
    module procedure la_saxpy
#endif
end interface

Licensing

LAPACK is a freely-available software package. It is available from netlib via anonymous ftp and the World Wide Web. Thus, it can be included in commercial software packages (and has been). Credit for the library should be given to the LAPACK authors. The license used for the software is the modified BSD license. According to the original license, we changed the name of the routines and commented the changes made to the original.

Acknowledgments

Part of this work was supported by the Sovereign Tech Fund.