linalg.f90

! linalg.f90
 
 
 
module linalg
    use, intrinsic :: iso_fortran_env, only : int32, real64
    use ferror, only : errors
    implicit none
 
    private
    public :: mtx_mult
    public :: rank1_update
    public :: diag_mtx_mult
    public :: trace
    public :: mtx_rank
    public :: det
    public :: swap
    public :: recip_mult_array
    public :: tri_mtx_mult
    public :: lu_factor
    public :: form_lu
    public :: qr_factor
    public :: form_qr
    public :: mult_qr
    public :: qr_rank1_update
    public :: cholesky_factor
    public :: cholesky_rank1_update
    public :: cholesky_rank1_downdate
    public :: rz_factor
    public :: mult_rz
    public :: svd
    public :: solve_triangular_system
    public :: solve_lu
    public :: solve_qr
    public :: solve_cholesky
    public :: mtx_inverse
    public :: mtx_pinverse
    public :: solve_least_squares
    public :: solve_least_squares_full
    public :: solve_least_squares_svd
    public :: eigen
    public :: sort
    public :: lq_factor
    public :: form_lq
    public :: mult_lq
    public :: solve_lq
    public :: band_mtx_mult
    public :: band_mtx_to_full_mtx
    public :: band_diag_mtx_mult
    public :: banded_to_dense
    public :: dense_to_banded
    public :: extract_diagonal
    public :: csr_matrix
    public :: msr_matrix
    public :: size
    public :: create_empty_csr_matrix
    public :: create_empty_msr_matrix
    public :: nonzero_count
    public :: dense_to_csr
    public :: diag_to_csr
    public :: banded_to_csr
    public :: csr_to_dense
    public :: csr_to_msr
    public :: msr_to_csr
    public :: dense_to_msr
    public :: msr_to_dense
    public :: create_csr_matrix
    public :: matmul
    public :: operator(+)
    public :: operator(-)
    public :: operator(*)
    public :: operator(/)
    public :: assignment(=)
    public :: transpose
    public :: sparse_direct_solve
    public :: pgmres_solver
    public :: la_no_operation
    public :: la_transpose
    public :: la_hermitian_transpose
    public :: la_no_error
    public :: la_invalid_input_error
    public :: la_array_size_error
    public :: la_singular_matrix_error
    public :: la_matrix_format_error
    public :: la_out_of_memory_error
    public :: la_convergence_error
    public :: la_invalid_operation_error
 
! ******************************************************************************
! TYPES
! ------------------------------------------------------------------------------
    type :: csr_matrix
        integer(int32), allocatable, dimension(:) :: row_indices
        integer(int32), allocatable, dimension(:) :: column_indices
        real(real64), allocatable, dimension(:) :: values
        integer(int32), private :: n = 0
    contains
        procedure, public :: get => csr_get_element
    end type
 
! ------------------------------------------------------------------------------
    type :: msr_matrix
        integer(int32), allocatable, dimension(:) :: indices
        real(real64), allocatable, dimension(:) :: values
        integer(int32), private :: m = 0
        integer(int32), private :: n = 0
        integer(int32), private :: nnz = 0
    end type
 
! ******************************************************************************
! CONSTANTS
! ------------------------------------------------------------------------------
    integer(int32), parameter :: la_no_operation = 0
    integer(int32), parameter :: la_transpose = 1
    integer(int32), parameter :: la_hermitian_transpose = 2
 
! ******************************************************************************
! ERROR FLAGS
! ------------------------------------------------------------------------------
    integer(int32), parameter :: la_no_error = 0
    integer(int32), parameter :: la_invalid_input_error = 101
    integer(int32), parameter :: la_array_size_error = 102
    integer(int32), parameter :: la_singular_matrix_error = 103
    integer(int32), parameter :: la_matrix_format_error = 104
    integer(int32), parameter :: la_out_of_memory_error = 105
    integer(int32), parameter :: la_convergence_error = 106
    integer(int32), parameter :: la_invalid_operation_error = 107
 
! ******************************************************************************
! INTERFACES
! ------------------------------------------------------------------------------
interface mtx_mult
    module procedure :: mtx_mult_mtx
    module procedure :: mtx_mult_vec
    module procedure :: cmtx_mult_mtx
    module procedure :: cmtx_mult_vec
end interface
 
! ------------------------------------------------------------------------------
interface rank1_update
    module procedure :: rank1_update_dbl
    module procedure :: rank1_update_cmplx
end interface
 
! ------------------------------------------------------------------------------
interface diag_mtx_mult
    module procedure :: diag_mtx_mult_mtx
    module procedure :: diag_mtx_mult_mtx2
    module procedure :: diag_mtx_mult_mtx3
    module procedure :: diag_mtx_mult_mtx4
    module procedure :: diag_mtx_mult_mtx_cmplx
    module procedure :: diag_mtx_mult_mtx2_cmplx
    module procedure :: diag_mtx_mult_mtx_mix
    module procedure :: diag_mtx_mult_mtx2_mix
    module procedure :: diag_mtx_sparse_mult
end interface
 
! ------------------------------------------------------------------------------
interface trace
    module procedure :: trace_dbl
    module procedure :: trace_cmplx
end interface
 
! ------------------------------------------------------------------------------
interface mtx_rank
    module procedure :: mtx_rank_dbl
    module procedure :: mtx_rank_cmplx
end interface
 
! ------------------------------------------------------------------------------
interface det
    module procedure :: det_dbl
    module procedure :: det_cmplx
end interface
 
! ------------------------------------------------------------------------------
interface swap
    module procedure :: swap_dbl
    module procedure :: swap_cmplx
end interface
 
! ------------------------------------------------------------------------------
interface recip_mult_array
    module procedure :: recip_mult_array_dbl
end interface
 
! ------------------------------------------------------------------------------
interface tri_mtx_mult
    module procedure :: tri_mtx_mult_dbl
    module procedure :: tri_mtx_mult_cmplx
end interface
 
! ------------------------------------------------------------------------------
interface lu_factor
    module procedure :: lu_factor_dbl
    module procedure :: lu_factor_cmplx
    module procedure :: csr_lu_factor
end interface
 
interface form_lu
    module procedure :: form_lu_all
    module procedure :: form_lu_all_cmplx
    module procedure :: form_lu_only
    module procedure :: form_lu_only_cmplx
end interface
 
! ------------------------------------------------------------------------------
interface qr_factor
    module procedure :: qr_factor_no_pivot
    module procedure :: qr_factor_no_pivot_cmplx
    module procedure :: qr_factor_pivot
    module procedure :: qr_factor_pivot_cmplx
end interface
 
! ------------------------------------------------------------------------------
interface form_qr
    module procedure :: form_qr_no_pivot
    module procedure :: form_qr_no_pivot_cmplx
    module procedure :: form_qr_pivot
    module procedure :: form_qr_pivot_cmplx
end interface
 
! ------------------------------------------------------------------------------
interface mult_qr
    module procedure :: mult_qr_mtx
    module procedure :: mult_qr_mtx_cmplx
    module procedure :: mult_qr_vec
    module procedure :: mult_qr_vec_cmplx
end interface
 
! ------------------------------------------------------------------------------
interface qr_rank1_update
    module procedure :: qr_rank1_update_dbl
    module procedure :: qr_rank1_update_cmplx
end interface
 
! ------------------------------------------------------------------------------
interface cholesky_factor
    module procedure :: cholesky_factor_dbl
    module procedure :: cholesky_factor_cmplx
end interface
 
! ------------------------------------------------------------------------------
interface cholesky_rank1_update
    module procedure :: cholesky_rank1_update_dbl
    module procedure :: cholesky_rank1_update_cmplx
end interface
 
