           !------------------------!
           ! Cubpack User Interface !
           !------------------------!
Module CUI

USE Precision_Model
USE internal_types


Implicit NONE

PRIVATE

PUBLIC :: CUBATR, CUBPACK_INFO
!-----------------------------------------------------------------------
!***BEGIN PROLOGUE CUBATR
!***DATE WRITTEN   901114   (YYMMDD)
!***REVISION DATE  970620   (YYMMDD)
!***REVISION DATE  980406   (YYMMDD) (MDIV removed)
!***REVISION DATE  000809   (YYMMDD)
!***REVISION DATE  010719   (YYMMDD)
!***REVISION DATE  020715   (YYMMDD) (CUBPACK_INFO added)
!***AUTHOR
!          Ronald Cools, Dept. of Computer Science,
!          Katholieke Universiteit Leuven, Celestijnenlaan 200A,
!          B-3001 Heverlee, Belgium
!          Email:  Ronald.Cools@cs.kuleuven.ac.be
!***REVISION DATE  150729   (YYMMDD) (Made threadsafe / re-entrant)
!***MODIFIER
!          Pat Scott, Dept. of Physics, Imperial College London
!          Email: p.scott@imperial
!
!***PURPOSE  Computation of integrals over a collection of regions.
!
!***DESCRIPTION
!          CUBATR is the driver routine for CUBPACK and the only
!          routine that a user has to deal with (at the moment).
!
!-----------------------------------------------------------------------
PRIVATE :: CUBATR_X, CUBATR_1, CUBATR_CLEAR

!
! Module variables
!
INTEGER, SAVE, PRIVATE :: PreJob=0
INTEGER, PRIVATE                :: BOTTIH,BOTTRH
INTEGER,         DIMENSION(:), PRIVATE, ALLOCATABLE  :: IWork
REAL(kind=stnd), DIMENSION(:), PRIVATE, ALLOCATABLE  :: RWork
TYPE(EPSALG_MEM), PRIVATE              :: M
!$omp threadprivate (PreJob, BOTTIH, BOTTRH, IWork, RWork, M)

INTERFACE CUBATR
   MODULE PROCEDURE CUBATR_X, CUBATR_1, CUBATR_CLEAR
END INTERFACE

CONTAINS

