- Generated by
1.13.2
|
fortran-lapack
|
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.
All procedures work with all types (real, complex) and kinds (32, 64, 128-bit floats).
c = chol(a [, lower] [, other_zeroed])
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.
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..c: A real or complex matrix of size \( [n,n] \), containing the Cholesky factors. The returned matrix is triangular (upper or lower, as selected).error stop.call cholesky(a [, c] [, lower] [, other_zeroed])
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.
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..error stop.call eig(a [, b] [, lambda] [, right] [, left] [, overwrite_a] [, overwrite_b] [, err])
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.
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.err is not provided, exceptions will trigger an error stop.A or B can improve performance but destroys the original matrix data.call eigh(a, lambda [, vectors] [, upper_a] [, overwrite_a] [, err])
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.
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.err is not provided, exceptions will trigger an error stop.A can improve performance but destroys the original matrix data.lambda = eigvals(a [, b] [, err])
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.
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.lambda: A complex array of eigenvalues, computed from the input matrix \(A\) (and \(B\) if in the generalized case).err is not provided, exceptions will trigger an error stop.lambda = eigvalsh(a [, upper_a] [, err])
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\).
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.
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.lambda: A real array containing the computed eigenvalues of the matrix \(A\).err is not provided, execution will stop on errors.real array, with the same kind as the input matrix a.x = solve(a, b [, overwrite_a] [, err])
Solve linear systems - one (b(:)) or many (b(:,:)).
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.For a full-rank matrix, returns an array value that represents the solution to the linear system of equations.
err is not present, exceptions trigger an error stop.call solve_lu(a, b, x [, pivot] [, overwrite_a] [, err])
Solve linear systems - one (b(:)) or many (b(:,:)) - writing the result into the caller's array x instead of returning a new one. Storage for the pivot indices may be supplied as well: when x and pivot are both provided and overwrite_a=.true., the call performs no internal allocation, which makes it suited to a loop over many systems of the same size. The routine is pure.
a: A real or complex coefficient matrix of size \( [n,n] \). 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).x: An array of the same shape and kind as b. On output it holds the solution.pivot (optional): An integer(ilp) array of size n that receives the diagonal pivot indices of the LU factorization.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.a, b, x or pivot have invalid/incompatible sizes.err is not present, exceptions trigger an error stop.x.call solve_chol(a, b, x [, lower] [, overwrite_a] [, err])
Factorize a real symmetric or complex Hermitian positive definite matrix and solve \( A x = b \) in one call, for one (b(:)) or many (b(:,:)) right-hand sides. Only the triangle lower selects is read. The result is written into the caller's array x. The routine is pure.
a: A real symmetric or complex Hermitian positive definite matrix of size \( [n,n] \). If overwrite_a=.true., it is overwritten with its Cholesky factor.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).x: An array of the same shape and kind as b. On output it holds the solution.lower (optional, default = .true.): If .true., the lower triangle of a is read and the factorization is \( A = L L^H \); otherwise the upper triangle is read and the factorization is \( A = U^H U \).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.a is not positive definite.a, b or x have invalid/incompatible sizes.err is not present, exceptions trigger an error stop.call solve_lower_chol(l, b, x [, err])
call solve_upper_chol(u, b, x [, err])
Solve \( A x = b \) for one or many right-hand sides from a Cholesky factor computed earlier, without factorizing again. Each call costs two triangular solves. Both routines are pure.
l / u: The lower or upper Cholesky factor of size \( [n,n] \), as returned by cholesky with lower=.true. or lower=.false..b: A rank-1 (one system) or rank-2 (many systems) array of the same kind as the factor, containing the right-hand-side vector(s).x: An array of the same shape and kind as b. On output it holds the solution.err (optional): A type(la_state) variable.b or x have invalid/incompatible sizes.err is not present, exceptions trigger an error stop.x = lstsq(a, b [, cond] [, overwrite_a] [, rank] [, err])
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).
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.Returns the solution array \( x \) with size \( n \) (for a single right-hand side) or \( n \times nrhs \) (for multiple right-hand sides).
a and the right-hand side b have incompatible sizes.err is not provided, the function stops execution on error.overwrite_a is enabled, the original contents of a and b may be lost.call lstsq_space(a, b, lrwork, liwork) for real data call lstsq_space(a, b, lrwork, liwork, lcwork) for complex data
This subroutine returns the sizes of the working arrays that solve_lstsq needs for a problem of the shape of a and b, so that a repeated solve of problems of the same size performs no internal allocation.
a: A real or complex matrix of size \( [m,n] \). Only its shape is read.b: A right-hand-side vector of size \(m\) or matrix of size \( [m,nrhs] \). Only its shape is read.lrwork: An integer returning the size of the real working array.liwork: An integer returning the size of the integer working array.lcwork: An integer returning the size of the complex working array. Complex data only.The three (two, for real data) workspace sizes are returned in the output arguments.
pure and cannot fail: it reads the shapes of its inputs and nothing else.call solve_lstsq(a, b, x [, real_storage] [, int_storage] [, cmpl_storage] [, cond] [, singvals] [, overwrite_a] [, rank] [, err])
This subroutine computes the least-squares solution of \( A \cdot x = b \) and writes it into the caller's array. Given the optional working arrays, it allocates nothing, so a repeated solve of problems of the same size runs without any memory traffic.
a: A real or complex matrix of size \( [m,n] \). It is an intent(inout) argument, and is destroyed if overwrite_a is true.b: The right-hand-side vector of size \(m\) or matrix of size \( [m,nrhs] \). It is an intent(in) argument.x: The solution vector of size \( \ge n \) or matrix of size \( [\ge n,nrhs] \). It is an intent(inout) contiguous argument.real_storage (optional): A real working array of size at least the lrwork returned by lstsq_space.int_storage (optional): An integer working array of size at least liwork.cmpl_storage (optional): A complex working array of size at least lcwork. Complex data only.cond (optional): The cutoff for rank evaluation: singular values \( s_i \le \text{cond} \cdot \max(s) \) are treated as zero.singvals (optional): A real array of size at least \( \min(m,n) \) returning the singular values in decreasing order.overwrite_a (optional): If true, a may be overwritten and destroyed. Default is false.rank (optional): An integer returning the rank of a.err (optional): A state return flag of type(la_state). If an error occurs and err is not provided, the subroutine will stop execution.The solution is written into x.
x is too small, or if a working array is too small.err is not provided, the subroutine will stop execution on errors.lstsq is this subroutine with the solution allocated for the caller.x = weighted_lstsq(w, a, b [, cond] [, overwrite_a] [, rank] [, err])
This function minimizes \( \|D (b - A \cdot x)\| \) with \( D = \mathrm{diag}(\sqrt{w}) \): the \(i\)-th equation carries the weight \( w_i \). Both sides are scaled and the resulting ordinary least-squares problem is solved with GELSD.
w: A real vector of size \(m\). The weights are always real and must all be positive.a: A real or complex matrix of size \( [m,n] \). It is an intent(inout) argument, and is destroyed if overwrite_a is true.b: The right-hand-side vector of size \(m\). It is an intent(in) argument.cond (optional): The cutoff for rank evaluation.overwrite_a (optional): If true, a may be overwritten and destroyed. Default is false.rank (optional): An integer returning the rank of a.err (optional): A state return flag of type(la_state). If an error occurs and err is not provided, the function will stop execution.The function returns the solution vector of size \(n\).
w or b does not have one entry per row of a, or if any weight is not positive.err is not provided, the function will stop execution on errors.lstsq.call solve_weighted_lstsq(w, a, b, x [, cond] [, overwrite_a] [, rank] [, err])
This subroutine is the subroutine form of weighted_lstsq: it writes the solution into the caller's array instead of allocating it.
w: A real vector of size \(m\) of positive weights.a: A real or complex matrix of size \( [m,n] \). It is an intent(inout) argument, and is destroyed if overwrite_a is true.b: The right-hand-side vector of size \(m\). It is an intent(in) argument.x: The solution vector of size \(n\). It is an intent(inout) contiguous argument.cond (optional): The cutoff for rank evaluation.overwrite_a (optional): If true, a may be overwritten and destroyed. Default is false.rank (optional): An integer returning the rank of a.err (optional): A state return flag of type(la_state). If an error occurs and err is not provided, the subroutine will stop execution.The solution is written into x.
w or b does not have one entry per row of a, or if any weight is not positive.err is not provided, the subroutine will stop execution on errors.a unless overwrite_a is set, so a is unchanged by default.x = constrained_lstsq(a, b, c, d [, overwrite_matrices] [, err])
This function minimizes \( \|b - A \cdot x\| \) subject to \( C \cdot x = d \), with \( A \) of size \( [m,n] \) and \( C \) of size \( [p,n] \). The problem has a unique solution when \( p \le n \le m+p \), \( \mathrm{rank}(C) = p \) and the stacked matrix has rank \( n \).
a: A real or complex matrix of size \( [m,n] \), the least-squares cost. It is an intent(inout) argument, and is destroyed if overwrite_matrices is true.b: The least-squares right-hand-side vector of size \(m\). It is an intent(inout) argument.c: The constraint matrix of size \( [p,n] \). It is an intent(inout) argument.d: The constraint right-hand-side vector of size \(p\). It is an intent(inout) argument.overwrite_matrices (optional): If true, a, b, c and d may be overwritten and destroyed. Default is false.err (optional): A state return flag of type(la_state). If an error occurs and err is not provided, the function will stop execution.The function returns the solution vector of size \(n\).
c is rank deficient or if the stacked matrix is.err is not provided, the function will stop execution on errors.call solve_constrained_lstsq(a, b, c, d, x [, storage] [, overwrite_matrices] [, err])
This subroutine is the subroutine form of constrained_lstsq: it writes the solution into the caller's array and can reuse a caller-provided workspace.
a: A real or complex matrix of size \( [m,n] \), the least-squares cost. It is an intent(inout) argument.b: The least-squares right-hand-side vector of size \(m\). It is an intent(inout) argument.c: The constraint matrix of size \( [p,n] \). It is an intent(inout) argument.d: The constraint right-hand-side vector of size \(p\). It is an intent(inout) argument.x: The solution vector of size \(n\). It is an intent(out) argument.storage (optional): A working array of size at least the lwork returned by constrained_lstsq_space.overwrite_matrices (optional): If true, a, b, c and d may be overwritten and destroyed. Default is false.err (optional): A state return flag of type(la_state). If an error occurs and err is not provided, the subroutine will stop execution.The solution is written into x.
storage is too small.c is rank deficient or if the stacked matrix is.err is not provided, the subroutine will stop execution on errors.overwrite_matrices, all four inputs are copied, so none of them is changed.call constrained_lstsq_space(a, c, lwork [, err])
This subroutine asks LAPACK for the optimal size of the workspace array that solve_constrained_lstsq needs for a problem of the shape of a and c.
a: A real or complex matrix of size \( [m,n] \). Only its shape is read.c: The constraint matrix of size \( [p,n] \). Only its shape is read.lwork: An integer returning the size of the workspace array.err (optional): A state return flag of type(la_state). If an error occurs and err is not provided, the subroutine will stop execution.err is not provided, the subroutine will stop execution on errors.d = det(a [, overwrite_a] [, err])
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).
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.The function returns a real scalar value representing the determinant of the input matrix \( A \), with the same kind as \( A \).
a is not square.err is not provided, the function will stop execution on errors.overwrite_a is enabled, the input matrix a will be destroyed during the computation process.This operator computes the determinant of a square real or complex matrix \( A \) from its LU factorization, in the same way as det. It is pure, so it can be used inside pure procedures and do concurrent blocks; it takes no overwrite_a flag and never modifies its operand, which is copied internally.
A: A real or complex square matrix of size \( [n,n] \).A scalar of the same type and kind as A, holding its determinant.
err argument, so a non-square or singular matrix triggers an error stop.err argument instead.inv_a = inv(a [, err])
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.
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.inv_a: A real or complex square matrix of size \( [n,n] \), representing the inverse of a.a is singular or has invalid size.err is not provided, exceptions will trigger an error stop.call invert(a [, pivot] [, err])
call invert(a, inva [, pivot] [, err])
This subroutine computes the inverse \( A^{-1} \) of a real or complex square matrix \( A \). The first form works in-place, modifying a directly; the second writes the inverse into a second matrix inva of the same shape and leaves a untouched. Both use 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:
The inverse is then obtained by solving \( A X = I \) using the LU factors.
a: A real or complex square matrix of size \( [n,n] \). In the in-place form it is replaced with its inverse \( A^{-1} \) on output; in the split form it is read only.inva (split form only): A matrix of the same shape and kind as a, which receives the inverse \( A^{-1} \).pivot (optional): An integer(ilp) array of size at least n that receives the diagonal pivot indices of the LU factorization. Supplying it avoids the internal allocation of the pivot array.err (optional): A type(la_state) variable that returns the error state. If not provided, the function will stop execution on error.a has invalid size, if inva does not match the shape of a, or if pivot is shorter than n.err is not provided, exceptions will trigger an error stop.a. If the original matrix needs to be preserved, use the split form or inv instead.a can be computed before inversion using det to check for singularity.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:
A: A real or complex square matrix of size \( [n,n] \), representing the input matrix to be inverted.invA: A real or complex square matrix of size \( [n,n] \), and same kind as A representing its inverse. A is singular or the inversion fails, an empty matrix (size \( [0,0] \)) is returned instead of raising an error.A is singular or an error occurs during inversion, the function returns an empty matrix (size \( [0,0] \)) instead of raising an exception.pinva = pinv(a [, rtol] [, err])
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.
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.pinva: A real or complex matrix of size \( [n,m] \), representing the pseudo-inverse of a.err is not provided, exceptions will trigger an error stop.*GESVD.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.pseudoinvert.pinva = .pinv. a
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 \).
a: A real or complex matrix of size \( [m,n] \), representing the input matrix to be inverted.pinva: A real or complex matrix of size \( [n,m] \), representing the pseudo-inverse of a.pinv(a).call svd(a, s [, u] [, vt] [, overwrite_a] [, full_matrices] [, err])
This subroutine computes the Singular Value Decomposition (SVD) of a matrix \( A \):
\[A = U \cdot S \cdot V^T \]
where:
The singular values are returned in the array \( S \), and optionally, the matrices \( U \) and \( V^T \) are computed and returned.
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.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 \).
err is not provided, exceptions will trigger an error stop.overwrite_a is enabled, the input matrix a may be overwritten during computation.s = svdvals(a [, err])
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: 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.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) \).err is not provided, exceptions will trigger an error stop.d = diag(n, source [, err]) for scalar input d = diag(source(:) [, err]) for array input d = diag(source(:), k) for array input placed on the k-th diagonal v = diag(a(:,:) [, k]) to extract a diagonal of a matrix
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. Given the offset k, the array is placed on the k-th superdiagonal ( \(k>0\)) or subdiagonal ( \(k<0\)) instead, and the matrix grows to \((n+|k|) \times (n+|k|)\). Given a matrix instead of a vector, the function extracts the requested diagonal and returns it as a vector.
n: The size of the square matrix (only used if a scalar is provided for the diagonal).source:a: A matrix whose k-th diagonal is returned as a vector.k (optional): The index of the diagonal: 0 is the main diagonal, \(k>0\) the k-th superdiagonal, \(k<0\) the k-th subdiagonal.err (optional): A state return flag of type(la_state). If an error occurs and err is not provided, the function will stop execution.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. With the offset k, the matrix is of size \((n+|k|) \times (n+|k|)\). Given a matrix, the function returns the requested diagonal as a vector; the vector has size zero when the requested diagonal lies outside the matrix.
err is not provided, the function will stop execution on errors.err parameter is provided, the error state of the function will be returned.eye = eye(m [, n] [, mold] [, err])
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.
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.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.
err is not provided, the function will stop execution on errors.real(real64) if no type is specified.mold scalar is used to provide a function return type.err parameter is provided, the error state of the function will be returned.t = trace(a)
This function returns the sum of the main diagonal elements of a matrix. The matrix does not need to be square: for a \( m \times n \) matrix, the first \( \min(m,n) \) diagonal elements are summed.
a: A real or complex matrix of size \( [m,n] \). It is an intent(in) argument.The function returns a scalar of the same type and kind as a, equal to \( \sum_i a_{ii} \).
pure and cannot fail: an empty matrix returns zero.trace(a) is equivalent to sum(diag(a)), computed without building the intermediate vector.c = outer_product(u, v)
This function returns the outer product \( u \otimes v \) of two vectors, the matrix whose \((i,j)\) element is \( u_i v_j \).
u: A real or complex vector of size \(m\). It is an intent(in) argument.v: A vector of the same type and kind as u, of size \(n\). It is an intent(in) argument.The function returns the \( m \times n \) matrix \( u \otimes v \), of the same type and kind as the inputs.
pure and cannot fail: the two vectors may have any lengths.v. For the Hermitian outer product, pass conjg(v).c = cross_product(a, b)
This function returns the cross product \( a \times b \) of two vectors of size 3, the vector orthogonal to both inputs.
a: A real or complex vector of size 3. It is an intent(in) argument.b: A vector of size 3, of the same type and kind as a. It is an intent(in) argument.The function returns a vector of size 3, of the same type and kind as the inputs.
pure and cannot fail: both arguments are declared of fixed size 3, so a wrong length is a compile-time error.c = kronecker_product(a, b)
This function returns the Kronecker product \( A \otimes B \): given \( A \) of size \( m_1 \times n_1 \) and \( B \) of size \( m_2 \times n_2 \), the result is the \( (m_1 m_2) \times (n_1 n_2) \) block matrix whose \((i,j)\) block is \( A_{ij} B \).
a: A real or complex matrix of size \( [m_1,n_1] \). It is an intent(in) argument.b: A matrix of size \( [m_2,n_2] \), of the same type and kind as a. It is an intent(in) argument.The function returns the \( (m_1 m_2) \times (n_1 n_2) \) Kronecker product matrix, of the same type and kind as the inputs.
pure and cannot fail: the two matrices may have any shapes.kronecker_product(a, b) and kronecker_product(b, a) differ by a permutation of rows and columns.ah = hermitian(a)
This function returns the Hermitian transpose of a matrix: conjg(transpose(a)) for a complex matrix, transpose(a) for a real one.
a: A real or complex matrix of size \( [m,n] \). It is an intent(in) argument.The function returns the \( n \times m \) matrix \( a^H \), of the same type and kind as a.
pure and cannot fail.l = is_square(a)
This function returns .true. if the input matrix has as many rows as columns.
a: A real or complex matrix of size \( [m,n] \). It is an intent(in) argument.The function returns a logical flag, .true. if \( m = n \).
pure and cannot fail.l = is_diagonal(a)
This function returns .true. if every entry of the input matrix outside the main diagonal is exactly zero. The matrix does not need to be square.
a: A real or complex matrix of size \( [m,n] \). It is an intent(in) argument.The function returns a logical flag, .true. if \( a_{ij} = 0 \) for all \( i \neq j \).
pure and cannot fail.l = is_symmetric(a)
This function returns .true. if the input matrix equals its own transpose.
a: A real or complex matrix of size \( [m,n] \). It is an intent(in) argument.The function returns a logical flag, .true. if \( a = a^T \). A non-square matrix returns .false..
pure and cannot fail.complex matrix this is the transpose without conjugation. Use is_hermitian for the conjugate test.l = is_skew_symmetric(a)
This function returns .true. if the input matrix equals the negative of its own transpose.
a: A real or complex matrix of size \( [m,n] \). It is an intent(in) argument.The function returns a logical flag, .true. if \( a = -a^T \). A non-square matrix returns .false..
pure and cannot fail.l = is_hermitian(a)
This function returns .true. if the input matrix equals its own conjugate transpose.
a: A real or complex matrix of size \( [m,n] \). It is an intent(in) argument.The function returns a logical flag, .true. if \( a = a^H \). A non-square matrix returns .false..
pure and cannot fail.real matrix this is the same test as is_symmetric.l = is_triangular(a, uplo [, err])
This function returns .true. if every entry of the input matrix below (‘uplo = 'U’) or above (uplo = 'L'`) the main diagonal is exactly zero. The matrix does not need to be square.
a: A real or complex matrix of size \( [m,n] \). It is an intent(in) argument.uplo: A character flag selecting the triangle to test, ‘'U’for upper or'L'for lower. It is anintent(in)argument. -err<tt>(optional): A state return flag of [type(la_state)](@ref la_state_type::la_state). If an error occurs anderr` is not provided, the function will stop execution.
The function returns a logical flag, .true. if a is triangular of the requested type.
uplo is neither ‘'U’nor'L', and returns.false..Iferr` is not provided, the function will stop execution on errors.err is pure.l = is_hessenberg(a, uplo [, err])
This function returns .true. if every entry of the input matrix more than one row below (‘uplo = 'U’) or more than one row above (uplo = 'L'`) the main diagonal is exactly zero. The matrix does not need to be square.
a: A real or complex matrix of size \( [m,n] \). It is an intent(in) argument.uplo: A character flag selecting the Hessenberg form to test, ‘'U’for upper or'L'for lower. It is anintent(in)argument. -err<tt>(optional): A state return flag of [type(la_state)](@ref la_state_type::la_state). If an error occurs anderr` is not provided, the function will stop execution.
The function returns a logical flag, .true. if a is Hessenberg of the requested type.
uplo is neither ‘'U’nor'L', and returns.false..Iferr` is not provided, the function will stop execution on errors.call qr(a, q, r [, overwrite_a] [, storage] [, err])
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:
shape(Q) == [m,m], shape(R) == [m,n]), the full problem is solved.shape(Q) == [m,k], shape(R) == [k,n]), the reduced problem is solved.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.The QR factorization matrices \( Q \) and \( R \) are returned in the corresponding arguments.
err is not provided, exceptions will trigger an error stop.*GEQRF.overwrite_a is enabled, the input matrix a will be modified during computation.call qr_space(a, lwork [, err])
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.
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.The workspace size \( lwork \) that should be allocated before calling the QR factorization routine is returned.
err is not provided, exceptions will trigger an error stop.call schur(a, t, z [, eigvals] [, overwrite_a] [, storage] [, err])
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:
complex matrices and quasi-upper-triangular for real matrices (with possible \( 2 \times 2 \) blocks on the diagonal).If a pre-allocated workspace is provided, no internal memory allocations take place.
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.The Schur decomposition matrices \( T \) and optionally \( Z \) are returned in the corresponding arguments.
err is not provided, exceptions will trigger an error stop.overwrite_a is enabled, the input matrix a will be modified during computation.call schur_space(a, lwork [, err])
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.
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.The workspace size \( lwork \) that should be allocated before calling the Schur decomposition routine is returned.
err is not provided, exceptions will trigger an error stop.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:
stdlib-compatible names)implicit none(type, external) everywherepure procedures where possibleintent added to all procedure argumentsDO 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 presentstdlib_, currently).parameters removed, and numeric constants moved to the top of each module.0.0, 0.d0, (1.0,0.0)) removedGrouping every precision of a topic in one module hopefully allows for cross-procedural inlining which is otherwise impossible without link-time optimization.
An automated build is currently available via the Fortran Package Manager. To add fortran-lapack to your project, simply add it as a dependency:
fortran-lapack is compatible with the LAPACK API. If high-performance external BLAS/LAPACK libraries are available, it is sufficient to define macros
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:
The Fortran under src/ and test/ is generated from the kind-templated fypp sources under fypp/src/ and fypp/test/, with the shared kind algebra in include/. The generated files are committed, so building or installing the package never needs fypp. After editing a template, regenerate with
--check is what continuous integration runs; it prints the templates it does not own yet and the reason for each. la_blas and la_lapack are umbrella modules that re-export 57 topic modules, ten for BLAS and 47 for LAPACK; the generic interfaces they publish are data tables under include/, regenerated with python3 scripts/templatize.py --blas-interfaces and python3 scripts/templatize.py --lapack-interfaces. Two further scripts support that layout: scripts/templatize.py converts per-kind Fortran into one template per topic, driven by scripts/la_modules.tsv, and scripts/check_generated.py compares the regenerated tree against a git reference routine by routine.
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.
Part of this work was supported by the Sovereign Tech Fund.
1.13.2