! ------------------------------------------------------------------------------
interface cholesky_rank1_downdate
    module procedure :: cholesky_rank1_downdate_dbl
    module procedure :: cholesky_rank1_downdate_cmplx
end interface
 
! ------------------------------------------------------------------------------
interface rz_factor
    module procedure :: rz_factor_dbl
    module procedure :: rz_factor_cmplx
end interface
 
! ------------------------------------------------------------------------------
interface mult_rz
    module procedure :: mult_rz_mtx
    module procedure :: mult_rz_mtx_cmplx
    module procedure :: mult_rz_vec
    module procedure :: mult_rz_vec_cmplx
end interface
 
! ------------------------------------------------------------------------------
interface svd
    module procedure :: svd_dbl
    module procedure :: svd_cmplx
end interface
 
! ------------------------------------------------------------------------------
interface solve_triangular_system
    module procedure :: solve_tri_mtx
    module procedure :: solve_tri_mtx_cmplx
    module procedure :: solve_tri_vec
    module procedure :: solve_tri_vec_cmplx
end interface
 
! ------------------------------------------------------------------------------
interface solve_lu
    module procedure :: solve_lu_mtx
    module procedure :: solve_lu_mtx_cmplx
    module procedure :: solve_lu_vec
    module procedure :: solve_lu_vec_cmplx
    module procedure :: csr_lu_solve
end interface
 
! ------------------------------------------------------------------------------
interface solve_qr
    module procedure :: solve_qr_no_pivot_mtx
    module procedure :: solve_qr_no_pivot_mtx_cmplx
    module procedure :: solve_qr_no_pivot_vec
    module procedure :: solve_qr_no_pivot_vec_cmplx
    module procedure :: solve_qr_pivot_mtx
    module procedure :: solve_qr_pivot_mtx_cmplx
    module procedure :: solve_qr_pivot_vec
    module procedure :: solve_qr_pivot_vec_cmplx
end interface
 
! ------------------------------------------------------------------------------
interface solve_cholesky
    module procedure :: solve_cholesky_mtx
    module procedure :: solve_cholesky_mtx_cmplx
    module procedure :: solve_cholesky_vec
    module procedure :: solve_cholesky_vec_cmplx
end interface
 
! ------------------------------------------------------------------------------
interface solve_least_squares
    module procedure :: solve_least_squares_mtx
    module procedure :: solve_least_squares_mtx_cmplx
    module procedure :: solve_least_squares_vec
    module procedure :: solve_least_squares_vec_cmplx
end interface
 
! ------------------------------------------------------------------------------
interface solve_least_squares_full
    module procedure :: solve_least_squares_mtx_pvt
    module procedure :: solve_least_squares_mtx_pvt_cmplx
    module procedure :: solve_least_squares_vec_pvt
    module procedure :: solve_least_squares_vec_pvt_cmplx
end interface
 
! ------------------------------------------------------------------------------
interface solve_least_squares_svd
    module procedure :: solve_least_squares_mtx_svd
    module procedure :: solve_least_squares_vec_svd
end interface
 
! ------------------------------------------------------------------------------
interface mtx_inverse
    module procedure :: mtx_inverse_dbl
    module procedure :: mtx_inverse_cmplx
end interface
 
! ------------------------------------------------------------------------------
interface mtx_pinverse
    module procedure :: mtx_pinverse_dbl
    module procedure :: mtx_pinverse_cmplx
end interface
 
! ------------------------------------------------------------------------------
interface eigen
    module procedure :: eigen_symm
    module procedure :: eigen_asymm
    module procedure :: eigen_gen
    module procedure :: eigen_cmplx
end interface
 
! ------------------------------------------------------------------------------
interface sort
    module procedure :: sort_dbl_array
    module procedure :: sort_dbl_array_ind
    module procedure :: sort_cmplx_array
    module procedure :: sort_cmplx_array_ind
    module procedure :: sort_eigen_cmplx
    module procedure :: sort_eigen_dbl
    module procedure :: sort_int32_array
    module procedure :: sort_int32_array_ind
end interface
 
interface lq_factor
    module procedure :: lq_factor_no_pivot
    module procedure :: lq_factor_no_pivot_cmplx
end interface
 
interface form_lq
    module procedure :: form_lq_no_pivot
    module procedure :: form_lq_no_pivot_cmplx
end interface
 
interface mult_lq
    module procedure :: mult_lq_mtx
    module procedure :: mult_lq_mtx_cmplx
    module procedure :: mult_lq_vec
    module procedure :: mult_lq_vec_cmplx
end interface
 
interface solve_lq
    module procedure :: solve_lq_mtx
    module procedure :: solve_lq_mtx_cmplx
    module procedure :: solve_lq_vec
    module procedure :: solve_lq_vec_cmplx
end interface
 
! ------------------------------------------------------------------------------
interface band_mtx_mult
    module procedure :: band_mtx_vec_mult_dbl
    module procedure :: band_mtx_vec_mult_cmplx
end interface
 
interface band_mtx_to_full_mtx
    module procedure :: band_to_full_mtx_dbl
    module procedure :: band_to_full_mtx_cmplx
end interface
 
interface band_diag_mtx_mult
    module procedure :: band_diag_mtx_mult_dbl
    module procedure :: band_diag_mtx_mult_cmplx
end interface
 
interface banded_to_dense
    module procedure :: banded_to_dense_dbl
    module procedure :: banded_to_dense_cmplx
end interface
 
interface dense_to_banded
    module procedure :: dense_to_banded_dbl
    module procedure :: dense_to_banded_cmplx
end interface
 
interface extract_diagonal
    module procedure :: extract_diagonal_dbl
    module procedure :: extract_diagonal_cmplx
    module procedure :: extract_diagonal_csr
end interface
 