SUBROUTINE CUBATR_X     &
     (DIMENS,NumFun,Integrand,NumRgn,Vertices,RgType,Value,AbsErr,    &
!   and optional parameters
IFAIL,Neval,EpsAbs,EpsRel,Restart,MinPts,MaxPts,Key,Job,Tune)
!-----------------------------------------------------------------------
!   Input parameters
!   ----------------
!
!   DIMENS Integer.
!          The dimension of the region of integration.
!
!   NumFun Integer.
!          Number of components of the integrand.
!
!   Integrand
!          Externally declared function for computing all components
!          of the integrand at the given evaluation point.
!          It must have input parameter X:
!              X(1)   The x-coordinate of the evaluation point.
!              X(2)   The y-coordinate of the evaluation point.
!              ...
!              X(DIMENS) The z-coordinate of the evaluation point.
!         and NumFun, the number of components of the integrand.
!         It must be compatible with the following interface:
!           INTERFACE 
!              FUNCTION Integrand(NUMFUN,X) RESULT(Value)
!                USE Precision_Model
!                INTEGER, INTENT(IN) :: NUMFUN
!                REAL(kind=stnd), DIMENSION(:), INTENT(IN) :: X
!                REAL(kind=stnd), DIMENSION(NUMFUN) :: Value
!              END FUNCTION Integrand
!           END INTERFACE
!
!   NumRgn Integer.
!          The number of given regions.
!
!   Vertices
!          Real array of dimension (DIMENS,DIMENS+1,NumRgn).
!          Vertices(1:DIMENS,K,L) are the x, y, ... coordinates
!          respectively of vertex K of region L, where
!          K = 1,...,DIMENS+1 and L = 1,...,NumRgn.
!
!   RgType Integer array of dimension (NumRgn).
!          RgType(L) describes the type of region L.
!
!   Value  Real array of dimension NumFun.
!          Approximations to all components of the integral if
!          the procedure is restarted.
!
!   AbsErr Real array of dimension NumFun.
!          Estimates of absolute errors if the procedure is restarted.
!
!   IFAIL  Optional integer argument.
!          This follows the NAG convention:
!          IFAIL = 1 : soft silent error
!                      Control returned to calling program.
!          IFAIL = -1: soft noisy error
!                      Error message is printed.
!                      Control returned to calling program.
!          IFAIL = 0 : hard noisy error
!                      Error message is printed and program is stopped.
!          Default IFAIL = -1.
!
!   EpsAbs Optional real argument.
!          Requested absolute error.
!          Default  EpsAbs = 0.
!
!   EpsRel Optional real argument.
!          Requested relative error.
!          Default EpsRel = sqrt(machine precision).
!
!   Restart Optional boolean argument.
!          If Restart = FALSE, this is the first attempt to compute
!                              the integral.
!          If Restart = TRUE, then we restart a previous attempt.
!          In this case the only parameters for CUBATR that may
!          be changed (with respect to the previous call of CUBATR)
!          are MinPts, MaxPts, EpsAbs, EpsRel, Key and Restart.
!          Default Restart = FALSE.
!
!   MinPts Optional integer argument.
!          The minimum allowed number of integrand evaluations.
!          Default MinPts = 0.
!
!   MaxPts Optional integer argument.
!          The maximum allowed number of integrand evaluations.
!          Default MaxPts = enough to do 500 subdivisions.
!
!   Key    Optional integer argument.
!          Can be used by Rule_General to choose between several
!          local integration rules.
!          Default Key = 2 if Dimension=1 and extrapolation is used 
!                                        (This corresponds to QAGS)
!          Default Key = 0 otherwise
!
!   Job    Optional integer argument.
!          If |Job| = 0, then nothing will be done except freeing all
!                        allocated memory.
!                        This is usefull after a call of CUBATR if no
!                        Restart will be done later and memory usage
!                        might become an issue later.
!                        Equivalently, one can call CUBATR()
!                        without any arguments.
!                   = 1, the global adaptive algorithm is called
!                   = 2, extrapolation using the epsilon algorithm is used.
!                   = 11, a region will be divided in 2**DIMENS subregions
!                        and the global adaptive algorithm is called.
!                        In combination with Key=0, this resembles DUCTRI and DCUTET.
!                   = 12, a region will be divided in 2 subregions
!                        and the global adaptive algorithm is called.
!                        In combination with Key=3 or 4, this resembles DCUHRE.
!          If Job < 0, then an overview of the Region Collection is dumped.
!          This will create the files tmp_integerstore and tmp_realstore.
!          Default Job = 1.
!
!   Tune   Optional real argument.
!          Can be used by Global_Adapt or the local error estimators
!          to influence the reliability. 0 <= Tune <= 1.
!          Tune = 1 is the most reliable available.
!          Default Tune = 1.
!          Note that this is an experimental and controversial parameter.
!          In this version, only Tune = 1 is supported for all regions.
!
!   Output parameters
!   -----------------
!
!   Value  Real array of dimension NumFun.
!          Approximations to all components of the integral
!
!   AbsErr Real array of dimension NumFun.
!          Estimates of absolute errors.
!
!   NEval  Optional Integer.
!          Number of integrand evaluations used by CUBATR for this call.
!
!   IFAIL  Optional Integer.
!          IFAIL = 0 for normal exit.
!
!            AbsErr(K) <=  EpsAbs or
!            AbsErr(K) <=  ABS(Value(K))*EpsRel with MaxPts or less
!            function evaluations for all values of K,
!            1 <= K <= NumFun .
!
!          IFAIL = 1 if MaxPts was too small to obtain the required
!            accuracy. In this case Global_Adapt returns values of
!            Value with estimated absolute errors AbsErr.
!
!          IFAIL > 1 in more serious case of trouble.
!-----------------------------------------------------------------------
! MODULES USED
USE Check_Input
USE Error_Handling
USE DS_ROUTINES, ONLY: DSCOPY, DSINIT, DSUSED, DSSTAT, DSSUM, DSPINT
USE CubatureRule_General, ONLY: Rule_Cost
USE Global_Adaptive_Algorithm
!***END PROLOGUE CUBATR
!-----------------------------------------------------------------------
!
! Global variables
!
INTERFACE 
   FUNCTION Integrand(NUMFUN,X) RESULT(Value)
     USE Precision_Model
     INTEGER, INTENT(IN) :: NUMFUN
     REAL(kind=stnd), DIMENSION(:), INTENT(IN) :: X
     REAL(kind=stnd), DIMENSION(NUMFUN) :: Value
   END FUNCTION Integrand
