pure subroutine gesv_rel(A, b, x, info)
use external_interfaces_solver
! inputs:
real(rk), dimension(:, :), contiguous, intent(in) :: A ! input matrix A
real(rk), dimension(:), contiguous, intent(in) :: b ! right-hand side matrix b
! outputs:
real(rk), dimension(max(1, size(A, 2))), intent(out) :: x ! solution matrix x
integer, intent(out) :: info ! result info
! local variables
integer :: n, lda, ldb, nrhs
integer, dimension(size(A, 2)) :: ipiv
real(rk), dimension(:,:), allocatable :: a_copy
real(rk), dimension(:,:), allocatable :: b_copy
! get dimensions
nrhs = 1 ! size(b, 2)
n = size(A, 2)
lda = max(1, n)
ldb = max(1, n)
! copy the input matrices
a_copy = a
allocate(b_copy(ldb, nrhs))
b_copy(:, 1) = b
! call gels subroutine
call gesv(n, nrhs, a_copy, lda, ipiv, b_copy, ldb, info)
! copy the solution matrix
if (info == 0) then
x = b_copy(1:ldb, 1) ! nrhs = 1
else
error stop 'gesv failed'
end if
end subroutine gesv_rel