      SUBROUTINE initial(ns,m,emax,x,hx,hpx,lb,xlb,ub,xub,ifault,
     1             iwv,rwv,eval)
C                                                                       
C This subroutine takes as input the number of starting values m        
C  and the starting values x(i),hx(i),hpx(i)  i=1,m                     
C As output we have pointer ipt along with ilow and ihigh and the lower 
C and upper hulls defined  by z,hz,scum,cu,hulb,huub stored in working 
C vectors iwv and rwv                   
C Ifault detects wrong starting points or non-concavity                 
C                                                                       
      INTEGER ns,nn,m,ilow,ihigh,ifault,i,iwv(*)
      INTEGER iipt,iz,ihuz,iscum,ix,ihx,ihpx
      LOGICAL ub,lb,horiz
      DOUBLE PRECISION xlb,xub,emax,x(*),hx(*),hpx(*),rwv(*),expon
      DOUBLE PRECISION hulb,huub,eps,cu,alcu,huzmax,zlog,zmax
      DOUBLE PRECISION a,b
      external eval
      zlog(a)=dlog(a)
      zmax(a,b)=dmax1(a,b)
C
C DESCRIPTION OF PARAMETERS and place of storage
C
C     lb   iwv(5) : boolean indicating if there is a lower bound to the
C                    domain
C     ub   iwv(6) : boolean indicating if there is an upper bound
C     xlb  rwv(8) : value of the lower bound
C     xub  rwv(9) : value of the upper bound
C     emax rwv(3) : large value for which it is possible to compute
C                   an exponential, eps=exp(-emax) is taken as a small
C                   value used to test for numerical unstability
C     m    iwv(4) : number of starting points
C     ns   iwv(3) : maximum number of points defining the hulls
C     x    rwv(ix+1)  : vector containing the abscissae of the starting
C                       points
C     hx   rwv(ihx+1) : vector containing the ordinates 
C     hpx  rwv(ihpx+1): vector containing the derivatives
C     ifault      : diagnostic
C     iwv,rwv     : integer and real working vectors
C                                                                       
      eps=expon(-emax,emax)
      ifault=0                                                          
      ilow=1                                                            
      ihigh=1                                                           
      nn=ns+1
C at least one starting point
      IF (m.LT.1) ifault=1
      huzmax=hx(1)                                                      
      IF (.NOT.ub) xub=0.0
      IF (.NOT.lb) xlb=0.0
      hulb=(xlb-x(1))*hpx(1)+hx(1)                                      
      huub=(xub-x(1))*hpx(1)+hx(1)                                      
C if bounded on both sides                                              
      IF ((ub).AND.(lb)) THEN                               
        huzmax=zmax(huub,hulb)                                         
        horiz=(abs(hpx(1)).LT.eps)
        IF (horiz) THEN
          cu=expon((huub+hulb)*0.5-huzmax,emax)*(xub-xlb)
        ELSE
          cu=expon(huub-huzmax,emax)*(1-expon(hulb-huub,emax))/hpx(1)
        ENDIF
      ELSE IF ((ub).AND.(.NOT.lb))THEN                           
C if bounded on the right and unbounded on the left 
        huzmax=huub                                                     
        cu=1.0/hpx(1)
      ELSE IF ((.NOT.ub).AND.(lb))THEN                           
C if bounded on the left and unbounded on the right                     
        huzmax=hulb                                                     
        cu=-1.0/hpx(1)
C if unbounded at least 2 starting points
      ELSE
        cu=0.0
        IF (m.LT.2) ifault=1
      ENDIF                                                             
      IF (cu.GT.0.0) alcu=zlog(cu)
C set pointers 
      iipt=6
      iz=9
      ihuz=nn+iz
      iscum=nn+ihuz
      ix=nn+iscum
      ihx=nn+ix
      ihpx=nn+ihx