END INTERFACE
INTEGER, INTENT(IN)               :: DIMENS,NumFun,NumRgn
INTEGER, DIMENSION(:), INTENT(IN) :: RgType
LOGICAL, INTENT(IN),    OPTIONAL  :: Restart
INTEGER, INTENT(OUT),   OPTIONAL  :: NEval
INTEGER, INTENT(IN),    OPTIONAL  :: Job,Key,MaxPts,MinPts
INTEGER, INTENT(IN OUT), OPTIONAL  :: IFAIL
REAL(kind=stnd), INTENT(IN), OPTIONAL :: Tune,EpsAbs,EpsRel
REAL(kind=stnd), INTENT(IN), DIMENSION(:,:,:) :: Vertices
REAL(kind=stnd), INTENT(IN OUT), DIMENSION(:)  :: AbsErr, Value
!
! Named constants
!
INTEGER, PARAMETER ::  NRINFO=1, NIINFO=5
!
!  Local variables
!
INTEGER :: BLOCK,i,Inform,Leval,LJob,LMaxPts,LMinPts,MinCost, &
           NRVERT,NrSub,NRVACA,MAXRGN,RULCLS,STATUS,Tmp
LOGICAL                                     :: EpsAlg,LRestart
REAL(kind=stnd)                             :: LEpsAbs, LEpsRel
REAL(kind=stnd), DIMENSION(:), ALLOCATABLE  :: TmpRWork
INTEGER,    DIMENSION(:), ALLOCATABLE       :: TmpIWork
TYPE(INTEGRATOR_INFO)                       :: CINFO
TYPE(USER_INFO)                             :: UINFO
!-----------------------------------------------------------------------
!
! Check array sizes
!      Array size mismatch results in hard error.
   Inform = 0
   IF (size(rgtype) < numrgn) THEN
       write(unit=*,fmt=*) "Error: size(rgtype) < numrgn"
       Inform = Inform + 1
   END IF
   IF (size(abserr) < numfun) THEN
       write(unit=*,fmt=*) "Error: size(abserr) < numfun"
       Inform = Inform + 1
   END IF
   IF (size(Value) < numfun) THEN
       write(unit=*,fmt=*) "Error: size(Value) < numfun"
       Inform = Inform + 1
   END IF
   IF ((size(vertices,1) /= dimens) .or. (size(vertices,2) /= dimens+1) &
          .or. (size(vertices,3) < numrgn)) THEN
       Inform = Inform + 1
       write(unit=*,fmt=*)"Error: size(vertices) /= (/dimens,dimens+1,numrgn/)"
   END IF
   IF (Inform /= 0) THEN
      WRITE(unit=*,fmt=*)  "Array size mismatch results in hard error."
      STOP   ! "Array size mismatch results in hard error."
   END IF

