# `stdlib_linalg` ```{contents} :local: :depth: 2 ``` The `stdlib` linear algebra library provides high-level APIs for dealing with common linear algebra operations. ## BLAS and LAPACK ### Status Experimental ### Description `BLAS` and `LAPACK` backends provide efficient low level implementations of many linear algebra algorithms, and are employed for non-trivial operators. A Modern Fortran version of the [Reference-LAPACK 3.10.1](http://github.com/reference-LAPACK) implementation is provided as a backend. Modern Fortran modules with full explicit typing features are provided after an [automated conversion](https://github.com/perazz/fortran-lapack/blob/main/scripts/modularize_blas.py) of the legacy codes: - [stdlib_linalg_blas(module)], [stdlib_linalg_lapack(module)] provide kind-agnostic interfaces to all functions. - Both libraries are available for 32- (`sp`), 64- (`dp`) and 128-bit (`qp`) `real` and `complex` numbers (the latter if available in the current build) - Free format, lower-case style - `implicit none(type, external)` applied to all procedures and modules - `intent` added and all `pure` procedures where possible - `stdlib` provides all procedures in two different flavors: (a) original BLAS/LAPACK names with a prefix `stdlib_?` (ex: `stdlib_dgemv`, `stdlib_sgemv`); (b) A generic, kind agnostic ``, i.e. `gemv`. - F77-style `parameter`s removed, and all numeric constants have been generalized with KIND-dependent Fortran intrinsics. - preprocessor-based OpenMP directives retained. The single-source module structure hopefully allows for cross-procedural inlining which is otherwise impossible without link-time optimization. When available, highly optimized libraries that take advantage of specialized processor instructions should be preferred over the `stdlib` implementation. Examples of such libraries are: OpenBLAS, MKL (TM), Accelerate, and ATLAS. In order to enable their usage, simply ensure that the following pre-processor macros are defined: - `STDLIB_EXTERNAL_BLAS` wraps all BLAS procedures (except for the 128-bit ones) to an external library - `STDLIB_EXTERNAL_LAPACK` wraps all LAPACK procedures (except for the 128-bit ones) to an external library These can be enabled during the build process. For example, with CMake, one can enable these preprocessor directives using `add_compile_definitions(STDLIB_EXTERNAL_BLAS STDLIB_EXTERNAL_LAPACK)`. The same is possible from the `fpm` branch, where the `cpp` preprocessor is enabled by default. For example, the macros can be added to the project's manifest: ```toml # Link against appropriate external BLAS and LAPACK libraries, if necessary [build] link = ["blas", "lapack"] [dependencies] stdlib="*" # Macros are only needed if using an external library [preprocess] [preprocess.cpp] macros = ["STDLIB_EXTERNAL_BLAS", "STDLIB_EXTERNAL_LAPACK"] ``` or directly via compiler flags: `fpm build --flag "-DSTDLIB_EXTERNAL_BLAS -DSTDLIB_EXTERNAL_LAPACK -lblas -llapack"`. ### Syntax All procedures in the `BLAS` and `LAPACK` backends follow the standard interfaces from the [Reference LAPACK](https://www.netlib.org/lapack/). So, the online [Users Guide](https://www.netlib.org/lapack/explore-html/) should be consulted for the full API and descriptions of procedure arguments and their usage. The `stdlib` implementation makes both kind-agnostic and specific procedure interfaces available via modules [stdlib_linalg_blas(module)] and [stdlib_linalg_lapack(module)]. Because all procedures start with a letter [that indicates the base datatype](https://www.netlib.org/lapack/lug/node24.html), the `stdlib` generic interface drops the heading letter and contains all kind-dependent implementations. For example, the generic interface to the `axpy` function looks like: ```fortran !> AXPY: constant times a vector plus a vector. interface axpy module procedure stdlib_saxpy module procedure stdlib_daxpy module procedure stdlib_qaxpy module procedure stdlib_caxpy module procedure stdlib_zaxpy module procedure stdlib_waxpy end interface axpy ``` The generic interface is the endpoint for using an external library. Whenever the latter is used, references to the internal `module procedure`s are replaced with interfaces to the external library, for example: ```fortran !> AXPY: constant times a vector plus a vector. interface axpy pure subroutine caxpy(n,ca,cx,incx,cy,incy) import sp,dp,qp,ilp,lk implicit none complex(sp), intent(in) :: ca,cx(*) integer(ilp), intent(in) :: incx,incy,n complex(sp), intent(inout) :: cy(*) end subroutine caxpy ! [....] module procedure stdlib_qaxpy end interface axpy ``` Note that the 128-bit functions are only provided by `stdlib` and always point to the internal implementation. Because 128-bit precision is identified as [stdlib_kinds(module):qp], initials for 128-bit procedures were labelled as `q` (quadruple-precision reals) and `w` ("wide" or quadruple-precision complex numbers). Extended precision ([stdlib_kinds(module):xdp]) calculations are labelled as `x` (extended-precision reals). and `y` (extended-precision complex numbers). ### Example ```{literalinclude} ../../example/linalg/example_blas_gemv.f90 :language: fortran ``` ```{literalinclude} ../../example/linalg/example_lapack_getrf.f90 :language: fortran ``` ### Licensing The Fortran Standard Library is distributed under the MIT License. `LAPACK` and its contained `BLAS` are a freely-available software package. They are available from [netlib](https://www.netlib.org/lapack/) via anonymous ftp and the World Wide Web. Thus, they can be included in commercial software packages (and have been). The license used for the `BLAS` and `LAPACK` backends is the [modified BSD license](https://www.netlib.org/lapack/LICENSE.txt). The header of the `LICENSE.txt` file has as its licensing requirements: Copyright (c) 1992-2013 The University of Tennessee and The University of Tennessee Research Foundation. All rights reserved. Copyright (c) 2000-2013 The University of California Berkeley. All rights reserved. Copyright (c) 2006-2013 The University of Colorado Denver. All rights reserved. $COPYRIGHT$ Additional copyrights may follow $HEADER$ Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: - Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. - Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer listed in this license in the documentation and/or other materials provided with the distribution. - Neither the name of the copyright holders nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. The copyright holders provide no reassurances that the source code provided does not infringe any patent, copyright, or any other intellectual property rights of third parties. The copyright holders disclaim any liability to any recipient for claims brought against recipient by any third party for infringement of that parties intellectual property rights. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. So the license for the `LICENSE.txt` code is compatible with the use of modified versions of the code in the Fortran Standard Library under the MIT license. Credit for the `BLAS`, `LAPACK` libraries should be given to the [LAPACK authors](https://www.netlib.org/lapack/contributor-list.html). According to the original license, we also changed the name of the routines and commented the changes made to the original. ## `diag` - Create a diagonal array or extract the diagonal elements of an array ### Status Stable ### Class Pure function. ### Description Create a diagonal array or extract the diagonal elements of an array ### Syntax `d = ` `diag` `(a [, k])` ### Arguments `a`: Shall be a rank-1 or or rank-2 array. If `a` is a rank-1 array (i.e. a vector) then `diag` returns a rank-2 array with the elements of `a` on the diagonal. If `a` is a rank-2 array (i.e. a matrix) then `diag` returns a rank-1 array of the diagonal elements. `k` (optional): Shall be a scalar of type `integer` and specifies the diagonal. The default `k = 0` represents the main diagonal, `k > 0` are diagonals above the main diagonal, `k < 0` are diagonals below the main diagonal. ### Return value Returns a diagonal array or a vector with the extracted diagonal elements. ### Example ```{literalinclude} ../../example/linalg/example_diag1.f90 :language: fortran ``` ```{literalinclude} ../../example/linalg/example_diag2.f90 :language: fortran ``` ```{literalinclude} ../../example/linalg/example_diag3.f90 :language: fortran ``` ```{literalinclude} ../../example/linalg/example_diag4.f90 :language: fortran ``` ```{literalinclude} ../../example/linalg/example_diag5.f90 :language: fortran ``` ## `eye` - Construct the identity matrix ### Status Stable ### Class Pure function. ### Description Constructs the identity matrix. ### Syntax `I = ` `eye` `(dim1 [, dim2] [, mold])` ### Arguments - `dim1`: A scalar of type `integer`. This is an `intent(in)` argument and specifies the number of rows. - `dim2`: A scalar of type `integer`. This is an optional `intent(in)` argument specifying the number of columns. If not provided, the matrix is square (`dim1 = dim2`). - `mold`: A scalar of any supported `integer`, `real`, or `complex` type. This is an optional `intent(in)` argument. If provided, the returned identity matrix will have the same type and kind as `mold`. If not provided, the matrix will be of type `real(real64)` by default. ### Return value Returns the identity matrix, with ones on the main diagonal and zeros elsewhere. - By default, the return value is of type `real(real64)`, which is recommended for arithmetic safety. - If the `mold` argument is provided, the return value will match the type and kind of `mold`, allowing for arbitrary `integer`, `real`, or `complex` return types. ### Example ```fortran !> Return default type (real64) A = eye(2,2)/2 !! A == diag([0.5_dp, 0.5_dp]) !> Return 32-bit complex A = eye(2,2, mold=(0.0,0.0))/2 !! A == diag([(0.5,0.5), (0.5,0.5)]) ``` ```{literalinclude} ../../example/linalg/example_eye1.f90 :language: fortran ``` ```{literalinclude} ../../example/linalg/example_eye2.f90 :language: fortran ``` ## `trace` - Trace of a matrix ### Status Stable ### Description Trace of a matrix (rank-2 array) ### Syntax `result = ` `trace` `(A)` ### Arguments `A`: Shall be a rank-2 array. If `A` is not square, then `trace(A)` will return the sum of diagonal values from the square sub-section of `A`. ### Return value Returns the trace of the matrix, i.e. the sum of diagonal elements. ### Example ```{literalinclude} ../../example/linalg/example_trace.f90 :language: fortran ``` ## `outer_product` - Computes the outer product of two vectors ### Status Experimental ### Description Computes the outer product of two vectors ### Syntax `d = ` `outer_product` `(u, v)` ### Arguments `u`: Shall be a rank-1 array `v`: Shall be a rank-1 array ### Return value Returns a rank-2 array equal to `u v^T` (where `u, v` are considered column vectors). The shape of the returned array is `[size(u), size(v)]`. ### Example ```{literalinclude} ../../example/linalg/example_outer_product.f90 :language: fortran ``` ## `kronecker_product` - Computes the Kronecker product of two rank-2 arrays ### Status Experimental ### Description Computes the Kronecker product of two rank-2 arrays ### Syntax `C = ` `kronecker_product` `(A, B)` ### Arguments `A`: Shall be a rank-2 array with dimensions M1, N1 `B`: Shall be a rank-2 array with dimensions M2, N2 ### Return value Returns a rank-2 array equal to `A \otimes B`. The shape of the returned array is `[M1*M2, N1*N2]`. ### Example ```{literalinclude} ../../example/linalg/example_kronecker_product.f90 :language: fortran ``` ## `cross_product` - Computes the cross product of two vectors ### Status Experimental ### Description Computes the cross product of two vectors ### Syntax `c = ` `cross_product` `(a, b)` ### Arguments `a`: Shall be a rank-1 and size-3 array `b`: Shall be a rank-1 and size-3 array ### Return value Returns a rank-1 and size-3 array which is perpendicular to both `a` and `b`. ### Example ```{literalinclude} ../../example/linalg/example_cross_product.f90 :language: fortran ``` ## `is_square` - Checks if a matrix is square ### Status Experimental ### Description Checks if a matrix is square ### Syntax `d = ` `is_square` `(A)` ### Arguments `A`: Shall be a rank-2 array ### Return value Returns a `logical` scalar that is `.true.` if the input matrix is square, and `.false.` otherwise. ### Example ```{literalinclude} ../../example/linalg/example_is_square.f90 :language: fortran ``` ## `is_diagonal` - Checks if a matrix is diagonal ### Status Experimental ### Description Checks if a matrix is diagonal ### Syntax `d = ` `is_diagonal` `(A)` ### Arguments `A`: Shall be a rank-2 array ### Return value Returns a `logical` scalar that is `.true.` if the input matrix is diagonal, and `.false.` otherwise. Note that nonsquare matrices may be diagonal, so long as `a_ij = 0` when `i /= j`. ### Example ```{literalinclude} ../../example/linalg/example_is_diagonal.f90 :language: fortran ``` ## `is_symmetric` - Checks if a matrix is symmetric ### Status Experimental ### Description Checks if a matrix is symmetric ### Syntax `d = ` `is_symmetric` `(A)` ### Arguments `A`: Shall be a rank-2 array ### Return value Returns a `logical` scalar that is `.true.` if the input matrix is symmetric, and `.false.` otherwise. ### Example ```{literalinclude} ../../example/linalg/example_is_symmetric.f90 :language: fortran ``` ## `is_skew_symmetric` - Checks if a matrix is skew-symmetric ### Status Experimental ### Description Checks if a matrix is skew-symmetric ### Syntax `d = ` `is_skew_symmetric` `(A)` ### Arguments `A`: Shall be a rank-2 array ### Return value Returns a `logical` scalar that is `.true.` if the input matrix is skew-symmetric, and `.false.` otherwise. ### Example ```{literalinclude} ../../example/linalg/example_is_skew_symmetric.f90 :language: fortran ``` ## `hermitian` - Compute the Hermitian version of a rank-2 matrix ### Status Experimental ### Description Compute the Hermitian version of a rank-2 matrix. For `complex` matrices, the function returns the conjugate transpose (`conjg(transpose(a))`). For `real` or `integer` matrices, the function returns the transpose (`transpose(a)`). ### Syntax `h = ` `hermitian` `(a)` ### Arguments `a`: Shall be a rank-2 array of type `integer`, `real`, or `complex`. The input matrix `a` is not modified. ### Return value Returns a rank-2 array of the same shape and type as `a`. If `a` is of type `complex`, the Hermitian matrix is computed as `conjg(transpose(a))`. For `real` or `integer` types, it is equivalent to the intrinsic `transpose(a)`. ### Example ```{literalinclude} ../../example/linalg/example_hermitian.f90 :language: fortran ``` ## `is_hermitian` - Checks if a matrix is Hermitian ### Status Experimental ### Description Checks if a matrix is Hermitian ### Syntax `d = ` `is_hermitian` `(A)` ### Arguments `A`: Shall be a rank-2 array ### Return value Returns a `logical` scalar that is `.true.` if the input matrix is Hermitian, and `.false.` otherwise. ### Example ```{literalinclude} ../../example/linalg/example_is_hermitian.f90 :language: fortran ``` ## `is_triangular` - Checks if a matrix is triangular ### Status Experimental ### Description Checks if a matrix is triangular ### Syntax `d = ` `is_triangular` `(A,uplo)` ### Arguments `A`: Shall be a rank-2 array `uplo`: Shall be a single character from `{'u','U','l','L'}` ### Return value Returns a `logical` scalar that is `.true.` if the input matrix is the type of triangular specified by `uplo` (upper or lower), and `.false.` otherwise. Note that the definition of triangular used in this implementation allows nonsquare matrices to be triangular. Specifically, upper triangular matrices satisfy `a_ij = 0` when `j < i`, and lower triangular matrices satisfy `a_ij = 0` when `j > i`. ### Example ```{literalinclude} ../../example/linalg/example_is_triangular.f90 :language: fortran ``` ## `is_hessenberg` - Checks if a matrix is hessenberg ### Status Experimental ### Description Checks if a matrix is Hessenberg ### Syntax `d = ` `is_hessenberg` `(A,uplo)` ### Arguments `A`: Shall be a rank-2 array `uplo`: Shall be a single character from `{'u','U','l','L'}` ### Return value Returns a `logical` scalar that is `.true.` if the input matrix is the type of Hessenberg specified by `uplo` (upper or lower), and `.false.` otherwise. Note that the definition of Hessenberg used in this implementation allows nonsquare matrices to be Hessenberg. Specifically, upper Hessenberg matrices satisfy `a_ij = 0` when `j < i-1`, and lower Hessenberg matrices satisfy `a_ij = 0` when `j > i+1`. ### Example ```{literalinclude} ../../example/linalg/example_is_hessenberg.f90 :language: fortran ``` ## `solve` - Solves a linear matrix equation or a linear system of equations. ### Status Stable ### Description This function computes the solution to a linear matrix equation \( A \cdot x = b \), where \( A \) is a square, full-rank, `real` or `complex` matrix. Result vector or array `x` returns the exact solution to within numerical precision, provided that the matrix is not ill-conditioned. An error is returned if the matrix is rank-deficient or singular to working precision. The solver is based on LAPACK's `*GESV` backends. ### Syntax `Pure` interface: `x = ` `solve` `(a, b)` Expert interface: `x = ` `solve` `(a, b [, overwrite_a], err)` ### Arguments `a`: Shall be a rank-2 `real` or `complex` square array containing the coefficient matrix. It is normally an `intent(in)` argument. If `overwrite_a=.true.`, it is an `intent(inout)` argument and is destroyed by the call. `b`: Shall be a rank-1 or rank-2 array of the same kind as `a`, containing the right-hand-side vector(s). It is an `intent(in)` argument. `overwrite_a` (optional): Shall be an input logical flag. if `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument. `err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. The function is not `pure` if this argument is provided. ### Return value For a full-rank matrix, returns an array value that represents the solution to the linear system of equations. 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`. ### Example ```{literalinclude} ../../example/linalg/example_solve1.f90 :language: fortran ``` ```{literalinclude} ../../example/linalg/example_solve2.f90 :language: fortran ``` ## `solve_lu` - Solves a linear matrix equation or a linear system of equations (subroutine interface). ### Status Stable ### Description This subroutine computes the solution to a linear matrix equation \( A \cdot x = b \), where \( A \) is a square, full-rank, `real` or `complex` matrix. Result vector or array `x` returns the exact solution to within numerical precision, provided that the matrix is not ill-conditioned. An error is returned if the matrix is rank-deficient or singular to working precision. If all optional arrays are provided by the user, no internal allocations take place. The solver is based on LAPACK's `*GESV` backends. ### Syntax Simple (`Pure`) interface: `call ` `solve_lu` `(a, b, x)` Expert (`Pure`) interface: `call ` `solve_lu` `(a, b, x [, pivot, overwrite_a, err])` ### Arguments `a`: Shall be a rank-2 `real` or `complex` square array containing the coefficient matrix. It is normally an `intent(in)` argument. If `overwrite_a=.true.`, it is an `intent(inout)` argument and is destroyed by the call. `b`: Shall be a rank-1 or rank-2 array of the same kind as `a`, containing the right-hand-side vector(s). It is an `intent(in)` argument. `x`: Shall be a rank-1 or rank-2 array of the same kind and size as `b`, that returns the solution(s) to the system. It is an `intent(inout)` argument, and must have the `contiguous` property. `pivot` (optional): Shall be a rank-1 array of the same kind and matrix dimension as `a`, providing storage for the diagonal pivot indices. It is an `intent(inout)` arguments, and returns the diagonal pivot indices. `overwrite_a` (optional): Shall be an input logical flag. if `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument. `err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. ### Return value For a full-rank matrix, returns an array value that represents the solution to the linear system of equations. 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`. ### Example ```{literalinclude} ../../example/linalg/example_solve3.f90 :language: fortran ``` ## `solve_chol` - Solves a linear system using Cholesky factorization (one-shot interface). ### Status Experimental ### Description This subroutine computes the solution to a linear matrix equation \( A \cdot x = b \), where \( A \) is a symmetric (or Hermitian) positive definite matrix. It combines Cholesky factorization and the solve step in a single call. Result vector or array `x` returns the exact solution to within numerical precision, provided that the matrix is positive definite. An error is returned if the matrix is not positive definite or has incompatible sizes with the right-hand-side. Use this routine for one-time solves. For repeated solves with the same matrix but different right-hand sides, use `cholesky` followed by `solve_lower_chol`/`solve_upper_chol` for better performance. The solver is based on LAPACK's `*POSV` backends. ### Syntax Simple (`Pure`) interface: `call ` `solve_chol` `(a, b, x)` Expert (`Pure`) interface: `call ` `solve_chol` `(a, b, x [, lower, overwrite_a, err])` ### Arguments `a`: Shall be a rank-2 `real` or `complex` square array containing the coefficient matrix. It is an `intent(inout)` argument. By default (`overwrite_a=.false.`) its contents are preserved; if `overwrite_a=.true.`, it is used as temporary storage and its contents are destroyed by the call. `b`: Shall be a rank-1 or rank-2 array of the same kind as `a`, containing the right-hand-side vector(s). It is an `intent(in)` argument. `x`: Shall be a rank-1 or rank-2 array of the same kind and size as `b`, that returns the solution(s) to the system. It is an `intent(inout)` argument, and must have the `contiguous` property. `lower` (optional): Shall be an input `logical` flag. If `.true.` (default), the lower triangular Cholesky factorization is computed. If `.false.`, the upper triangular factorization is computed. It is an `intent(in)` argument. `overwrite_a` (optional): Shall be an input `logical` flag. If `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument. `err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. ### Return value For a positive definite matrix, returns an array value that represents the solution to the linear system of equations. Raises `LINALG_ERROR` if the matrix is not positive definite. Raises `LINALG_VALUE_ERROR` if the matrix and rhs vectors have invalid/incompatible sizes. If `err` is not present, exceptions trigger an `error stop`. ### Example ```{literalinclude} ../../example/linalg/example_solve_chol.f90 :language: fortran ``` ## `solve_lower_chol` - Solves a linear system using pre-computed lower Cholesky factor. ### Status Experimental ### Description This subroutine computes the solution to a linear matrix equation \( A \cdot x = b \), where \( A \) is a symmetric (or Hermitian) positive definite matrix that has been **previously factorized** using the Cholesky decomposition (via `cholesky` with `lower=.true.`). Result vector or array `x` returns the exact solution to within numerical precision, provided that the factorization is correct. An error is returned if the matrix and right-hand-side have incompatible sizes. The solver is based on LAPACK's `*POTRS` backends. ### Syntax `call ` `solve_lower_chol` `(l, b, x [, err])` ### Arguments `l`: Shall be a rank-2 `real` or `complex` square array containing the **lower** Cholesky factor `L` (output of `cholesky(..., lower=.true.)`). It is an `intent(in)` argument. `b`: Shall be a rank-1 or rank-2 array of the same kind as `l`, containing the right-hand-side vector(s). It is an `intent(in)` argument. `x`: Shall be a rank-1 or rank-2 array of the same kind and size as `b`, that returns the solution(s) to the system. It is an `intent(inout)` argument, and must have the `contiguous` property. `err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. ### Return value For a correctly factorized matrix, returns an array value that represents the solution to the linear system of equations. Raises `LINALG_VALUE_ERROR` if the matrix and rhs vectors have invalid/incompatible sizes. If `err` is not present, exceptions trigger an `error stop`. ### Example ```{literalinclude} ../../example/linalg/example_solve_lower_chol.f90 :language: fortran ``` ## `solve_upper_chol` - Solves a linear system using pre-computed upper Cholesky factor. ### Status Experimental ### Description This subroutine computes the solution to a linear matrix equation \( A \cdot x = b \), where \( A \) is a symmetric (or Hermitian) positive definite matrix that has been **previously factorized** using the Cholesky decomposition (via `cholesky` with `lower=.false.`). Result vector or array `x` returns the exact solution to within numerical precision, provided that the factorization is correct. An error is returned if the matrix and right-hand-side have incompatible sizes. The solver is based on LAPACK's `*POTRS` backends. ### Syntax `call ` `solve_upper_chol` `(u, b, x [, err])` ### Arguments `u`: Shall be a rank-2 `real` or `complex` square array containing the **upper** Cholesky factor `U` (output of `cholesky(..., lower=.false.)`). It is an `intent(in)` argument. `b`: Shall be a rank-1 or rank-2 array of the same kind as `u`, containing the right-hand-side vector(s). It is an `intent(in)` argument. `x`: Shall be a rank-1 or rank-2 array of the same kind and size as `b`, that returns the solution(s) to the system. It is an `intent(inout)` argument, and must have the `contiguous` property. `err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. ### Return value For a correctly factorized matrix, returns an array value that represents the solution to the linear system of equations. Raises `LINALG_VALUE_ERROR` if the matrix and rhs vectors have invalid/incompatible sizes. If `err` is not present, exceptions trigger an `error stop`. ### Example ```{literalinclude} ../../example/linalg/example_solve_upper_chol.f90 :language: fortran ``` ## `lstsq` - Computes the least squares solution to a linear matrix equation. ### Status Stable ### Description This function computes the least-squares solution to a linear matrix equation \( A \cdot x = b \). Result vector `x` returns the approximate solution that minimizes the 2-norm \( || A \cdot x - b ||_2 \), i.e., it contains the least-squares solution to the problem. Matrix `A` may be full-rank, over-determined, or under-determined. The solver is based on LAPACK's `*GELSD` backends. ### Syntax `x = ` `lstsq` `(a, b, [, cond, overwrite_a, rank, err])` ### Arguments `a`: Shall be a rank-2 `real` or `complex` array containing the coefficient matrix. It is an `intent(inout)` argument. `b`: Shall be a rank-1 or rank-2 array of the same kind as `a`, containing one or more right-hand-side vector(s), each in its leading dimension. It is an `intent(in)` argument. `cond` (optional): Shall be a scalar `real` value cut-off threshold for rank evaluation: `s_i >= cond*maxval(s), i=1:rank`. Shall be a scalar, `intent(in)` argument. `overwrite_a` (optional): Shall be an input `logical` flag. If `.true.`, input matrix `A` will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument. `rank` (optional): Shall be an `integer` scalar value, that contains the rank of input matrix `A`. This is an `intent(out)` argument. `err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. ### Return value Returns an array value of the same kind and rank as `b`, containing the solution(s) to the least squares system. Raises `LINALG_ERROR` if the underlying Singular Value Decomposition process did not converge. Raises `LINALG_VALUE_ERROR` if the matrix and right-hand-side vector have invalid/incompatible sizes. Exceptions trigger an `error stop`. ### Example ```{literalinclude} ../../example/linalg/example_lstsq1.f90 :language: fortran ``` ## `solve_lstsq` - Compute the least squares solution to a linear matrix equation (subroutine interface). ### Status Stable ### Description This subroutine computes the least-squares solution to a linear matrix equation \( A \cdot x = b \). Result vector `x` returns the approximate solution that minimizes the 2-norm \( || A \cdot x - b ||_2 \), i.e., it contains the least-squares solution to the problem. Matrix `A` may be full-rank, over-determined, or under-determined. The solver is based on LAPACK's `*GELSD` backends. ### Syntax `call ` `solve_lstsq` `(a, b, x, [, real_storage, int_storage, [cmpl_storage, ] cond, singvals, overwrite_a, rank, err])` ### Arguments `a`: Shall be a rank-2 `real` or `complex` array containing the coefficient matrix. It is an `intent(inout)` argument. `b`: Shall be a rank-1 or rank-2 array of the same kind as `a`, containing one or more right-hand-side vector(s), each in its leading dimension. It is an `intent(in)` argument. `x`: Shall be an array of same kind and rank as `b`, and leading dimension of at least `n`, containing the solution(s) to the least squares system. It is an `intent(inout)` argument. `real_storage` (optional): Shall be a `real` rank-1 array of the same kind `a`, providing working storage for the solver. It minimum size can be determined with a call to `lstsq_space`. It is an `intent(inout)` argument. `int_storage` (optional): Shall be an `integer` rank-1 array, providing working storage for the solver. It minimum size can be determined with a call to `lstsq_space`. It is an `intent(inout)` argument. `cmpl_storage` (optional): For `complex` systems, it shall be a `complex` rank-1 array, providing working storage for the solver. It minimum size can be determined with a call to `lstsq_space`. It is an `intent(inout)` argument. `cond` (optional): Shall be a scalar `real` value cut-off threshold for rank evaluation: `s_i >= cond*maxval(s), i=1:rank`. Shall be a scalar, `intent(in)` argument. `singvals` (optional): Shall be a `real` rank-1 array of the same kind `a` and size at least `min(m,n)`, returning the list of singular values `s(i)>=cond*maxval(s)` from the internal SVD, in descending order of magnitude. It is an `intent(out)` argument. `overwrite_a` (optional): Shall be an input `logical` flag. If `.true.`, input matrix `A` will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument. `rank` (optional): Shall be an `integer` scalar value, that contains the rank of input matrix `A`. This is an `intent(out)` argument. `err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. ### Return value Returns an array value that represents the solution to the least squares system. Raises `LINALG_ERROR` if the underlying Singular Value Decomposition process did not converge. Raises `LINALG_VALUE_ERROR` if the matrix and right-hand-side vector have invalid/incompatible sizes. Exceptions trigger an `error stop`. ### Example ```{literalinclude} ../../example/linalg/example_lstsq2.f90 :language: fortran ``` ## `lstsq_space` - Compute internal working space requirements for the least squares solver. ### Status Stable ### Description This subroutine computes the internal working space requirements for the least-squares solver, `solve_lstsq` . ### Syntax `call ` `lstsq_space` `(a, b, lrwork, liwork [, lcwork])` ### Arguments `a`: Shall be a rank-2 `real` or `complex` array containing the linear system coefficient matrix. It is an `intent(in)` argument. `b`: Shall be a rank-1 or rank-2 array of the same kind as `a`, containing the system's right-hand-side vector(s). It is an `intent(in)` argument. `lrwork`: Shall be an `integer` scalar, that returns the minimum array size required for the `real` working storage to this system. `liwork`: Shall be an `integer` scalar, that returns the minimum array size required for the `integer` working storage to this system. `lcwork` (`complex` `a`, `b`): For a `complex` system, shall be an `integer` scalar, that returns the minimum array size required for the `complex` working storage to this system. ## `constrained_lstsq` - Compute the solution of the equality-constrained least-squares problem {#constrained-lstsq} ### Status Experimental ### Description This function computes the solution \(x\) of the equality-constrained linear least-squares problem $$ \begin{aligned} \mathrm{minimize} & \quad \| Ax - b \|^2 \\ \mathrm{subject~to} & \quad Cx = d, \end{aligned} $$ where \(A\) is an \( m \times n \) matrix (with \(m \geq n\)) and \(C\) a \( p \times n\) matrix (with \(p \leq n\)). The solver is based on LAPACK's `*GLSE` backends. ### Syntax `x = ` `constrained_lstsq` `(A, b, C, d[, overwrite_matrices, err])` ### Arguments `a`: Shall be a rank-2 `real` or `complex` array used in the definition of the least-squares cost. It is an `intent(inout)` argument. `b`: Shall be a rank-1 array of the same kind as `a` appearing in the definition of the least-squares cost. It is an `intent(inout)` argument. `c`: Shall be a rank-2 `real` or `complex` array of the same kind as `a` defining the linear equality constraints. It is an `intent(inout)` argument. `d`: Shall be a rank-1 array of the same kind as `a` appearing in the definition of the linear equality constraints. `overwrite_matrices` (optional): Shall be an input `logical` flag. If `.true.`, the input matrices and vectors will be overwritten during the computation of the solution. It is an `intent(in)` argument. `err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. ### Return value Returns an array of the same kind as `a` containing the solution of the equality constrained least-squares problem. Raises `LINALG_ERROR` if the underlying constrained least-squares solver did not converge. Raises `LINALG_VALUE_ERROR` if the matrices and vectors have invalid/incompatible dimensions. Exceptions trigger an `error stop`. ### Example ```{literalinclude} ../../example/linalg/example_constrained_lstsq1.f90 :language: fortran ``` ## `solve_constrained_lstsq` - Compute the solution of the equality-constrained least squares problem (subroutine interface) {#solve-constrained-lstsq} ### Status Experimental ### Description This subroutine computes the solution \(x\) of the equality-constrained linear least-squares problem $$ \begin{aligned} \mathrm{minimize} & \quad \| Ax - b \|^2 \\ \mathrm{subject~to} & \quad Cx = d, \end{aligned} $$ where \(A\) is an \( m \times n \) matrix (with \(m \geq n\)) and \(C\) a \( p \times n\) matrix (with \(p \leq n\)). The solver is based on LAPACK's `*GLSE` backends. ### Syntax `call ` `solve_constrained_lstsq` `(a, b, c, d, x [, storage, overwrite_matrices, err])` ### Arguments `a`: Shall be a rank-2 `real` or `complex` array used in the definition of the least-squares cost. It is an `intent(inout)` argument. `b`: Shall be a rank-1 array of the same kind as `a` appearing in the definition of the least-squares cost. It is an `intent(inout)` argument. `c`: Shall be a rank-2 `real` or `complex` array of the same kind as `a` defining the linear equality constraints. It is an `intent(inout)` argument. `d`: Shall be a rank-1 array of the same kind as `a` appearing in the definition of the linear equality constraints. `x`: Shall be a rank-1 array of the same kind as `a`. On exit, it contains the solution of the constrained least-squares problem. It is an `intent(out)` argument. `storage` (optional): Shall be a rank-1 array of the same kind as `a` providing working storage for the solver. Its minimum size can be determined with a call to [stdlib_linalg(module):constrained_lstsq_space(interface)]. It is an `intent(out)` argument. `overwrite_matrices` (optional): Shall be an input `logical` flag. If `.true.`, the input matrices and vectors will be overwritten during the computation of the solution. It is an `intent(in)` argument. `err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. ### Example ```{literalinclude} ../../example/linalg/example_constrained_lstsq2.f90 :language: fortran ``` ## `constrained_lstsq_space` - Compute internal workspace requirements for the constrained least-square solver {#constrained-lstsq-space} ### Status Experimental ### Description This subroutine computes the internal workspace requirements for the constrained least-squares solver, [stdlib_linalg(module):solve_constrained_lstsq(interface)]. ### Syntax call [stdlib_linalg(module):constrained_lstsq_space(interface)]`(a, c, lwork [, err])` ### Arguments `a`: Shall be a rank-2 `real` or `complex` array used in the definition of the least-squares cost. It is an `intent(in)` argument. `c`: Shall be a rank-2 `real` or `complex` array of the same kind as `a` defining the linear equality constraints. It is an `intent(in)` argument. `lwork`: Shall be an `integer` scalar returning the optimal size required for the workspace array to solve the constrained least-squares problem. ## `weighted_lstsq` - Computes the weighted least squares solution to a linear matrix equation {#weighted-lstsq} ### Status Experimental ### Description This function computes the weighted least-squares solution to a linear matrix equation \( A \cdot x = b \) where each observation has a different weight. The solver minimizes the weighted 2-norm \( \| D(Ax - b) \|_2^2 \) where \( D = \mathrm{diag}(\sqrt{w}) \) is a diagonal matrix formed from the square roots of the weight vector. This is equivalent to transforming the problem to ordinary least-squares through row scaling. The solver is based on LAPACK's `*GELSD` backends. ### Syntax `x = ` `weighted_lstsq` `(w, a, b [, cond, overwrite_a, rank, err])` ### Arguments `w`: Shall be a rank-1 `real` array containing the weight vector. For complex `a` and `b`, `w` shall use the same real kind as the components of `a`. All weights must be positive. It is an `intent(in)` argument. `a`: Shall be a rank-2 `real` or `complex` array containing the coefficient matrix. It is an `intent(inout)` argument. `b`: Shall be a rank-1 array of the same kind as `a`, containing the right-hand-side vector. It is an `intent(in)` argument. `cond` (optional): Shall be a scalar `real` value cut-off threshold for rank evaluation: singular values with `s_i <= cond*maxval(s)` are treated as zero, and the rank counts values with `s_i > cond*maxval(s), i=1:rank`. Shall be a scalar, `intent(in)` argument. `overwrite_a` (optional): Shall be an input `logical` flag. If `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument. `rank` (optional): Shall be an `integer` scalar value, that contains the rank of input matrix `a`. This is an `intent(out)` argument. `err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. ### Return value Returns an array value of the same kind as `a`, containing the weighted least-squares solution. Raises `LINALG_VALUE_ERROR` if any weight is non-positive or if the matrix and weight/right-hand-side have incompatible sizes. Raises `LINALG_ERROR` if the underlying Singular Value Decomposition process did not converge. Exceptions trigger an `error stop`. ### Example ```{literalinclude} ../../example/linalg/example_weighted_lstsq.f90 :language: fortran ``` ## `solve_weighted_lstsq` - Compute the weighted least squares solution to a linear matrix equation (subroutine interface). {#solve-weighted-lstsq} ### Status Experimental ### Description This subroutine computes the weighted least-squares solution to a linear matrix equation \( A \cdot x = b \) where each observation has a different weight. The solver minimizes the weighted 2-norm \( \| D(Ax - b) \|_2^2 \) where \( D = \mathrm{diag}(\sqrt{w}) \) is a diagonal matrix formed from the square roots of the weight vector. This is equivalent to transforming the problem to ordinary least-squares through row scaling. The solver is based on LAPACK's `*GELSD` backends. ### Syntax `call ` `solve_weighted_lstsq` `(w, a, b, x [, cond, overwrite_a, rank, err])` ### Arguments `w`: Shall be a rank-1 `real` array containing the weight vector. For complex `a` and `b`, `w` shall use the same real kind as the components of `a`. All weights must be positive. It is an `intent(in)` argument. `a`: Shall be a rank-2 `real` or `complex` array containing the coefficient matrix. It is an `intent(inout)` argument. `b`: Shall be a rank-1 array of the same kind as `a`, containing the right-hand-side vector. It is an `intent(in)` argument. `x`: Shall be a rank-1 array of the same kind as `a`, and size of at least `n`, containing the solution to the weighted least squares system. It is an `intent(inout)` argument. `cond` (optional): Shall be a scalar `real` value cut-off threshold for rank evaluation: singular values with `s_i <= cond*maxval(s)` are treated as zero, and the rank counts values with `s_i > cond*maxval(s), i=1:rank`. Shall be a scalar, `intent(in)` argument. `overwrite_a` (optional): Shall be an input `logical` flag. If `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument. `rank` (optional): Shall be an `integer` scalar value, that contains the rank of input matrix `a`. This is an `intent(out)` argument. `err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. ### Return value Returns an array value of the same kind as `a`, containing the weighted least-squares solution. Raises `LINALG_VALUE_ERROR` if any weight is non-positive or if the matrix and weight/right-hand-side have incompatible sizes. Raises `LINALG_ERROR` if the underlying Singular Value Decomposition process did not converge. Exceptions trigger an `error stop`. ### Example ```{literalinclude} ../../example/linalg/example_weighted_lstsq.f90 :language: fortran ``` ## `det` - Computes the determinant of a square matrix ### Status Stable ### Description This function computes the determinant of a `real` or `complex` square matrix. This interface comes with a `pure` version `det(a)`, and a non-pure version `det(a,overwrite_a,err)` that allows for more expert control. ### Syntax `c = ` `det` `(a [, overwrite_a, err])` ### Arguments `a`: Shall be a rank-2 square array `overwrite_a` (optional): Shall be an input `logical` flag. if `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument. `err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. ### Return value Returns a `real` scalar value of the same kind of `a` that represents the determinant of the matrix. Raises `LINALG_ERROR` if the matrix is singular. Raises `LINALG_VALUE_ERROR` if the matrix is non-square. Exceptions are returned to the `err` argument if provided; an `error stop` is triggered otherwise. ### Example ```{literalinclude} ../../example/linalg/example_determinant.f90 :language: fortran ``` ## `.det.` - Determinant operator of a square matrix ### Status Stable ### Description This operator returns the determinant of a real square matrix. This interface is equivalent to the `pure` version of determinant `det`. ### Syntax `c = ` `operator(.det.)` `a` ### Arguments `a`: Shall be a rank-2 square array of any `real` or `complex` kinds. It is an `intent(in)` argument. ### Return value Returns a real scalar value that represents the determinant of the matrix. Raises `LINALG_ERROR` if the matrix is singular. Raises `LINALG_VALUE_ERROR` if the matrix is non-square. Exceptions trigger an `error stop`. ### Example ```{literalinclude} ../../example/linalg/example_determinant2.f90 :language: fortran ``` ## `qr` - Compute the QR factorization of a matrix ### Status Experimental ### Description This subroutine computes the QR factorization of a `real` or `complex` matrix: \( A = Q R \) where \( Q \) is orthonormal and \( R \) is upper-triangular. Matrix \( A \) has size `[m,n]`, with \( m \ge n \). The results are returned in output matrices \( Q \) and \(R \), that have the same type and kind as \( A \). Given `k = min(m,n)`, one can write \( A = \( Q_1 Q_2 \) \cdot \( \frac{R_1}{0}\) \). Because the lower rows of \( R \) are zeros, a reduced problem \( A = Q_1 R_1 \) may be solved. The size of the input arguments determines what problem is solved: on full matrices (`shape(Q)==[m,m]`, `shape(R)==[m,n]`), the full problem is solved. On reduced matrices (`shape(Q)==[m,k]`, `shape(R)==[k,n]`), the reduced problem is solved. ### Syntax `call ` `qr` `(a, q, r [, pivots] [, overwrite_a] [, storage] [, err])` ### Arguments `a`: Shall be a rank-2 `real` or `complex` array containing the coefficient matrix of size `[m,n]`. It is an `intent(in)` argument, if `overwrite_a=.false.`. Otherwise, it is an `intent(inout)` argument, and is destroyed upon return. `q`: Shall be a rank-2 array of the same kind as `a`, containing the orthonormal matrix `q`. It is an `intent(out)` argument. It should have a shape equal to either `[m,m]` or `[m,k]`, whether the full or the reduced problem is sought for. `r`: Shall be a rank-2 array of the same kind as `a`, containing the upper triangular matrix `r`. It is an `intent(out)` argument. It should have a shape equal to either `[m,n]` or `[k,n]`, whether the full or the reduced problem is sought for. `pivots` (optional): Shall be an `integer` array of size `n`. If provided, QR factorization with column-pivoting is being computed. It is an `intent(out)` argument. `overwrite_a` (optional): Shall be an input `logical` flag (default: `.false.`). If `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation. It is an `intent(in)` argument. `storage` (optional): Shall be a rank-1 array of the same type and kind as `a`, providing working storage for the solver. Its minimum size can be determined with a call to `qr_space`. It is an `intent(out)` argument. `err` (optional): Shall be a `type(linalg_state_type)` value. It is an `intent(out)` argument. ### Return value Returns the QR factorization matrices into the \( Q \) and \( R \) arguments and the optional pivots in `pivots`. Raises `LINALG_VALUE_ERROR` if any of the matrices has invalid or unsuitable size for the full/reduced problem. Raises `LINALG_ERROR` on insufficient user storage space. If the state argument `err` is not present, exceptions trigger an `error stop`. ### Example ```{literalinclude} ../../example/linalg/example_qr.f90 :language: fortran ``` ```{literalinclude} ../../example/linalg/example_pivoting_qr.f90 :language: fortran ``` ## `qr_space` - Compute internal working space requirements for the QR factorization. ### Status Experimental ### Description This subroutine computes the internal working space requirements for the QR factorization, `qr` . ### Syntax `call ` `qr_space` `(a, lwork, [, pivoting] [, err])` ### Arguments `a`: Shall be a rank-2 `real` or `complex` array containing the coefficient matrix. It is an `intent(in)` argument. `lwork`: Shall be an `integer` scalar, that returns the minimum array size required for the working storage in `qr` to factorize `a`. `pivoting` (optional): Shall a `logical` flag (default: `.false.`). If `.true.`, on exit `lwork` is the optimal workspace size for the QR factorization with column pivoting. If `.false.`, `lwork` is the optimal workspace size for the standard QR factorization. It is an `intent(in)` argument. `err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. ### Example ```{literalinclude} ../../example/linalg/example_qr_space.f90 :language: fortran ``` ```{literalinclude} ../../example/linalg/example_pivoting_qr_space.f90 :language: fortran ``` ## `schur` - Compute the Schur decomposition of a matrix ### Status Experimental ### Description This subroutine computes the Schur decomposition of a `real` or `complex` matrix: \( A = Z T Z^H \), where \( Z \) is unitary (orthonormal) and \( T \) is upper-triangular (for complex matrices) or quasi-upper-triangular (for real matrices, with possible \( 2 \times 2 \) blocks on the diagonal). Matrix \( A \) has size `[n,n]`. The results are returned in output matrices \( T \) and \( Z \). Matrix \( T \) is the Schur form, and matrix \( Z \) is the unitary transformation matrix such that \( A = Z T Z^H \). If requested, the eigenvalues of \( T \) can also be returned as a `complex` array of size `[n]`. ### Syntax `call ` `schur` `(a, t [, z,] [, eigvals] [, overwrite_a] [, storage] [, err])` ### Arguments - `a`: Shall be a rank-2 `real` or `complex` array containing the matrix to be decomposed. It is an `intent(inout)` argument if `overwrite_a = .true.`; otherwise, it is an `intent(in)` argument. - `t`: Shall be a rank-2 array of the same kind as `a`, containing the Schur form \( T \) of the matrix. It is an `intent(out)` argument and should have a shape of `[n,n]`. - `z`: Shall be a rank-2 array of the same kind as `a`, containing the unitary matrix \( Z \). It is an `intent(out)` argument and is optional. If provided, it should have the shape `[n,n]`. - `eigvals` (optional): Shall be a rank-1 `complex` or `real` array of the same kind as `a`, containing the eigenvalues of \( A \) (the diagonal elements of \( T \)), or their `real` component only. The array must be of size `[n]`. If not provided, the eigenvalues are not returned. It is an `intent(out)` argument. - `overwrite_a` (optional): Shall be a `logical` flag (default: `.false.`). If `.true.`, the input matrix `a` will be overwritten and destroyed upon return, avoiding internal data allocation. It is an `intent(in)` argument. - `storage` (optional): Shall be a rank-1 array of the same type and kind as `a`, providing working storage for the solver. Its minimum size can be determined with a call to `schur_space`. It is an `intent(inout)` argument. - `err` (optional): Shall be a `type(linalg_state_type)` value. It is an `intent(out)` argument. If not provided, exceptions trigger an `error stop`. ### Return value Returns the Schur decomposition matrices into the \( T \) and \( Z \) arguments. If `eigvals` is provided, it will also return the eigenvalues of the matrix \( A \). Raises `LINALG_VALUE_ERROR` if any of the matrices have invalid or unsuitable size for the decomposition. Raises `LINALG_VALUE_ERROR` if the `real` component is only requested, but the eigenvalues have non-trivial imaginary parts. Raises `LINALG_ERROR` on insufficient user storage space. If the state argument `err` is not present, exceptions trigger an `error stop`. ### Example ```{literalinclude} ../../example/linalg/example_schur_eigvals.f90 :language: fortran ``` --- ## `schur_space` - Compute internal working space requirements for the Schur decomposition ### Status Experimental ### Description This subroutine computes the internal working space requirements for the Schur decomposition, `schur`. ### Syntax `call ` `schur_space` `(a, lwork, [, err])` ### Arguments - `a`: Shall be a rank-2 `real` or `complex` array containing the matrix to be decomposed. It is an `intent(in)` argument. - `lwork`: Shall be an `integer` scalar that returns the minimum array size required for the working storage in `schur` to decompose `a`. It is an `intent(out)` argument. - `err` (optional): Shall be a `type(linalg_state_type)` value. It is an `intent(out)` argument. ### Return value Returns the required working storage size `lwork` for the Schur decomposition. This value can be used to pre-allocate a workspace array in case multiple Schur decompositions of the same matrix size are needed. If pre-allocated working arrays are provided, no internal memory allocations will take place during the decomposition. ## `eig` - Eigenvalues and Eigenvectors of a Square Matrix ### Status Stable ### Description This subroutine computes the solution to the eigenproblem \( A \cdot \bar{v} - \lambda \cdot \bar{v} \), where \( A \) is a square, full-rank, `real` or `complex` matrix, or to the generalized eigenproblem \( A \cdot \bar{v} - \lambda \cdot B \cdot \bar{v} \), where \( B \) is a square matrix with the same type, kind and size as \( A \). Result array `lambda` returns the eigenvalues of \( A \). The user can request eigenvectors to be returned: if provided, on output `left` will contain the left eigenvectors, `right` the right eigenvectors of \( A \). Both `left` and `right` are rank-2 arrays, where eigenvectors are stored as columns. The solver is based on LAPACK's `*GEEV` (standard eigenproblem) and `*GGEV` (generalized eigenproblem) backends. ### Syntax For the standard eigenproblem: `call ` `eig` `(a, lambda [, right] [,left] [,overwrite_a] [,err])` For the generalized eigenproblem: `call ` `eig` `(a, b, lambda [, right] [, left] [, overwrite_a] [, overwrite_b] [, err]) ### Arguments `a` : `real` or `complex` square array containing the coefficient matrix. If `overwrite_a=.false.`, it is an `intent(in)` argument. Otherwise, it is an `intent(inout)` argument and is destroyed by the call. `b`: `real` or `complex` square array containing the second coefficient matrix. If `overwrite_b=.false.`, it is an `intent(in)` argument. Otherwise, it is an `intent(inout)` argument and is destroyed by the call. `lambda`: Shall be a `complex` or `real` rank-1 array of the same kind as `a`, containing the eigenvalues, or their `real` component only. It is an `intent(out)` argument. `right` (optional): Shall be a `complex` rank-2 array of the same size and kind as `a`, containing the right eigenvectors of `a`. It is an `intent(out)` argument. `left` (optional): Shall be a `complex` rank-2 array of the same size and kind as `a`, containing the left eigenvectors of `a`. It is an `intent(out)` argument. `overwrite_a` (optional): Shall be an input logical flag. if `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument. `overwrite_b` (optional): Shall be an input logical flag. If `.true.`, input matrix `b` will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument. `err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. ### Return value Raises `LINALG_ERROR` if the calculation did not converge. Raises `LINALG_VALUE_ERROR` if any matrix or arrays have invalid/incompatible sizes. Raises `LINALG_VALUE_ERROR` if the `real` component is only requested, but the eigenvalues have non-trivial imaginary parts. If `err` is not present, exceptions trigger an `error stop`. ### Example ```{literalinclude} ../../example/linalg/example_eig.f90 :language: fortran ``` ## `eigh` - Eigenvalues and Eigenvectors of a Real symmetric or Complex Hermitian Square Matrix ### Status Stable ### Description This subroutine computes the solution to the eigendecomposition \( A \cdot \bar{v} - \lambda \cdot \bar{v} \), where \( A \) is a square, full-rank, `real` symmetric \( A = A^T \) or `complex` Hermitian \( A = A^H \) matrix. Result array `lambda` returns the `real` eigenvalues of \( A \). The user can request the orthogonal eigenvectors to be returned: on output `vectors` may contain the matrix of eigenvectors, returned as a column. Normally, only the lower triangular part of \( A \) is accessed. On input, `logical` flag `upper_a` allows the user to request what triangular part of the matrix should be used. The solver is based on LAPACK's `*SYEV` and `*HEEV` backends. ### Syntax `call ` `eigh` `(a, lambda [, vectors] [, upper_a] [, overwrite_a] [,err])` ### Arguments `a` : `real` or `complex` square array containing the coefficient matrix. It is an `intent(in)` argument. If `overwrite_a=.true.`, it is an `intent(inout)` argument and is destroyed by the call. `lambda`: Shall be a `complex` rank-1 array of the same precision as `a`, containing the eigenvalues. It is an `intent(out)` argument. `vectors` (optional): Shall be a rank-2 array of the same type, size and kind as `a`, containing the eigenvectors of `a`. It is an `intent(out)` argument. `upper_a` (optional): Shall be an input `logical` flag. If `.true.`, the upper triangular part of `a` will be accessed. Otherwise, the lower triangular part will be accessed. It is an `intent(in)` argument. `overwrite_a` (optional): Shall be an input `logical` flag. If `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument. `err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. ### Return value Raises `LINALG_ERROR` if the calculation did not converge. Raises `LINALG_VALUE_ERROR` if any matrix or arrays have invalid/incompatible sizes. If `err` is not present, exceptions trigger an `error stop`. ### Example ```{literalinclude} ../../example/linalg/example_eigh.f90 :language: fortran ``` ## `eigvals` - Eigenvalues of a Square Matrix ### Status Stable ### Description This function computes the eigenvalues for either a standard or generalized eigenproblem: - **Standard eigenproblem**: \( A \cdot \bar{v} - \lambda \cdot \bar{v} \), where \( A \) is a square, full-rank `real` or `complex` matrix. - **Generalized eigenproblem**: \( A \cdot \bar{v} - \lambda \cdot B \cdot \bar{v} \), where \( B \) is a square matrix with the same type and kind as \( A \). The eigenvalues are stored in the result array `lambda`, which is `complex` (even for real input matrices). The solver uses LAPACK's `*GEEV` and `*GGEV` backends for the standard and generalized problems, respectively. ### Syntax For the standard eigenproblem: `lambda = ` `eigvals` `(a [, err])` For the generalized eigenproblem: `lambda = ` `eigvals` `(a, b [, err])` ### Arguments `a`: Shall be a `real` or `complex` square array containing the coefficient matrix. It is an `intent(in)` argument. `b` (optional): Shall be a `real` or `complex` square array containing the second coefficient matrix for the generalized problem. It is an `intent(in)` argument. `err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. ### Return Value Returns a `complex` rank-1 array containing the eigenvalues of the problem. Raises `LINALG_ERROR` if the calculation did not converge. Raises `LINALG_VALUE_ERROR` if any matrix or arrays have invalid/incompatible sizes. If `err` is not present, exceptions trigger an `error stop`. ### Example ```{literalinclude} ../../example/linalg/example_eigvals.f90 :language: fortran ``` ## `eigvalsh` - Eigenvalues of a Real Symmetric or Complex Hermitian Square Matrix ### Status Stable ### Description This function returns the eigenvalues to matrix \( A \): a where \( A \) is a square, full-rank, `real` symmetric \( A = A^T \) or `complex` Hermitian \( A = A^H \) matrix. The eigenvalues are solutions to the eigenproblem \( A \cdot \bar{v} - \lambda \cdot \bar{v} \). Result array `lambda` is `real`, and returns the eigenvalues of \( A \). The solver is based on LAPACK's `*SYEV` and `*HEEV` backends. ### Syntax `lambda = ` `eigvalsh` `(a, [, upper_a] [,err])` ### Arguments `a` : `real` or `complex` square array containing the coefficient matrix. It is an `intent(in)` argument. `upper_a` (optional): Shall be an input logical flag. If `.true.`, the upper triangular part of `a` will be used accessed. Otherwise, the lower triangular part will be accessed (default). It is an `intent(in)` argument. `err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. ### Return value Returns a `real` array containing the eigenvalues of `a`. Raises `LINALG_ERROR` if the calculation did not converge. Raises `LINALG_VALUE_ERROR` if any matrix or arrays have invalid/incompatible sizes. If `err` is not present, exceptions trigger an `error stop`. ### Example ```{literalinclude} ../../example/linalg/example_eigvalsh.f90 :language: fortran ``` ## `svd` - Compute the singular value decomposition of a rank-2 array (matrix). ### Status Stable ### Description This subroutine computes the singular value decomposition of a `real` or `complex` rank-2 array (matrix) \( A = U \cdot S \cdot \V^T \). The solver is based on LAPACK's `*GESDD` backends. Result vector `s` returns the array of singular values on the diagonal of \( S \). If requested, `u` contains the left singular vectors, as columns of \( U \). If requested, `vt` contains the right singular vectors, as rows of \( V^T \). ### Syntax `call ` `svd` `(a, s, [, u, vt, overwrite_a, full_matrices, err])` ### Class Subroutine ### Arguments `a`: Shall be a rank-2 `real` or `complex` array containing the coefficient matrix of size `[m,n]`. It is an `intent(inout)` argument, but returns unchanged unless `overwrite_a=.true.`. `s`: Shall be a rank-1 `real` array, returning the list of `k = min(m,n)` singular values. It is an `intent(out)` argument. `u` (optional): Shall be a rank-2 array of same kind as `a`, returning the left singular vectors of `a` as columns. Its size should be `[m,m]` unless `full_matrices=.false.`, in which case, it can be `[m,min(m,n)]`. It is an `intent(out)` argument. `vt` (optional): Shall be a rank-2 array of same kind as `a`, returning the right singular vectors of `a` as rows. Its size should be `[n,n]` unless `full_matrices=.false.`, in which case, it can be `[min(m,n),n]`. It is an `intent(out)` argument. `overwrite_a` (optional): Shall be an input `logical` flag. If `.true.`, input matrix `A` will be used as temporary storage and overwritten. This avoids internal data allocation. By default, `overwrite_a=.false.`. It is an `intent(in)` argument. `full_matrices` (optional): Shall be an input `logical` flag. If `.true.` (default), matrices `u` and `vt` shall be full-sized. Otherwise, their secondary dimension can be resized to `min(m,n)`. See `u`, `v` for details. `err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. ### Return values Returns an array `s` that contains the list of singular values of matrix `a`. If requested, returns a rank-2 array `u` that contains the left singular vectors of `a` along its columns. If requested, returns a rank-2 array `vt` that contains the right singular vectors of `a` along its rows. Raises `LINALG_ERROR` if the underlying Singular Value Decomposition process did not converge. Raises `LINALG_VALUE_ERROR` if the matrix or any of the output arrays invalid/incompatible sizes. Exceptions trigger an `error stop`, unless argument `err` is present. ### Example ```{literalinclude} ../../example/linalg/example_svd.f90 :language: fortran :encoding: iso-8859-1 ``` ## `svdvals` - Compute the singular values of a rank-2 array (matrix). ### Status Stable ### Description This subroutine computes the singular values of a `real` or `complex` rank-2 array (matrix) from its singular value decomposition \( A = U \cdot S \cdot \V^T \). The solver is based on LAPACK's `*GESDD` backends. Result vector `s` returns the array of singular values on the diagonal of \( S \). ### Syntax `s = ` `svdvals` `(a [, err])` ### Arguments `a`: Shall be a rank-2 `real` or `complex` array containing the coefficient matrix of size `[m,n]`. It is an `intent(in)` argument. `err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. ### Return values Returns an array `s` that contains the list of singular values of matrix `a`. Raises `LINALG_ERROR` if the underlying Singular Value Decomposition process did not converge. Raises `LINALG_VALUE_ERROR` if the matrix or any of the output arrays invalid/incompatible sizes. Exceptions trigger an `error stop`, unless argument `err` is present. ### Example ```{literalinclude} ../../example/linalg/example_svdvals.f90 :language: fortran ``` ## `cholesky` - Compute the Cholesky factorization of a rank-2 square array (matrix) ### Status Experimental ### Description This subroutine computes the Cholesky factorization of a `real` or `complex` rank-2 square array (matrix), \( A = L \cdot L^T \), or \( A = U^T \cdot U \). \( A \) is symmetric or complex Hermitian, and \( L \), \( U \) are lower- or upper-triangular, respectively. The solver is based on LAPACK's `*POTRF` backends. ### Syntax `call ` `cholesky` `(a, c, lower [, other_zeroed] [, err])` ### Class Subroutine ### Arguments `a`: Shall be a rank-2 square `real` or `complex` array containing the coefficient matrix of size `[n,n]`. It is an `intent(inout)` argument, but returns unchanged if the argument `c` is present. `c` (optional): Shall be a rank-2 square `real` or `complex` of the same size and kind as `a`. It is an `intent(out)` argument, that returns the triangular Cholesky matrix `L` or `U`. `lower`: Shall be an input `logical` flag. If `.true.`, the lower triangular decomposition \( A = L \cdot L^T \) will be performed. If `.false.`, the upper decomposition \( A = U^T \cdot U \) will be performed. `other_zeroed` (optional): Shall be an input `logical` flag. If `.true.`, the unused part of the output matrix will contain zeroes. Otherwise, it will not be accessed. This saves cpu time. By default, `other_zeroed=.true.`. It is an `intent(in)` argument. `err` (optional): Shall be a `type(linalg_state_type)` value. It is an `intent(out)` argument. ### Return values The factorized matrix is returned in-place overwriting `a` if no other arguments are provided. Otherwise, it can be provided as a second argument `c`. In this case, `a` is not overwritten. The `logical` flag `lower` determines whether the lower- or the upper-triangular factorization should be performed. Results are returned on the applicable triangular region of the output matrix, while the unused triangular region is filled by zeroes by default. Optional argument `other_zeroed`, if `.false.` allows the expert user to avoid zeroing the unused part; however, in this case, the unused region of the matrix is not accessed and will usually contain invalid values. Raises `LINALG_ERROR` if the underlying process did not converge. Raises `LINALG_VALUE_ERROR` if the matrix or any of the output arrays invalid/incompatible sizes. Exceptions trigger an `error stop`, unless argument `err` is present. ### Example ```{literalinclude} ../../example/linalg/example_cholesky.f90 :language: fortran ``` ## `chol` - Compute the Cholesky factorization of a rank-2 square array (matrix) ### Status Experimental ### Description This is a `pure` functional interface for the Cholesky factorization of a `real` or `complex` rank-2 square array (matrix) computed as \( A = L \cdot L^T \), or \( A = U^T \cdot U \). \( A \) is symmetric or complex Hermitian, and \( L \), \( U \) are lower- or upper-triangular, respectively. The solver is based on LAPACK's `*POTRF` backends. Result matrix `c` has the same size and kind as `a`, and returns the lower or upper triangular factor. ### Syntax `c = ` `chol` `(a, lower [, other_zeroed])` ### Arguments `a`: Shall be a rank-2 square `real` or `complex` array containing the coefficient matrix of size `[n,n]`. It is an `intent(inout)` argument, but returns unchanged if argument `c` is present. `lower`: Shall be an input `logical` flag. If `.true.`, the lower triangular decomposition \( A = L \cdot L^T \) will be performed. If `.false.`, the upper decomposition \( A = U^T \cdot U \) will be performed. `other_zeroed` (optional): Shall be an input `logical` flag. If `.true.`, the unused part of the output matrix will contain zeroes. Otherwise, it will not be accessed. This saves cpu time. By default, `other_zeroed=.true.`. It is an `intent(in)` argument. ### Return values Returns a rank-2 array `c` of the same size and kind as `a`, that contains the triangular Cholesky matrix `L` or `U`. Raises `LINALG_ERROR` if the underlying process did not converge. Raises `LINALG_VALUE_ERROR` if the matrix or any of the output arrays invalid/incompatible sizes. Exceptions trigger an `error stop`, unless argument `err` is present. ### Example ```{literalinclude} ../../example/linalg/example_chol.f90 :language: fortran ``` ## `.inv.` - Inverse operator of a square matrix ### Status Stable ### Description This operator returns the inverse of a `real` or `complex` square matrix \( A \). The inverse \( A^{-1} \) is defined such that \( A \cdot A^{-1} = A^{-1} \cdot A = I_n \). This interface is equivalent to the function `inv`. ### Syntax `b = ` `operator(.inv.)` `a` ### Arguments `a`: Shall be a rank-2 square array of any `real` or `complex` kinds. It is an `intent(in)` argument. ### Return value Returns a rank-2 square array with the same type, kind and rank as `a`, that contains the inverse of `a`. If an exception occurred on input errors, or singular matrix, `NaN`s will be returned. For fine-grained error control in case of singular matrices prefer the `subroutine` and the `function` interfaces. ### Example ```{literalinclude} ../../example/linalg/example_inverse_operator.f90 :language: fortran ``` ## `invert` - Inversion of a square matrix ### Status Stable ### Description This subroutine inverts a square `real` or `complex` matrix in-place. The inverse \( A^{-1} \) is defined such that \( A \cdot A^{-1} = A^{-1} \cdot A = I_n \). On return, the input matrix `a` is replaced by its inverse. The solver is based on LAPACK's `*GETRF` and `*GETRI` backends. ### Syntax `call ` `invert` `(a, [,inva] [, pivot] [, err])` ### Arguments `a`: Shall be a rank-2, square, `real` or `complex` array containing the coefficient matrix. If `inva` is provided, it is an `intent(in)` argument. If `inva` is not provided, it is an `intent(inout)` argument: on output, it is replaced by the inverse of `a`. `inva` (optional): Shall be a rank-2, square, `real` or `complex` array with the same size, and kind as `a`. On output, it contains the inverse of `a`. `pivot` (optional): Shall be a rank-1 array of the same kind and matrix dimension as `a`, that contains the diagonal pivot indices on return. It is an `intent(inout)` argument. `err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. ### Return value Computes the inverse of the matrix \( A \), \(A^{-1}\, and returns it either in \( A \) or in another matrix. Raises `LINALG_ERROR` if the matrix is singular or has invalid size. Raises `LINALG_VALUE_ERROR` if `inva` and `a` do not have the same size. If `err` is not present, exceptions trigger an `error stop`. ### Example ```{literalinclude} ../../example/linalg/example_inverse_inplace.f90 :language: fortran ``` ```{literalinclude} ../../example/linalg/example_inverse_subroutine.f90 :language: fortran ``` ## `inv` - Inverse of a square matrix. ### Status Stable ### Description This function returns the inverse of a square `real` or `complex` matrix in-place. The inverse, \( A^{-1} \), is defined such that \( A \cdot A^{-1} = A^{-1} \cdot A = I_n \). The solver is based on LAPACK's `*GETRF` and `*GETRI` backends. ### Syntax `b ` `inv` `(a, [, err])` ### Arguments `a`: Shall be a rank-2, square, `real` or `complex` array containing the coefficient matrix. It is an `intent(inout)` argument. `err` (optional): Shall be a `type(linalg_state_type)` value. It is an `intent(out)` argument. ### Return value Returns an array value of the same type, kind and rank as `a`, that contains the inverse matrix \(A^{-1}\). Raises `LINALG_ERROR` if the matrix is singular or has invalid size. If `err` is not present, exceptions trigger an `error stop`. ### Example ```{literalinclude} ../../example/linalg/example_inverse_function.f90 :language: fortran ``` ## `pinv` - Moore-Penrose pseudo-inverse of a matrix ### Status Experimental ### Description This function computes the Moore-Penrose pseudo-inverse of a `real` or `complex` matrix. The pseudo-inverse, \( A^{+} \), generalizes the matrix inverse and satisfies the conditions: - \( A \cdot A^{+} \cdot A = A \) - \( A^{+} \cdot A \cdot A^{+} = A^{+} \) - \( (A \cdot A^{+})^T = A \cdot A^{+} \) - \( (A^{+} \cdot A)^T = A^{+} \cdot A \) The computation is based on singular value decomposition (SVD). Singular values below a relative tolerance threshold \( \text{rtol} \cdot \sigma_{\max} \), where \( \sigma_{\max} \) is the largest singular value, are treated as zero. ### Syntax `b =` `pinv` `(a, [, rtol, err])` ### Arguments `a`: Shall be a rank-2, `real` or `complex` array of shape `[m, n]` containing the coefficient matrix. It is an `intent(in)` argument. `rtol` (optional): Shall be a scalar `real` value specifying the relative tolerance for singular value cutoff. If `rtol` is not provided, the default relative tolerance is \( \text{rtol} = \text{max}(m, n) \cdot \epsilon \), where \( \epsilon \) is the machine precision for the element type of `a`. It is an `intent(in)` argument. `err` (optional): Shall be a `type(linalg_state_type)` value. It is an `intent(out)` argument. ### Return value Returns an array value of the same type, kind, and rank as `a` with shape `[n, m]`, that contains the pseudo-inverse matrix \( A^{+} \). Raises `LINALG_ERROR` if the underlying SVD did not converge. Raises `LINALG_VALUE_ERROR` if `a` has invalid size. If `err` is not present, exceptions trigger an `error stop`. ### Example ```{literalinclude} ../../example/linalg/example_pseudoinverse.f90 :language: fortran ``` ## `pseudoinvert` - Moore-Penrose pseudo-inverse of a matrix ### Status Experimental ### Description This subroutine computes the Moore-Penrose pseudo-inverse of a `real` or `complex` matrix. The pseudo-inverse \( A^{+} \) is a generalization of the matrix inverse and satisfies the following properties: - \( A \cdot A^{+} \cdot A = A \) - \( A^{+} \cdot A \cdot A^{+} = A^{+} \) - \( (A \cdot A^{+})^T = A \cdot A^{+} \) - \( (A^{+} \cdot A)^T = A^{+} \cdot A \) The computation is based on singular value decomposition (SVD). Singular values below a relative tolerance threshold \( \text{rtol} \cdot \sigma_{\max} \), where \( \sigma_{\max} \) is the largest singular value, are treated as zero. On return, matrix `pinva` `[n, m]` will store the pseudo-inverse of `a` `[m, n]`. ### Syntax `call ` `pseudoinvert` `(a, pinva [, rtol] [, err])` ### Arguments `a`: Shall be a rank-2, `real` or `complex` array containing the coefficient matrix. It is an `intent(in)` argument. `pinva`: Shall be a rank-2 array of the same kind as `a`, and size equal to that of `transpose(a)`. On output, it contains the Moore-Penrose pseudo-inverse of `a`. `rtol` (optional): Shall be a scalar `real` value specifying the relative tolerance for singular value cutoff. If not provided, the default threshold is \( \text{max}(m, n) \cdot \epsilon \), where \( \epsilon \) is the machine precision for the element type of `a`. `err` (optional): Shall be a `type(linalg_state_type)` value. It is an `intent(out)` argument. ### Return value Computes the Moore-Penrose pseudo-inverse of the matrix \( A \), \( A^{+} \), and returns it in matrix `pinva`. Raises `LINALG_ERROR` if the underlying SVD did not converge. Raises `LINALG_VALUE_ERROR` if `pinva` and `a` have degenerate or incompatible sizes. If `err` is not present, exceptions trigger an `error stop`. ### Example ```{literalinclude} ../../example/linalg/example_pseudoinverse.f90 :language: fortran ``` ## `.pinv.` - Moore-Penrose Pseudo-Inverse operator ### Status Experimental ### Description This operator returns the Moore-Penrose pseudo-inverse of a `real` or `complex` matrix \( A \). The pseudo-inverse \( A^{+} \) is computed using Singular Value Decomposition (SVD), and singular values below a given threshold are treated as zero. This interface is equivalent to the function `pinv`. ### Syntax `b = ` `operator(.pinv.)` `a` ### Arguments `a`: Shall be a rank-2 array of any `real` or `complex` kinds, with arbitrary dimensions \( m \times n \). It is an `intent(in)` argument. ### Return value Returns a rank-2 array with the same type, kind, and rank as `a`, that contains the Moore-Penrose pseudo-inverse of `a`. If an exception occurs, or if the input matrix is degenerate (e.g., rank-deficient), the returned matrix will contain `NaN`s. For more detailed error handling, it is recommended to use the `subroutine` or `function` interfaces. ### Example ```{literalinclude} ../../example/linalg/example_pseudoinverse.f90 :language: fortran ``` ## `get_norm` - Computes the vector norm of a generic-rank array. ### Status Experimental ### Description This `pure subroutine` interface computes one of several vector norms of `real` or `complex` array \( A \), depending on the value of the `order` input argument. \( A \) may be an array of any rank. Result `nrm` returns a `real`, scalar norm value for the whole array; if `dim` is specified, `nrm` is a rank n-1 array with the same shape as \(A \) and dimension `dim` dropped, containing all norms evaluated along `dim`. ### Syntax `call ` `get_norm` `(a, nrm, order, [, dim, err])` ### Arguments `a`: Shall be a rank-n `real` or `complex` array containing the data. It is an `intent(in)` argument. `nrm`: if `dim` is absent, shall be a scalar with the norm evaluated over all the elements of the array. Otherwise, an array of rank `n-1`, and a shape similar to that of `a` with dimension `dim` dropped. `order`: Shall be an `integer` value or a `character` flag that specifies the norm type, as follows. It is an `intent(in)` argument. | Integer input | Character Input | Norm type | |------------------|------------------|---------------------------------------------------------| | `-huge(0)` | `'-inf', '-Inf'` | Minimum absolute value \( \min_i{ \left|a_i\right| } \) | | `1` | `'1'` | 1-norm \( \sum_i{ \left|a_i\right| } \) | | `2` | `'2'` | Euclidean norm \( \sqrt{\sum_i{ a_i^2 }} \) | | `>=3` | `'3','4',...` | p-norm \( \left( \sum_i{ \left|a_i\right|^p }\right) ^{1/p} \) | | `huge(0)` | `'inf', 'Inf'` | Maximum absolute value \( \max_i{ \left|a_i\right| } \) | `dim` (optional): Shall be a scalar `integer` value with a value in the range from `1` to `n`, where `n` is the rank of the array. It is an `intent(in)` argument. `err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. If `err` is not present, the function is `pure`. ### Return value By default, the return value `nrm` is a scalar, and contains the norm as evaluated over all elements of the generic-rank array \( A \). If the optional `dim` argument is present, `nrm` is a rank `n-1` array with the same shape as \( A \) except for dimension `dim`, that is collapsed. Each element of `nrm` contains the 1D norm of the elements of \( A \), evaluated along dimension `dim` only. Raises `LINALG_ERROR` if the requested norm type is invalid. Raises `LINALG_VALUE_ERROR` if any of the arguments has an invalid size. If `err` is not present, exceptions trigger an `error stop`. ### Example ```{literalinclude} ../../example/linalg/example_get_norm.f90 :language: fortran ``` ## `norm` - Computes the vector norm of a generic-rank array. ### Status Experimental ### Description This function computes one of several vector norms of `real` or `complex` array \( A \), depending on the value of the `order` input argument. \( A \) may be an array of any rank. ### Syntax `x = ` `norm` `(a, order, [, dim, err])` ### Arguments `a`: Shall be a rank-n `real` or `complex` array containing the data. It is an `intent(in)` argument. `order`: Shall be an `integer` value or a `character` flag that specifies the norm type, as follows. It is an `intent(in)` argument. | Integer input | Character Input | Norm type | |------------------|------------------|---------------------------------------------------------| | `-huge(0)` | `'-inf', '-Inf'` | Minimum absolute value \( \min_i{ \left|a_i\right| } \) | | `1` | `'1'` | 1-norm \( \sum_i{ \left|a_i\right| } \) | | `2` | `'2'` | Euclidean norm \( \sqrt{\sum_i{ a_i^2 }} \) | | `>=3` | `'3','4',...` | p-norm \( \left( \sum_i{ \left|a_i\right|^p }\right) ^{1/p} \) | | `huge(0)` | `'inf', 'Inf'` | Maximum absolute value \( \max_i{ \left|a_i\right| } \) | `dim` (optional): Shall be a scalar `integer` value with a value in the range from `1` to `n`, where `n` is the rank of the array. It is an `intent(in)` argument. `err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. If `err` is not present, the function is `pure`. ### Return value By default, the return value `x` is a scalar, and contains the norm as evaluated over all elements of the generic-rank array \( A \). If the optional `dim` argument is present, `x` is a rank `n-1` array with the same shape as \( A \) except for dimension `dim`, that is dropped. Each element of `x` contains the 1D norm of the elements of \( A \), evaluated along dimension `dim` only. Raises `LINALG_ERROR` if the requested norm type is invalid. Raises `LINALG_VALUE_ERROR` if any of the arguments has an invalid size. If `err` is not present, exceptions trigger an `error stop`. ### Example ```{literalinclude} ../../example/linalg/example_norm.f90 :language: fortran ``` ## `mnorm` - Computes the matrix norm of a generic-rank array. ### Status Experimental ### Description This function computes one of several matrix norms of `real` or `complex` array \( A \), depending on the value of the `order` input argument. \( A \) must be an array of rank 2 or higher. For arrays of rank > 2, matrix norms are computed over specified dimensions. ### Syntax `x = ` `mnorm` `(a [, order, dim, err])` ### Arguments `a`: Shall be a rank-n `real` or `complex` array containing the data, where n >= 2. It is an `intent(in)` argument. `order` (optional): Shall be an `integer` value or a `character` flag that specifies the norm type, as follows. It is an `intent(in)` argument. | Integer input | Character Input | Norm type | |------------------|---------------------------------|-----------------------------------------------------------------------------| | `1` | `'1'` | 1-norm (maximum column sum) \( \max_j \sum_i{ \left|a_{i,j}\right| } \) | | `2` | `'2'` | 2-norm (largest singular value) | | (not prov.) | `'Euclidean','Frobenius','Fro'` | Frobenius norm \( \sqrt{\sum_{i,j}{ \left|a_{i,j}\right|^2 }} \) | | `huge(0)` | `'inf', 'Inf', 'INF'` | Infinity norm (maximum row sum) \( \max_i \sum_j{ \left|a_{i,j}\right| } \) | `dim` (optional): For arrays of rank > 2, shall be an integer array of size 2 specifying the dimensions over which to compute the matrix norm. Default value is `[1,2]`. It is an `intent(in)` argument. `err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. ### Return value For rank-2 input arrays, the return value `x` is a scalar containing the matrix norm. For arrays of rank > 2, if the optional `dim` argument is present, `x` is a rank `n-2` array with the same shape as \( A \) except for dimensions `dim(1)` and `dim(2)`, which are dropped. Each element of `x` contains the matrix norm of the corresponding submatrix of \( A \), evaluated over the specified dimensions only, with the given order. If an invalid norm type is provided, defaults to 1-norm and raises `LINALG_ERROR`. Raises `LINALG_VALUE_ERROR` if any of the arguments has an invalid size. If `err` is not present, exceptions trigger an `error stop`. ### Example ```{literalinclude} ../../example/linalg/example_mnorm.f90 :language: fortran ``` ## `expm` - Computes the matrix exponential {#expm} ### Status Experimental ### Description Given a matrix \(A\), this function computes its matrix exponential \(E = \exp(A)\) using a Pade approximation. ### Syntax `E = ` `expm` `(a [, order])` ### Arguments `a`: Shall be a rank-2 `real` or `complex` array containing the data. It is an `intent(in)` argument. `order` (optional): Shall be a non-negative `integer` value specifying the order of the Pade approximation. By default `order=10`. It is an `intent(in)` argument. ### Return value The returned array `E` contains the Pade approximation of \(\exp(A)\). If `A` is non-square or `order` is negative, it raises a `LINALG_VALUE_ERROR`. ### Example ```{literalinclude} ../../example/linalg/example_expm.f90 :language: fortran ``` ## `matrix_exp` - Computes the matrix exponential {#matrix_exp} ### Status Experimental ### Description Given a matrix \(A\), this function computes its matrix exponential \(E = \exp(A)\) using a Pade approximation. ### Syntax `call ` `matrix_exp` `(a [, e, order, err])` ### Arguments `a`: Shall be a rank-2 `real` or `complex` array containing the data. If `e` is not passed, it is an `intent(inout)` argument and is overwritten on exit by the matrix exponential. If `e` is passed, it is an `intent(in)` argument and is left unchanged. `e` (optional): Shall be a rank-2 `real` or `complex` array with the same dimensions as `a`. It is an `intent(out)` argument. On exit, it contains the matrix exponential of `a`. `order` (optional): Shall be a non-negative `integer` value specifying the order of the Pade approximation. By default `order=10`. It is an `intent(in)` argument. `err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. ### Return value The returned array `A` (in-place) or `E` (out-of-place) contains the Pade approximation of \(\exp(A)\). If `A` is non-square or `order` is negative, it raises a `LINALG_VALUE_ERROR`. If `err` is not present, exceptions trigger an `error stop`. ### Example ```{literalinclude} ../../example/linalg/example_matrix_exp.f90 :language: fortran ```