C store values in working vectors
      iwv(1)=ilow
      iwv(2)=ihigh
      iwv(3)=ns
      iwv(4)=1
      IF (lb) THEN
        iwv(5)=1
      ELSE
        iwv(5)=0
      ENDIF 
      IF (ub) THEN
        iwv(6)=1
      ELSE
        iwv(6)=0
      ENDIF 
      IF (ns.LT.m) ifault=2
      iwv(iipt+1)=0                                                     
      rwv(1)=hulb
      rwv(2)=huub
      rwv(3)=emax
      rwv(4)=eps
      rwv(5)=cu
      rwv(6)=alcu
      rwv(7)=huzmax
      rwv(8)=xlb
      rwv(9)=xub
      rwv(iscum+1)=1.0
      DO 9 i=1,m
        rwv(ix+i)=x(i)
        rwv(ihx+i)=hx(i)
        rwv(ihpx+i)=hpx(i)
  9   CONTINUE
C create lower and upper hulls                                          
      i=1                                                               
 10   IF (i.LT.m)THEN                                                   
         CALL update(iwv(4),iwv(1),iwv(2),iwv(iipt+1),rwv(iscum+1),
     +           rwv(5),rwv(ix+1),rwv(ihx+1),rwv(ihpx+1),rwv(iz+1),
     +           rwv(ihuz+1),rwv(7),rwv(3),lb,rwv(8),rwv(1),ub,rwv(9),
     +           rwv(2),ifault,rwv(4),rwv(6),eval) 
         i=iwv(4)
         IF (ifault.NE.0) RETURN
      GOTO 10                                                           
      ENDIF                                                             
C test for wrong starting points                                        
      IF ((.NOT.lb).AND.(hpx(iwv(1)).LT.eps)) ifault=3
      IF ((.NOT.ub).AND.(hpx(iwv(2)).GT.-eps)) ifault=4
      RETURN
      END
C                                                                       
C***********************************************************************
C                                                                       
      SUBROUTINE sample(iwv,rwv,eval,beta,ifault,iseed)
      DOUBLE PRECISION beta,rwv(*)
      INTEGER iipt,iz,ns,nn,ihuz,iscum,ix,ihx,ihpx,ifault,iwv(*)
      INTEGER iseed
      LOGICAL ub,lb
C
C     set pointers
C
      iipt=6
      iz=9
      ns=iwv(3)
      nn=ns+1
      ihuz=nn+iz
      iscum=nn+ihuz
      ix=nn+iscum
      ihx=nn+ix
      ihpx=nn+ihx
      lb=.FALSE.
      ub=.FALSE.
      IF (iwv(5).EQ.1) lb=.TRUE.
      IF (iwv(6).EQ.1) ub=.TRUE.
C
C     call sampling subroutine
C
      CALL spl1(ns,iwv(4),iwv(1),iwv(2),iwv(iipt+1),rwv(iscum+1),
     +  rwv(5),rwv(ix+1),rwv(ihx+1),rwv(ihpx+1),rwv(iz+1),rwv(ihuz+1),  
     +  rwv(7),lb,rwv(8),rwv(1),ub,rwv(9),rwv(2),eval,beta,ifault,
     +  rwv(3),rwv(4),rwv(6),iseed)    
      RETURN
      END
C                                                                       
C***********************************************************************
C                                                                       
      SUBROUTINE spl1(ns,n,ilow,ihigh,ipt,scum,cu,x,hx,hpx,z,huz,    
     +  huzmax,lb,xlb,hulb,ub,xub,huub,eval,beta,ifault,emax,eps,alcu,
     +  iseed)  
C                                                                       
C this subroutine performs the adaptive rejection sampling, it calls 
C subroutine splhull to sample from the upper hull ,if the sampling 
C involves a function evaluation it calls the updating subroutine
C ifault is a diagnostic of any problem: non concavity, 0 random number 
C    or numerical imprecision
C                                                                       
      INTEGER ns,n,ilow,ihigh,ifault,i,j,ipt(*),n1,l,gencnt
      INTEGER iseed
      LOGICAL ub,lb,sampld
      DOUBLE PRECISION z(*),huz(*),x(*),hx(*),hpx(*),scum(*)
      DOUBLE PRECISION xlb,xub,emax,u1,u2,beta,alu1,fx
      DOUBLE PRECISION hulb,huub,eps,cu,alcu,huzmax,alhl,alhu
      DOUBLE PRECISION zlog,a
      DOUBLE PRECISION DRNUNF
      EXTERNAL eval,DRNUNF
      zlog(a)=dlog(a)