!-----------------------------------------------------------------------
IF (PRESENT(NEval)) THEN
    NEval = 0
END IF
!-----------------------------------------------------------------------
IF (PRESENT(Job)) THEN
    LJob = Job
    IF (Job == 0) THEN
       CALL CUBATR_CLEAR()
       RETURN
    END IF
 ELSE
    LJob = 1 
END IF
!-----------------------------------------------------------------------
!
! Set optional arguments
!
IF ( PRESENT(Restart)) THEN 
                              LRestart = Restart
                       ELSE 
                              LRestart = .FALSE.
END IF
IF ( PRESENT(Key)) THEN 
                        CINFO%Key = Key
                   ELSE 
                       IF ((ABS(LJob) == 2) .AND. (DIMENS == 1)) THEN
                           CINFO%Key = 2   ! simulate QAGS
                       ELSE
                           CINFO%Key = 0
                       END IF
END IF
!-----------------------------------------------------------------------
!
!  Check input parameters
!
IF ( .NOT. LRestart) THEN
         CALL CHECK(DIMENS,NumFun,NumRgn,RgType,Inform, &
                    LJob,IFAIL,EpsAbs,EpsRel,MinPts,MaxPts,Tune)
ELSE IF ( ALLOCATED(IWork) ) THEN
         CALL CHECK(DIMENS,NumFun,BOTTRH,BOTTIH,IWork,  &
                    Inform,LJob,IFAIL,EpsAbs,EpsRel,MinPts,MaxPts,Tune)
ELSE
   Inform = 4096  ! There is nothing to restart from
END IF
IF (Inform /= 0) THEN 
    CALL Handle_Error(Inform,IFAIL)
    RETURN
END IF
!-----------------------------------------------------------------------
RULCLS = Rule_Cost( DIMENS, RgType(1), CINFO%Key )
MinCost = RULCLS
DO i = 2,NumRgn
   Tmp = Rule_Cost( DIMENS, RgType(i), CINFO%Key)
   RULCLS = max(RULCLS,Tmp)
   MinCost = MinCost + Tmp
END DO
!-----------------------------------------------------------------------
!
! Set optional arguments
!
   
IF ( PRESENT(MinPts)) THEN 
                              LMinPts = MinPts
                      ELSE 
                              LMinPts = 0
END IF
IF ( PRESENT(MaxPts)) THEN 
                              LMaxPts = MaxPts
                      ELSE 
                              LMaxPts = 500*RULCLS
END IF
IF ( PRESENT(Tune)) THEN 
                              CINFO%Tune = Tune
                    ELSE 
                              CINFO%Tune = 1
END IF
IF ( PRESENT(EpsAbs)) THEN 
                              LEpsAbs = EpsAbs
                      ELSE 
                              LEpsAbs = 0
END IF
IF ( PRESENT(EpsRel)) THEN 
                              LEpsRel = EpsRel
                      ELSE 
                              LEpsRel = SQRT(EPSILON(LEpsRel))
END IF
!-----------------------------------------------------------------------
!
! Set other parameters of the Global Adaptive algorithm
!      
!
!  NrSub is an upper limit for the number of subregions after subdivision.
!  This influence memory managment, so don't exagerate here.
IF (DIMENS <= 3) THEN 
                              NrSub = 2**DIMENS
                 ELSE 
                              NrSub = 4
END IF
!
EpsAlg = ( ABS(LJob) == 2 )
CINFO%UNIFORM_SUBDIV = EpsAlg
CINFO%NrSub = NrSub

IF (( ABS(LJob) == 11) .AND. (DIMENS <= 3)) THEN
   ! simulate dcutri and dcutet ; NrSub = 2**DIMENS ; EpsAlg = .FALSE.
   CINFO%UNIFORM_SUBDIV = .TRUE.
END IF

IF ( ABS(LJob) == 12 ) THEN
   ! simulate dcuhre; NrSub = 2 ; EpsAlg = .FALSE.
   CINFO%NrSub = 2