! ******************************************************************************
! LINALG_BASIC.F90
! ------------------------------------------------------------------------------
interface
    module subroutine mtx_mult_mtx(transa, transb, alpha, a, b, beta, c, err)
        logical, intent(in) :: transa, transb
        real(real64), intent(in) :: alpha, beta
        real(real64), intent(in), dimension(:,:) :: a, b
        real(real64), intent(inout), dimension(:,:) :: c
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine mtx_mult_vec(trans, alpha, a, b, beta, c, err)
        logical, intent(in) :: trans
        real(real64), intent(in) :: alpha, beta
        real(real64), intent(in), dimension(:,:) :: a
        real(real64), intent(in), dimension(:) :: b
        real(real64), intent(inout), dimension(:) :: c
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine cmtx_mult_mtx(opa, opb, alpha, a, b, beta, c, err)
        integer(int32), intent(in) :: opa, opb
        complex(real64), intent(in) :: alpha, beta
        complex(real64), intent(in), dimension(:,:) :: a, b
        complex(real64), intent(inout), dimension(:,:) :: c
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine cmtx_mult_vec(opa, alpha, a, b, beta, c, err)
        integer(int32), intent(in) :: opa
        complex(real64), intent(in) :: alpha, beta
        complex(real64), intent(in), dimension(:,:) :: a
        complex(real64), intent(in), dimension(:) :: b
        complex(real64), intent(inout), dimension(:) :: c
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine rank1_update_dbl(alpha, x, y, a, err)
        real(real64), intent(in) :: alpha
        real(real64), intent(in), dimension(:) :: x, y
        real(real64), intent(inout), dimension(:,:) :: a
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine rank1_update_cmplx(alpha, x, y, a, err)
        complex(real64), intent(in) :: alpha
        complex(real64), intent(in), dimension(:) :: x, y
        complex(real64), intent(inout), dimension(:,:) :: a
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine diag_mtx_mult_mtx(lside, trans, alpha, a, b, beta, c, err)
        logical, intent(in) :: lside, trans
        real(real64) :: alpha, beta
        real(real64), intent(in), dimension(:) :: a
        real(real64), intent(in), dimension(:,:) :: b
        real(real64), intent(inout), dimension(:,:) :: c
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine diag_mtx_mult_mtx2(lside, alpha, a, b, err)
        logical, intent(in) :: lside
        real(real64), intent(in) :: alpha
        real(real64), intent(in), dimension(:) :: a
        real(real64), intent(inout), dimension(:,:) :: b
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine diag_mtx_mult_mtx3(lside, trans, alpha, a, b, beta, c, err)
        logical, intent(in) :: lside, trans
        real(real64) :: alpha, beta
        complex(real64), intent(in), dimension(:) :: a
        real(real64), intent(in), dimension(:,:) :: b
        complex(real64), intent(inout), dimension(:,:) :: c
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine diag_mtx_mult_mtx4(lside, opb, alpha, a, b, beta, c, err)
        logical, intent(in) :: lside
        integer(int32), intent(in) :: opb
        real(real64) :: alpha, beta
        complex(real64), intent(in), dimension(:) :: a
        complex(real64), intent(in), dimension(:,:) :: b
        complex(real64), intent(inout), dimension(:,:) :: c
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine diag_mtx_mult_mtx_cmplx(lside, opb, alpha, a, b, beta, c, err)
        logical, intent(in) :: lside
        integer(int32), intent(in) :: opb
        complex(real64) :: alpha, beta
        complex(real64), intent(in), dimension(:) :: a
        complex(real64), intent(in), dimension(:,:) :: b
        complex(real64), intent(inout), dimension(:,:) :: c
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine diag_mtx_mult_mtx2_cmplx(lside, alpha, a, b, err)
        logical, intent(in) :: lside
        complex(real64), intent(in) :: alpha
        complex(real64), intent(in), dimension(:) :: a
        complex(real64), intent(inout), dimension(:,:) :: b
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine diag_mtx_mult_mtx_mix(lside, opb, alpha, a, b, beta, c, err)
        logical, intent(in) :: lside
        integer(int32), intent(in) :: opb
        complex(real64) :: alpha, beta
        real(real64), intent(in), dimension(:) :: a
        complex(real64), intent(in), dimension(:,:) :: b
        complex(real64), intent(inout), dimension(:,:) :: c
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine diag_mtx_mult_mtx2_mix(lside, alpha, a, b, err)
        logical, intent(in) :: lside
        complex(real64), intent(in) :: alpha
        real(real64), intent(in), dimension(:) :: a
        complex(real64), intent(inout), dimension(:,:) :: b
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine diag_mtx_sparse_mult(lside, alpha, a, b, err)
        logical, intent(in) :: lside
        real(real64), intent(in) :: alpha
        real(real64), intent(in), dimension(:) :: a
        type(csr_matrix), intent(inout) :: b
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    pure module function trace_dbl(x) result(y)
        real(real64), intent(in), dimension(:,:) :: x
        real(real64) :: y
    end function
    
    pure module function trace_cmplx(x) result(y)
        complex(real64), intent(in), dimension(:,:) :: x
        complex(real64) :: y
    end function
    
    module function mtx_rank_dbl(a, tol, work, olwork, err) result(rnk)
        real(real64), intent(inout), dimension(:,:) :: a
        real(real64), intent(in), optional :: tol
        real(real64), intent(out), target, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
        integer(int32) :: rnk
    end function
 
    module function mtx_rank_cmplx(a, tol, work, olwork, rwork, err) result(rnk)
        complex(real64), intent(inout), dimension(:,:) :: a
        real(real64), intent(in), optional :: tol
        complex(real64), intent(out), target, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        real(real64), intent(out), target, optional, dimension(:) :: rwork
        class(errors), intent(inout), optional, target :: err
        integer(int32) :: rnk
    end function
    
    module function det_dbl(a, iwork, err) result(x)
        real(real64), intent(inout), dimension(:,:) :: a
        integer(int32), intent(out), target, optional, dimension(:) :: iwork
        class(errors), intent(inout), optional, target :: err
        real(real64) :: x
    end function
    
    module function det_cmplx(a, iwork, err) result(x)
        complex(real64), intent(inout), dimension(:,:) :: a
        integer(int32), intent(out), target, optional, dimension(:) :: iwork
        class(errors), intent(inout), optional, target :: err
        complex(real64) :: x
    end function
    
    module subroutine swap_dbl(x, y, err)
        real(real64), intent(inout), dimension(:) :: x, y
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine swap_cmplx(x, y, err)
        complex(real64), intent(inout), dimension(:) :: x, y
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine recip_mult_array_dbl(a, x)
        real(real64), intent(in) :: a
        real(real64), intent(inout), dimension(:) :: x
    end subroutine
    
    module subroutine tri_mtx_mult_dbl(upper, alpha, a, beta, b, err)
        logical, intent(in) :: upper
        real(real64), intent(in) :: alpha, beta
        real(real64), intent(in), dimension(:,:) :: a
        real(real64), intent(inout), dimension(:,:) :: b
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine tri_mtx_mult_cmplx(upper, alpha, a, beta, b, err)
        logical, intent(in) :: upper
        complex(real64), intent(in) :: alpha, beta
        complex(real64), intent(in), dimension(:,:) :: a
        complex(real64), intent(inout), dimension(:,:) :: b
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine band_mtx_vec_mult_dbl(trans, kl, ku, alpha, a, x, beta, &
        y, err)
        logical, intent(in) :: trans
        integer(int32), intent(in) :: kl, ku
        real(real64), intent(in) :: alpha, beta
        real(real64), intent(in), dimension(:,:) :: a
        real(real64), intent(in), dimension(:) :: x
        real(real64), intent(inout), dimension(:) :: y
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine band_mtx_vec_mult_cmplx(trans, kl, ku, alpha, a, x, &
        beta, y, err)
        integer(int32), intent(in) :: trans
        integer(int32), intent(in) :: kl, ku
        complex(real64), intent(in) :: alpha, beta
        complex(real64), intent(in), dimension(:,:) :: a
        complex(real64), intent(in), dimension(:) :: x
        complex(real64), intent(inout), dimension(:) :: y
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine band_to_full_mtx_dbl(kl, ku, b, f, err)
        integer(int32), intent(in) :: kl, ku
        real(real64), intent(in), dimension(:,:) :: b
        real(real64), intent(out), dimension(:,:) :: f
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine band_to_full_mtx_cmplx(kl, ku, b, f, err)
        integer(int32), intent(in) :: kl, ku
        complex(real64), intent(in), dimension(:,:) :: b
        complex(real64), intent(out), dimension(:,:) :: f
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine band_diag_mtx_mult_dbl(left, m, kl, ku, alpha, a, b, err)
        logical, intent(in) :: left
        integer(int32), intent(in) :: m, kl, ku
        real(real64), intent(in) :: alpha
        real(real64), intent(inout), dimension(:,:) :: a
        real(real64), intent(in), dimension(:) :: b
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine band_diag_mtx_mult_cmplx(left, m, kl, ku, alpha, a, b, err)
        logical, intent(in) :: left
        integer(int32), intent(in) :: m, kl, ku
        complex(real64), intent(in) :: alpha
        complex(real64), intent(inout), dimension(:,:) :: a
        complex(real64), intent(in), dimension(:) :: b
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine banded_to_dense_dbl(m, kl, ku, a, x, err)
        integer(int32), intent(in) :: m, kl, ku
        real(real64), intent(in), dimension(:,:) :: a
        real(real64), intent(out), dimension(:,:) :: x
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine banded_to_dense_cmplx(m, kl, ku, a, x, err)
        integer(int32), intent(in) :: m, kl, ku
        complex(real64), intent(in), dimension(:,:) :: a
        complex(real64), intent(out), dimension(:,:) :: x
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine dense_to_banded_dbl(a, kl, ku, x, err)
        real(real64), intent(in), dimension(:,:) :: a
        integer(int32), intent(in) :: kl, ku
        real(real64), intent(out), dimension(:,:) :: x
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine dense_to_banded_cmplx(a, kl, ku, x, err)
        complex(real64), intent(in), dimension(:,:) :: a
        integer(int32), intent(in) :: kl, ku
        complex(real64), intent(out), dimension(:,:) :: x
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine extract_diagonal_dbl(a, diag, err)
        ! Arguments
        real(real64), intent(in), dimension(:,:) :: a
        real(real64), intent(out), dimension(:) :: diag
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine extract_diagonal_cmplx(a, diag, err)
        ! Arguments
        complex(real64), intent(in), dimension(:,:) :: a
        complex(real64), intent(out), dimension(:) :: diag
        class(errors), intent(inout), optional, target :: err
    end subroutine