C                                                                       
      sampld=.FALSE.
      gencnt=0.0                                                    
10    IF (.NOT.sampld) THEN
        gencnt=gencnt+1
C        PRINT*, 'Number of generations =',gencnt
        u2=DRNUNF() 
C test for zero random number                                         
        IF (u2.EQ.0.0) THEN                                             
          ifault=6
          RETURN
        ENDIF                                                           
      CALL splhull(u2,ipt,ilow,lb,xlb,hulb,huzmax,alcu,x,hx,hpx,
     +                              z,huz,scum,eps,emax,beta,i,j)
C      PRINT*, 'Exit from splhull'
C sample u1 to compute rejection                                        
        u1=DRNUNF()                                                    
        IF (u1.EQ.0.0) ifault=6                                         
        alu1=zlog(u1)                                                   
C compute alhu: upper hull at point u1                                  
        alhu=hpx(i)*(beta-x(i))+hx(i)-huzmax                          
        IF ((beta.GT.x(ilow)).AND.(beta.LT.x(ihigh))) THEN              
C compute alhl: value of the lower hull at point u1                     
          IF (beta.GT.x(i)) THEN                                        
            j=i
            i=ipt(i)
          ENDIF
          alhl=hx(i)+(beta-x(i))*(hx(i)-hx(j))/(x(i)-x(j))-huzmax 
C squeezing test                                                        
          IF ((alhl-alhu).GT.alu1) THEN                                 
             sampld=.TRUE.                                              
          ENDIF                                                         
        ENDIF                                                           
C if not sampled evaluate the function ,do the rejection test and update
        IF (.NOT.sampld) THEN                                           
          n1=n+1                                                        
          x(n1)=beta                                                    
          CALL eval(x(n1),hx(n1),hpx(n1))                               
          fx=hx(n1)-huzmax                                              
          IF (alu1.lt.(fx-alhu)) sampld=.TRUE.                          
C update while the number of points defining the hulls is lower than ns
          IF (n.LT.ns)  THEN
C            PRINT*, 'Entering update'
            CALL update(n,ilow,ihigh,ipt,scum,cu,x,hx,hpx,z,huz,huzmax, 
     +           emax,lb,xlb,hulb,ub,xub,huub,ifault,eps,alcu,eval) 
C            PRINT*, 'Exit from update'
          ENDIF
          IF (ifault.NE.0) RETURN                   
        ENDIF                                                           
        GOTO 10                                                         
      ENDIF            
C      PRINT*, 'Exit from spl1'                                             
      RETURN                                                            
      END                                                               
C                                                                       
C***********************************************************************
C                                                                       
      SUBROUTINE splhull(u2,ipt,ilow,lb,xlb,hulb,huzmax,alcu,x,hx,hpx,
     +                              z,huz,scum,eps,emax,beta,i,j)
C
C this subroutine samples beta from the normalised upper hull
C
      DOUBLE PRECISION z(*),huz(*),x(*),hx(*),hpx(*),scum(*)
      DOUBLE PRECISION expon,xlb,emax,eh,beta,u2,a,logdu,logtg,sign
      DOUBLE PRECISION hulb,eps,alcu,huzmax,zlog
      INTEGER ilow,i,j,ipt(*)
      LOGICAL lb,horiz
      zlog(a)=dlog(a)
C
        i=ilow                                                          
C
C find from which exponential piece you sample                        
 20   IF (u2.gt.scum(i))THEN                                          
        j=i                                                           
        i=ipt(i)                                                      
        GOTO 20                                                       
      ENDIF                                                           
      IF (i.eq.ilow) THEN                                             
C sample below z(ilow),depending on the existence of a lower bound  
        IF (lb) THEN                                            
          eh=hulb-huzmax-alcu                                         
          horiz=(abs(hpx(ilow)).LT.eps)
          IF (horiz) THEN
            beta=xlb+u2*expon(-eh,emax)
          ELSE
           sign=abs(hpx(i))/hpx(i)
           logtg=zlog(abs(hpx(i)))
           logdu=zlog(u2)
           eh=logdu+logtg-eh
           IF (eh.LT.emax) THEN
            beta=xlb+zlog(1.0+sign*expon(eh,emax))/hpx(i)
           ELSE
            beta=xlb+eh/hpx(i)
           ENDIF                                                       
          ENDIF
        ELSE                                                          
