dprepi Subroutine

subroutine dprepi(Neq, Y, S, Yh, Savr, Ewt, Rtem, Ia, Ja, Ic, Jc, Wk, Iwk, Ipper, res, jac, adda)

Uses

  • proc~~dprepi~2~~UsesGraph proc~dprepi~2 dprepi.inc::dprepi module~m_odepack M_odepack proc~dprepi~2->module~m_odepack

This routine performs preprocessing related to the sparse linear systems that must be solved.

The operations that are performed here are:

  • compute sparseness structure of the iteration matrix P = A - con*J according to MOSS,
  • compute grouping of column indices (MITER = 2),
  • compute a new ordering of rows and columns of the matrix,
  • reorder JA corresponding to the new ordering,
  • perform a symbolic LU factorization of the matrix, and
  • set pointers for segments of the IWK/WK array.

In addition to variables described previously, DPREPI uses the following for communication:

YH

the history array. Only the first column, containing the current Y vector, is used. Used only if MOSS .ne. 0.

S

array of length NEQ, identical to YDOTI in the driver, used only if MOSS .ne. 0.

SAVR

a work array of length NEQ, used only if MOSS .ne. 0.

EWT

array of length NEQ containing (inverted) error weights. Used only if MOSS = 2 or 4 or if ISTATE = MOSS = 1.

RTEM

a work array of length NEQ, identical to ACOR in the driver, used only if MOSS = 2 or 4.

WK

a real work array of length LENWK, identical to WM in the driver.

IWK

integer work array, assumed to occupy the same space as WK.

LENWK

the length of the work arrays WK and IWK.

ISTATC

a copy of the driver input argument ISTATE (= 1 on the first call, = 3 on a continuation call).

IYS

flag value from ODRV or CDRV.

IPPER

output error flag, with the following values and meanings: = 0 no error. = -1 insufficient storage for internal structure pointers. = -2 insufficient storage for JGROUP. = -3 insufficient storage for ODRV. = -4 other error flag from ODRV (should never occur). = -5 insufficient storage for CDRV. = -6 other error flag from CDRV. = -7 if the RES routine returned error flag IRES = IER = 2. = -8 if the RES routine returned error flag IRES = IER = 3.

Arguments

Type IntentOptional Attributes Name
integer :: Neq(*)
real(kind=dp), intent(inout) :: Y(*)
real(kind=dp) :: S(*)
real(kind=dp), intent(in) :: Yh(*)
real(kind=dp), intent(inout) :: Savr(*)
real(kind=dp), intent(in) :: Ewt(*)
real(kind=dp), intent(inout) :: Rtem(*)
integer, intent(in) :: Ia(*)
integer, intent(in) :: Ja(*)
integer, intent(in) :: Ic(*)
integer, intent(in) :: Jc(*)
real(kind=dp), intent(inout) :: Wk(*)
integer, intent(inout) :: Iwk(*)
integer, intent(out) :: Ipper
real :: res
integer :: jac
real :: adda

Calls

proc~~dprepi~2~~CallsGraph proc~dprepi~2 dprepi.inc::dprepi proc~adjlr~2 M_odepack::adjlr proc~dprepi~2->proc~adjlr~2 proc~cdrv M_odepack::cdrv proc~dprepi~2->proc~cdrv proc~cntnzu~2 M_odepack::cntnzu proc~dprepi~2->proc~cntnzu~2 proc~jgroup M_odepack::jgroup proc~dprepi~2->proc~jgroup proc~odrv~2 M_odepack::odrv proc~dprepi~2->proc~odrv~2

Variables

Type Visibility Attributes Name Initial
integer, public, parameter :: dp = kind(0.0d0)
real(kind=dp), public :: erwt
real(kind=dp), public :: fac
integer, public :: i
integer, public :: ibr
integer, public :: ier
integer, public :: ipil
integer, public :: ipiu
integer, public :: iptt1
integer, public :: iptt2
integer, public :: j
integer, public :: k
integer, public :: kamax
integer, public :: kamin
integer, public :: kcmax
integer, public :: kcmin
integer, public :: knew
integer, public :: ldif
integer, public :: lenigp
integer, public :: lenwk1
integer, public :: liwk
integer, public :: ljfo
integer, public :: maxg
integer, public :: np1
integer, public :: nzsut
real(kind=dp), public :: yj