end interface
 
! ******************************************************************************
! LINALG_FACTOR.F90
! ------------------------------------------------------------------------------
interface
    module subroutine lu_factor_dbl(a, ipvt, err)
        real(real64), intent(inout), dimension(:,:) :: a
        integer(int32), intent(out), dimension(:) :: ipvt
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine lu_factor_cmplx(a, ipvt, err)
        complex(real64), intent(inout), dimension(:,:) :: a
        integer(int32), intent(out), dimension(:) :: ipvt
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine form_lu_all(lu, ipvt, u, p, err)
        real(real64), intent(inout), dimension(:,:) :: lu
        integer(int32), intent(in), dimension(:) :: ipvt
        real(real64), intent(out), dimension(:,:) :: u, p
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine form_lu_all_cmplx(lu, ipvt, u, p, err)
        complex(real64), intent(inout), dimension(:,:) :: lu
        integer(int32), intent(in), dimension(:) :: ipvt
        complex(real64), intent(out), dimension(:,:) :: u
        real(real64), intent(out), dimension(:,:) :: p
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine form_lu_only(lu, u, err)
        real(real64), intent(inout), dimension(:,:) :: lu
        real(real64), intent(out), dimension(:,:) :: u
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine form_lu_only_cmplx(lu, u, err)
        complex(real64), intent(inout), dimension(:,:) :: lu
        complex(real64), intent(out), dimension(:,:) :: u
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine qr_factor_no_pivot(a, tau, work, olwork, err)
        real(real64), intent(inout), dimension(:,:) :: a
        real(real64), intent(out), dimension(:) :: tau
        real(real64), intent(out), target, dimension(:), optional :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine qr_factor_no_pivot_cmplx(a, tau, work, olwork, err)
        complex(real64), intent(inout), dimension(:,:) :: a
        complex(real64), intent(out), dimension(:) :: tau
        complex(real64), intent(out), target, dimension(:), optional :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine qr_factor_pivot(a, tau, jpvt, work, olwork, err)
        real(real64), intent(inout), dimension(:,:) :: a
        real(real64), intent(out), dimension(:) :: tau
        integer(int32), intent(inout), dimension(:) :: jpvt
        real(real64), intent(out), target, dimension(:), optional :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine qr_factor_pivot_cmplx(a, tau, jpvt, work, olwork, rwork, &
            err)
        complex(real64), intent(inout), dimension(:,:) :: a
        complex(real64), intent(out), dimension(:) :: tau
        integer(int32), intent(inout), dimension(:) :: jpvt
        complex(real64), intent(out), target, dimension(:), optional :: work
        integer(int32), intent(out), optional :: olwork
        real(real64), intent(out), target, dimension(:), optional :: rwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine form_qr_no_pivot(r, tau, q, work, olwork, err)
        real(real64), intent(inout), dimension(:,:) :: r
        real(real64), intent(in), dimension(:) :: tau
        real(real64), intent(out), dimension(:,:) :: q
        real(real64), intent(out), target, dimension(:), optional :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine form_qr_no_pivot_cmplx(r, tau, q, work, olwork, err)
        complex(real64), intent(inout), dimension(:,:) :: r
        complex(real64), intent(in), dimension(:) :: tau
        complex(real64), intent(out), dimension(:,:) :: q
        complex(real64), intent(out), target, dimension(:), optional :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine form_qr_pivot(r, tau, pvt, q, p, work, olwork, err)
        real(real64), intent(inout), dimension(:,:) :: r
        real(real64), intent(in), dimension(:) :: tau
        integer(int32), intent(in), dimension(:) :: pvt
        real(real64), intent(out), dimension(:,:) :: q, p
        real(real64), intent(out), target, dimension(:), optional :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine form_qr_pivot_cmplx(r, tau, pvt, q, p, work, olwork, err)
        complex(real64), intent(inout), dimension(:,:) :: r
        complex(real64), intent(in), dimension(:) :: tau
        integer(int32), intent(in), dimension(:) :: pvt
        complex(real64), intent(out), dimension(:,:) :: q, p
        complex(real64), intent(out), target, dimension(:), optional :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine mult_qr_mtx(lside, trans, a, tau, c, work, olwork, err)
        logical, intent(in) :: lside, trans
        real(real64), intent(in), dimension(:) :: tau
        real(real64), intent(inout), dimension(:,:) :: a, c
        real(real64), intent(out), target, dimension(:), optional :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine mult_qr_mtx_cmplx(lside, trans, a, tau, c, work, olwork, err)
        logical, intent(in) :: lside, trans
        complex(real64), intent(in), dimension(:) :: tau
        complex(real64), intent(inout), dimension(:,:) :: a, c
        complex(real64), intent(out), target, dimension(:), optional :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine mult_qr_vec(trans, a, tau, c, work, olwork, err)
        logical, intent(in) :: trans
        real(real64), intent(inout), dimension(:,:) :: a
        real(real64), intent(in), dimension(:) :: tau
        real(real64), intent(inout), dimension(:) :: c
        real(real64), intent(out), target, dimension(:), optional :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine mult_qr_vec_cmplx(trans, a, tau, c, work, olwork, err)
        logical, intent(in) :: trans
        complex(real64), intent(inout), dimension(:,:) :: a
        complex(real64), intent(in), dimension(:) :: tau
        complex(real64), intent(inout), dimension(:) :: c
        complex(real64), intent(out), target, dimension(:), optional :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine qr_rank1_update_dbl(q, r, u, v, work, err)
        real(real64), intent(inout), dimension(:,:) :: q, r
        real(real64), intent(inout), dimension(:) :: u, v
        real(real64), intent(out), target, optional, dimension(:) :: work
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine qr_rank1_update_cmplx(q, r, u, v, work, rwork, err)
        complex(real64), intent(inout), dimension(:,:) :: q, r
        complex(real64), intent(inout), dimension(:) :: u, v
        complex(real64), intent(out), target, optional, dimension(:) :: work
        real(real64), intent(out), target, optional, dimension(:) :: rwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine cholesky_factor_dbl(a, upper, err)
        real(real64), intent(inout), dimension(:,:) :: a
        logical, intent(in), optional :: upper
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine cholesky_factor_cmplx(a, upper, err)
        complex(real64), intent(inout), dimension(:,:) :: a
        logical, intent(in), optional :: upper
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine cholesky_rank1_update_dbl(r, u, work, err)
        real(real64), intent(inout), dimension(:,:) :: r
        real(real64), intent(inout), dimension(:) :: u
        real(real64), intent(out), target, optional, dimension(:) :: work
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine cholesky_rank1_update_cmplx(r, u, work, err)
        complex(real64), intent(inout), dimension(:,:) :: r
        complex(real64), intent(inout), dimension(:) :: u
        real(real64), intent(out), target, optional, dimension(:) :: work
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine cholesky_rank1_downdate_dbl(r, u, work, err)
        real(real64), intent(inout), dimension(:,:) :: r
        real(real64), intent(inout), dimension(:) :: u
        real(real64), intent(out), target, optional, dimension(:) :: work
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine cholesky_rank1_downdate_cmplx(r, u, work, err)
        complex(real64), intent(inout), dimension(:,:) :: r
        complex(real64), intent(inout), dimension(:) :: u
        real(real64), intent(out), target, optional, dimension(:) :: work
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine rz_factor_dbl(a, tau, work, olwork, err)
        real(real64), intent(inout), dimension(:,:) :: a
        real(real64), intent(out), dimension(:) :: tau
        real(real64), intent(out), target, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine rz_factor_cmplx(a, tau, work, olwork, err)
        complex(real64), intent(inout), dimension(:,:) :: a
        complex(real64), intent(out), dimension(:) :: tau
        complex(real64), intent(out), target, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine mult_rz_mtx(lside, trans, l, a, tau, c, work, olwork, err)
        logical, intent(in) :: lside, trans
        integer(int32), intent(in) :: l
        real(real64), intent(inout), dimension(:,:) :: a, c
        real(real64), intent(in), dimension(:) :: tau
        real(real64), intent(out), target, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine mult_rz_mtx_cmplx(lside, trans, l, a, tau, c, work, olwork, err)
        logical, intent(in) :: lside, trans
        integer(int32), intent(in) :: l
        complex(real64), intent(inout), dimension(:,:) :: a, c
        complex(real64), intent(in), dimension(:) :: tau
        complex(real64), intent(out), target, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine mult_rz_vec(trans, l, a, tau, c, work, olwork, err)
        logical, intent(in) :: trans
        integer(int32), intent(in) :: l
        real(real64), intent(inout), dimension(:,:) :: a
        real(real64), intent(in), dimension(:) :: tau
        real(real64), intent(inout), dimension(:) :: c
        real(real64), intent(out), target, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine mult_rz_vec_cmplx(trans, l, a, tau, c, work, olwork, err)
        logical, intent(in) :: trans
        integer(int32), intent(in) :: l
        complex(real64), intent(inout), dimension(:,:) :: a
        complex(real64), intent(in), dimension(:) :: tau
        complex(real64), intent(inout), dimension(:) :: c
        complex(real64), intent(out), target, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine svd_dbl(a, s, u, vt, work, olwork, err)
        real(real64), intent(inout), dimension(:,:) :: a
        real(real64), intent(out), dimension(:) :: s
        real(real64), intent(out), optional, dimension(:,:) :: u, vt
        real(real64), intent(out), target, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine svd_cmplx(a, s, u, vt, work, olwork, rwork, err)
        complex(real64), intent(inout), dimension(:,:) :: a
        real(real64), intent(out), dimension(:) :: s
        complex(real64), intent(out), optional, dimension(:,:) :: u, vt
        complex(real64), intent(out), target, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        real(real64), intent(out), target, optional, dimension(:) :: rwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine lq_factor_no_pivot(a, tau, work, olwork, err)
        real(real64), intent(inout), dimension(:,:) :: a
        real(real64), intent(out), dimension(:) :: tau
        real(real64), intent(out), target, dimension(:), optional :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine lq_factor_no_pivot_cmplx(a, tau, work, olwork, err)
        complex(real64), intent(inout), dimension(:,:) :: a
        complex(real64), intent(out), dimension(:) :: tau
        complex(real64), intent(out), target, dimension(:), optional :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine form_lq_no_pivot(l, tau, q, work, olwork, err)
        real(real64), intent(inout), dimension(:,:) :: l
        real(real64), intent(in), dimension(:) :: tau
        real(real64), intent(out), dimension(:,:) :: q
        real(real64), intent(out), target, dimension(:), optional :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine form_lq_no_pivot_cmplx(l, tau, q, work, olwork, err)
        complex(real64), intent(inout), dimension(:,:) :: l
        complex(real64), intent(in), dimension(:) :: tau
        complex(real64), intent(out), dimension(:,:) :: q
        complex(real64), intent(out), target, dimension(:), optional :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine mult_lq_mtx(lside, trans, a, tau, c, work, olwork, err)
        logical, intent(in) :: lside, trans
        real(real64), intent(in), dimension(:,:) :: a
        real(real64), intent(in), dimension(:) :: tau
        real(real64), intent(inout), dimension(:,:) :: c
        real(real64), intent(out), target, dimension(:), optional :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine mult_lq_mtx_cmplx(lside, trans, a, tau, c, work, olwork, err)
        logical, intent(in) :: lside, trans
        complex(real64), intent(in), dimension(:,:) :: a
        complex(real64), intent(in), dimension(:) :: tau
        complex(real64), intent(inout), dimension(:,:) :: c
        complex(real64), intent(out), target, dimension(:), optional :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine mult_lq_vec(trans, a, tau, c, work, olwork, err)
        logical, intent(in) :: trans
        real(real64), intent(in), dimension(:,:) :: a
        real(real64), intent(in), dimension(:) :: tau
        real(real64), intent(inout), dimension(:) :: c
        real(real64), intent(out), target, dimension(:), optional :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine mult_lq_vec_cmplx(trans, a, tau, c, work, olwork, err)
        logical, intent(in) :: trans
        complex(real64), intent(in), dimension(:,:) :: a
        complex(real64), intent(in), dimension(:) :: tau
        complex(real64), intent(inout), dimension(:) :: c
        complex(real64), intent(out), target, dimension(:), optional :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