C     hpx(i) must be positive , x(ilow) is left of the mode
          beta=(zlog(hpx(i)*u2)+alcu-hx(i)+x(i)*hpx(i)+huzmax)/hpx(i)   
        ENDIF                                                         
      ELSE                                                            
C   sample above(j)                                                   
        eh=huz(j)-huzmax-alcu                                         
        horiz=(abs(hpx(i)).LT.eps)
        IF (horiz) THEN
         beta=z(j)+(u2-scum(j))*expon(-eh,emax)
        ELSE
         sign=abs(hpx(i))/hpx(i)
         logtg=zlog(abs(hpx(i)))
         logdu=zlog(u2-scum(j))
         eh=logdu+logtg-eh
         IF (eh.LT.emax) THEN
          beta=z(j)+(zlog(1.0+sign*expon(eh,emax)))/hpx(i)  
         ELSE                                                        
          beta=z(j)+eh/hpx(i)
         ENDIF
        ENDIF
      ENDIF                                                           
      RETURN
      END
C                                                                       
C***********************************************************************
      SUBROUTINE intersection(x1,y1,yp1,x2,y2,yp2,z1,hz1,eps,ifault)    
C                                                                       
C computes the intersection (z1,hz1) between 2 tangents defined by
C   x1,y1,yp1 and x2,y2,yp2
C
      DOUBLE PRECISION x1,y1,yp1,x2,y2,yp2,z1,hz1,eps,y12,y21,dh
      INTEGER ifault
C
C first test for non-concavity                                          
      y12=y1+yp1*(x2-x1)
      y21=y2+yp2*(x1-x2)
      IF ((y21.lt.y1).OR.(y12.LT.y2)) THEN                                 
         ifault=5                                                       
         RETURN                                                         
      ENDIF                                                             
      dh=yp2-yp1                                                        
C  IF the lines are nearly parallel,
C  the intersection is taken at the midpoint
      IF (abs(dh).LE.eps)THEN                                           
        z1=0.5*(x1+x2)                                                  
        hz1=0.5*(y1+y2)                                                 
C  Else compute from the left or the right for greater numerical
C       precision
      ELSE IF (abs(yp1).LT.abs(yp2)) THEN                               
        z1=x2+(y1-y2+yp1*(x2-x1))/dh                                    
        hz1=yp1*(z1-x1)+y1                                              
      ELSE                                                              
        z1=x1+(y1-y2+yp2*(x2-x1))/dh                                    
        hz1=yp2*(z1-x2)+y2                                              
      ENDIF                                                             
C  test for misbehaviour due to numerical imprecision
      IF ((z1.LT.x1).OR.(z1.GT.x2)) ifault=7
      RETURN                                                            
      END                                                               
C                                                                       
C***********************************************************************
C                                                                       
      SUBROUTINE update(n,ilow,ihigh,ipt,scum,cu,x,hx,hpx,z,huz,        
     +  huzmax,emax,lb,xlb,hulb,ub,xub,huub,ifault,eps,alcu,eval)
C                                                                       
C this subroutine increments n and updates all the parameters which
C define the lower and the upper hull
C                                                                       
      INTEGER n,ilow,ihigh,ifault,i,j,ipt(*)
      LOGICAL ub,lb,horiz
      DOUBLE PRECISION xlb,xub,x(*),hx(*),hpx(*),z(*),huz(*),scum(*)
      DOUBLE PRECISION hulb,huub,eps,cu,alcu,huzmax,emax,expon,dh,u
      DOUBLE PRECISION a,b,zlog,zmax
      EXTERNAL eval
      zlog(a)=dlog(a)
      zmax(a,b)=dmax1(a,b)
