This subroutine solves the linear system for the sparse preconditioner , given a vector , using the incomplete LU decomposition produced by jac_ilupre. The solution is carried out by lusol from SPARSKIT.
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
integer, | intent(in) | :: | neq |
Problem size. |
||
real(kind=rk), | intent(in) | :: | t |
Current independent variable (not used). |
||
real(kind=rk), | intent(in) | :: | y(*) |
Current dependent variables (not used). |
||
real(kind=rk), | intent(in) | :: | ydot(*) |
Current derivatives of dependent variables (not used). |
||
real(kind=rk), | intent(in) | :: | r0(*) |
Current residual evaluated at |
||
real(kind=rk), | intent(in) | :: | wk(*) |
Real work space available to this subroutine. |
||
real(kind=rk), | intent(in) | :: | cj |
Scalar used in forming the system Jacobian. |
||
real(kind=rk), | intent(in) | :: | wght(*) |
Error weights for computing norms. |
||
real(kind=rk), | intent(inout) | :: | rwp(*) |
Matrix elements of ILU. |
||
integer, | intent(inout) | :: | iwp(*) |
Array indices for elements of ILU |
||
real(kind=rk), | intent(inout) | :: | b(*) |
Right-hand side vector on input; solution on output. |
||
real(kind=rk), | intent(in) | :: | epslin |
Tolerance for linear system (not used). |
||
integer, | intent(inout) | :: | ierr |
Error flag (not used). |
||
real(kind=rk), | intent(inout) | :: | rpar(*) |
Real array used for communication between the calling program and user routines (not used). |
||
integer, | intent(inout) | :: | ipar(*) |
Integer array used for communication between the calling program and user routines. |
subroutine psol_ilupre( & neq, t, y, ydot, r0, wk, cj, wght, rwp, iwp, b, epslin, ierr, rpar, ipar) !! This subroutine solves the linear system \(P x = b\) for the sparse preconditioner \(P\), !! given a vector \(b\), using the incomplete LU decomposition produced by [[jac_ilupre]]. !! The solution is carried out by [[lusol]] from SPARSKIT. use dsparskit, only: lusol, dvperm integer, intent(in) :: neq !! Problem size. real(rk), intent(in) :: t !! Current independent variable (not used). real(rk), intent(in) :: y(*) !! Current dependent variables (not used). real(rk), intent(in) :: ydot(*) !! Current derivatives of dependent variables (not used). real(rk), intent(in) :: r0(*) !! Current residual evaluated at `(t, y, ydot)` (not used). real(rk), intent(in) :: wk(*) !! Real work space available to this subroutine. real(rk), intent(in) :: cj !! Scalar used in forming the system Jacobian. real(rk), intent(in) :: wght(*) !! Error weights for computing norms. real(rk), intent(inout) :: rwp(*) !! Matrix elements of ILU. integer, intent(inout) :: iwp(*) !! Array indices for elements of ILU real(rk), intent(inout) :: b(*) !! Right-hand side vector on input; solution on output. real(rk), intent(in) :: epslin !! Tolerance for linear system (not used). integer, intent(inout) :: ierr !! Error flag (not used). real(rk), intent(inout) :: rpar(*) !! Real array used for communication between the calling program and user routines !! (not used). integer, intent(inout) :: ipar(*) !! Integer array used for communication between the calling program and user routines. integer :: i, lplu, lju, ljlu, lrownrms, lperm, lqperm, ireorder, isrnorm, ipremeth, & jscalcol ! Load IPREMETH, IREORDER and ISRNORM values from IPAR. ipremeth = ipar(5) ireorder = ipar(7) isrnorm = ipar(8) jscalcol = ipar(11) ! Load pointers into RWP and IWP arrays. lplu = ipar(21) lju = ipar(22) ljlu = ipar(23) lrownrms = ipar(24) lperm = ipar(25) lqperm = ipar(26) ! Scale c by multiplying by row-normalization factors, if used. if (isrnorm == 1) then do i = 1, neq b(i) = b(i)*rwp(lrownrms + i - 1) end do end if ! Solve P*x=b for a preconditioner stored as a sparse matrix in compressed sparse row ! format. If rows and columns of P were reordered (permuted), permute b, then use inverse ! permutation on x. if (ipremeth == 1 .or. ipremeth == 2) then if (ireorder == 1) call dvperm(neq, b, iwp(lperm)) call lusol(neq, b, wk, rwp(lplu), iwp(ljlu), iwp(lju)) if (ireorder == 1) call dvperm(neq, wk, iwp(lqperm)) end if ! Unscale x by dividing by column scaling vector WGHT. if (jscalcol == 1) then b(1:neq) = wk(1:neq)/wght(1:neq) else b(1:neq) = wk(1:neq) end if ierr = 0 end subroutine psol_ilupre