stdlib_linalg.fypp

#:include "common.fypp"
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
#:set RCI_KINDS_TYPES = RC_KINDS_TYPES + INT_KINDS_TYPES
#:set RHS_SUFFIX = ["one","many"]
#:set RHS_SYMBOL = [ranksuffix(r) for r in [1,2]]
#:set RHS_EMPTY  = [emptyranksuffix(r) for r in [1,2]]
#:set ALL_RHS    = list(zip(RHS_SYMBOL,RHS_SUFFIX,RHS_EMPTY))
#:set EIG_PROBLEM  = ["standard", "generalized"]
#:set EIG_FUNCTION = ["geev","ggev"]
#:set EIG_PROBLEM_LIST = list(zip(EIG_PROBLEM, EIG_FUNCTION))
module stdlib_linalg
  !!Provides a support for various linear algebra procedures
  !! ([Specification](../page/specs/stdlib_linalg.html))
  use stdlib_kinds, only: xdp, int8, int16, int32, int64
  use stdlib_linalg_constants, only: sp, dp, qp, lk, ilp
  use stdlib_error, only: error_stop
  use stdlib_optval, only: optval
  use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling
  implicit none
  private

  public :: chol
  public :: cholesky
  public :: det
  public :: operator(.det.)
  public :: diag
  public :: eig
  public :: eigh
  public :: eigvals
  public :: eigvalsh
  public :: expm
  public :: eye
  public :: inv
  public :: invert
  public :: operator(.inv.)
  public :: pinv
  public :: pseudoinvert
  public :: operator(.pinv.)
  public :: lstsq
  public :: lstsq_space
  public :: norm
  public :: mnorm
  public :: get_norm
  public :: solve
  public :: solve_lu  
  public :: solve_lstsq
  public :: trace
  public :: svd
  public :: svdvals
  public :: outer_product
  public :: kronecker_product
  public :: cross_product
  public :: qr
  public :: qr_space
  public :: schur
  public :: schur_space
  public :: is_square
  public :: is_diagonal
  public :: is_symmetric
  public :: is_skew_symmetric
  public :: hermitian
  public :: is_hermitian
  public :: is_triangular
  public :: is_hessenberg
  
  ! Export linalg error handling
  public :: linalg_state_type, linalg_error_handling

  interface chol
    !! version: experimental 
    !!
    !! Computes the Cholesky factorization \( A = L \cdot L^T \), or \( A = U^T \cdot U \). 
    !! ([Specification](../page/specs/stdlib_linalg.html#chol-compute-the-cholesky-factorization-of-a-rank-2-square-array-matrix))
    !! 
    !!### Summary 
    !! Pure function interface for computing the Cholesky triangular factors. 
    !!
    !!### Description
    !! 
    !! This interface provides methods for computing the lower- or upper- triangular matrix from the 
    !! Cholesky factorization of a `real` symmetric or `complex` Hermitian matrix.
    !! Supported data types include `real` and `complex`.    
    !! 
    !!@note The solution is based on LAPACK's `*POTRF` methods.
    !!     
     #:for rk,rt,ri in RC_KINDS_TYPES
     pure module function stdlib_linalg_${ri}$_cholesky_fun(a,lower,other_zeroed) result(c)
         !> Input matrix a[m,n]
         ${rt}$, intent(in) :: a(:,:)
         !> [optional] is the lower or upper triangular factor required? Default = lower
         logical(lk), optional, intent(in) :: lower
         !> [optional] should the unused half of the return matrix be zeroed out? Default: yes
         logical(lk), optional, intent(in) :: other_zeroed
         !> Output matrix with Cholesky factors c[n,n]
         ${rt}$ :: c(size(a,1),size(a,2))                  
     end function stdlib_linalg_${ri}$_cholesky_fun
     #:endfor    
  end interface chol  

  interface cholesky
    !! version: experimental 
    !!
    !! Computes the Cholesky factorization \( A = L \cdot L^T \), or \( A = U^T \cdot U \). 
    !! ([Specification](../page/specs/stdlib_linalg.html#cholesky-compute-the-cholesky-factorization-of-a-rank-2-square-array-matrix))
    !! 
    !!### Summary 
    !! Pure subroutine interface for computing the Cholesky triangular factors. 
    !!
    !!### Description
    !! 
    !! This interface provides methods for computing the lower- or upper- triangular matrix from the 
    !! Cholesky factorization of a `real` symmetric or `complex` Hermitian matrix.
    !! Supported data types include `real` and `complex`.    
    !! The factorization is computed in-place if only one matrix argument is present; or returned into 
    !! a second matrix argument, if present. The `lower` `logical` flag allows to select between upper or 
    !! lower factorization; the `other_zeroed` optional `logical` flag allows to choose whether the unused
    !! part of the triangular matrix should be filled with zeroes.
    !! 
    !!@note The solution is based on LAPACK's `*POTRF` methods.
    !!         
     #:for rk,rt,ri in RC_KINDS_TYPES
     pure module subroutine stdlib_linalg_${ri}$_cholesky_inplace(a,lower,other_zeroed,err) 
         !> Input matrix a[m,n]
         ${rt}$, intent(inout), target :: a(:,:)
         !> [optional] is the lower or upper triangular factor required? Default = lower
         logical(lk), optional, intent(in) :: lower
         !> [optional] should the unused half of the return matrix be zeroed out? Default: yes
         logical(lk), optional, intent(in) :: other_zeroed
         !> [optional] state return flag. On error if not requested, the code will stop
         type(linalg_state_type), optional, intent(out) :: err
     end subroutine stdlib_linalg_${ri}$_cholesky_inplace
     
     pure module subroutine stdlib_linalg_${ri}$_cholesky(a,c,lower,other_zeroed,err) 
         !> Input matrix a[n,n]
         ${rt}$, intent(in) :: a(:,:)
         !> Output matrix with Cholesky factors c[n,n]
         ${rt}$, intent(out) :: c(:,:)         
         !> [optional] is the lower or upper triangular factor required? Default = lower
         logical(lk), optional, intent(in) :: lower
         !> [optional] should the unused half of the return matrix be zeroed out? Default: yes
         logical(lk), optional, intent(in) :: other_zeroed
         !> [optional] state return flag. On error if not requested, the code will stop
         type(linalg_state_type), optional, intent(out) :: err     
     end subroutine stdlib_linalg_${ri}$_cholesky
     #:endfor
  end interface cholesky
     
  interface diag
    !! version: experimental
    !!
    !! Creates a diagonal array or extract the diagonal elements of an array
    !! ([Specification](../page/specs/stdlib_linalg.html#
    !! diag-create-a-diagonal-array-or-extract-the-diagonal-elements-of-an-array))
      !
      ! Vector to matrix
      !
    #:for k1, t1 in RCI_KINDS_TYPES
      module function diag_${t1[0]}$${k1}$(v) result(res)
        ${t1}$, intent(in) :: v(:)
        ${t1}$ :: res(size(v),size(v))
      end function diag_${t1[0]}$${k1}$
    #:endfor
    #:for k1, t1 in RCI_KINDS_TYPES
      module function diag_${t1[0]}$${k1}$_k(v,k) result(res)
        ${t1}$, intent(in) :: v(:)
        integer, intent(in) :: k
        ${t1}$ :: res(size(v)+abs(k),size(v)+abs(k))
      end function diag_${t1[0]}$${k1}$_k
    #:endfor

      !
      ! Matrix to vector
      !
    #:for k1, t1 in RCI_KINDS_TYPES
      module function diag_${t1[0]}$${k1}$_mat(A) result(res)
        ${t1}$, intent(in) :: A(:,:)
        ${t1}$ :: res(minval(shape(A)))
      end function diag_${t1[0]}$${k1}$_mat
    #:endfor
    #:for k1, t1 in RCI_KINDS_TYPES
      module function diag_${t1[0]}$${k1}$_mat_k(A,k) result(res)
        ${t1}$, intent(in) :: A(:,:)
        integer, intent(in) :: k
        ${t1}$ :: res(minval(shape(A))-abs(k))
      end function diag_${t1[0]}$${k1}$_mat_k
    #:endfor
  end interface


  ! Matrix trace
  interface trace
    !! version: experimental
    !!
    !! Computes the trace of a matrix
    !! ([Specification](../page/specs/stdlib_linalg.html#
    !! trace-trace-of-a-matrix))
    #:for k1, t1 in RCI_KINDS_TYPES
      module procedure trace_${t1[0]}$${k1}$
    #:endfor
  end interface

  ! Identity matrix 
  interface eye
    !! version: experimental
    !!
    !! Constructs the identity matrix
    !! ([Specification](../page/specs/stdlib_linalg.html#eye-construct-the-identity-matrix))    
    #:for k1, t1 in RCI_KINDS_TYPES
      module procedure eye_${t1[0]}$${k1}$
    #:endfor
  end interface eye

  ! Outer product (of two vectors)
  interface outer_product
    !! version: experimental
    !!
    !! Computes the outer product of two vectors, returning a rank-2 array
    !! ([Specification](../page/specs/stdlib_linalg.html#
    !! outer_product-computes-the-outer-product-of-two-vectors))
    #:for k1, t1 in RCI_KINDS_TYPES
      pure module function outer_product_${t1[0]}$${k1}$(u, v) result(res)
        ${t1}$, intent(in) :: u(:), v(:)
        ${t1}$ :: res(size(u),size(v))
      end function outer_product_${t1[0]}$${k1}$
    #:endfor
  end interface outer_product

  interface kronecker_product
    !! version: experimental
    !!
    !! Computes the Kronecker product of two arrays of size M1xN1, and of M2xN2, returning an (M1*M2)x(N1*N2) array
    !! ([Specification](../page/specs/stdlib_linalg.html#
    !! kronecker_product-computes-the-kronecker-product-of-two-matrices))
    #:for k1, t1 in RCI_KINDS_TYPES
      pure module function kronecker_product_${t1[0]}$${k1}$(A, B) result(C)
        ${t1}$, intent(in) :: A(:,:), B(:,:)
        ${t1}$ :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
      end function kronecker_product_${t1[0]}$${k1}$
    #:endfor
  end interface kronecker_product


  ! Cross product (of two vectors)
  interface cross_product
    !! version: experimental
    !!
    !! Computes the cross product of two vectors, returning a rank-1 and size-3 array
    !! ([Specification](../page/specs/stdlib_linalg.html#cross_product-computes-the-cross-product-of-two-3-d-vectors))
    #:for k1, t1 in RCI_KINDS_TYPES
      pure module function cross_product_${t1[0]}$${k1}$(a, b) result(res)
        ${t1}$, intent(in) :: a(3), b(3)
        ${t1}$ :: res(3)
      end function cross_product_${t1[0]}$${k1}$
    #:endfor
  end interface cross_product


  ! Check for squareness
  interface is_square
    !! version: experimental
    !!
    !! Checks if a matrix (rank-2 array) is square
    !! ([Specification](../page/specs/stdlib_linalg.html#
    !! is_square-checks-if-a-matrix-is-square))
    #:for k1, t1 in RCI_KINDS_TYPES
      module procedure is_square_${t1[0]}$${k1}$
    #:endfor
  end interface is_square


  ! Check for diagonality
  interface is_diagonal
    !! version: experimental
    !!
    !! Checks if a matrix (rank-2 array) is diagonal
    !! ([Specification](../page/specs/stdlib_linalg.html#
    !! is_diagonal-checks-if-a-matrix-is-diagonal))
    #:for k1, t1 in RCI_KINDS_TYPES
      module procedure is_diagonal_${t1[0]}$${k1}$
    #:endfor
  end interface is_diagonal


  ! Check for symmetry
  interface is_symmetric
    !! version: experimental
    !!
    !! Checks if a matrix (rank-2 array) is symmetric
    !! ([Specification](../page/specs/stdlib_linalg.html#
    !! is_symmetric-checks-if-a-matrix-is-symmetric))
    #:for k1, t1 in RCI_KINDS_TYPES
      module procedure is_symmetric_${t1[0]}$${k1}$
    #:endfor
  end interface is_symmetric


  ! Check for skew-symmetry
  interface is_skew_symmetric
    !! version: experimental
    !!
    !! Checks if a matrix (rank-2 array) is skew-symmetric
    !! ([Specification](../page/specs/stdlib_linalg.html#
    !! is_skew_symmetric-checks-if-a-matrix-is-skew-symmetric))
    #:for k1, t1 in RCI_KINDS_TYPES
      module procedure is_skew_symmetric_${t1[0]}$${k1}$
    #:endfor
  end interface is_skew_symmetric


  ! Check for Hermiticity
  interface is_hermitian
    !! version: experimental
    !!
    !! Checks if a matrix (rank-2 array) is Hermitian
    !! ([Specification](../page/specs/stdlib_linalg.html#
    !! is_hermitian-checks-if-a-matrix-is-hermitian))
    #:for k1, t1 in RCI_KINDS_TYPES
      module procedure is_hermitian_${t1[0]}$${k1}$
    #:endfor
  end interface is_hermitian

  interface hermitian
    !! version: experimental
    !!
    !! Computes the Hermitian version of a rank-2 matrix.
    !! For complex matrices, this returns `conjg(transpose(a))`.
    !! For real or integer matrices, this returns `transpose(a)`.
    !!
    !! Usage:
    !! ```
    !! A  = reshape([(1, 2), (3, 4), (5, 6), (7, 8)], [2, 2])
    !! AH = hermitian(A)
    !! ```
    !!
    !! [Specification](../page/specs/stdlib_linalg.html#hermitian-compute-the-hermitian-version-of-a-rank-2-matrix)
    !!

    #:for k1, t1 in RCI_KINDS_TYPES
    pure module function hermitian_${t1[0]}$${k1}$(a) result(ah)
        ${t1}$, intent(in) :: a(:,:)
        ${t1}$ :: ah(size(a, 2), size(a, 1))
    end function hermitian_${t1[0]}$${k1}$
    #:endfor

  end interface hermitian


  ! Check for triangularity
  interface is_triangular
    !! version: experimental
    !!
    !! Checks if a matrix (rank-2 array) is triangular
    !! ([Specification](../page/specs/stdlib_linalg.html#
    !! is_triangular-checks-if-a-matrix-is-triangular))
    #:for k1, t1 in RCI_KINDS_TYPES
      module procedure is_triangular_${t1[0]}$${k1}$
    #:endfor
  end interface is_triangular
  

  ! Check for matrix being Hessenberg
  interface is_hessenberg
    !! version: experimental
    !!
    !! Checks if a matrix (rank-2 array) is Hessenberg
    !! ([Specification](../page/specs/stdlib_linalg.html#
    !! is_hessenberg-checks-if-a-matrix-is-hessenberg))
    #:for k1, t1 in RCI_KINDS_TYPES
      module procedure is_Hessenberg_${t1[0]}$${k1}$
    #:endfor
  end interface is_hessenberg

  ! Solve linear system system Ax=b.
  interface solve
     !! version: experimental 
     !!
     !! Solves the linear system \( A \cdot x = b \) for the unknown vector \( x \) from a square matrix \( A \). 
     !! ([Specification](../page/specs/stdlib_linalg.html#solve-solves-a-linear-matrix-equation-or-a-linear-system-of-equations))
     !!
     !!### Summary 
     !! Interface for solving a linear system arising from a general matrix.
     !!
     !!### Description
     !! 
     !! This interface provides methods for computing the solution of a linear matrix system.
     !! Supported data types include `real` and `complex`. No assumption is made on the matrix 
     !! structure. 
     !! The function can solve simultaneously either one (from a 1-d right-hand-side vector `b(:)`) 
     !! or several (from a 2-d right-hand-side vector `b(:,:)`) systems.
     !! 
     !!@note The solution is based on LAPACK's generic LU decomposition based solvers `*GESV`.
     !!    
     #:for nd,ndsuf,nde in ALL_RHS
     #:for rk,rt,ri in RC_KINDS_TYPES
     module function stdlib_linalg_${ri}$_solve_${ndsuf}$(a,b,overwrite_a,err) result(x)
         !> Input matrix a[n,n]
         ${rt}$, intent(inout), target :: a(:,:)
         !> Right hand side vector or array, b[n] or b[n,nrhs]
         ${rt}$, intent(in) :: b${nd}$     
         !> [optional] Can A data be overwritten and destroyed?
         logical(lk), optional, intent(in) :: overwrite_a
         !> [optional] state return flag. On error if not requested, the code will stop
         type(linalg_state_type), intent(out) :: err
         !> Result array/matrix x[n] or x[n,nrhs]
         ${rt}$, allocatable, target :: x${nd}$
     end function stdlib_linalg_${ri}$_solve_${ndsuf}$
     pure module function stdlib_linalg_${ri}$_pure_solve_${ndsuf}$(a,b) result(x)
         !> Input matrix a[n,n]
         ${rt}$, intent(in) :: a(:,:)
         !> Right hand side vector or array, b[n] or b[n,nrhs]
         ${rt}$, intent(in) :: b${nd}$
         !> Result array/matrix x[n] or x[n,nrhs]
         ${rt}$, allocatable, target :: x${nd}$
     end function stdlib_linalg_${ri}$_pure_solve_${ndsuf}$
     #:endfor
     #:endfor    
  end interface solve

  ! Solve linear system Ax = b using LU decomposition (subroutine interface).
  interface solve_lu
     !! version: experimental 
     !!
     !! Solves the linear system \( A \cdot x = b \) for the unknown vector \( x \) from a square matrix \( A \). 
     !! ([Specification](../page/specs/stdlib_linalg.html#solve-lu-solves-a-linear-matrix-equation-or-a-linear-system-of-equations-subroutine-interface))
     !!
     !!### Summary 
     !! Subroutine interface for solving a linear system using LU decomposition.
     !!
     !!### Description
     !! 
     !! This interface provides methods for computing the solution of a linear matrix system using
     !! a subroutine. Supported data types include `real` and `complex`. No assumption is made on the matrix 
     !! structure. Preallocated space for the solution vector `x` is user-provided, and it may be provided
     !! for the array of pivot indices, `pivot`. If all pre-allocated work spaces are provided, no internal 
     !! memory allocations take place when using this interface.     
     !! The function can solve simultaneously either one (from a 1-d right-hand-side vector `b(:)`) 
     !! or several (from a 2-d right-hand-side vector `b(:,:)`) systems.
     !! 
     !!@note The solution is based on LAPACK's generic LU decomposition based solvers `*GESV`.
     !!        
     #:for nd,ndsuf,nde in ALL_RHS
     #:for rk,rt,ri in RC_KINDS_TYPES
     pure module subroutine stdlib_linalg_${ri}$_solve_lu_${ndsuf}$(a,b,x,pivot,overwrite_a,err)     
         !> Input matrix a[n,n]
         ${rt}$, intent(inout), target :: a(:,:)
         !> Right hand side vector or array, b[n] or b[n,nrhs]
         ${rt}$, intent(in) :: b${nd}$
         !> Result array/matrix x[n] or x[n,nrhs]     
         ${rt}$, intent(inout), contiguous, target :: x${nd}$
         !> [optional] Storage array for the diagonal pivot indices
         integer(ilp), optional, intent(inout), target :: pivot(:)
         !> [optional] Can A data be overwritten and destroyed?
         logical(lk), optional, intent(in) :: overwrite_a
         !> [optional] state return flag. On error if not requested, the code will stop
         type(linalg_state_type), optional, intent(out) :: err
     end subroutine stdlib_linalg_${ri}$_solve_lu_${ndsuf}$
     #:endfor
     #:endfor    
  end interface solve_lu     
     
  ! Least squares solution to system Ax=b, i.e. such that the 2-norm abs(b-Ax) is minimized.
  interface lstsq
    !! version: experimental 
    !!
    !! Computes the squares solution to system \( A \cdot x = b \). 
    !! ([Specification](../page/specs/stdlib_linalg.html#lstsq-computes-the-least-squares-solution-to-a-linear-matrix-equation))
    !! 
    !!### Summary 
    !! Interface for computing least squares, i.e. the 2-norm \( || (b-A \cdot x ||_2 \) minimizing solution.
    !!
    !!### Description
    !! 
    !! This interface provides methods for computing the least squares of a linear matrix system.
    !! Supported data types include `real` and `complex`.
    !! 
    !!@note The solution is based on LAPACK's singular value decomposition `*GELSD` methods.
    !! 
    #:for nd,ndsuf,nde in ALL_RHS
    #:for rk,rt,ri in RC_KINDS_TYPES
      module function stdlib_linalg_${ri}$_lstsq_${ndsuf}$(a,b,cond,overwrite_a,rank,err) result(x)
         !> Input matrix a[n,n]
         ${rt}$, intent(inout), target :: a(:,:)
         !> Right hand side vector or array, b[n] or b[n,nrhs]
         ${rt}$, intent(in) :: b${nd}$
         !> [optional] cutoff for rank evaluation: singular values s(i)<=cond*maxval(s) are considered 0.
         real(${rk}$), optional, intent(in) :: cond
         !> [optional] Can A,b data be overwritten and destroyed?
         logical(lk), optional, intent(in) :: overwrite_a
         !> [optional] Return rank of A
         integer(ilp), optional, intent(out) :: rank
         !> [optional] state return flag. On error if not requested, the code will stop
         type(linalg_state_type), optional, intent(out) :: err
         !> Result array/matrix x[n] or x[n,nrhs]
         ${rt}$, allocatable, target :: x${nd}$
      end function stdlib_linalg_${ri}$_lstsq_${ndsuf}$
    #:endfor
    #:endfor
  end interface lstsq

  ! Least squares solution to system Ax=b, i.e. such that the 2-norm abs(b-Ax) is minimized.
  interface solve_lstsq
    !! version: experimental 
    !!
    !! Computes the squares solution to system \( A \cdot x = b \). 
    !! ([Specification](../page/specs/stdlib_linalg.html#solve-lstsq-compute-the-least-squares-solution-to-a-linear-matrix-equation-subroutine-interface))
    !! 
    !!### Summary 
    !! Subroutine interface for computing least squares, i.e. the 2-norm \( || (b-A \cdot x ||_2 \) minimizing solution.
    !!
    !!### Description
    !! 
    !! This interface provides methods for computing the least squares of a linear matrix system using 
    !! a subroutine. Supported data types include `real` and `complex`. If pre-allocated work spaces 
    !! are provided, no internal memory allocations take place when using this interface.
    !! 
    !!@note The solution is based on LAPACK's singular value decomposition `*GELSD` methods.
    !! 
    #:for nd,ndsuf,nde in ALL_RHS
    #:for rk,rt,ri in RC_KINDS_TYPES
      module subroutine stdlib_linalg_${ri}$_solve_lstsq_${ndsuf}$(a,b,x,real_storage,int_storage,&
                        #{if rt.startswith('c')}#cmpl_storage,#{endif}#cond,singvals,overwrite_a,rank,err) 
         !> Input matrix a[n,n]
         ${rt}$, intent(inout), target :: a(:,:)
         !> Right hand side vector or array, b[n] or b[n,nrhs]
         ${rt}$, intent(in) :: b${nd}$
         !> Result array/matrix x[n] or x[n,nrhs]
         ${rt}$, intent(inout), contiguous, target :: x${nd}$     
         !> [optional] real working storage space
         real(${rk}$), optional, intent(inout), target :: real_storage(:)
         !> [optional] integer working storage space
         integer(ilp), optional, intent(inout), target :: int_storage(:)
         #:if rt.startswith('complex')   
         !> [optional] complex working storage space
         ${rt}$, optional, intent(inout), target :: cmpl_storage(:)                  
         #:endif
         !> [optional] cutoff for rank evaluation: singular values s(i)<=cond*maxval(s) are considered 0.
         real(${rk}$), optional, intent(in) :: cond
         !> [optional] list of singular values [min(m,n)], in descending magnitude order, returned by the SVD
         real(${rk}$), optional, intent(out), target :: singvals(:)                  
         !> [optional] Can A,b data be overwritten and destroyed?
         logical(lk), optional, intent(in) :: overwrite_a
         !> [optional] Return rank of A
         integer(ilp), optional, intent(out) :: rank
         !> [optional] state return flag. On error if not requested, the code will stop
         type(linalg_state_type), optional, intent(out) :: err
      end subroutine stdlib_linalg_${ri}$_solve_lstsq_${ndsuf}$
    #:endfor
    #:endfor
  end interface solve_lstsq

  ! Return the working array space required by the least squares solver
  interface lstsq_space
    !! version: experimental 
    !!
    !! Computes the integer, real [, complex] working space required by the least-squares solver
    !! ([Specification](../page/specs/stdlib_linalg.html#lstsq-space-compute-internal-working-space-requirements-for-the-least-squares-solver))
    !! 
    !!### Description
    !! 
    !! This interface provides sizes of integer, real [, complex] working spaces required by the 
    !! least-squares solver. These sizes can be used to pre-allocated working arrays in case several 
    !! repeated least-squares solutions to a same system are sought. If pre-allocated working arrays 
    !! are provided, no internal allocations will take place.
    !! 
    #:for nd,ndsuf,nde in ALL_RHS
    #:for rk,rt,ri in RC_KINDS_TYPES
      pure module subroutine stdlib_linalg_${ri}$_lstsq_space_${ndsuf}$(a,b,lrwork,liwork#{if rt.startswith('c')}#,lcwork#{endif}#)
         !> Input matrix a[m,n]
         ${rt}$, intent(in), target :: a(:,:)
         !> Right hand side vector or array, b[n] or b[n,nrhs]
         ${rt}$, intent(in) :: b${nd}$
         !> Size of the working space arrays                
         integer(ilp), intent(out) :: lrwork,liwork#{if rt.startswith('c')}#,lcwork#{endif}#         
      end subroutine stdlib_linalg_${ri}$_lstsq_space_${ndsuf}$
    #:endfor
    #:endfor
  end interface lstsq_space

  ! QR factorization of rank-2 array A
  interface qr
    !! version: experimental 
    !!
    !! Computes the QR factorization of matrix \( A = Q R \). 
    !! ([Specification](../page/specs/stdlib_linalg.html#qr-compute-the-qr-factorization-of-a-matrix))
    !! 
    !!### Summary 
    !! Compute 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 \). 
    !!
    !!### Description
    !! 
    !! This interface provides methods for computing the QR factorization of a matrix. 
    !! Supported data types include `real` and `complex`. If a pre-allocated work space 
    !! is provided, no internal memory allocations take place when using this interface.
    !!
    !! Given `k = min(m,n)`, one can write \( A = \( Q_1  Q_2 \) \cdot \( \frac{R_1}{0}\) \). 
    !! The user may want the full problem (provide `shape(Q)==[m,m]`, `shape(R)==[m,n]`) or the reduced  
    !! problem only: \( A = Q_1 R_1 \) (provide `shape(Q)==[m,k]`, `shape(R)==[k,n]`).
    !! 
    !!@note The solution is based on LAPACK's QR factorization (`*GEQRF`) and ordered matrix output (`*ORGQR`, `*UNGQR`). 
    !!     
    #:for rk,rt,ri in RC_KINDS_TYPES
      pure module subroutine stdlib_linalg_${ri}$_qr(a,q,r,overwrite_a,storage,err) 
         !> Input matrix a[m,n]
         ${rt}$, intent(inout), target :: a(:,:)
         !> Orthogonal matrix Q ([m,m], or [m,k] if reduced)
         ${rt}$, intent(out), contiguous, target :: q(:,:)
         !> Upper triangular matrix R ([m,n], or [k,n] if reduced)
         ${rt}$, intent(out), contiguous, target :: r(:,:)
         !> [optional] Can A data be overwritten and destroyed?
         logical(lk), optional, intent(in) :: overwrite_a
         !> [optional] Provide pre-allocated workspace, size to be checked with qr_space
         ${rt}$, intent(out), optional, target :: storage(:)
         !> [optional] state return flag. On error if not requested, the code will stop
         type(linalg_state_type), optional, intent(out) :: err
      end subroutine stdlib_linalg_${ri}$_qr
    #:endfor    
  end interface qr

  ! Return the working array space required by the QR factorization solver
  interface qr_space
    !! version: experimental 
    !!
    !! Computes the working array space required by the QR factorization solver
    !! ([Specification](../page/specs/stdlib_linalg.html#qr-space-compute-internal-working-space-requirements-for-the-qr-factorization))
    !! 
    !!### Description
    !! 
    !! This interface returns the size of the `real` or `complex` working storage required by the 
    !! QR factorization solver. The working size only depends on the kind (`real` or `complex`) and size of
    !! the matrix being factorized. Storage size can be used to pre-allocate a working array in case several 
    !! repeated QR factorizations to a same-size matrix are sought. If pre-allocated working arrays 
    !! are provided, no internal allocations will take place during the factorization.
    !!     
    #:for rk,rt,ri in RC_KINDS_TYPES
      pure module subroutine get_qr_${ri}$_workspace(a,lwork,err)
         !> Input matrix a[m,n]
         ${rt}$, intent(in), target :: a(:,:)
         !> Minimum workspace size for both operations
         integer(ilp), intent(out) :: lwork
         !> State return flag. Returns an error if the query failed
         type(linalg_state_type), optional, intent(out) :: err
      end subroutine get_qr_${ri}$_workspace
    #:endfor
  end interface qr_space
 
  ! Schur decomposition of rank-2 array A
  interface schur
    !! version: experimental
    !!
    !! Computes the Schur decomposition of matrix \( A = Z T Z^H \).
    !! ([Specification](../page/specs/stdlib_linalg.html#schur-compute-the-schur-decomposition-of-a-matrix))
    !!
    !!### Summary
    !! Compute the Schur decomposition of a `real` or `complex` matrix: \( A = Z T Z^H \), where \( Z \) is
    !! orthonormal/unitary and \( T \) is upper-triangular or quasi-upper-triangular. Matrix \( A \) has size `[m,m]`.
    !!
    !!### Description
    !! 
    !! This interface provides methods for computing the Schur decomposition of a matrix.
    !! Supported data types include `real` and `complex`. If a pre-allocated workspace is provided, no internal 
    !! memory allocations take place when using this interface.
    !!
    !! The output matrix \( T \) is upper-triangular for `complex` input matrices and quasi-upper-triangular
    !! for `real` input matrices (with possible \( 2 \times 2 \) blocks on the diagonal).
    !!
    !!@note The solution is based on LAPACK's Schur decomposition routines (`*GEES`). 
    !! 
    #:for rk,rt,ri in RC_KINDS_TYPES
      module subroutine stdlib_linalg_${ri}$_schur(a, t, z, eigvals, overwrite_a, storage, err)
         !> Input matrix a[m,m]
         ${rt}$, intent(inout), target :: a(:,:)
         !> Schur form of A: upper-triangular or quasi-upper-triangular matrix T
         ${rt}$, intent(out), contiguous, target :: t(:,:)
         !> Unitary/orthonormal transformation matrix Z
         ${rt}$, optional, intent(out), contiguous, target :: z(:,:)
         !> [optional] Output eigenvalues that appear on the diagonal of T
         complex(${rk}$), optional, intent(out), contiguous, target :: eigvals(:)
         !> [optional] Can A data be overwritten and destroyed?
         logical(lk), optional, intent(in) :: overwrite_a                 
         !> [optional] Provide pre-allocated workspace, size to be checked with schur_space         
         ${rt}$, optional, intent(inout), target :: storage(:)
         !> [optional] State return flag. On error if not requested, the code will stop
         type(linalg_state_type), optional, intent(out) :: err
      end subroutine stdlib_linalg_${ri}$_schur
      
      ! Schur decomposition subroutine: real eigenvalue interface
      module subroutine stdlib_linalg_real_eig_${ri}$_schur(a,t,z,eigvals,overwrite_a,storage,err)
         !> Input matrix a[m,m]
         ${rt}$, intent(inout), target :: a(:,:)
         !> Schur form of A: upper-triangular or quasi-upper-triangular matrix T
         ${rt}$, intent(out), contiguous, target :: t(:,:)
         !> Unitary/orthonormal transformation matrix Z
         ${rt}$, optional, intent(out), contiguous, target :: z(:,:)
         !> Output real eigenvalues that appear on the diagonal of T
         real(${rk}$), intent(out), contiguous, target :: eigvals(:)
         !> [optional] Provide pre-allocated workspace, size to be checked with schur_space
         ${rt}$, optional, intent(inout), target :: storage(:)
         !> [optional] Can A data be overwritten and destroyed?
         logical(lk), optional, intent(in) :: overwrite_a        
         !> [optional] State return flag. On error if not requested, the code will stop
         type(linalg_state_type), optional, intent(out) :: err        
      end subroutine stdlib_linalg_real_eig_${ri}$_schur 
    #:endfor
  end interface schur

  ! Return the working array space required by the Schur decomposition solver
  interface schur_space
    !! version: experimental
    !!
    !! Computes the working array space required by the Schur decomposition solver
    !! ([Specification](../page/specs/stdlib_linalg.html#schur-space-compute-internal-working-space-requirements-for-the-schur-decomposition))
    !!
    !!### Description
    !! 
    !! This interface returns the size of the `real` or `complex` working storage required by the 
    !! Schur decomposition solver. The working size only depends on the kind (`real` or `complex`) and size of
    !! the matrix being decomposed. Storage size can be used to pre-allocate a working array in case several 
    !! repeated Schur decompositions of same-size matrices are sought. If pre-allocated working arrays 
    !! are provided, no internal allocations will take place during the decomposition.
    !!     
    #:for rk,rt,ri in RC_KINDS_TYPES
      module subroutine get_schur_${ri}$_workspace(a,lwork,err)
         !> Input matrix a[m,m]
         ${rt}$, intent(in), target :: a(:,:)
         !> Minimum workspace size for the decomposition operation
         integer(ilp), intent(out) :: lwork
         !> State return flag. Returns an error if the query failed
         type(linalg_state_type), optional, intent(out) :: err
      end subroutine get_schur_${ri}$_workspace
    #:endfor
  end interface schur_space

  interface det
    !! version: experimental 
    !!
    !! Computes the determinant of a square matrix
    !! ([Specification](../page/specs/stdlib_linalg.html#det-computes-the-determinant-of-a-square-matrix))
    !! 
    !!### Summary 
    !! Interface for computing matrix determinant.
    !!
    !!### Description
    !! 
    !! This interface provides methods for computing the determinant of a matrix.
    !! Supported data types include `real` and `complex`.
    !! 
    !!@note The provided functions are intended for square matrices only.          
    !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
    !! 
    !!### Example
    !!
    !!```fortran
    !!
    !!    real(sp) :: a(3,3), d
    !!    type(linalg_state_type) :: state  
    !!    a = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
    !!
    !!    ! ...
    !!    d = det(a,err=state)
    !!    if (state%ok()) then 
    !!       print *, 'Success! det=',d
    !!    else
    !!       print *, state%print()
    !!    endif
    !!    ! ...
    !!```
    !!     
    #:for rk,rt in RC_KINDS_TYPES
    #:if rk!="xdp" 
    module procedure stdlib_linalg_${rt[0]}$${rk}$determinant
    module procedure stdlib_linalg_pure_${rt[0]}$${rk}$determinant
    #:endif
    #:endfor
  end interface det

  interface operator(.det.)
    !! version: experimental 
    !!
    !! Determinant operator of a square matrix
    !! ([Specification](../page/specs/stdlib_linalg.html#det-determinant-operator-of-a-square-matrix))
    !!
    !!### Summary
    !! Pure operator interface for computing matrix determinant.
    !!
    !!### Description
    !! 
    !! This pure operator interface provides a convenient way to compute the determinant of a matrix.
    !! Supported data types include real and complex.
    !!
    !!@note The provided functions are intended for square matrices.
    !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
    !!
    !!### Example
    !!
    !!```fortran
    !!
    !!    ! ...
    !!    real(sp) :: matrix(3,3), d
    !!    matrix = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
    !!    d = .det.matrix
    !!    ! ...
    !! 
    !!```
    !     
    #:for rk,rt in RC_KINDS_TYPES
    #:if rk!="xdp"
    module procedure stdlib_linalg_pure_${rt[0]}$${rk}$determinant
    #:endif
    #:endfor        
  end interface operator(.det.)

  interface
    #:for rk,rt in RC_KINDS_TYPES
    #:if rk!="xdp" 
    module function stdlib_linalg_${rt[0]}$${rk}$determinant(a,overwrite_a,err) result(det)
        !> Input matrix a[m,n]
        ${rt}$, intent(inout), target :: a(:,:)
        !> [optional] Can A data be overwritten and destroyed?
        logical(lk), optional, intent(in) :: overwrite_a
        !> State return flag. 
        type(linalg_state_type), intent(out) :: err
        !> Matrix determinant
        ${rt}$ :: det
    end function stdlib_linalg_${rt[0]}$${rk}$determinant
    pure module function stdlib_linalg_pure_${rt[0]}$${rk}$determinant(a) result(det)
        !> Input matrix a[m,n]
        ${rt}$, intent(in) :: a(:,:)
        !> Matrix determinant
        ${rt}$ :: det                
    end function stdlib_linalg_pure_${rt[0]}$${rk}$determinant
    #:endif
    #:endfor
  end interface  

  ! Matrix Inverse: Function interface
  interface inv
    !! version: experimental 
    !!
    !! Inverse of a square matrix
    !! ([Specification](../page/specs/stdlib_linalg.html#inv-inverse-of-a-square-matrix))
    !!
    !!### Summary
    !! This interface provides methods for computing the inverse of a square `real` or `complex` matrix.
    !! The inverse \( A^{-1} \) is defined such that \( A \cdot A^{-1} = A^{-1} \cdot A = I_n \).
    !!
    !!### Description
    !!     
    !! This function interface provides methods that return the inverse of a square matrix.    
    !! Supported data types include `real` and `complex`. 
    !! The inverse matrix \( A^{-1} \) is returned as a function result. 
    !! Exceptions are raised in case of singular matrix or invalid size, and trigger an `error stop` if 
    !! the state flag `err` is not provided. 
    !!
    !!@note The provided functions are intended for square matrices.
    !!       
    #:for rk,rt,ri in RC_KINDS_TYPES
    module function stdlib_linalg_inverse_${ri}$(a,err) result(inva)
         !> Input matrix a[n,n]
         ${rt}$, intent(in) :: a(:,:)
         !> Output matrix inverse
         ${rt}$, allocatable :: inva(:,:)
         !> [optional] state return flag. On error if not requested, the code will stop
         type(linalg_state_type), optional, intent(out) :: err
    end function stdlib_linalg_inverse_${ri}$
    #:endfor
  end interface inv

  ! Matrix Inverse: Subroutine interface - in-place inversion
  interface invert
    !! version: experimental 
    !!
    !! Inversion of a square matrix
    !! ([Specification](../page/specs/stdlib_linalg.html#invert-inversion-of-a-square-matrix))
    !!
    !!### Summary
    !! This interface provides methods for inverting 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 \).
    !!
    !!### Description
    !!     
    !! This subroutine interface provides a way to compute the inverse of a matrix.    
    !! Supported data types include `real` and `complex`. 
    !! The user may provide a unique matrix argument `a`. In this case, `a` is replaced by the inverse matrix.
    !! on output. Otherwise, one may provide two separate arguments: an input matrix `a` and an output matrix 
    !! `inva`. In this case, `a` will not be modified, and the inverse is returned in `inva`.
    !! Pre-allocated storage may be provided for the array of pivot indices, `pivot`. If all pre-allocated 
    !! work spaces are provided, no internal memory allocations take place when using this interface.             
    !!
    !!@note The provided subroutines are intended for square matrices.
    !!      
    #:for rk,rt,ri in RC_KINDS_TYPES
    module subroutine stdlib_linalg_invert_inplace_${ri}$(a,pivot,err)
         !> Input matrix a[n,n]
         ${rt}$, intent(inout) :: a(:,:)
         !> [optional] Storage array for the diagonal pivot indices
         integer(ilp), optional, intent(inout), target :: pivot(:)         
         !> [optional] state return flag. On error if not requested, the code will stop
         type(linalg_state_type), optional, intent(out) :: err
    end subroutine stdlib_linalg_invert_inplace_${ri}$
    ! Compute the square matrix inverse of a
    module subroutine stdlib_linalg_invert_split_${ri}$(a,inva,pivot,err)
         !> Input matrix a[n,n].
         ${rt}$, intent(in) :: a(:,:)
         !> Inverse matrix a[n,n]. 
         ${rt}$, intent(out) :: inva(:,:)         
         !> [optional] Storage array for the diagonal pivot indices
         integer(ilp), optional, intent(inout), target :: pivot(:)
         !> [optional] state return flag. On error if not requested, the code will stop
         type(linalg_state_type), optional, intent(out) :: err    
    end subroutine stdlib_linalg_invert_split_${ri}$
    #:endfor
  end interface invert

  ! Matrix Inverse: Operator interface
  interface operator(.inv.)
    !! version: experimental 
    !!
    !! Inverse operator of a square matrix
    !! ([Specification](../page/specs/stdlib_linalg.html#inv-inverse-operator-of-a-square-matrix))
    !!
    !!### Summary
    !! Operator interface for computing the inverse of a square `real` or `complex` matrix.
    !!
    !!### Description
    !! 
    !! This operator interface provides a convenient way to compute the inverse of a matrix.
    !! Supported data types include `real` and `complex`. On input errors or singular matrix, 
    !! NaNs will be returned.
    !!
    !!@note The provided functions are intended for square matrices.
    !!
    #:for rk,rt,ri in RC_KINDS_TYPES
    module function stdlib_linalg_inverse_${ri}$_operator(a) result(inva)
         !> Input matrix a[n,n]
         ${rt}$, intent(in) :: a(:,:)
         !> Result matrix
         ${rt}$, allocatable :: inva(:,:)        
    end function stdlib_linalg_inverse_${ri}$_operator
    #:endfor
  end interface operator(.inv.)


  ! Moore-Penrose Pseudo-Inverse: Function interface
  interface pinv
    !! version: experimental 
    !!
    !! Pseudo-inverse of a matrix
    !! ([Specification](../page/specs/stdlib_linalg.html#pinv-moore-penrose-pseudo-inverse-of-a-matrix))
    !!
    !!### Summary
    !! This interface provides methods for computing the Moore-Penrose pseudo-inverse of a matrix.
    !! The pseudo-inverse \( A^{+} \) is a generalization of the matrix inverse, computed for square, singular, 
    !! or rectangular matrices. It is defined such that it 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 \)
    !!
    !!### Description
    !!     
    !! This function interface provides methods that return the Moore-Penrose pseudo-inverse of a matrix.    
    !! Supported data types include `real` and `complex`. 
    !! The pseudo-inverse \( A^{+} \) is returned as a function result. The computation is based on the 
    !! singular value decomposition (SVD). An optional relative tolerance `rtol` is provided to control the 
    !! inclusion of singular values during inversion. Singular values below \( \text{rtol} \cdot \sigma_{\max} \) 
    !! are treated as zero, where \( \sigma_{\max} \) is the largest singular value. If `rtol` is not provided, 
    !! a default threshold is applied.
    !! 
    !! Exceptions are raised in case of computational errors or invalid input, and trigger an `error stop` 
    !! if the state flag `err` is not provided. 
    !!
    !!@note The provided functions are intended for both rectangular and square matrices.
    !!       
    #:for rk,rt,ri in RC_KINDS_TYPES
    module function stdlib_linalg_pseudoinverse_${ri}$(a,rtol,err) result(pinva)
        !> Input matrix a[m,n]
        ${rt}$, intent(in), target :: a(:,:)
        !> [optional] Relative tolerance for singular value cutoff
        real(${rk}$), optional, intent(in) :: rtol         
        !> [optional] State return flag. On error if not requested, the code will stop
        type(linalg_state_type), optional, intent(out) :: err
        !> Output matrix pseudo-inverse [n,m]
        ${rt}$ :: pinva(size(a,2,kind=ilp),size(a,1,kind=ilp))         
     end function stdlib_linalg_pseudoinverse_${ri}$
    #:endfor
  end interface pinv

  ! Moore-Penrose Pseudo-Inverse: Subroutine interface 
  interface pseudoinvert
    !! version: experimental 
    !!
    !! Computation of the Moore-Penrose pseudo-inverse
    !! ([Specification](../page/specs/stdlib_linalg.html#pseudoinvert-moore-penrose-pseudo-inverse-of-a-matrix))
    !!
    !!### Summary
    !! This interface provides methods for computing the Moore-Penrose pseudo-inverse of a rectangular 
    !! or square `real` or `complex` matrix.
    !! The pseudo-inverse \( A^{+} \) generalizes the matrix inverse and satisfies the 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 \)
    !!
    !!### Description
    !!     
    !! This subroutine interface provides a way to compute the Moore-Penrose pseudo-inverse of a matrix.    
    !! Supported data types include `real` and `complex`. 
    !! Users must provide two matrices: the input matrix `a` [m,n] and the output pseudo-inverse `pinva` [n,m]. 
    !! The input matrix `a` is used to compute the pseudo-inverse and is not modified. The computed 
    !! pseudo-inverse is stored in `pinva`. The computation is based on the singular value decomposition (SVD).
    !! 
    !! An optional relative tolerance `rtol` is used to control the inclusion of singular values in the 
    !! computation. Singular values below \( \text{rtol} \cdot \sigma_{\max} \) are treated as zero, 
    !! where \( \sigma_{\max} \) is the largest singular value. If `rtol` is not provided, a default 
    !! threshold is applied. 
    !! 
    !! Exceptions are raised in case of computational errors or invalid input, and trigger an `error stop` 
    !! if the state flag `err` is not provided.
    !!
    !!@note The provided subroutines are intended for both rectangular and square matrices.
    !!       
    #:for rk,rt,ri in RC_KINDS_TYPES
    module subroutine stdlib_linalg_pseudoinvert_${ri}$(a,pinva,rtol,err)
        !> Input matrix a[m,n]
        ${rt}$, intent(inout) :: a(:,:)
        !> Output pseudo-inverse matrix [n,m]
        ${rt}$, intent(out) :: pinva(:,:)
        !> [optional] Relative tolerance for singular value cutoff
        real(${rk}$), optional, intent(in) :: rtol
        !> [optional] State return flag. On error if not requested, the code will stop
        type(linalg_state_type), optional, intent(out) :: err
    end subroutine stdlib_linalg_pseudoinvert_${ri}$
    #:endfor
  end interface pseudoinvert

  ! Moore-Penrose Pseudo-Inverse: Operator interface
  interface operator(.pinv.)
    !! version: experimental 
    !!
    !! Pseudo-inverse operator of a matrix
    !! ([Specification](../page/specs/stdlib_linalg.html#pinv-moore-penrose-pseudo-inverse-operator))
    !!
    !!### Summary
    !! Operator interface for computing the Moore-Penrose pseudo-inverse of a `real` or `complex` matrix.
    !!
    !!### Description
    !! 
    !! This operator interface provides a convenient way to compute the Moore-Penrose pseudo-inverse 
    !! of a matrix. Supported data types include `real` and `complex`. The pseudo-inverse \( A^{+} \) 
    !! is computed using singular value decomposition (SVD), with singular values below an internal 
    !! threshold treated as zero.
    !! 
    !! For computational errors or invalid input, the function may return a matrix filled with NaNs.
    !!
    !!@note The provided functions are intended for both rectangular and square matrices.
    !!
    #:for rk,rt,ri in RC_KINDS_TYPES
    module function stdlib_linalg_pinv_${ri}$_operator(a) result(pinva)
         !> Input matrix a[m,n]
         ${rt}$, intent(in), target :: a(:,:)
         !> Result pseudo-inverse matrix
         ${rt}$ :: pinva(size(a,2,kind=ilp),size(a,1,kind=ilp))
    end function stdlib_linalg_pinv_${ri}$_operator
    #:endfor
  end interface operator(.pinv.)


  ! Eigendecomposition of a square matrix: eigenvalues, and optionally eigenvectors
  interface eig    
     !! version: experimental 
     !!
     !! Solves the eigendecomposition \( A \cdot \bar{v} - \lambda \cdot \bar{v} \) for square matrix \( A \). 
     !! ([Specification](../page/specs/stdlib_linalg.html#eig-eigenvalues-and-eigenvectors-of-a-square-matrix))
     !!
     !!### Summary 
     !! Subroutine interface for computing eigenvalues and eigenvectors of a square matrix.
     !!
     !!### Description
     !! 
     !! This interface provides methods for computing the eigenvalues, and optionally eigenvectors, 
     !! of a general square matrix. Supported data types include `real` and `complex`, and no assumption is 
     !! made on the matrix structure. The user may request either left, right, or both 
     !! eigenvectors to be returned. They are returned as columns of a square matrix with the same size as `A`. 
     !! Preallocated space for both eigenvalues `lambda` and the eigenvector matrices must be user-provided.      
     !! 
     !!@note The solution is based on LAPACK's general eigenproblem solvers `*GEEV`.
     !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
     !!       
    #:for rk,rt,ri in RC_KINDS_TYPES
    #:for ep,ei in EIG_PROBLEM_LIST
    module subroutine stdlib_linalg_eig_${ep}$_${ri}$(a#{if ei=='ggev'}#,b#{endif}#,lambda,right,left, &
                                                      overwrite_a#{if ei=='ggev'}#,overwrite_b#{endif}#,err)
     !! Eigendecomposition of matrix A returning an array `lambda` of eigenvalues, 
     !! and optionally right or left eigenvectors.        
         !> Input matrix A[m,n]
         ${rt}$, intent(inout), target :: a(:,:)
         #:if ei=='ggev'
         !> Generalized problem matrix B[n,n]
         ${rt}$, intent(inout), target :: b(:,:)
         #:endif         
         !> Array of eigenvalues
         complex(${rk}$), intent(out) :: lambda(:)
         !> The columns of RIGHT contain the right eigenvectors of A
         complex(${rk}$), optional, intent(out), target :: right(:,:)
         !> The columns of LEFT contain the left eigenvectors of A
         complex(${rk}$), optional, intent(out), target :: left(:,:)
         !> [optional] Can A data be overwritten and destroyed?
         logical(lk), optional, intent(in) :: overwrite_a
         #:if ei=='ggev'
         !> [optional] Can B data be overwritten and destroyed? (default: no)
         logical(lk), optional, intent(in) :: overwrite_b
         #:endif         
         !> [optional] state return flag. On error if not requested, the code will stop
         type(linalg_state_type), optional, intent(out) :: err
    end subroutine stdlib_linalg_eig_${ep}$_${ri}$

    module subroutine stdlib_linalg_real_eig_${ep}$_${ri}$(a#{if ei=='ggev'}#,b#{endif}#,lambda,right,left, &
                                                           overwrite_a#{if ei=='ggev'}#,overwrite_b#{endif}#,err)
     !! Eigendecomposition of matrix A returning an array `lambda` of real eigenvalues, 
     !! and optionally right or left eigenvectors. Returns an error if the eigenvalues had
     !! non-trivial imaginary parts.
         !> Input matrix A[m,n]
         ${rt}$, intent(inout), target :: a(:,:)
         #:if ei=='ggev'
         !> Generalized problem matrix B[n,n]
         ${rt}$, intent(inout), target :: b(:,:)  
         #:endif
         !> Array of real eigenvalues
         real(${rk}$), intent(out) :: lambda(:)
         !> The columns of RIGHT contain the right eigenvectors of A
         complex(${rk}$), optional, intent(out), target :: right(:,:)
         !> The columns of LEFT contain the left eigenvectors of A
         complex(${rk}$), optional, intent(out), target :: left(:,:)
         !> [optional] Can A data be overwritten and destroyed?
         logical(lk), optional, intent(in) :: overwrite_a
         #:if ei=='ggev'
         !> [optional] Can B data be overwritten and destroyed? (default: no)
         logical(lk), optional, intent(in) :: overwrite_b
         #:endif         
         !> [optional] state return flag. On error if not requested, the code will stop
         type(linalg_state_type), optional, intent(out) :: err
    end subroutine stdlib_linalg_real_eig_${ep}$_${ri}$
    #:endfor    
    #:endfor
  end interface eig

  ! Eigenvalues of a square matrix
  interface eigvals
     !! version: experimental 
     !!
     !! Returns the eigenvalues \( lambda \), \( A \cdot \bar{v} - \lambda \cdot \bar{v} \), for square matrix \( A \). 
     !! ([Specification](../page/specs/stdlib_linalg.html#eigvals-eigenvalues-of-a-square-matrix))
     !!
     !!### Summary 
     !! Function interface for computing the eigenvalues of a square matrix.
     !!
     !!### Description
     !! 
     !! This interface provides functions for returning the eigenvalues of a general square matrix. 
     !! Supported data types include `real` and `complex`, and no assumption is made on the matrix structure. 
     !! An `error stop` is thrown in case of failure; otherwise, error information can be returned 
     !! as an optional `type(linalg_state_type)` output flag. 
     !! 
     !!@note The solution is based on LAPACK's general eigenproblem solvers `*GEEV`.
     !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
     !!       
    #:for rk,rt,ri in RC_KINDS_TYPES
    #:for ep,ei in EIG_PROBLEM_LIST
    module function stdlib_linalg_eigvals_${ep}$_${ri}$(a#{if ei=='ggev'}#,b#{endif}#,err) result(lambda)
     !! Return an array of eigenvalues of matrix A.
         !> Input matrix A[m,n]
         ${rt}$, intent(in), dimension(:,:), target :: a 
         #:if ei=='ggev'
         !> Generalized problem matrix B[n,n]
         ${rt}$, intent(inout), dimension(:,:), target :: b         
         #:endif                  
         !> [optional] state return flag. On error if not requested, the code will stop
         type(linalg_state_type), intent(out) :: err
         !> Array of singular values
         complex(${rk}$), allocatable :: lambda(:)        
    end function stdlib_linalg_eigvals_${ep}$_${ri}$
    
    module function stdlib_linalg_eigvals_noerr_${ep}$_${ri}$(a#{if ei=='ggev'}#,b#{endif}#) result(lambda)
     !! Return an array of eigenvalues of matrix A.
         !> Input matrix A[m,n]
         ${rt}$, intent(in), dimension(:,:), target :: a 
         #:if ei=='ggev'
         !> Generalized problem matrix B[n,n]
         ${rt}$, intent(inout), dimension(:,:), target :: b         
         #:endif                  
         !> Array of singular values
         complex(${rk}$), allocatable :: lambda(:)
    end function stdlib_linalg_eigvals_noerr_${ep}$_${ri}$
    #:endfor
    #:endfor
  end interface eigvals
     
  ! Eigendecomposition of a real symmetric or complex hermitian matrix
  interface eigh
     !! version: experimental 
     !!
     !! Solves the eigendecomposition \( A \cdot \bar{v} - \lambda \cdot \bar{v} \) for a real symmetric 
     !! \( A = A^T \) or complex Hermitian \( A = A^H \) square matrix. 
     !! ([Specification](../page/specs/stdlib_linalg.html#eigh-eigenvalues-and-eigenvectors-of-a-real-symmetric-or-complex-hermitian-square-matrix))
     !!
     !!### Summary 
     !! Subroutine interface for computing eigenvalues and eigenvectors of a real symmetric or complex Hermitian square matrix.
     !!
     !!### Description
     !! 
     !! This interface provides methods for computing the eigenvalues, and optionally eigenvectors, 
     !! of a real symmetric or complex Hermitian square matrix. Supported data types include `real` and `complex`. 
     !! The matrix must be symmetric (if `real`) or Hermitian (if `complex`). Only the lower or upper 
     !! half of the matrix is accessed, and the user can select which using the optional `upper_a` 
     !! flag (default: use lower half). The vectors are orthogonal, and may be returned as columns of an optional 
     !! matrix `vectors` with the same kind and size as `A`. 
     !! Preallocated space for both eigenvalues `lambda` and the eigenvector matrix must be user-provided.      
     !! 
     !!@note The solution is based on LAPACK's eigenproblem solvers `*SYEV`/`*HEEV`.
     !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
     !!      
    #:for rk,rt,ri in RC_KINDS_TYPES
    module subroutine stdlib_linalg_eigh_${ri}$(a,lambda,vectors,upper_a,overwrite_a,err)
     !! Eigendecomposition of a real symmetric or complex Hermitian matrix A returning an array `lambda` 
     !! of eigenvalues, and optionally right or left eigenvectors.        
         !> Input matrix A[m,n]
         ${rt}$, intent(inout), target :: a(:,:)
         !> Array of eigenvalues
         real(${rk}$), intent(out) :: lambda(:)
         !> The columns of vectors contain the orthonormal eigenvectors of A
         ${rt}$, optional, intent(out), target :: vectors(:,:)
         !> [optional] Can A data be overwritten and destroyed?
         logical(lk), optional, intent(in) :: overwrite_a
         !> [optional] Should the upper/lower half of A be used? Default: lower
         logical(lk), optional, intent(in) :: upper_a
         !> [optional] state return flag. On error if not requested, the code will stop
         type(linalg_state_type), optional, intent(out) :: err
    end subroutine stdlib_linalg_eigh_${ri}$
    #:endfor        
  end interface eigh
     
  ! Eigenvalues of a real symmetric or complex hermitian matrix
  interface eigvalsh
     !! version: experimental 
     !!
     !! Returns the eigenvalues \( lambda \), \( A \cdot \bar{v} - \lambda \cdot \bar{v} \), for a real
     !! symmetric \( A = A^T \) or complex Hermitian \( A = A^H \) square matrix. 
     !! ([Specification](../page/specs/stdlib_linalg.html#eigvalsh-eigenvalues-of-a-real-symmetric-or-complex-hermitian-square-matrix))
     !!
     !!### Summary 
     !! Function interface for computing the eigenvalues of a real symmetric or complex hermitian square matrix.
     !!
     !!### Description
     !! 
     !! This interface provides functions for returning the eigenvalues of a real symmetric or complex Hermitian
     !! square matrix. Supported data types include `real` and `complex`. The matrix must be symmetric 
     !! (if `real`) or Hermitian (if `complex`). Only the lower or upper half of the matrix is accessed, 
     !! and the user can select which using the optional `upper_a` flag (default: use lower half). 
     !! An `error stop` is thrown in case of failure; otherwise, error information can be returned 
     !! as an optional `type(linalg_state_type)` output flag. 
     !! 
     !!@note The solution is based on LAPACK's eigenproblem solvers `*SYEV`/`*HEEV`.
     !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
     !!         
    #:for rk,rt,ri in RC_KINDS_TYPES
    module function stdlib_linalg_eigvalsh_${ri}$(a,upper_a,err) result(lambda)
     !! Return an array of eigenvalues of real symmetric / complex hermitian A
         !> Input matrix A[m,n]
         ${rt}$, intent(in), target :: a(:,:)
         !> [optional] Should the upper/lower half of A be used? Default: lower
         logical(lk), optional, intent(in) :: upper_a         
         !> [optional] state return flag. On error if not requested, the code will stop
         type(linalg_state_type), intent(out) :: err
         !> Array of singular values
         real(${rk}$), allocatable :: lambda(:)    
    end function stdlib_linalg_eigvalsh_${ri}$
    
    module function stdlib_linalg_eigvalsh_noerr_${ri}$(a,upper_a) result(lambda)
     !! Return an array of eigenvalues of real symmetric / complex hermitian A        
         !> Input matrix A[m,n]
         ${rt}$, intent(in), target :: a(:,:)
         !> [optional] Should the upper/lower half of A be used? Default: lower
         logical(lk), optional, intent(in) :: upper_a         
         !> Array of singular values
         real(${rk}$), allocatable :: lambda(:)
    end function stdlib_linalg_eigvalsh_noerr_${ri}$
    #:endfor                
  end interface eigvalsh

  ! Singular value decomposition  
  interface svd 
    !! version: experimental 
    !!
    !! Computes the singular value decomposition of a `real` or `complex` 2d matrix.
    !! ([Specification](../page/specs/stdlib_linalg.html#svd-compute-the-singular-value-decomposition-of-a-rank-2-array-matrix))
    !! 
    !!### Summary 
    !! Interface for computing the singular value decomposition of a `real` or `complex` 2d matrix.
    !!
    !!### Description
    !! 
    !! This interface provides methods for computing the singular value decomposition of a matrix.
    !! Supported data types include `real` and `complex`. The subroutine returns a `real` array of 
    !! singular values, and optionally, left- and right- singular vector matrices, `U` and `V`. 
    !! For a matrix `A` with size [m,n], full matrix storage for `U` and `V` should be [m,m] and [n,n]. 
    !! It is possible to use partial storage [m,k] and [k,n], `k=min(m,n)`, choosing `full_matrices=.false.`.
    !! 
    !!@note The solution is based on LAPACK's singular value decomposition `*GESDD` methods.
    !! 
    !!### Example
    !!
    !!```fortran
    !!    real(sp) :: a(2,3), s(2), u(2,2), vt(3,3) 
    !!    a = reshape([3,2, 2,3, 2,-2],[2,3])
    !!
    !!    call svd(A,s,u,v)
    !!    print *, 'singular values = ',s
    !!```
    !!         
    #:for rk,rt,ri in RC_KINDS_TYPES
    module subroutine stdlib_linalg_svd_${ri}$(a,s,u,vt,overwrite_a,full_matrices,err)
     !!### Summary
     !! Compute singular value decomposition of a matrix \( A = U \cdot S \cdot \V^T \)
     !!
     !!### Description
     !!
     !! This function computes the singular value decomposition of a `real` or `complex` matrix \( A \), 
     !! and returns the array of singular values, and optionally the left matrix \( U \) containing the 
     !! left unitary singular vectors, and the right matrix \( V^T \), containing the right unitary 
     !! singular vectors.
     !!
     !! param: a Input matrix of size [m,n].
     !! param: s Output `real` array of size [min(m,n)] returning a list of singular values.
     !! param: u [optional] Output left singular matrix of size [m,m] or [m,min(m,n)] (.not.full_matrices). Contains singular vectors as columns.
     !! param: vt [optional] Output right singular matrix of size [n,n] or [min(m,n),n] (.not.full_matrices). Contains singular vectors as rows.     
     !! param: overwrite_a [optional] Flag indicating if the input matrix can be overwritten.
     !! param: full_matrices [optional] If `.true.` (default), matrices \( U \) and \( V^T \) have size [m,m], [n,n]. Otherwise, they are [m,k], [k,n] with `k=min(m,n)`.
     !! param: err [optional] State return flag.      
     !!            
         !> Input matrix A[m,n]
         ${rt}$, intent(inout), target :: a(:,:)
         !> Array of singular values
         real(${rk}$), intent(out) :: s(:)
         !> The columns of U contain the left singular vectors 
         ${rt}$, optional, intent(out), target :: u(:,:)
         !> The rows of V^T contain the right singular vectors
         ${rt}$, optional, intent(out), target :: vt(:,:)
         !> [optional] Can A data be overwritten and destroyed?
         logical(lk), optional, intent(in) :: overwrite_a
         !> [optional] full matrices have shape(u)==[m,m], shape(vh)==[n,n] (default); otherwise
         !> they are shape(u)==[m,k] and shape(vh)==[k,n] with k=min(m,n)
         logical(lk), optional, intent(in) :: full_matrices
         !> [optional] state return flag. On error if not requested, the code will stop
         type(linalg_state_type), optional, intent(out) :: err
      end subroutine stdlib_linalg_svd_${ri}$
    #:endfor     
  end interface svd

  ! Singular values
  interface svdvals
    !! version: experimental 
    !!
    !! Computes the singular values of a `real` or `complex` 2d matrix.
    !! ([Specification](../page/specs/stdlib_linalg.html#svdvals-compute-the-singular-values-of-a-rank-2-array-matrix))
    !! 
    !!### Summary 
    !! 
    !! Function interface for computing the array of singular values from the singular value decomposition 
    !! of a `real` or `complex` 2d matrix.
    !!
    !!### Description
    !! 
    !! This interface provides methods for computing the singular values a 2d matrix.
    !! Supported data types include `real` and `complex`. The function returns a `real` array of 
    !! singular values, with size [min(m,n)]. 
    !! 
    !!@note The solution is based on LAPACK's singular value decomposition `*GESDD` methods.
    !! 
    !!### Example
    !!
    !!```fortran
    !!    real(sp) :: a(2,3), s(2)
    !!    a = reshape([3,2, 2,3, 2,-2],[2,3])
    !!
    !!    s = svdvals(A)
    !!    print *, 'singular values = ',s
    !!```
    !!   
    #:for rk,rt,ri in RC_KINDS_TYPES
    module function stdlib_linalg_svdvals_${ri}$(a,err) result(s)
     !!### Summary
     !! Compute singular values \(S \) from the singular-value decomposition of a matrix \( A = U \cdot S \cdot \V^T \).
     !!
     !!### Description
     !!
     !! This function returns the array of singular values from the singular value decomposition of a `real` 
     !! or `complex` matrix \( A = U \cdot S \cdot V^T \).
     !!
     !! param: a Input matrix of size [m,n].
     !! param: err [optional] State return flag.      
     !!
     !!### Return value
     !! 
     !! param: s `real` array of size [min(m,n)] returning a list of singular values.
     !!           
         !> Input matrix A[m,n]
         ${rt}$, intent(in), target :: a(:,:)
         !> [optional] state return flag. On error if not requested, the code will stop
         type(linalg_state_type), optional, intent(out) :: err
         !> Array of singular values
         real(${rk}$), allocatable :: s(:)
      end function stdlib_linalg_svdvals_${ri}$
    #:endfor             
  end interface svdvals  


  #! Allow for integer or character norm input: i.e., norm(a,2) or norm(a, '2')
  #:set NORM_INPUT_TYPE    = ["character(len=*)","integer(ilp)"]
  #:set NORM_INPUT_SHORT   = ["char","int"]
  #:set NORM_INPUT_OPTIONS = list(zip(NORM_INPUT_TYPE,NORM_INPUT_SHORT))
  ! Vector norms: function interface
  interface norm
     !! version: experimental 
     !!
     !! Computes the vector norm of a generic-rank array \( A \). 
     !! ([Specification](../page/specs/stdlib_linalg.html#norm-computes-the-vector-norm-of-a-generic-rank-array))
     !! 
     !!### Summary 
     !! Return one of several scalar norm metrics of a `real` or `complex` input array \( A \), 
     !! that can have any rank. For generic rank-n arrays, the scalar norm over the whole 
     !! array is returned by default. If `n>=2` and the optional input dimension `dim` is specified, 
     !! a rank `n-1` array is returned with dimension `dim` collapsed, containing all 1D array norms 
     !! evaluated along dimension `dim` only.
     !! 
     !!
     !!### Description
     !! 
     !! This interface provides methods for computing the vector norm(s) of an array.  
     !! Supported data types include `real` and `complex`. 
     !! Input arrays may have generic rank from 1 to ${MAXRANK}$.
     !!
     !! Norm type input is mandatory, and it is provided via the `order` argument. 
     !! This can be provided as either an `integer` value or a `character` string. 
     !! Allowed metrics are: 
     !! - 1-norm \( \sum_i{ \left|a_i\right| } \): `order` = 1 or '1'    
     !! - Euclidean norm \( \sqrt{\sum_i{ a_i^2 }} \): `order` = 2 or '2'
     !! - p-norm \( \left( \sum_i{ \left|a_i\right|^p }\right) ^{1/p} \): `integer` `order`, order>=3
     !! - Infinity norm \( \max_i{ \left|a_i\right| } \): order = huge(0) or 'inf'
     !! - Minus-infinity norm \( \min_i{ \left|a_i\right| } \): order = -huge(0) or '-inf'
     !! 
     !!### Example
     !!
     !!```fortran
     !!
     !!    real(sp) :: a(3,3), na, rown(3)
     !!    a = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
     !!
     !!    ! L2 norm: whole matrix
     !!    na = norm(a, 2)
     !!   
     !!    ! Infinity norm of each row
     !!    rown = norm(a, 'inf', dim=2)
     !!     
     !!```     
     !!
     #:for rk,rt,ri in RC_KINDS_TYPES
     #:for it,ii in NORM_INPUT_OPTIONS
     !> Scalar norms: ${rt}$
     #:for rank in range(1, MAXRANK + 1)        
     pure module function stdlib_linalg_norm_${rank}$D_order_${ii}$_${ri}$(a, order) result(nrm)
        !> Input ${rank}$-d matrix a${ranksuffix(rank)}$
        ${rt}$, intent(in) :: a${ranksuffix(rank)}$
        !> Order of the matrix norm being computed.
        ${it}$, intent(in) :: order
        !> Norm of the matrix.
        real(${rk}$) :: nrm       
     end function stdlib_linalg_norm_${rank}$D_order_${ii}$_${ri}$
     module function stdlib_linalg_norm_${rank}$D_order_err_${ii}$_${ri}$(a, order, err) result(nrm)
        !> Input ${rank}$-d matrix a${ranksuffix(rank)}$
        ${rt}$, intent(in) :: a${ranksuffix(rank)}$
        !> Order of the matrix norm being computed.
        ${it}$, intent(in) :: order
        !> Output state return flag. 
        type(linalg_state_type), intent(out) :: err                         
        !> Norm of the matrix.
        real(${rk}$) :: nrm               
     end function stdlib_linalg_norm_${rank}$D_order_err_${ii}$_${ri}$
     #:endfor
     !> Array norms: ${rt}$
     #:for rank in range(2, MAXRANK + 1)
     pure module function stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$(a, order, dim) result(nrm)
        !> Input matrix a[..]
        ${rt}$, intent(in), target :: a${ranksuffix(rank)}$
        !> Order of the matrix norm being computed.
        ${it}$, intent(in) :: order
        !> Dimension the norm is computed along 
        integer(ilp), intent(in) :: dim
        !> Norm of the matrix. (Same shape as `a`, with `dim` dropped).
        real(${rk}$) :: nrm${reduced_shape('a', rank, 'dim')}$   
     end function stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$
     module function stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_err_${ii}$_${ri}$(a, order, dim, err) result(nrm)
        !> Input matrix a[..]
        ${rt}$, intent(in), target :: a${ranksuffix(rank)}$
        !> Order of the matrix norm being computed.
        ${it}$, intent(in) :: order
        !> Dimension the norm is computed along
        integer(ilp), intent(in) :: dim
        !> Output state return flag. 
        type(linalg_state_type), intent(out) :: err                                 
        !> Norm of the matrix. (Same shape as `a`, with `dim` dropped).
        real(${rk}$) :: nrm${reduced_shape('a', rank, 'dim')}$     
     end function stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_err_${ii}$_${ri}$
     #:endfor
     #:endfor        
     #:endfor            
  end interface norm  

  !> Vector norm: subroutine interface
  interface get_norm
     !! version: experimental 
     !!
     !! Computes the vector norm of a generic-rank array \( A \). 
     !! ([Specification](../page/specs/stdlib_linalg.html#get-norm-computes-the-vector-norm-of-a-generic-rank-array))
     !! 
     !!### Summary 
     !! Subroutine interface that returns one of several scalar norm metrics of a `real` or `complex` 
     !! input array \( A \), that can have any rank. For generic rank-n arrays, the scalar norm over 
     !! the whole array is returned by default. If `n>=2` and the optional input dimension `dim` is 
     !! specified, a rank `n-1` array is returned with dimension `dim` collapsed, containing all 1D 
     !! array norms evaluated along dimension `dim` only.
     !! 
     !!
     !!### Description
     !! 
     !! This `pure subroutine `interface provides methods for computing the vector norm(s) of an array.  
     !! Supported data types include `real` and `complex`. 
     !! Input arrays may have generic rank from 1 to ${MAXRANK}$.
     !!
     !! Norm type input is mandatory, and it is provided via the `order` argument. 
     !! This can be provided as either an `integer` value or a `character` string. 
     !! Allowed metrics are: 
     !! - 1-norm \( \sum_i{ \left|a_i\right| } \): `order` = 1 or '1'    
     !! - Euclidean norm \( \sqrt{\sum_i{ a_i^2 }} \): `order` = 2 or '2'
     !! - p-norm \( \left( \sum_i{ \left|a_i\right|^p }\right) ^{1/p} \): `integer` `order`, order>=3
     !! - Infinity norm \( \max_i{ \left|a_i\right| } \): order = huge(0) or 'inf'
     !! - Minus-infinity norm \( \min_i{ \left|a_i\right| } \): order = -huge(0) or '-inf'
     !!     
     !!### Example
     !!
     !!```fortran
     !!
     !!    real(sp) :: a(3,3), na, rown(3)
     !!    type(linalg_state_type) :: err
     !!    a = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
     !!
     !!    ! L2 norm: whole matrix
     !!    call get_norm(a, na, 2)
     !!   
     !!    ! Infinity norms of each row, with error control
     !!    call get_norm(a, rown, 'inf', dim=2, err=err)     
     !!     
     !!```     
     !!     
     #:for rk,rt,ri in RC_KINDS_TYPES
     #:for it,ii in NORM_INPUT_OPTIONS
        !> Scalar norms: ${rt}$
        #:for rank in range(1, MAXRANK + 1)                    
        pure module subroutine norm_${rank}$D_${ii}$_${ri}$(a, nrm, order, err)
           !> Input ${rank}$-d matrix a${ranksuffix(rank)}$
           ${rt}$, intent(in), target :: a${ranksuffix(rank)}$
           !> Norm of the matrix.
           real(${rk}$), intent(out) :: nrm
           !> Order of the matrix norm being computed.
           ${it}$, intent(in) :: order
           !> [optional] state return flag. On error if not requested, the code will stop
           type(linalg_state_type), intent(out), optional :: err   
        end subroutine norm_${rank}$D_${ii}$_${ri}$
        #:endfor
        !> Array norms: ${rt}$
        #:for rank in range(2, MAXRANK + 1)
        pure module subroutine norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$(a, nrm, order, dim, err)
           !> Input matrix a[..]
           ${rt}$, intent(in) :: a${ranksuffix(rank)}$
           !> Dimension the norm is computed along
           integer(ilp), intent(in) :: dim        
           !> Norm of the matrix. (Same shape as `a`, with `dim` dropped).
           real(${rk}$), intent(out) :: nrm${reduced_shape('a', rank, 'dim')}$     
           !> Order of the matrix norm being computed.
           ${it}$, intent(in) :: order
           !> [optional] state return flag. On error if not requested, the code will stop
           type(linalg_state_type), intent(out), optional :: err           
        end subroutine  norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$
        #:endfor
     #:endfor
     #:endfor
  end interface get_norm

  !> Matrix norms: function interface
  interface mnorm
     !! version: experimental 
     !!
     !! Computes the matrix norm of a generic-rank array \( A \). 
     !! ([Specification](../page/specs/stdlib_linalg.html#mnorm-computes-the-matrix-norm-of-a-generic-rank-array))
     !! 
     !!### Summary 
     !! Return one of several matrix norm metrics of a `real` or `complex` input array \( A \), 
     !! that can have rank 2 or higher. For rank-2 arrays, the matrix norm is returned.
     !! If rank>2 and the optional input dimensions `dim` are specified, 
     !! a rank `n-2` array is returned with dimensions `dim(1),dim(2)` collapsed, containing all 
     !! matrix norms evaluated over the specified dimensions only. `dim==[1,2]` are assumed as default
     !! dimensions if not specified.
     !! 
     !!### Description
     !! 
     !! This interface provides methods for computing the matrix norm(s) of an array.  
     !! Supported data types include `real` and `complex`. 
     !! Input arrays must have rank >= 2.
     !!
     !! Norm type input is optional, and it is provided via the `order` argument. 
     !! This can be provided as either an `integer` value or a `character` string. 
     !! Allowed metrics are: 
     !! - 1-norm: `order` = 1 or '1'    
     !! - 2-norm: `order` = 2 or '2'
     !! - Euclidean/Frobenius: `order` = 'Euclidean','Frobenius', or argument not specified
     !! - Infinity norm: `order` = huge(0) or 'Inf'
     !! 
     !! If an invalid norm type is provided, the routine returns an error state.
     !!
     !!### Example
     !!
     !!```fortran
     !!    real(sp) :: a(3,3), na
     !!    real(sp) :: b(3,3,4), nb(4)  ! Array of 4 3x3 matrices
     !!    a = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
     !!    
     !!    ! Euclidean/Frobenius norm of single matrix
     !!    na = mnorm(a)
     !!    na = mnorm(a, 'Euclidean')
     !!   
     !!    ! 1-norm of each 3x3 matrix in b
     !!    nb = mnorm(b, 1, dim=[1,2])
     !!     
     !!    ! Infinity-norm 
     !!    na = mnorm(b, 'inf', dim=[3,2])
     !!```     
     !!
      #:for rk,rt,ri in RC_KINDS_TYPES
      #:for it,ii in NORM_INPUT_OPTIONS   
      
      !> Matrix norms: ${rt}$ rank-2 arrays
      module function matrix_norm_${ii}$_${ri}$(a, order, err) result(nrm)
        !> Input matrix a(m,n)
        ${rt}$, intent(in), target :: a(:,:)
        !> Norm of the matrix.        
        real(${rk}$) :: nrm
        !> Order of the matrix norm being computed.
        ${it}$, #{if 'integer' in it}#optional, #{endif}#intent(in) :: order
        !> [optional] state return flag. On error if not requested, the code will stop
        type(linalg_state_type), intent(out), optional :: err      
      end function matrix_norm_${ii}$_${ri}$
      
      !> Matrix norms: ${rt}$ higher rank arrays
      #:for rank in range(3, MAXRANK + 1)
      module function matrix_norm_${rank}$D_to_${rank-2}$D_${ii}$_${ri}$(a, order, dim, err) result(nrm)
          !> Input matrix a(m,n)
          ${rt}$, intent(in), contiguous, target :: a${ranksuffix(rank)}$
          !> Norm of the matrix.        
          real(${rk}$), allocatable :: nrm${ranksuffix(rank-2)}$
          !> Order of the matrix norm being computed.
          ${it}$, #{if 'integer' in it}#optional, #{endif}#intent(in) :: order
          !> [optional] dimensions of the sub-matrices the norms should be evaluated at (default = [1,2])
          integer(ilp), optional, intent(in) :: dim(2)
          !> [optional] state return flag. On error if not requested, the code will stop
          type(linalg_state_type), intent(out), optional :: err        
      end function matrix_norm_${rank}$D_to_${rank-2}$D_${ii}$_${ri}$
      #:endfor            
      #:endfor
      #:endfor            
  end interface mnorm

  !> Matrix exponential: function interface
  interface expm
    !! version : experimental
    !!
    !! Computes the exponential of a matrix using a rational Pade approximation.
    !! ([Specification](../page/specs/stdlib_linalg.html#expm))
    !!
    !! ### Description
    !!
    !! This interface provides methods for computing the exponential of a matrix
    !! represented as a standard Fortran rank-2 array. Supported data types include
    !! `real` and `complex`.
    !!
    !! By default, the order of the Pade approximation is set to 10. It can be changed
    !! via the `order` argument which must be non-negative.
    !!
    !! If the input matrix is non-square or the order of the Pade approximation is
    !! negative, the function returns an error state.
    !!
    !! ### Example
    !!
    !! ```fortran
    !!  real(dp) :: A(3, 3), E(3, 3)
    !!
    !!  A = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
    !!
    !!  ! Default Pade approximation of the matrix exponential.
    !!  E = expm(A)
    !!
    !!  ! Pade approximation with specified order.
    !!  E = expm(A, order=12)
    !! ```
    !!
    #:for rk,rt,ri in RC_KINDS_TYPES
    module function stdlib_expm_${ri}$(A, order, err) result(E)
        !> Input matrix a(n, n).
        ${rt}$, intent(in) :: A(:, :)
        !> [optional] Order of the Pade approximation (default `order=10`)
        integer(ilp), optional, intent(in) :: order
        !> [optional] State return flag. On error, if not requested, the code will stop.
        type(linalg_state_type), optional, intent(out) :: err
        !> Exponential of the input matrix E = exp(A).
        ${rt}$, allocatable :: E(:, :)
    end function stdlib_expm_${ri}$
    #:endfor
  end interface expm

contains


    !> Version: experimental
    !>
    !> Constructs the identity matrix.
    !> ([Specification](../page/specs/stdlib_linalg.html#eye-construct-the-identity-matrix))
    #:for k1, t1 in RCI_KINDS_TYPES
    pure function eye_${t1[0]}$${k1}$(dim1, dim2, mold) result(result)

        integer, intent(in) :: dim1
        integer, intent(in), optional :: dim2
        ${t1}$, intent(in) #{if t1 == 'real(dp)'}#, optional #{endif}#:: mold        
        ${t1}$, allocatable :: result(:, :)

        integer :: dim2_
        integer :: i

        dim2_ = optval(dim2, dim1)
        allocate(result(dim1, dim2_))
        
        result = 0
        do i = 1, min(dim1, dim2_)
            result(i, i) = 1
        end do

    end function eye_${t1[0]}$${k1}$
    #:endfor

    #:for k1, t1 in RCI_KINDS_TYPES
      function trace_${t1[0]}$${k1}$(A) result(res)
        ${t1}$, intent(in) :: A(:,:)
        ${t1}$ :: res
        integer :: i
        res = 0
        do i = 1, minval(shape(A))
          res = res + A(i,i)
        end do
      end function trace_${t1[0]}$${k1}$
    #:endfor


    #:for k1, t1 in RCI_KINDS_TYPES
      pure function is_square_${t1[0]}$${k1}$(A) result(res)
        ${t1}$, intent(in) :: A(:,:)
        logical :: res
        res = (size(A,1) == size(A,2))
      end function is_square_${t1[0]}$${k1}$
    #:endfor


    #:for k1, t1 in RCI_KINDS_TYPES
      pure function is_diagonal_${t1[0]}$${k1}$(A) result(res)
        ${t1}$, intent(in) :: A(:,:)
        logical :: res
        ${t1}$, parameter :: zero = 0 !zero of relevant type
        integer :: m, n, o, i, j
        m = size(A,1)
        n = size(A,2)
        do j = 1, n !loop over all columns
            o = min(j-1,m) !index of row above diagonal (or last row)
            do i = 1, o !loop over rows above diagonal
                if (A(i,j) /= zero) then
                  res = .false.
                  return
                end if
            end do
            do i = o+2, m !loop over rows below diagonal
                if (A(i,j) /= zero) then
                  res = .false.
                  return
                end if
            end do
        end do
        res = .true. !otherwise A is diagonal
      end function is_diagonal_${t1[0]}$${k1}$
    #:endfor


    #:for k1, t1 in RCI_KINDS_TYPES
      pure function is_symmetric_${t1[0]}$${k1}$(A) result(res)
        ${t1}$, intent(in) :: A(:,:)
        logical :: res
        integer :: n, i, j
        if (.not. is_square(A)) then
           res = .false.
           return !nonsquare matrices cannot be symmetric
        end if
        n = size(A,1) !symmetric dimension of A
        do j = 1, n !loop over all columns
            do i = 1, j-1 !loop over all rows above diagonal
                if (A(i,j) /= A(j,i)) then
                  res = .false.
                  return
                end if
            end do
        end do
        res = .true. !otherwise A is symmetric
      end function is_symmetric_${t1[0]}$${k1}$
    #:endfor


    #:for k1, t1 in RCI_KINDS_TYPES
      pure function is_skew_symmetric_${t1[0]}$${k1}$(A) result(res)
        ${t1}$, intent(in) :: A(:,:)
        logical :: res
        integer :: n, i, j
        if (.not. is_square(A)) then
           res = .false.
           return !nonsquare matrices cannot be skew-symmetric
        end if
        n = size(A,1) !symmetric dimension of A
        do j = 1, n !loop over all columns
            do i = 1, j !loop over all rows above diagonal (and diagonal)
                if (A(i,j) /= -A(j,i)) then
                  res = .false.
                  return
                end if
            end do
        end do
        res = .true. !otherwise A is skew-symmetric
      end function is_skew_symmetric_${t1[0]}$${k1}$
    #:endfor


    #:for k1, t1 in (REAL_KINDS_TYPES + INT_KINDS_TYPES)
      pure function is_hermitian_${t1[0]}$${k1}$(A) result(res)
        ${t1}$, intent(in) :: A(:,:)
        logical :: res
        res = is_symmetric(A) !symmetry and Hermiticity are equivalent for real matrices
      end function is_hermitian_${t1[0]}$${k1}$
    #:endfor
    #:for k1, t1 in CMPLX_KINDS_TYPES
      pure function is_hermitian_${t1[0]}$${k1}$(A) result(res)
        ${t1}$, intent(in) :: A(:,:)
        logical :: res
        integer :: n, i, j
        if (.not. is_square(A)) then
           res = .false.
           return !nonsquare matrices cannot be Hermitian
        end if
        n = size(A,1) !symmetric dimension of A
        do j = 1, n !loop over all columns
            do i = 1, j !loop over all rows above diagonal (and diagonal)
                if (A(i,j) /= conjg(A(j,i))) then
                  res = .false.
                  return
                end if
            end do
        end do
        res = .true. !otherwise A is Hermitian
      end function is_hermitian_${t1[0]}$${k1}$
    #:endfor

    #:for k1, t1 in RCI_KINDS_TYPES
      pure module function hermitian_${t1[0]}$${k1}$(a) result(ah)
        ${t1}$, intent(in) :: a(:,:)
        ${t1}$ :: ah(size(a, 2), size(a, 1))
        #:if t1.startswith('complex')
        ah = conjg(transpose(a))
        #:else
        ah = transpose(a)
        #:endif
      end function hermitian_${t1[0]}$${k1}$
    #:endfor

    #:for k1, t1 in RCI_KINDS_TYPES
      function is_triangular_${t1[0]}$${k1}$(A,uplo) result(res)
        ${t1}$, intent(in) :: A(:,:)
        character, intent(in) :: uplo
        logical :: res
        ${t1}$, parameter :: zero = 0 !zero of relevant type
        integer :: m, n, o, i, j
        m = size(A,1)
        n = size(A,2)
        if ((uplo == 'u') .or. (uplo == 'U')) then !check for upper triangularity
          do j = 1, n !loop over all columns
              o = min(j-1,m) !index of row above diagonal (or last row)
              do i = o+2, m !loop over rows below diagonal
                  if (A(i,j) /= zero) then
                    res = .false.
                    return
                  end if
              end do
          end do
        else if ((uplo == 'l') .or. (uplo == 'L')) then !check for lower triangularity
          do j=1,n !loop over all columns
              o = min(j-1,m) !index of row above diagonal (or last row)
              do i=1,o !loop over rows above diagonal
                  if (A(i,j) /= zero) then
                    res = .false.
                    return
                  end if
              end do
           end do
        else
           call error_stop("ERROR (is_triangular): second argument must be one of {'u','U','l','L'}")
        end if
     
        res = .true. !otherwise A is triangular of the requested type
      end function is_triangular_${t1[0]}$${k1}$
    #:endfor


    #:for k1, t1 in RCI_KINDS_TYPES
      function is_hessenberg_${t1[0]}$${k1}$(A,uplo) result(res)
        ${t1}$, intent(in) :: A(:,:)
        character, intent(in) :: uplo
        logical :: res
        ${t1}$, parameter :: zero = 0 !zero of relevant type
        integer :: m, n, o, i, j
        m = size(A,1)
        n = size(A,2)
        if ((uplo == 'u') .or. (uplo == 'U')) then !check for upper Hessenberg
          do j = 1, n !loop over all columns
              o = min(j-2,m) !index of row two above diagonal (or last row)
              do i = o+4, m !loop over rows two or more below main diagonal
                  if (A(i,j) /= zero) then
                    res = .false.
                    return
                  end if
              end do
          end do
        else if ((uplo == 'l') .or. (uplo == 'L')) then !check for lower Hessenberg
          do j = 1, n !loop over all columns
              o = min(j-2,m) !index of row two above diagonal (or last row)
              do i = 1, o !loop over rows one or more above main diagonal
                  if (A(i,j) /= zero) then
                    res = .false.
                    return
                  end if
              end do
           end do
        else
           call error_stop("ERROR (is_hessenberg): second argument must be one of {'u','U','l','L'}")
        end if
        res = .true. !otherwise A is Hessenberg of the requested type
      end function is_hessenberg_${t1[0]}$${k1}$
    #:endfor
    
end module stdlib_linalg