end interface
 
! ******************************************************************************
! LINALG_SOLVE.F90
! ------------------------------------------------------------------------------
interface
    module subroutine solve_tri_mtx(lside, upper, trans, nounit, alpha, a, b, err)
        logical, intent(in) :: lside, upper, trans, nounit
        real(real64), intent(in) :: alpha
        real(real64), intent(in), dimension(:,:) :: a
        real(real64), intent(inout), dimension(:,:) :: b
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine solve_tri_mtx_cmplx(lside, upper, trans, nounit, alpha, a, b, err)
        logical, intent(in) :: lside, upper, trans, nounit
        complex(real64), intent(in) :: alpha
        complex(real64), intent(in), dimension(:,:) :: a
        complex(real64), intent(inout), dimension(:,:) :: b
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine solve_tri_vec(upper, trans, nounit, a, x, err)
        logical, intent(in) :: upper, trans, nounit
        real(real64), intent(in), dimension(:,:) :: a
        real(real64), intent(inout), dimension(:) :: x
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine solve_tri_vec_cmplx(upper, trans, nounit, a, x, err)
        logical, intent(in) :: upper, trans, nounit
        complex(real64), intent(in), dimension(:,:) :: a
        complex(real64), intent(inout), dimension(:) :: x
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine solve_lu_mtx(a, ipvt, b, err)
        real(real64), intent(in), dimension(:,:) :: a
        integer(int32), intent(in), dimension(:) :: ipvt
        real(real64), intent(inout), dimension(:,:) :: b
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine solve_lu_mtx_cmplx(a, ipvt, b, err)
        complex(real64), intent(in), dimension(:,:) :: a
        integer(int32), intent(in), dimension(:) :: ipvt
        complex(real64), intent(inout), dimension(:,:) :: b
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine solve_lu_vec(a, ipvt, b, err)
        real(real64), intent(in), dimension(:,:) :: a
        integer(int32), intent(in), dimension(:) :: ipvt
        real(real64), intent(inout), dimension(:) :: b
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine solve_lu_vec_cmplx(a, ipvt, b, err)
        complex(real64), intent(in), dimension(:,:) :: a
        integer(int32), intent(in), dimension(:) :: ipvt
        complex(real64), intent(inout), dimension(:) :: b
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine solve_qr_no_pivot_mtx(a, tau, b, work, olwork, err)
        real(real64), intent(inout), dimension(:,:) :: a, b
        real(real64), intent(in), dimension(:) :: tau
        real(real64), intent(out), target, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine solve_qr_no_pivot_mtx_cmplx(a, tau, b, work, olwork, err)
        complex(real64), intent(inout), dimension(:,:) :: a, b
        complex(real64), intent(in), dimension(:) :: tau
        complex(real64), intent(out), target, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine solve_qr_no_pivot_vec(a, tau, b, work, olwork, err)
        real(real64), intent(inout), dimension(:,:) :: a
        real(real64), intent(in), dimension(:) :: tau
        real(real64), intent(inout), dimension(:) :: b
        real(real64), intent(out), target, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine solve_qr_no_pivot_vec_cmplx(a, tau, b, work, olwork, err)
        complex(real64), intent(inout), dimension(:,:) :: a
        complex(real64), intent(in), dimension(:) :: tau
        complex(real64), intent(inout), dimension(:) :: b
        complex(real64), intent(out), target, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine solve_qr_pivot_mtx(a, tau, jpvt, b, work, olwork, err)
        real(real64), intent(inout), dimension(:,:) :: a
        real(real64), intent(in), dimension(:) :: tau
        integer(int32), intent(in), dimension(:) :: jpvt
        real(real64), intent(inout), dimension(:,:) :: b
        real(real64), intent(out), target, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine solve_qr_pivot_mtx_cmplx(a, tau, jpvt, b, work, olwork, err)
        complex(real64), intent(inout), dimension(:,:) :: a
        complex(real64), intent(in), dimension(:) :: tau
        integer(int32), intent(in), dimension(:) :: jpvt
        complex(real64), intent(inout), dimension(:,:) :: b
        complex(real64), intent(out), target, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine solve_qr_pivot_vec(a, tau, jpvt, b, work, olwork, err)
        real(real64), intent(inout), dimension(:,:) :: a
        real(real64), intent(in), dimension(:) :: tau
        integer(int32), intent(in), dimension(:) :: jpvt
        real(real64), intent(inout), dimension(:) :: b
        real(real64), intent(out), target, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine solve_qr_pivot_vec_cmplx(a, tau, jpvt, b, work, olwork, err)
        complex(real64), intent(inout), dimension(:,:) :: a
        complex(real64), intent(in), dimension(:) :: tau
        integer(int32), intent(in), dimension(:) :: jpvt
        complex(real64), intent(inout), dimension(:) :: b
        complex(real64), intent(out), target, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine solve_cholesky_mtx(upper, a, b, err)
        logical, intent(in) :: upper
        real(real64), intent(in), dimension(:,:) :: a
        real(real64), intent(inout), dimension(:,:) :: b
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine solve_cholesky_mtx_cmplx(upper, a, b, err)
        logical, intent(in) :: upper
        complex(real64), intent(in), dimension(:,:) :: a
        complex(real64), intent(inout), dimension(:,:) :: b
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine solve_cholesky_vec(upper, a, b, err)
        logical, intent(in) :: upper
        real(real64), intent(in), dimension(:,:) :: a
        real(real64), intent(inout), dimension(:) :: b
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine solve_cholesky_vec_cmplx(upper, a, b, err)
        logical, intent(in) :: upper
        complex(real64), intent(in), dimension(:,:) :: a
        complex(real64), intent(inout), dimension(:) :: b
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine solve_least_squares_mtx(a, b, work, olwork, err)
        real(real64), intent(inout), dimension(:,:) :: a, b
        real(real64), intent(out), target, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine solve_least_squares_mtx_cmplx(a, b, work, olwork, err)
        complex(real64), intent(inout), dimension(:,:) :: a, b
        complex(real64), intent(out), target, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine solve_least_squares_vec(a, b, work, olwork, err)
        real(real64), intent(inout), dimension(:,:) :: a
        real(real64), intent(inout), dimension(:) :: b
        real(real64), intent(out), target, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine solve_least_squares_vec_cmplx(a, b, work, olwork, err)
        complex(real64), intent(inout), dimension(:,:) :: a
        complex(real64), intent(inout), dimension(:) :: b
        complex(real64), intent(out), target, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine solve_least_squares_mtx_pvt(a, b, ipvt, arnk, work, olwork, err)
        real(real64), intent(inout), dimension(:,:) :: a, b
        integer(int32), intent(inout), target, optional, dimension(:) :: ipvt
        integer(int32), intent(out), optional :: arnk
        real(real64), intent(out), target, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine solve_least_squares_mtx_pvt_cmplx(a, b, ipvt, arnk, &
            work, olwork, rwork, err)
        complex(real64), intent(inout), dimension(:,:) :: a, b
        integer(int32), intent(inout), target, optional, dimension(:) :: ipvt
        integer(int32), intent(out), optional :: arnk
        complex(real64), intent(out), target, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        real(real64), intent(out), target, optional, dimension(:) :: rwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine solve_least_squares_vec_pvt(a, b, ipvt, arnk, work, olwork, err)
        real(real64), intent(inout), dimension(:,:) :: a
        real(real64), intent(inout), dimension(:) :: b
        integer(int32), intent(inout), target, optional, dimension(:) :: ipvt
        integer(int32), intent(out), optional :: arnk
        real(real64), intent(out), target, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine solve_least_squares_vec_pvt_cmplx(a, b, ipvt, arnk, &
            work, olwork, rwork, err)
        complex(real64), intent(inout), dimension(:,:) :: a
        complex(real64), intent(inout), dimension(:) :: b
        integer(int32), intent(inout), target, optional, dimension(:) :: ipvt
        integer(int32), intent(out), optional :: arnk
        complex(real64), intent(out), target, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        real(real64), intent(out), target, optional, dimension(:) :: rwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine solve_least_squares_mtx_svd(a, b, s, arnk, work, olwork, err)
        real(real64), intent(inout), dimension(:,:) :: a, b
        integer(int32), intent(out), optional :: arnk
        real(real64), intent(out), target, optional, dimension(:) :: work, s
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine solve_least_squares_mtx_svd_cmplx(a, b, s, arnk, work, &
            olwork, rwork, err)
        complex(real64), intent(inout), dimension(:,:) :: a, b
        integer(int32), intent(out), optional :: arnk
        complex(real64), intent(out), target, optional, dimension(:) :: work
        real(real64), intent(out), target, optional, dimension(:) :: rwork, s
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine solve_least_squares_vec_svd(a, b, s, arnk, work, olwork, err)
        real(real64), intent(inout), dimension(:,:) :: a
        real(real64), intent(inout), dimension(:) :: b
        integer(int32), intent(out), optional :: arnk
        real(real64), intent(out), target, optional, dimension(:) :: work, s
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine solve_least_squares_vec_svd_cmplx(a, b, s, arnk, work, &
            olwork, rwork, err)
        complex(real64), intent(inout), dimension(:,:) :: a
        complex(real64), intent(inout), dimension(:) :: b
        integer(int32), intent(out), optional :: arnk
        complex(real64), intent(out), target, optional, dimension(:) :: work
        real(real64), intent(out), target, optional, dimension(:) :: rwork, s
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine mtx_inverse_dbl(a, iwork, work, olwork, err)
        real(real64), intent(inout), dimension(:,:) :: a
        integer(int32), intent(out), target, optional, dimension(:) :: iwork
        real(real64), intent(out), target, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine mtx_inverse_cmplx(a, iwork, work, olwork, err)
        complex(real64), intent(inout), dimension(:,:) :: a
        integer(int32), intent(out), target, optional, dimension(:) :: iwork
        complex(real64), intent(out), target, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine mtx_pinverse_dbl(a, ainv, tol, work, olwork, err)
        real(real64), intent(inout), dimension(:,:) :: a
        real(real64), intent(out), dimension(:,:) :: ainv
        real(real64), intent(in), optional :: tol
        real(real64), intent(out), target, dimension(:), optional :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine mtx_pinverse_cmplx(a, ainv, tol, work, olwork, rwork, err)
        complex(real64), intent(inout), dimension(:,:) :: a
        complex(real64), intent(out), dimension(:,:) :: ainv
        real(real64), intent(in), optional :: tol
        complex(real64), intent(out), target, dimension(:), optional :: work
        integer(int32), intent(out), optional :: olwork
        real(real64), intent(out), target, dimension(:), optional :: rwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine solve_lq_mtx(a, tau, b, work, olwork, err)
        real(real64), intent(in), dimension(:,:) :: a
        real(real64), intent(in), dimension(:) :: tau
        real(real64), intent(inout), dimension(:,:) :: b
        real(real64), intent(out), target, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine solve_lq_mtx_cmplx(a, tau, b, work, olwork, err)
        complex(real64), intent(in), dimension(:,:) :: a
        complex(real64), intent(in), dimension(:) :: tau
        complex(real64), intent(inout), dimension(:,:) :: b
        complex(real64), intent(out), target, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine solve_lq_vec(a, tau, b, work, olwork, err)
        real(real64), intent(in), dimension(:,:) :: a
        real(real64), intent(in), dimension(:) :: tau
        real(real64), intent(inout), dimension(:) :: b
        real(real64), intent(out), target, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine solve_lq_vec_cmplx(a, tau, b, work, olwork, err)
        complex(real64), intent(in), dimension(:,:) :: a
        complex(real64), intent(in), dimension(:) :: tau
        complex(real64), intent(inout), dimension(:) :: b
        complex(real64), intent(out), target, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