Source Code

subroutine dprepi(Neq,Y,S,Yh,Savr,Ewt,Rtem,Ia,Ja,Ic,Jc,Wk,Iwk,Ipper,res,jac,adda)
Use M_odepack
implicit none
integer,parameter :: dp=kind(0.0d0)
!
integer                     :: Neq(*)
real(kind=dp),intent(inout) :: Y(*)
real(kind=dp)               :: S(*)
real(kind=dp),intent(in)    :: Yh(*)
real(kind=dp),intent(inout) :: Savr(*)
real(kind=dp),intent(in)    :: Ewt(*)
real(kind=dp),intent(inout) :: Rtem(*)
integer,intent(in)          :: Ia(*)
integer,intent(in)          :: Ja(*)
integer,intent(in)          :: Ic(*)
integer,intent(in)          :: Jc(*)
real(kind=dp),intent(inout) :: Wk(*)
integer,intent(inout)       :: Iwk(*)
integer,intent(out)         :: Ipper
external res
external jac
external adda

real(kind=dp) :: erwt, fac, yj
integer :: i, ibr, ier, ipil, ipiu, iptt1, iptt2, j, k, kamax, kamin, kcmax, kcmin, knew, ldif, lenigp, lenwk1, &
 & liwk, ljfo, maxg, np1, nzsut

   dlss%ibian = dlss%lrat*2
   dlss%ipian = dlss%ibian + 1
   np1 = dls1%n + 1
   dlss%ipjan = dlss%ipian + np1
   dlss%ibjan = dlss%ipjan - 1
   lenwk1 = dlss%lenwk - dls1%n
   liwk = dlss%lenwk*dlss%lrat
   if ( dlss%moss==0 ) liwk = liwk - dls1%n
   if ( dlss%moss==1 .or. dlss%moss==2 ) liwk = lenwk1*dlss%lrat
   if ( dlss%ipjan+dls1%n-1>liwk ) goto 600
   if ( dlss%moss/=0 ) then
      !
      if ( dlss%istatc/=3 ) then
         !  ISTATE = 1 and MOSS .ne. 0.  Perturb Y for structure determination.
         !  Initialize S with random nonzero elements for structure determination.
         do i = 1, dls1%n
            erwt = 1.0D0/Ewt(i)
            fac = 1.0D0 + 1.0D0/(i+1.0D0)
            Y(i) = Y(i) + fac*sign(erwt,Y(i))
            S(i) = 1.0D0 + fac*erwt
         enddo
         select case (dlss%moss)
         case (1)
            goto 100
         case (2)
            goto 200
         case (3)
            goto 300
         case (4)
            goto 400
         case default
         endselect
      endif
      !
      !  ISTATE = 3 and MOSS .ne. 0. Load Y from YH(*,1) and S from YH(*,2). --
      do i = 1, dls1%n
         Y(i) = Yh(i)
         S(i) = Yh(dls1%n+i)
      enddo
      select case (dlss%moss)
      case (1)
         goto 100
      case (2)
         goto 200
      case (3)
         goto 300
      case (4)
         goto 400
      case default
      endselect
   endif
   !
   !  MOSS = 0. Process user's IA,JA and IC,JC. ----------------------------
   knew = dlss%ipjan
   kamin = Ia(1)
   kcmin = Ic(1)
   Iwk(dlss%ipian) = 1
   do j = 1, dls1%n
      do i = 1, dls1%n
         Iwk(liwk+i) = 0
      enddo
      kamax = Ia(j+1) - 1
      if ( kamin<=kamax ) then
         do k = kamin, kamax
            i = Ja(k)
            Iwk(liwk+i) = 1
            if ( knew>liwk ) goto 600
            Iwk(knew) = i
            knew = knew + 1
         enddo
      endif
      kamin = kamax + 1
      kcmax = Ic(j+1) - 1
      if ( kcmin<=kcmax ) then
         do k = kcmin, kcmax
            i = Jc(k)
            if ( Iwk(liwk+i)==0 ) then
               if ( knew>liwk ) goto 600
               Iwk(knew) = i
               knew = knew + 1
            endif
         enddo
      endif
      Iwk(dlss%ipian+j) = knew + 1 - dlss%ipjan
      kcmin = kcmax + 1
   enddo
   goto 500
   !
   !  MOSS = 1. Compute structure from user-supplied Jacobian routine JAC. -
   !  A dummy call to RES allows user to create temporaries for use in JAC.
    100  continue
   ier = 1
   call res(Neq,dls1%tn,Y,S,Savr,ier)
   if ( ier>1 ) goto 1000
   do i = 1, dls1%n
      Savr(i) = 0.0D0
      Wk(lenwk1+i) = 0.0D0
   enddo
   k = dlss%ipjan
   Iwk(dlss%ipian) = 1
   do j = 1, dls1%n
      call adda(Neq,dls1%tn,Y,j,Iwk(dlss%ipian),Iwk(dlss%ipjan),Wk(lenwk1+1))
      call jac(Neq,dls1%tn,Y,S,j,Iwk(dlss%ipian),Iwk(dlss%ipjan),Savr)
      do i = 1, dls1%n
         ljfo = lenwk1 + i
         if ( Wk(ljfo)==0.0D0 ) then
            if ( Savr(i)==0.0D0 ) cycle
            Savr(i) = 0.0D0
         else
            Wk(ljfo) = 0.0D0
            Savr(i) = 0.0D0
         endif
         if ( k>liwk ) goto 600
         Iwk(k) = i
         k = k + 1
      enddo
      Iwk(dlss%ipian+j) = k + 1 - dlss%ipjan
   enddo
   goto 500
   !
   !  MOSS = 2. Compute structure from results of N + 1 calls to RES. ------
    200  continue
   do i = 1, dls1%n
      Wk(lenwk1+i) = 0.0D0
   enddo
   k = dlss%ipjan
   Iwk(dlss%ipian) = 1
   ier = -1
   if ( dls1%miter==1 ) ier = 1
   call res(Neq,dls1%tn,Y,S,Savr,ier)
   if ( ier>1 ) goto 1000
   do j = 1, dls1%n
      call adda(Neq,dls1%tn,Y,j,Iwk(dlss%ipian),Iwk(dlss%ipjan),Wk(lenwk1+1))
      yj = Y(j)
      erwt = 1.0D0/Ewt(j)
      Y(j) = yj + sign(erwt,yj)
      call res(Neq,dls1%tn,Y,S,Rtem,ier)
      if ( ier>1 ) return
      Y(j) = yj
      do i = 1, dls1%n
         ljfo = lenwk1 + i
         if ( Wk(ljfo)/=0.0D0 ) then
            Wk(ljfo) = 0.0D0
         elseif ( Rtem(i)==Savr(i) ) then
            cycle
         endif
         if ( k>liwk ) goto 600
         Iwk(k) = i
         k = k + 1
      enddo
      Iwk(dlss%ipian+j) = k + 1 - dlss%ipjan
   enddo
   goto 500
   !
   !  MOSS = 3. Compute structure from the user's IA/JA and JAC routine. ---
   !  A dummy call to RES allows user to create temporaries for use in JAC.
    300  continue
   ier = 1
   call res(Neq,dls1%tn,Y,S,Savr,ier)
   if ( ier>1 ) goto 1000
   do i = 1, dls1%n
      Savr(i) = 0.0D0
   enddo
   knew = dlss%ipjan
   kamin = Ia(1)
   Iwk(dlss%ipian) = 1
   do j = 1, dls1%n
      call jac(Neq,dls1%tn,Y,S,j,Iwk(dlss%ipian),Iwk(dlss%ipjan),Savr)
      kamax = Ia(j+1) - 1
      if ( kamin<=kamax ) then
         do k = kamin, kamax
            i = Ja(k)
            Savr(i) = 0.0D0
            if ( knew>liwk ) goto 600
            Iwk(knew) = i
            knew = knew + 1
         enddo
      endif
      kamin = kamax + 1
      do i = 1, dls1%n
         if ( Savr(i)/=0.0D0 ) then
            Savr(i) = 0.0D0
            if ( knew>liwk ) goto 600
            Iwk(knew) = i
            knew = knew + 1
         endif
      enddo
      Iwk(dlss%ipian+j) = knew + 1 - dlss%ipjan
   enddo
   goto 500
   !
   !  MOSS = 4. Compute structure from user's IA/JA and N + 1 RES calls. ---
    400  continue
   knew = dlss%ipjan
   kamin = Ia(1)
   Iwk(dlss%ipian) = 1
   ier = -1
   if ( dls1%miter==1 ) ier = 1
   call res(Neq,dls1%tn,Y,S,Savr,ier)
   if ( ier>1 ) goto 1000
   do j = 1, dls1%n
      yj = Y(j)
      erwt = 1.0D0/Ewt(j)
      Y(j) = yj + sign(erwt,yj)
      call res(Neq,dls1%tn,Y,S,Rtem,ier)
      if ( ier>1 ) return
      Y(j) = yj
      kamax = Ia(j+1) - 1
      if ( kamin<=kamax ) then
         do k = kamin, kamax
            i = Ja(k)
            Rtem(i) = Savr(i)
            if ( knew>liwk ) goto 600
            Iwk(knew) = i
            knew = knew + 1
         enddo
      endif
      kamin = kamax + 1
      do i = 1, dls1%n
         if ( Rtem(i)/=Savr(i) ) then
            if ( knew>liwk ) goto 600
            Iwk(knew) = i
            knew = knew + 1
         endif
      enddo
      Iwk(dlss%ipian+j) = knew + 1 - dlss%ipjan
   enddo
   !
    500  continue
   if ( dlss%moss/=0 .and. dlss%istatc/=3 ) then
   !  If ISTATE = 0 or 1 and MOSS .ne. 0, restore Y from YH. ---------------
      do i = 1, dls1%n
         Y(i) = Yh(i)
      enddo
   endif
   dlss%nnz = Iwk(dlss%ipian+dls1%n) - 1
   Ipper = 0
   dlss%ngp = 0
   lenigp = 0
   dlss%ipigp = dlss%ipjan + dlss%nnz
   if ( dls1%miter==2 ) then
   !
   !  Compute grouping of column indices (MITER = 2). ----------------------
   !
      maxg = np1
      dlss%ipjgp = dlss%ipjan + dlss%nnz
      dlss%ibjgp = dlss%ipjgp - 1
      dlss%ipigp = dlss%ipjgp + dls1%n
      iptt1 = dlss%ipigp + np1
      iptt2 = iptt1 + dls1%n
      dlss%lreq = iptt2 + dls1%n - 1
      if ( dlss%lreq>liwk ) goto 700
      call jgroup(dls1%n,Iwk(dlss%ipian),Iwk(dlss%ipjan),maxg,dlss%ngp,Iwk(dlss%ipigp),Iwk(dlss%ipjgp),Iwk(iptt1),Iwk(iptt2),ier)
      if ( ier/=0 ) goto 700
      lenigp = dlss%ngp + 1
   endif
   !
   !  Compute new ordering of rows/columns of Jacobian. --------------------
   dlss%ipr = dlss%ipigp + lenigp
   dlss%ipc = dlss%ipr
   dlss%ipic = dlss%ipc + dls1%n
   dlss%ipisp = dlss%ipic + dls1%n
   dlss%iprsp = (dlss%ipisp-2)/dlss%lrat + 2
   dlss%iesp = dlss%lenwk + 1 - dlss%iprsp
   if ( dlss%iesp<0 ) goto 800
   ibr = dlss%ipr - 1
   do i = 1, dls1%n
      Iwk(ibr+i) = i
   enddo
   dlss%nsp = liwk + 1 - dlss%ipisp
   call odrv(dls1%n,Iwk(dlss%ipian),Iwk(dlss%ipjan),Wk,Iwk(dlss%ipr),Iwk(dlss%ipic),dlss%nsp,Iwk(dlss%ipisp),1,dlss%iys)
   if ( dlss%iys==11*dls1%n+1 ) then
   !
      Ipper = -4
      return
   else
      if ( dlss%iys/=0 ) goto 800
   !
   !  Reorder JAN and do symbolic LU factorization of matrix. --------------
      dlss%ipa = dlss%lenwk + 1 - dlss%nnz
      dlss%nsp = dlss%ipa - dlss%iprsp
      dlss%lreq = max(12*dls1%n/dlss%lrat,6*dls1%n/dlss%lrat+2*dls1%n+dlss%nnz) + 3
      dlss%lreq = dlss%lreq + dlss%iprsp - 1 + dlss%nnz
      if ( dlss%lreq>dlss%lenwk ) goto 900
      dlss%iba = dlss%ipa - 1
      do i = 1, dlss%nnz
         Wk(dlss%iba+i) = 0.0D0
      enddo
      dlss%ipisp = dlss%lrat*(dlss%iprsp-1) + 1
      call cdrv(dls1%n,Iwk(dlss%ipr),Iwk(dlss%ipc),Iwk(dlss%ipic),Iwk(dlss%ipian), &
       & Iwk(dlss%ipjan),Wk(dlss%ipa),Wk(dlss%ipa), &
       & Wk(dlss%ipa),dlss%nsp,Iwk(dlss%ipisp),Wk(dlss%iprsp),dlss%iesp,5,dlss%iys)
      dlss%lreq = dlss%lenwk - dlss%iesp
      if ( dlss%iys==10*dls1%n+1 ) goto 900
      if ( dlss%iys/=0 ) then
   !
         Ipper = -6
         dlss%lreq = dlss%lenwk
         return
      else
         ipil = dlss%ipisp
         ipiu = ipil + 2*dls1%n + 1
         dlss%nzu = Iwk(ipil+dls1%n) - Iwk(ipil)
         dlss%nzl = Iwk(ipiu+dls1%n) - Iwk(ipiu)
         if ( dlss%lrat<=1 ) then
            call adjlr(dls1%n,Iwk(dlss%ipisp),ldif)
            dlss%lreq = dlss%lreq + ldif
         endif
         if ( dlss%lrat==2 .and. dlss%nnz==dls1%n ) dlss%lreq = dlss%lreq + 1
         dlss%nsp = dlss%nsp + dlss%lreq - dlss%lenwk
         dlss%ipa = dlss%lreq + 1 - dlss%nnz
         dlss%iba = dlss%ipa - 1
         Ipper = 0
         return
      endif
   endif
   !
    600  continue
   Ipper = -1
   dlss%lreq = 2 + (2*dls1%n+1)/dlss%lrat
   dlss%lreq = max(dlss%lenwk+1,dlss%lreq)
   return
   !
    700  continue
   Ipper = -2
   dlss%lreq = (dlss%lreq-1)/dlss%lrat + 1
   return
   !
    800  continue
   Ipper = -3
   call cntnzu(dls1%n,Iwk(dlss%ipian),Iwk(dlss%ipjan),nzsut)
   dlss%lreq = dlss%lenwk - dlss%iesp + (3*dls1%n+4*nzsut-1)/dlss%lrat + 1
   return
   !
    900  continue
   Ipper = -5
   return
   !
    1000 continue
   Ipper = -ier - 5
   dlss%lreq = 2 + (2*dls1%n+1)/dlss%lrat
end subroutine dprepi