solves an overdetermined or underdetermined linear system using gels.
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
real(kind=rk), | intent(in), | dimension(:, :), contiguous | :: | A | ||
real(kind=rk), | intent(in), | dimension(:), contiguous | :: | b | ||
real(kind=rk), | intent(out), | dimension(max(1, size(A, 2))) | :: | x | ||
integer, | intent(out) | :: | info |
pure subroutine gels_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 character(1) :: trans integer :: m, n, lda, ldb, lwork, nrhs real(rk), allocatable :: work(:) real(rk) :: work1(1) real(rk), dimension(:,:), allocatable :: a_copy real(rk), dimension(:,:), allocatable :: b_copy ! trans = 'n' ! get dimensions nrhs = 1 ! size(b, 2) m = size(A, 1) n = size(A, 2) lda = max(1, m) ldb = max(1, max(m, n)) ! copy the input matrices a_copy = a allocate(b_copy(ldb, nrhs)) b_copy(:, 1) = b ! calculate the optimal size of the work array call gels(trans, m, n, nrhs, a_copy, lda, b_copy, ldb, work1, -1, info) ! allocate work array lwork = nint(work1(1)) allocate(work(lwork)) ! call gels subroutine call gels(trans, m, n, nrhs, a_copy, lda, b_copy, ldb, work, lwork, info) ! copy the solution matrix if (info == 0) then if (trans == 'n') x = b_copy(1:n, 1) ! nrhs = 1 if (trans == 't') x = b_copy(1:m, 1) ! nrhs = 1 else error stop 'gels failed' end if ! deallocate workspace deallocate(work) end subroutine gels_rel