END IF

NRVERT = DIMENS + 1   ! Only cubes and simplices are implemented here.

UINFO%NumFun = NumFun  
UINFO%NumRgn = NumRgn
UINFO%MinPts = LMinPts 
UINFO%MaxPts = LMaxPts
UINFO%EpsAbs = LEpsAbs 
UINFO%EpsRel = LEpsRel
!-----------------------------------------------------------------------
IF (LRestart) THEN
  ! This requires allocating larger arrays and copying
  ! the region collection.
  ALLOCATE(TmpRWork(SIZE(RWork)),STAT=status)
  IF (status /=  0) THEN
     WRITE(unit=*,fmt=*)  "Problem allocating real workspace."
     STOP  !  "Problem allocating real workspace."
  END IF
  ALLOCATE(TmpIWork(SIZE(IWork)),STAT=status)
  IF (status /=  0) THEN
     WRITE(unit=*,fmt=*) "Problem allocating integer workspace."
     STOP  ! "Problem allocating integer workspace."
  END IF
  MAXRGN = DSUSED(IWork)
  CALL DSCOPY(IWork,RWork,TmpIWork,TmpRWork)
ELSE
  MAXRGN = NumRgn
END IF
!-----------------------------------------------------------------------
! NRVACA is the number of regions the global adaptive algorithm
!        removes from the data structure for further processing.
!        In some routines for shared memory parallel machines
!        this is the variable MDIV
NRVACA = 1
! MAXRGN depends on the number of function evalutions
MAXRGN = MAXRGN + 1 + (NrSub-1)*(LMaxPts - RULCLS*NumRgn)/(RULCLS*NrSub)
!
! Compute length of workspace needed.
!
BOTTIH = MAXRGN*(1+NIINFO) + 15 + NRVACA
BLOCK = NRINFO+NRVERT*DIMENS+2*NumFun
IF (NumFun > 1) THEN
    BLOCK = BLOCK + 1
END IF
BOTTRH = MAXRGN*BLOCK
!
! Allocate space for the region collection
!
IF (ALLOCATED(RWork)) THEN
    DEALLOCATE(RWork)
END IF
ALLOCATE(RWork(BOTTRH),STAT=status)
IF (status /=  0) THEN
   WRITE(unit=*,fmt=*) "Problem allocating real workspace."
   STOP  ! "Problem allocating real workspace."
END IF
IF (ALLOCATED(IWork)) THEN
   DEALLOCATE(IWork)
END IF
ALLOCATE(IWork(BOTTIH),STAT=status)
IF (status  /=  0) THEN
   WRITE(unit=*,fmt=*) "Problem allocating integer workspace."
   STOP  ! "Problem allocating integer workspace."
END IF
!
! Initialise region collection
!
IF ( LRestart) THEN
    CALL DSCOPY(TmpIWork,TmpRWork,IWork,RWork)
    DEALLOCATE(TmpIWork,TmpRWork)
ELSE
   IF ( MinCost > LMaxPts ) THEN
       Inform = 128    ! Dit nummer werd al gebruikt !
   ELSE
       CALL DSINIT(DIMENS,NRVERT,NIINFO,NRINFO,NumFun,NRVACA,   &
                   BOTTIH,BOTTRH,IWork,Inform)
   END IF
   IF (Inform /= 0) THEN 
       CALL Handle_Error(Inform,IFAIL)
       RETURN
   END IF
END IF
!-----------------------------------------------------------------------
!
!  Call integration routine
!
If (EpsAlg) THEN
   IF ( LRestart .AND. (PreJob /= ABS(LJob))) THEN
      Inform = 3
   ELSE
      ! Observe that only relevant array sections are passed !
      CALL Global_Adapt_Extrap(DIMENS,CINFO,UINFO,NRVERT,NIINFO,    &
                     NRINFO, Vertices(1:DIMENS,1:NRVERT,1:NUMRGN),  &
                     RgType(1:NUMRGN),Integrand,LRestart,           &
                     Value(1:NUMFUN),AbsErr(1:NUMFUN),LEval,Inform, &
                     RWork,IWork,M)
   END IF