end interface
 
! ******************************************************************************
! LINALG_EIGEN.F90
! ------------------------------------------------------------------------------
interface
    module subroutine eigen_symm(vecs, a, vals, work, olwork, err)
        logical, intent(in) :: vecs
        real(real64), intent(inout), dimension(:,:) :: a
        real(real64), intent(out), dimension(:) :: vals
        real(real64), intent(out), pointer, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine eigen_asymm(a, vals, vecs, work, olwork, err)
        real(real64), intent(inout), dimension(:,:) :: a
        complex(real64), intent(out), dimension(:) :: vals
        complex(real64), intent(out), optional, dimension(:,:) :: vecs
        real(real64), intent(out), pointer, optional, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
    
    module subroutine eigen_gen(a, b, alpha, beta, vecs, work, olwork, err)
        real(real64), intent(inout), dimension(:,:) :: a, b
        complex(real64), intent(out), dimension(:) :: alpha
        real(real64), intent(out), optional, dimension(:) :: beta
        complex(real64), intent(out), optional, dimension(:,:) :: vecs
        real(real64), intent(out), optional, pointer, dimension(:) :: work
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine eigen_cmplx(a, vals, vecs, work, olwork, rwork, err)
        complex(real64), intent(inout), dimension(:,:) :: a
        complex(real64), intent(out), dimension(:) :: vals
        complex(real64), intent(out), optional, dimension(:,:) :: vecs
        complex(real64), intent(out), target, optional, dimension(:) :: work
        real(real64), intent(out), target, optional, dimension(:) :: rwork
        integer(int32), intent(out), optional :: olwork
        class(errors), intent(inout), optional, target :: err
    end subroutine
