DPRJIS is called to compute and process the matrix P = A - HEL(1)J, where J is an approximation to the Jacobian dr/dy, where r = g(t,y) - A(t,y)*s.
J is computed by columns, either by the user-supplied routine JAC if MITER = 1, or by finite differencing if MITER = 2.
J is stored in WK, rescaled, and ADDA is called to generate P.
The matrix P is subjected to LU decomposition in CDRV. P and its LU decomposition are stored separately in WK.
In addition to variables described previously, communication with DPRJIS uses the following:
array containing predicted values on entry.
work array of length N (ACOR in DSTODI).
array containing r evaluated at predicted y. On output it contains the residual evaluated at current values of t and y.
array containing predicted values of dy/dt (SAVF in DSTODI).
real work space for matrices. On output it contains P and its sparse LU decomposition. Storage of matrix elements starts at WK(3). WK also contains the following matrix-related data. WK(1) = SQRT(UROUND), used in numerical Jacobian increments.
integer work space for matrix-related data, assumed to be equivalenced to WK. In addition, WK(IPRSP) and IWK(IPISP) are assumed to have identical locations.
EL(1) (input).
output error flag (in COMMON). = 0 if no error. = 1 if zero pivot found in CDRV. = IRES (= 2 or 3) if RES returned IRES = 2 or 3. = -1 if insufficient storage for CDRV (should not occur). = -2 if other error found in CDRV (should not occur here).
output flag = 1 to indicate that the Jacobian matrix (or approximation) is now current.
This routine also uses other variables in global structures.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | dimension(*) | :: | Neq | |||
| real(kind=dp), | intent(inout), | dimension(*) | :: | Y | ||
| real(kind=dp), | intent(in), | dimension(Nyh,*) | :: | Yh | ||
| integer, | intent(in) | :: | Nyh | |||
| real(kind=dp), | intent(in), | dimension(*) | :: | Ewt | ||
| real(kind=dp), | intent(inout), | dimension(*) | :: | Rtem | ||
| real(kind=dp), | dimension(*) | :: | Savr | |||
| real(kind=dp), | dimension(*) | :: | S | |||
| real(kind=dp), | intent(inout), | dimension(*) | :: | Wk | ||
| integer, | dimension(*) | :: | Iwk | |||
| real | :: | res | ||||
| integer | :: | jac | ||||
| real | :: | adda |
| Type | Visibility | Attributes | Name | Initial | |||
|---|---|---|---|---|---|---|---|
| real(kind=dp), | public | :: | con | ||||
| real(kind=dp), | public | :: | fac | ||||
| real(kind=dp), | public | :: | hl0 | ||||
| integer, | public | :: | i | ||||
| integer, | public | :: | imul | ||||
| integer, | public | :: | ires | ||||
| integer, | public | :: | j | ||||
| integer, | public | :: | jj | ||||
| integer, | public | :: | jmax | ||||
| integer, | public | :: | jmin | ||||
| integer, | public | :: | k | ||||
| integer, | public | :: | kmax | ||||
| integer, | public | :: | kmin | ||||
| integer, | public | :: | ng | ||||
| real(kind=dp), | public | :: | r | ||||
| real(kind=dp), | public | :: | srur |
subroutine dprjis(Neq,Y,Yh,Nyh,Ewt,Rtem,Savr,S,Wk,Iwk,res,jac,adda) ! integer, dimension(*) :: Neq real(kind=dp), intent(inout), dimension(*) :: Y integer, intent(in) :: Nyh real(kind=dp), intent(in), dimension(Nyh,*) :: Yh real(kind=dp), intent(in), dimension(*) :: Ewt real(kind=dp), intent(inout), dimension(*) :: Rtem real(kind=dp), dimension(*) :: Savr real(kind=dp), dimension(*) :: S real(kind=dp), intent(inout), dimension(*) :: Wk integer, dimension(*) :: Iwk external :: res external :: jac external :: adda ! real(kind=dp) :: con, fac, hl0, r, srur integer :: i, imul, ires, j, jj, jmax, jmin, k, kmax, kmin, ng ! hl0 = dls1%h*dls1%el0 con = -hl0 dls1%jcur = 1 dls1%nje = dls1%nje + 1 if ( dls1%miter==2 ) then ! ! If MITER = 2, make NGP + 1 calls to RES to approximate J and P. ------ ires = -1 call res(Neq,dls1%tn,Y,S,Savr,ires) dls1%nfe = dls1%nfe + 1 if ( ires>1 ) then ! Error return for IRES = 2 or IRES = 3 return from RES. --------------- dls1%ierpj = ires return else srur = Wk(1) jmin = Iwk(dlss%ipigp) do ng = 1, dlss%ngp jmax = Iwk(dlss%ipigp+ng) - 1 do j = jmin, jmax jj = Iwk(dlss%ibjgp+j) r = max(srur*abs(Y(jj)),0.01D0/Ewt(jj)) Y(jj) = Y(jj) + r enddo call res(Neq,dls1%tn,Y,S,Rtem,ires) dls1%nfe = dls1%nfe + 1 if ( ires>1 ) then dls1%ierpj = ires return else do j = jmin, jmax jj = Iwk(dlss%ibjgp+j) Y(jj) = Yh(jj,1) r = max(srur*abs(Y(jj)),0.01D0/Ewt(jj)) fac = -hl0/r kmin = Iwk(dlss%ibian+jj) kmax = Iwk(dlss%ibian+jj+1) - 1 do k = kmin, kmax i = Iwk(dlss%ibjan+k) Rtem(i) = (Rtem(i)-Savr(i))*fac enddo call adda(Neq,dls1%tn,Y,jj,Iwk(dlss%ipian),Iwk(dlss%ipjan),Rtem) do k = kmin, kmax i = Iwk(dlss%ibjan+k) Wk(dlss%iba+k) = Rtem(i) enddo enddo jmin = jmax + 1 endif enddo ires = 1 call res(Neq,dls1%tn,Y,S,Savr,ires) dls1%nfe = dls1%nfe + 1 if ( ires>1 ) then dls1%ierpj = ires return endif endif else ! ! If MITER = 1, call RES, then call JAC and ADDA for each column. ------ ires = 1 call res(Neq,dls1%tn,Y,S,Savr,ires) dls1%nfe = dls1%nfe + 1 if ( ires>1 ) then dls1%ierpj = ires return else kmin = Iwk(dlss%ipian) do j = 1, dls1%n kmax = Iwk(dlss%ipian+j) - 1 do i = 1, dls1%n Rtem(i) = 0.0D0 enddo call jac(Neq,dls1%tn,Y,S,j,Iwk(dlss%ipian),Iwk(dlss%ipjan),Rtem) do i = 1, dls1%n Rtem(i) = Rtem(i)*con enddo call adda(Neq,dls1%tn,Y,j,Iwk(dlss%ipian),Iwk(dlss%ipjan),Rtem) do k = kmin, kmax i = Iwk(dlss%ibjan+k) Wk(dlss%iba+k) = Rtem(i) enddo kmin = kmax + 1 enddo endif endif ! ! Do numerical factorization of P matrix. ------------------------------ dlss%nlu = dlss%nlu + 1 dls1%ierpj = 0 do i = 1, dls1%n Rtem(i) = 0.0D0 enddo call cdrv(dls1%n,Iwk(dlss%ipr),Iwk(dlss%ipc),Iwk(dlss%ipic), & & Iwk(dlss%ipian),Iwk(dlss%ipjan),Wk(dlss%ipa),Rtem,Rtem,dlss%nsp, & & Iwk(dlss%ipisp),Wk(dlss%iprsp),dlss%iesp,2,dlss%iys) if ( dlss%iys==0 ) return imul = (dlss%iys-1)/dls1%n dls1%ierpj = -2 if ( imul==8 ) dls1%ierpj = 1 if ( imul==10 ) dls1%ierpj = -1 end subroutine dprjis