ELSE
   IF ( LRestart .AND. (PreJob /= LJob)) THEN
       IF ( ASSOCIATED(M%RESLA)) THEN
           DEALLOCATE(M%RESLA,M%ERLARG,M%RESULT1,M%ABSERR1,M%RCOPY)
       END IF
       CALL DSPINT(IWork,RWork)
       CALL DSSUM(Value,Abserr,IWork,RWork,Inform)
   END IF
   ! Observe that only relevant array sections are passed !
   CALL Global_Adapt(DIMENS,CINFO,UINFO,NRVERT,NIINFO,NRINFO,    &
                  Vertices(1:DIMENS,1:NRVERT,1:NUMRGN),          &
                  RgType(1:NUMRGN),Integrand,LRestart,           &
                  Value(1:NUMFUN),AbsErr(1:NUMFUN),LEval,Inform, &
                  RWork,IWork)
END IF
IF (PRESENT(NEval)) THEN
   NEval = LEval
END IF
!-----------------------------------------------------------------------
 IF (LJob < 0) THEN
    WRITE(unit=*,fmt=*) "Debug mode: dumping region collection overview."
    CALL DSSTAT(IWork(:),RWork(:))          ! For debugging.
 END IF
!-----------------------------------------------------------------------
! IF ((Inform >= 8) .or. (Inform == 3)) THEN
!    Something went wrong but the data structure remains untouched
!    and so this call can be ignored.
IF ((Inform < 8) .AND. (Inform /= 3)) THEN
   PreJob = ABS(LJob)
END IF
CALL Handle_Error(Inform,IFAIL)
RETURN
END SUBROUTINE CUBATR_X


SUBROUTINE CUBATR_1                                              &
     (DIMENS,Integrand,SVertices,SRgType,SValue,SAbsErr,         &
!   and optional parameters                                      & 
      IFAIL,Neval,EpsAbs,EpsRel,Restart,MaxPts,Key,Job)