end interface
 
! ******************************************************************************
! LINALG_SORTING.F90
! ------------------------------------------------------------------------------
interface
    module subroutine sort_dbl_array(x, ascend)
        real(real64), intent(inout), dimension(:) :: x
        logical, intent(in), optional :: ascend
    end subroutine
    
    module subroutine sort_dbl_array_ind(x, ind, ascend, err)
        real(real64), intent(inout), dimension(:) :: x
        integer(int32), intent(inout), dimension(:) :: ind
        logical, intent(in), optional :: ascend
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine sort_cmplx_array(x, ascend)
        complex(real64), intent(inout), dimension(:) :: x
        logical, intent(in), optional :: ascend
    end subroutine
 
    module subroutine sort_cmplx_array_ind(x, ind, ascend, err)
        complex(real64), intent(inout), dimension(:) :: x
        integer(int32), intent(inout), dimension(:) :: ind
        logical, intent(in), optional :: ascend
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine sort_eigen_cmplx(vals, vecs, ascend, err)
        complex(real64), intent(inout), dimension(:) :: vals
        complex(real64), intent(inout), dimension(:,:) :: vecs
        logical, intent(in), optional :: ascend
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine sort_eigen_dbl(vals, vecs, ascend, err)
        real(real64), intent(inout), dimension(:) :: vals
        real(real64), intent(inout), dimension(:,:) :: vecs
        logical, intent(in), optional :: ascend
        class(errors), intent(inout), optional, target :: err
    end subroutine
 
    module subroutine sort_int32_array(x, ascend)
        integer(int32), intent(inout), dimension(:) :: x
        logical, intent(in), optional :: ascend
    end subroutine
 
    module subroutine sort_int32_array_ind(x, ind, ascend, err)
        integer(int32), intent(inout), dimension(:) :: x
        integer(int32), intent(inout), dimension(:) :: ind
        logical, intent(in), optional :: ascend
        class(errors), intent(inout), optional, target :: err
    end subroutine
end interface
 