C
C DESCRIPTION OF PARAMETERS and place of storage
C
C     ilow iwv(1)    : index of the smallest x(i)
C     ihigh iwv(2)   : index of the largest x(i)
C     n    iwv(4)    : number of points defining the hulls
C     ipt  iwv(iipt) : pointer array:  ipt(i) is the index of the x(.) 
C                      immediately larger than x(i)
C     hulb rwv(1)    : value of the upper hull at xlb
C     huub rwv(2)    : value of the upper hull at xub
C     cu   rwv(5)    : integral of the exponentiated upper hull divided
C                      by exp(huzmax)
C     alcu rwv(6)    : logarithm of cu
C     huzmax rwv(7)  : maximum of huz(i); i=1,n
C     z    rwv(iz+1) : z(i) is the abscissa of the intersection between
C                      the tangents at x(i) and x(ipt(i))
C     huz  rwv(ihuz+1): huz(i) is the ordinate of the intersection
C                        defined above
C     scum rwv(iscum): scum(i) is the cumulative probability of the 
C                      normalised exponential of the upper hull 
C                      calculated at z(i)
C     eps  rwv(4)    : =exp(-emax) a very small number 
C
      n=n+1                                                             
C update z,huz and ipt                                                  
      IF (x(n).LT.x(ilow)) THEN                                         
C insert x(n) below x(ilow)                                             
C   test for non-concavity                                              
        IF (hpx(ilow).GT.hpx(n)) ifault=5 
        ipt(n)=ilow
        CALL intersection(x(n),hx(n),hpx(n),x(ilow),hx(ilow),hpx(ilow), 
     +         z(n),huz(n),eps,ifault)                               
        IF (ifault.NE.0) RETURN                     
        IF (lb) hulb=hpx(n)*(xlb-x(n))+hx(n)                      
        ilow=n                                                          
      ELSE                                                              
        i=ilow                                                          
        j=i                                                             
C find where to insert x(n)                                             
  10    IF ((x(n).GE.x(i)).AND.(ipt(i).ne.0))THEN                       
          j=i                                                           
          i=ipt(i)                                                      
          GOTO 10                                                       
         ENDIF                                                          
        IF (x(n).GE.x(i)) THEN                                          
C insert above x(ihigh)                                                 
C   test for non-concavity                                              
          IF (hpx(i).LT.hpx(n)) ifault=5
          ihigh=n                                                       
          ipt(i)=n                                                      
          ipt(n)=0                                                      
          CALL intersection(x(i),hx(i),hpx(i),x(n),hx(n),hpx(n),        
     +         z(i),huz(i),eps,ifault)                               
          IF (ifault.NE.0) RETURN                   
          huub=hpx(n)*(xub-x(n))+hx(n)                                  
          z(n)=0.0                                                      
          huz(n)=0.0
        ELSE                                                            
C insert x(n) between x(j) and x(i)                                     
C   test for non-concavity                                              
          IF ((hpx(j).LT.hpx(n)).OR.(hpx(i).GT.hpx(n))) ifault=5
          ipt(j)=n                                                      
          ipt(n)=i                                                      
C     insert z(j) between x(j) and x(n)                                 
          CALL intersection(x(j),hx(j),hpx(j),x(n),hx(n),hpx(n),        
     +         z(j),huz(j),eps,ifault)                               
          IF (ifault.NE.0) RETURN                   
C     insert z(n) between x(n) and x(i)                                 
          CALL intersection(x(n),hx(n),hpx(n),x(i),hx(i),hpx(i),        
     +         z(n),huz(n),eps,ifault)                               
          IF (ifault.NE.0) RETURN                   
        ENDIF                                                           
      ENDIF                                                             
C update huzmax                                                         
      j=ilow                                                            
      i=ipt(j)                                                          
      huzmax=huz(j)                                                     
 20   IF ((huz(j).LT.huz(i)).AND.(ipt(i).NE.0))THEN                     
        j=i                                                             
        i=ipt(i)                                                        
        huzmax=zmax(huzmax,huz(j))                                     
        GOTO 20                                                         
      ENDIF                                                             
      IF (lb) huzmax=zmax(huzmax,hulb)
      IF (ub) huzmax=zmax(huzmax,huub)