!-----------------------------------------------------------------------
!   Input parameters
!   ----------------
!
!   DIMENS Integer.
!          The dimension of the region of integration.
!
!   Integrand
!          Externally declared function for computing all components
!          of the integrand at the given evaluation point.
!          It must have input parameter X:
!              X(1)   The x-coordinate of the evaluation point.
!              X(2)   The y-coordinate of the evaluation point.
!              ...
!              X(DIMENS) The z-coordinate of the evaluation point.
!         and NumFun, the number of components of the integrand.
!         It must be compatible with the following interface:
!           INTERFACE 
!              FUNCTION Integrand(NUMFUN,X) RESULT(Value)
!                USE Precision_Model
!                INTEGER, INTENT(IN) :: NUMFUN
!                REAL(kind=stnd), DIMENSION(:), INTENT(IN) :: X
!                REAL(kind=stnd), DIMENSION(NUMFUN) :: Value
!              END FUNCTION Integrand
!           END INTERFACE
!
!   SVertices
!          Real array of dimension (DIMENS,DIMENS+1).
!          Vertices(1:DIMENS,K) are the x, y, ... coordinates
!          respectively of vertex K of the region, where
!          K = 1,...,DIMENS+1.
!
!   SRgType Integer.
!          RgType describes the type of region L.
!
!   SValue Real.
!          Approximation to the integral if the procedure is restarted.
!
!   SAbsErr Real.
!          Estimate of the absolute error if the procedure is restarted.
!
!   IFAIL  Optional integer argument.
!          This follows the NAG convention:
!          IFAIL = 1 : soft silent error
!                      Control returned to calling program.
!          IFAIL = -1: soft noisy error
!                      Error message is printed.
!                      Control returned to calling program.
!          IFAIL = 0 : hard noisy error
!                      Error message is printed and program is stopped.
!          Default IFAIL = -1.
!
!   EpsAbs Optional real argument.
!          Requested absolute error.
!          Default  EpsAbs = 0.
!
!   EpsRel Optional real argument.
!          Requested relative error.
!          Default EpsRel = sqrt(machine precision).
!
!   Restart Optional boolean argument.
!          If Restart = FALSE, this is the first attempt to compute
!                              the integral.
!          If Restart = TRUE, then we restart a previous attempt.
!          In this case the only parameters for CUBATR that may
!          be changed (with respect to the previous call of CUBATR)
!          are MinPts, MaxPts, EpsAbs, EpsRel, Key and Restart.
!          Default Restart = FALSE.
!
!   MaxPts Optional integer argument.
!          The maximum allowed number of integrand evaluations.
!          Default MaxPts = enough to do 500 subdivisions.
!
!   Key    Optional integer argument.
!          Can be used by Rule_General to choose between several
!          local integration rules.
!          Default Key = 2 if Dimension=1 and extrapolation is used 
!                                        (This corresponds to QAGS)
!          Default Key = 0 otherwise
!
!   Job    Optional integer argument.
!          If |Job| = 0, then nothing will be done except freeing all
!                        allocated memory.
!                        This is usefull after a call of CUBATR if no
!                        Restart will be done later and memory usage
!                        might become an issue later.
!                        Equivalently, one can call CUBATR()
!                        without any arguments.
!                   = 1, the global adaptive algorithm is called
!                   = 2, extrapolation using the epsilon algorithm is used.
!                   = 11, a region will be divided in 2**DIMENS subregions
!                        and the global adaptive algorithm is called.
!                        In combination with Key=0, this resembles DUCTRI and DCUTET.
!                   = 12, a region will be divided in 2 subregions
!                        and the global adaptive algorithm is called.
!                        In combination with Key=3 or 4, this resembles DCUHRE.
!          If Job < 0, then an overview of the Region Collection is dumped.
!          This will create the files tmp_integerstore and tmp_realstore.
!          Default Job = 1.
!
!   Output parameters
!   -----------------
!
!   SValue Real.
!          Approximation to the integral
!
!   AbsErr Real.
!          Estimate of the absolute error.
!
!   NEval  Optional Integer.
!          Number of integrand evaluations used by CUBATR for this call.
!
!   IFAIL  Optional Integer.
!          IFAIL = 0 for normal exit.
!
!            AbsErr(K) <=  EpsAbs or
!            AbsErr(K) <=  ABS(Value(K))*EpsRel with MaxPts or less
!            function evaluations for all values of K,
!            1 <= K <= NumFun .
!
!          IFAIL = 1 if MaxPts was too small to obtain the required
!            accuracy. In this case Global_Adapt returns values of
!            Value with estimated absolute errors AbsErr.
!
!          IFAIL > 1 in more serious case of trouble.
!-----------------------------------------------------------------------
!
! Global variables
!
INTERFACE 
   FUNCTION Integrand(NUMFUN,X) RESULT(Value)
      USE Precision_Model
      INTEGER, INTENT(IN) :: NUMFUN
      REAL(kind=stnd), DIMENSION(:), INTENT(IN) :: X
      REAL(kind=stnd), DIMENSION(NUMFUN) :: Value
   END FUNCTION Integrand