! ******************************************************************************
! LINALG_SPARSE.F90
! ------------------------------------------------------------------------------
    interface nonzero_count
        module procedure :: nonzero_count_csr
        module procedure :: nonzero_count_msr
    end interface
 
    interface size
        module procedure :: csr_size
        module procedure :: msr_size
    end interface
 
    interface matmul
        module procedure :: csr_mtx_mtx_mult
        module procedure :: csr_mtx_vec_mult
    end interface
 
    interface operator(+)
        module procedure :: csr_mtx_add
    end interface
 
    interface operator(-)
        module procedure :: csr_mtx_sub
    end interface
 
    interface operator(*)
        module procedure :: csr_mtx_mult_scalar_1
        module procedure :: csr_mtx_mult_scalar_2
    end interface
 
    interface operator(/)
        module procedure :: csr_mtx_divide_scalar_1
    end interface
 
    interface assignment(=)
        module procedure :: csr_assign_to_dense
        module procedure :: dense_assign_to_csr
        module procedure :: msr_assign_to_dense
        module procedure :: dense_assign_to_msr
        module procedure :: csr_assign_to_msr
        module procedure :: msr_assign_to_csr
    end interface
 
    interface transpose
        module procedure :: csr_transpose
    end interface
 
    interface sparse_direct_solve
        module procedure :: csr_solve_sparse_direct
    end interface
 
    interface pgmres_solver
        module procedure :: csr_pgmres_solver
    end interface
 
    interface
        module function csr_get_element(this, i, j) result(rst)
            class(csr_matrix), intent(in) :: this
            integer(int32), intent(in) :: i, j
            real(real64) :: rst
        end function
        pure module function csr_size(x, dim) result(rst)
            class(csr_matrix), intent(in) :: x
            integer(int32), intent(in) :: dim
            integer(int32) :: rst
        end function
 
        module function create_empty_csr_matrix(m, n, nnz, err) result(rst)
            integer(int32), intent(in) :: m, n, nnz
            class(errors), intent(inout), optional, target :: err
            type(csr_matrix) :: rst
        end function
 
        pure module function nonzero_count_csr(x) result(rst)
            class(csr_matrix), intent(in) :: x
            integer(int32) :: rst
        end function
 
        module function dense_to_csr(a, err) result(rst)
            real(real64), intent(in), dimension(:,:) :: a
            class(errors), intent(inout), optional, target :: err
            type(csr_matrix) :: rst
        end function
 
        module subroutine csr_to_dense(a, x, err)
            class(csr_matrix), intent(in) :: a
            real(real64), intent(out), dimension(:,:) :: x
            class(errors), intent(inout), optional, target :: err
        end subroutine
 
        module function csr_mtx_mtx_mult(a, b) result(rst)
            class(csr_matrix), intent(in) :: a, b
            type(csr_matrix) :: rst
        end function
 
        module function csr_mtx_vec_mult(a, b) result(rst)
            class(csr_matrix), intent(in) :: a
            real(real64), intent(in), dimension(:) :: b
            real(real64), allocatable, dimension(:) :: rst
        end function
 
        module function csr_mtx_add(a, b) result(rst)
            class(csr_matrix), intent(in) :: a, b
            type(csr_matrix) :: rst
        end function
 
        module function csr_mtx_sub(a, b) result(rst)
            class(csr_matrix), intent(in) :: a, b
            type(csr_matrix) :: rst
        end function
 
        module function csr_mtx_mult_scalar_1(a, b) result(rst)
            class(csr_matrix), intent(in) :: a
            real(real64), intent(in) :: b
            type(csr_matrix) :: rst
        end function
 
        module function csr_mtx_mult_scalar_2(a, b) result(rst)
            real(real64), intent(in) :: a
            class(csr_matrix), intent(in) :: b
            type(csr_matrix) :: rst
        end function
 
        module function csr_mtx_divide_scalar_1(a, b) result(rst)
            class(csr_matrix), intent(in) :: a
            real(real64), intent(in) :: b
            type(csr_matrix) :: rst
        end function
 
        module function csr_transpose(a) result(rst)
            class(csr_matrix), intent(in) :: a
            type(csr_matrix) :: rst
        end function
 
        module subroutine extract_diagonal_csr(a, diag, err)
            class(csr_matrix), intent(in) :: a
            real(real64), intent(out), dimension(:) :: diag
            class(errors), intent(inout), optional, target :: err
        end subroutine
 
        module subroutine csr_solve_sparse_direct(a, b, x, droptol, err)
            class(csr_matrix), intent(in) :: a
            real(real64), intent(in), dimension(:) :: b
            real(real64), intent(out), dimension(:) :: x
            real(real64), intent(in), optional :: droptol
            class(errors), intent(inout), optional, target :: err
        end subroutine
 
        module function diag_to_csr(a, err) result(rst)
            real(real64), intent(in), dimension(:) :: a
            class(errors), intent(inout), optional, target :: err
            type(csr_matrix) :: rst
        end function
 
        module function banded_to_csr(m, ml, mu, a, err) result(rst)
            integer(int32), intent(in) :: m, ml, mu
            real(real64), intent(in), dimension(:,:) :: a
            class(errors), intent(inout), optional, target :: err
            type(csr_matrix) :: rst
        end function
 
        module subroutine csr_assign_to_dense(dense, sparse)
            real(real64), intent(out), dimension(:,:) :: dense
            class(csr_matrix), intent(in) :: sparse
        end subroutine
 
        module subroutine dense_assign_to_csr(sparse, dense)
            type(csr_matrix), intent(out) :: sparse
            real(real64), intent(in), dimension(:,:) :: dense
        end subroutine
 
        pure module function msr_size(x, dim) result(rst)
            class(msr_matrix), intent(in) :: x
            integer(int32), intent(in) :: dim
            integer(int32) :: rst
        end function
 
        pure module function nonzero_count_msr(x) result(rst)
            class(msr_matrix), intent(in) :: x
            integer(int32) :: rst
        end function
 
        module function create_empty_msr_matrix(m, n, nnz, err) result(rst)
            integer(int32), intent(in) :: m, n, nnz
            class(errors), intent(inout), optional, target :: err
            type(msr_matrix) :: rst
        end function
 
        module function csr_to_msr(a, err) result(rst)
            class(csr_matrix), intent(in) :: a
            class(errors), intent(inout), optional, target :: err
            type(msr_matrix) :: rst
        end function
 
        module function msr_to_csr(a, err) result(rst)
            class(msr_matrix), intent(in) :: a
            class(errors), intent(inout), optional, target :: err
            type(csr_matrix) :: rst
        end function
 
        module function dense_to_msr(a, err) result(rst)
            real(real64), intent(in), dimension(:,:) :: a
            class(errors), intent(inout), optional, target :: err
            type(msr_matrix) :: rst
        end function
 
        module subroutine msr_to_dense(a, x, err)
            ! Arguments
            class(msr_matrix), intent(in) :: a
            real(real64), intent(out), dimension(:,:) :: x
            class(errors), intent(inout), optional, target :: err
        end subroutine
 
        module subroutine msr_assign_to_dense(dense, msr)
            real(real64), intent(out), dimension(:,:) :: dense
            class(msr_matrix), intent(in) :: msr
        end subroutine
 
        module subroutine dense_assign_to_msr(msr, dense)
            type(msr_matrix), intent(out) :: msr
            real(real64), intent(in), dimension(:,:) :: dense
        end subroutine
 
        module subroutine csr_assign_to_msr(msr, csr)
            type(msr_matrix), intent(out) :: msr
            class(csr_matrix), intent(in) :: csr
        end subroutine
 
        module subroutine msr_assign_to_csr(csr, msr)
            type(csr_matrix), intent(out) :: csr
            class(msr_matrix), intent(in) :: msr
        end subroutine
 
        module subroutine csr_lu_factor(a, lu, ju, droptol, err)
            class(csr_matrix), intent(in) :: a
            type(msr_matrix), intent(out) :: lu
            integer(int32), intent(out), dimension(:) :: ju
            real(real64), intent(in), optional :: droptol
            class(errors), intent(inout), optional, target :: err
        end subroutine
 
        module subroutine csr_lu_solve(lu, ju, b, x, err)
            class(msr_matrix), intent(in) :: lu
            integer(int32), intent(in), dimension(:) :: ju
            real(real64), intent(in), dimension(:) :: b
            real(real64), intent(out), dimension(:) :: x
            class(errors), intent(inout), optional, target :: err
        end subroutine
 
        module subroutine csr_pgmres_solver(a, lu, ju, b, x, im, tol, maxits, &
            iout, err)
            class(csr_matrix), intent(in) :: a
            class(msr_matrix), intent(in) :: lu
            integer(int32), intent(in), dimension(:) :: ju
            real(real64), intent(inout), dimension(:) :: b
            real(real64), intent(out), dimension(:) :: x
            integer(int32), intent(in), optional :: im, maxits, iout
            real(real64), intent(in), optional :: tol
            class(errors), intent(inout), optional, target :: err
        end subroutine
 
        module function create_csr_matrix(m, n, rows, cols, vals, err) &
            result(rst)
            integer(int32), intent(in) :: m, n
            integer(int32), intent(in), dimension(:) :: rows, cols
            real(real64), intent(in), dimension(:) :: vals
            class(errors), intent(inout), optional, target :: err
            type(csr_matrix) :: rst
        end function
    end interface
 
! ------------------------------------------------------------------------------
end module