C update cu                                                             
C  scum receives area below exponentiated upper hull left of z(i)       
      i=ilow                                                            
      horiz=(abs(hpx(ilow)).LT.eps)
      IF ((.NOT.lb).AND.(.NOT.horiz)) THEN
        cu=expon(huz(i)-huzmax,emax)/hpx(i)                             
      ELSE IF (lb.AND.horiz) THEN
        cu=(z(ilow)-xlb)*expon(hulb-huzmax,emax)
      ELSE IF (lb.AND.(.NOT.horiz)) THEN
        dh=hulb-huz(i)
        IF (dh.GT.emax) THEN
          cu=-expon(hulb-huzmax,emax)/hpx(i)
        ELSE
          cu=expon(huz(i)-huzmax,emax)*(1-expon(dh,emax))/hpx(i)
        ENDIF
      ELSE
        cu=0
      ENDIF                                                             
      scum(i)=cu                                                        
      j=i                                                               
      i=ipt(i)                                                          
 30   IF (ipt(i).NE.0)THEN                                              
        dh=huz(j)-huz(i)                                                
        horiz=(abs(hpx(i)).LT.eps)
        IF (horiz) THEN
          cu=cu+(z(i)-z(j))*expon((huz(i)+huz(j))*0.5-huzmax,emax)
        ELSE
          IF (dh.LT.emax) THEN                                          
            cu=cu+expon(huz(i)-huzmax,emax)*(1-expon(dh,emax))/hpx(i) 
          ELSE                                                          
            cu=cu-expon(huz(j)-huzmax,emax)/hpx(i)                      
          ENDIF                                                         
        ENDIF
        j=i                                                             
        i=ipt(i)                                                        
        scum(j)=cu                                                      
        GOTO 30                                                         
      ENDIF                                                             
      horiz=(abs(hpx(i)).LT.eps)
C if the derivative is very small the tangent is nearly horizontal
      IF (.NOT.(ub.OR.horiz)) THEN
        cu=cu-expon(huz(j)-huzmax,emax)/hpx(i)                          
      ELSE IF (ub.AND.horiz) THEN
        cu=cu+(xub-x(i))*expon((huub+hx(i))*0.5-huzmax,emax)
      ELSE IF (ub.AND.(.NOT.horiz)) THEN
        dh=huz(j)-huub
        IF (dh.GT.emax) THEN
          cu=cu-expon(huz(j)-huzmax,emax)/hpx(i)
        ELSE
          cu=cu+expon(huub-huzmax,emax)*(1-expon(dh,emax))/hpx(i)        
        ENDIF
      ENDIF                                                             
      scum(i)=cu                                                        
      IF (cu.GT.0) alcu=zlog(cu)
C normalize scum to obtain a cumulative probability while excluding     
C    unnecessary points                                                 
      i=ilow                                                            
      u=(cu-scum(i))/cu                                                 
      IF ((u.EQ.1.0).AND.(hpx(ipt(i)).GT.0)) THEN
        ilow=ipt(i)                                                     
        scum(i)=0.0                                                     
      ELSE                                                              
        scum(i)=1.0-u                                                   
      ENDIF                                                             
      j=i                                                               
      i=ipt(i)                                                          
 40   IF (ipt(i).NE.0) THEN                                             
        j=i                                                             
        i=ipt(i)                                                        
        u=(cu-scum(j))/cu                                               
        IF ((u.EQ.1.0).AND.(hpx(i).GT.0)) THEN 
          ilow=i                                                        
        ELSE                                                            
          scum(j)=1.0-u                                                 
        ENDIF                                                           
        GOTO 40                                                         
      ENDIF                                                             
      scum(i)=1.0
      IF (ub) huub=hpx(ihigh)*(xub-x(ihigh))+hx(ihigh)
      IF (lb) hulb=hpx(ilow)*(xlb-x(ilow))+hx(ilow)                     
C      do 100 ii=1,n
C      100   write(7,110) ii,z(ii),scum(ii)
C      110  format('i= ',i2,2x,2(2x,E10.4))
C       write(7,*) ' huzmax = ',huzmax
      RETURN                                                            
      END                                                               
C                                                                       
C***********************************************************************
C                                                                       
      DOUBLE PRECISION FUNCTION expon(x,emax)   
C                                                                       
C performs an exponential without underflow                             
C                                                                       
      DOUBLE PRECISION x,emax,zexp
      zexp(x)=dexp(x)
C
      IF (x.LT.-emax) THEN                                          
        expon=0.0                                                       
      ELSE                                                              
        expon=zexp(min(x,emax))       
      ENDIF                                                             
      RETURN                                                            
      END                                                               