END INTERFACE
LOGICAL, OPTIONAL, INTENT(IN)                 :: Restart
INTEGER, INTENT(IN)                           :: DIMENS,SRgType
INTEGER, INTENT(OUT), OPTIONAL                :: NEval
INTEGER, INTENT(IN), OPTIONAL                 :: Key,MaxPts,Job
INTEGER, INTENT(IN OUT), OPTIONAL             :: IFAIL
REAL(kind=stnd), INTENT(IN), OPTIONAL         :: EpsAbs,EpsRel
REAL(kind=stnd), DIMENSION(:,:), INTENT(IN)   :: SVertices
REAL(kind=stnd), INTENT(IN OUT)               :: SValue,SAbsErr
!
! Local variables
!
INTEGER, DIMENSION(1)                         :: RgType
REAL(kind=stnd), DIMENSION(1)                 :: Value, AbsErr
REAL(kind=stnd), DIMENSION(DIMENS,DIMENS+1,1) :: Vertices
!-------------------
RgType(1) = SRgType
Vertices(:,:,1) = SVertices
IF (PRESENT(Restart)) THEN
    IF ( Restart ) THEN
       Value(1) = SValue
       AbsErr(1) = SAbsErr
    END IF
END IF
CALL CUBATR             &
  (DIMENS,1,Integrand,1,Vertices,RgType,Value,AbsErr,        &
   !   and optional parameters                               &
   ifail=IFAIL,neval=Neval,epsabs=EpsAbs,epsrel=EpsRel,      &
   restart=Restart,maxpts=MaxPts,key=key,job=Job)
SValue = Value(1)
SAbsErr = AbsErr(1)
RETURN
END SUBROUTINE CUBATR_1

SUBROUTINE CUBATR_CLEAR()
  IF ( ALLOCATED(Iwork) ) THEN
     DEALLOCATE(RWork,IWork)
  END IF
  IF ( ASSOCIATED(M%RESLA)) THEN
     DEALLOCATE(M%RESLA,M%ERLARG,M%RESULT1,M%ABSERR1,M%RCOPY)
  END IF
  PreJob = 0
  RETURN
END SUBROUTINE CUBATR_CLEAR

SUBROUTINE CUBPACK_INFO()

REAL(kind=stnd) :: x=1.0e-30  ! lowest accuracy of cubature formula constants

print *," ---------------------------------------------------------------"
print *,"              CUBPACK information"
print *,"              -------------------"
print *," The model for real numbers in the current installed version,"
print *," obtained with the declaration REAL(KIND=stnd), has the"
print *," following characteristics:"
print *,"                base =  ",radix(x)
print *," digits in this base = ", digits(x)
! print *," highest exponent = ", maxexponent(x)
! print *," lowest exponent = "minexponent(x), "(normalized numbers)
print *," This implies:"
print *,"            machine epsilon = ",epsilon(x)
print *,"        largest real number = ", huge(x)
print *," smallest normalized number = ", tiny(x)
print *
print *," The lowest relative error that may be obtained with this"
print "(""  version is about "",G8.2)",max(50*epsilon(x),x)
print *," Asking for lower error will push the routine to use the"
print *," maximal number of function evaluations it is allowed."
print *
print *
print *," This version accepts a collection of hyper-rectangles"
print *,"   (and parallelepipeds) and simplices as integration regions."
print *," Extrapolation using the epsilon-algorithm is available"
print *,"      for dimensions 1, 2 and 3."
print *," The following values of KEY give different integration rules:"
print *," - finite interval: KEY = 1, 2, 3, 4, 5."
print *,"      KEY < 1 defaults to 1; KEY > 5 defaults to 5."
print *," - n-cube:          KEY = 3, 4       uses rule of degree 2*KEY+1"
print *,"           otherwise, uses for a square a rule of degree 13"
print *,"                                 3-cube a rule of degree 11"
print *,"                                        a rule of degree  7"
print *," - n-simplex:       KEY = 1, 2, 3, 4 uses rule of degree 2*KEY+1" 
print *,"         otherwise, uses for a triangle a rule of degree 13"
print *,"                            tetrahedron a rule of degree  8"
print *,"                                        a rule of degree  7"
print *," KEY = 0 corresponds to our preferred choice."
print *," ---------------------------------------------------------------"

RETURN
END SUBROUTINE CUBPACK_INFO

END Module CUI
