This routine manages the solution of the linear system arising in the Newton iteration.
Real matrix information and real temporary storage is stored in the array rwm
, while
integer matrix information is stored in the array iwm
.
For a dense matrix, the LINPACK routine dgesl is called. For a banded matrix, the
LINPACK routine dgbsl is called.
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
integer, | intent(in) | :: | neq |
Problem size. |
||
real(kind=rk), | intent(inout) | :: | delta(neq) |
On entry, the right-hand side vector of the linear system to be solved, and, on exit, the solution vector. |
||
real(kind=rk), | intent(inout) | :: | rwm(*) |
Real workspace for the linear system solver. |
||
integer, | intent(inout) | :: | iwm(*) |
Integer workspace for the linear system solver. |
Type | Visibility | Attributes | Name | Initial | |||
---|---|---|---|---|---|---|---|
integer, | public | :: | lipvt | ||||
integer, | public | :: | meband | ||||
integer, | public | :: | mtype |
subroutine dslvd(neq, delta, rwm, iwm) !! This routine manages the solution of the linear system arising in the Newton iteration. !! Real matrix information and real temporary storage is stored in the array `rwm`, while !! integer matrix information is stored in the array `iwm`. !! For a dense matrix, the LINPACK routine [[dgesl]] is called. For a banded matrix, the !! LINPACK routine [[dgbsl]] is called. use daskr_kinds, only: rk use daskr use linpack, only: dgesl, dgbsl implicit none integer, intent(in) :: neq !! Problem size. real(rk), intent(inout) :: delta(neq) !! On entry, the right-hand side vector of the linear system to be solved, !! and, on exit, the solution vector. real(rk), intent(inout) :: rwm(*) !! Real workspace for the linear system solver. integer, intent(inout) :: iwm(*) !! Integer workspace for the linear system solver. integer :: lipvt, meband, mtype lipvt = iwm(iwloc_liwp) mtype = iwm(iwloc_mtype) select case (mtype) case (1, 2) ! dense matrix. call dgesl(rwm, neq, neq, iwm(lipvt), delta, 0) case (3) ! dummy section for mtype=3. error stop "error: mtype=3 not implemented" case (4, 5) ! banded matrix. meband = 2*iwm(iwloc_ml) + iwm(iwloc_mu) + 1 call dgbsl(rwm, meband, neq, iwm(iwloc_ml), iwm(iwloc_mu), iwm(lipvt), delta, 0) case default error stop "error: unexpected value for mtype" end select end subroutine dslvd