c      This is the updated version of the power simulation program in Fortran. March, 1999
c      Please send you comments to    minggao@cuhk.edu.hk
c
c  To run the program, download and compile this file with a Fortran compiler, 
c  download the sample input file and then run the program using the sample input 
c  file as a standard input file (or $a.out<dsspe2.inp). You may change the 
c  parameters defined in the sample input file according to the simple rules 
c  explained in the sample input file. 
c 
c      The major difference between this program (dssp2.f) and the previous 
c      version (dssp.f) is that the distributions are piecewise exponential function 
c      rather than piecewise linear function specified by the values at the 
c      l grid points. 
c      You may also specify the alternative distribution function by specify the 
c      hazard ratios over some intervals. 
c      Four boundary generating
c      method are allowed incuding a user specified boundary. For more details
c      of the program see the paper by Gu and Lai (1998) whic can be downloaded 
c      at the Minggao Gu's Website. 
c
c      This program is an alternative to the program of Jerry Halpern
c      and Byron Wm. Byron, Jr. which simulates clinical trials.
c
c      k-----------Number of looks, upto nk (=15) looks are permitted.
c      ih----------Number of Statistics in the self's Beta family to be specified
c                     in one simulation. The upper bound is nih (=5)
c      accut------Accrual time. Uniform accrual rate is assumed.
c      accrat-----Accrual rate, no of patients in 1 unit.
c      tlook(nk)------Data monitoring times.
c      rho(nih,1)----- rho's in the Beta family of statistics. 
c      rho(nih,2)----- tau's in the Beta family of statistics. 
c      v(nk)-----Theoritical variance of different statistics at different looking
c                     times. Calculated by numerical integration with delta=h=.01.
c      b---------the b's in the first k-1 boundary points, Haybittle (1977).
c      binput----The initial input of b
c      stat(1:ih)------The value of ih  statistics.
c      x(m)----life-times. Upto m life-times is permitted.
c      y(m)----Entry times.
c      z(m)----covariate, 0 for baseline; 1 for alternative.
c      mm-----is the sample size.
c      md-----number of death in one simulation.
c      n---------Number of simulation.
c      stp(1:ih)----Standard deviation of above power.
c      iseed-----Random seed for simulation.
c      pn-------Type-I error specified by the user.
c
c      l------- Number of intervals to specify the baseline survival function <=nl(=20).
c      l2------- Number of intervals to specify the alternative survival function <=nl(=20).
c      t(0:l)-------The grid points to specify the baseline survival function
c      surt(0:l)---The value of the baseline survival function at the grid
c
c      t2(0:l)-------The grid points to specify the alternative survival function
c      surt2(0:l)---The value of the alternative survival function at the grid
c
c      lc1-----Number of intervals to specify the censoring distribution for group 1
c      tc1(0:lc1)-----The grid points for the censoring survival function of group 1
c      sc1(0:lc1)-----The value at the grid points
c
c      lc2, tc2, sc2 are similar values for specifying censoring 
c      distribution fro group 2
c
c      ihaz------Number of intevals to specify the hazard ratio. <=nl (=20
c      rt(1:ihaz)------ The end points of intervals of hazard ratios
c      rhaz(1:ihaz)---- The hazard ratios.
c
c      jj-------The number of choiceis of the method.
c
c      ld1-------The number of intervals to specify the drop-in rates
c      ld2-------The number of intervals to specify the drop-out rates
c      td1(1:ld1)----the actual time points for drop-in
c      td2(1:ld2)----the actual time points for drop-out
c      rdin(1:ld1)----the droup in rate(s)
c      rdout(1:ld2)----the droup out rate(s)
c
       implicit none
       integer m, nih, nl, nk 
       parameter(m=7100,nih=5,nl=20,nk=15)
       real tlook(nk), rho(nih,2), v(nk), rho1, rho2
       real e(0:nk),bb(nk),va(nk),bbv(nk),cc(nk)
       real stat(nih), var(nih), z(m), pn, futility
       real x(m),y(m),xx(m),zz(m), c(m)
       real sumtime(nih), v0(nih), v1(nih)
       real t(0:nl), surt(0:nl), rt(nl), rhaz(nl), t2(0:nl), surt2(0:nl)
       integer idelta(m), icount(nih), idone(nih)
       integer jj, ik, iseed
       integer n, k, j, i, nside, ld1, ld2
       integer ih, ihaz, l, md, ii, mm, ilook, isum, ms, l2
       character*12 infile, outfile
       character char
       real be(nih,nk), ve(nih,nk) 
       real b, binput, tmp, tmpp, accut, accrat, tlk, haznot
       real ran3, pp, rn, phi, td1(nl), td2(nl), rdin(nl), rdout(nl)
       integer lc1, lc2, i_usefunction, nbtoosmall
       real tc1(0:nl), tc2(0:nl), sc1(0:nl), sc2(0:nl)
c
       common l, t, surt, ihaz, rt, rhaz, l2, t2, surt2
       common lc1, lc2, tc1, tc2, sc1, sc2
       read(5,*)
       read(5,*)  
       read(5,*)  
       read(5,*)  
       read(5,*)  
       read(5,*)  
       read(5,*)  
       read(5,*)  
       read(5,*)  
       read(5,*)  
       read(5,*) jj
       if ( jj .eq. 0 ) then
	  jj = 2
	  write(6,*) ' The Lan-Demets method is chosen.'
       endif
       read(5,*)
       read(5,*)
       read(5,*) outfile
c       if ( outfile .eq. " ") then
c       outfile="out.log"
c       endif
       open(4,file=outfile)
       read(5,*)
       read(5,*)
       read(5,*) iseed
       if( iseed .eq. 0 ) then
	  iseed=-123457
	  write(6,*) ' The seed is set to -123457.'
       endif
       read(5,*)
       read(5,*)
       read(5,*) k
       if ( k .le. 1) then
	  k=1
	  write(6,*) ' You have specified the non-sequential design.'
       endif 
       read(5,*)
       read(5,*)
       read(5,*) (tlook(i),i=1,k)
       if( tlook(k) .le. 0 ) then
	  tlook(k)=10.0
	  write(6,*) ' The last look has been set at 10.'
       endif
       read(5,*)
       read(5,*)
       read(5,*) ih
       if( ih .eq. 0) then
	  ih=1
	  write(6,*) ' One test statistic is specified.'
       endif
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*) (rho(i,1),i=1,ih)
       read(5,*)
       read(5,*)
       read(5,*) (rho(i,2),i=1,ih)
c  Input the v0 and v1 for the smallest information and the largest information
c   to run tests
c  if v0(i) is 0, then later it will taken to be 1/2 of the null
c  variance at the first look. 
       read(5,*)
c       read(5,*)
c       read(5,*)
c       read(5,*)
c       read(5,*) (v0(i), i=1,ih)
c  if v1(i) is 0, then later it will taken to be  the null
c  variance at the  last look. 
c       read(5,*)
c       read(5,*)
c       read(5,*)
c       read(5,*) (v1(i), i=1,ih)
	do 3, i=1,ih
	  v0(i)=0.
	  v1(i)=0.
3	continue
       read(5,*)
       read(5,*)
       read(5,*) n
       if(n .le. 0) then
	  n=1000
	  write(6,*) ' The number of simulation is set to 1000.'
       endif
c Input the type-I error:'
       read(5,*)
       read(5,*)
       read(5,*) pn
       if( pn .le. .0000000001) then
	  pn=.05
          write(6,*) ' The type-I error is set to 0.05.'
       endif
c Input 1 or 2 sided test, type 1 or 2:'
       read(5,*)
       read(5,*)
       read(5,*) nside
       if( nside .eq. 0) then
	  nside=2
	  write(6,*) ' 2 sided test is spcified.'
       endif
       e(0)=0.
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*) (e(i),i=1,k)
       if(jj .eq. 1) then
	   if( k .eq. 1 ) then
	      e(1)=pn
           else
           do 551 j=1,k
	      e(j)=e(j)+e(j-1)
551        continue
             if( abs( e(k)-pn ) .gt. .000001) then
       write(6,*) " The sum of the e's is not equal to the",
     & " type-I error."
	     endif
           endif
       endif
c  jj = 3
c  Input the b in the repeated significance test, if 0, it will be set to 3    
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*) binput
         if( binput .le. 1.) then
	     b=3
	     write(6,*) 'Your input of b is less than 1, b is re-set to 3.'
	   else
	     b=binput
         endif
c  jj = 4
c Input the boundary values: if 0 it will be set to 1.96.
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*) (bb(i),i=1,k)
       do 768 i=1,k
768           if( bb(i) .lt. 0.00001 ) bb(i)=1.96
c   jj = 2
c   input the choice of use function here
c   Input the choice of the use function: type 1, 2, 3 or 4:
c  write(6,*) '     1-------OBrien and Fleming type.' 
c  write(6,*) '     2-------Pocock type.' 
c  write(6,*) '     3-------Linear use function.' 
c  write(6,*) '     4-------Haybittle-Pocock use function.' 
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*) i_usefunction
       if( i_usefunction .le. 0.01 ) then
	   i_usefunction=1
           write(6,*) ' The OBrien and Fleming use function is chosen.'
       endif
c  Input the number of intervals to specify the baseline survival function: 
c  Up to nl are allowed:
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*) l
       if( l .eq. 0 ) then
	   l=1
	   write(6,*) ' The number of intervals is set to 1.'
       endif
c Input the right cutpoint(s) of ',i2, ' interval(s):
       t(0)=0.
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*) (t(i),i=1,l)
       if (t(l) .le. 0.0000000001) then
	  do 626 i=1,l
626       t(i)=10.*float(i)/float(l)
          write(6,627) 
       endif
       surt(0)=1.
c   The survivial function at t(0) has been set to 1.
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*) (surt(i), i=1,l)
       if( surt(l) .le. 0.00001 ) then
	  surt(l)=0.00001
	  write(6, 901) t(l)
901     format(' Since your last survival value is <.00001.'/' The ',
     & 'survival value at ', f9.3,' is set to .00001. Hope this is OK.')
	 endif
c   Input accrual period and accrual rate the product of these must be less than m
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*) accut, accrat
       if( accut .eq. 0 ) accut=10
       if(accrat .eq. 0 ) accrat=10
c Input the number of intervals (>=1) to specify the alternative distribution 
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*) ihaz
	if (ihaz .le. 0) then
	   ihaz=1
	   write(6,*) ' The number of intervals is set to 1.'
	endif
c Input the right cutpoint(s) of the interval(s) for specifying distributions:
       read(5,*)
       read(5,*)
       read(5,*) (rt(i),i=1,ihaz)
c  if the last endpoint < .0000000001, then the endpoints were be set to
c  equal space apart between 0 and 10.
        if ( rt(ihaz) .le. 0.0000000001 ) then
	  do 926 i=1,ihaz
926       rt(i)=10.*float(i)/float(ihaz)
          write(6,627) 
627     format(' The right cutpoints are equally spaced. The',
     &  ' largest cutpoint is set to 10.')
	endif
c  Want hazard ratio for specifying alternative distribution or not
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*) haznot
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*)
      if( haznot .gt. 1.5 ) then
c Input the ',i2,' non-negative hazard ratio(s):')
         read(5,*) (rhaz(i), i=1,ihaz)
c  if anything < .0000000001 for hazard ratio was inputed, it will be set to 1
	 do 927 i=1,ihaz
         if ( rhaz(i) .le. 0.0000000001) then
           rhaz(i)=1.0
	   write(6,1927)  rt(i)
1927       format('The hazard ratio on [0,',f5.2,'] is set to 1.')
	 endif
927      continue
      else
          surt2(0) = 1.
          l2=ihaz
          do 929 i=1, l2
929       t2(i)=rt(i)
	  read(5,*) (surt2(i), i=1,l2)
	if( surt2(l2) .le. 0.00001 ) then
	  surt2(l2)=0.00001
	  write(6, 901) t2(l2)
	endif
       endif
c  Input the censoring distribution for group 1
c Input the number of intervals to specify the censoring survival function 
c for the control group, up to nl=20 are allowed
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*) lc1
       if( lc1 .le. 0) then
	  lc1=1
          write(6,*) 'The number of intervals is set to 1.'
       endif
       tc1(0)=0.0
       read(5,*)
       read(5,*)
       read(5,*) (tc1(i),i=1,lc1)
	if( tc1(lc1) .le. 0.00000001 ) then
	   do 826 i=1,lc1
826        tc1(i)=20.*float(i)/float(lc1)
           write(6,92) 
92        format(' The right cutpoints are equally spaced. The',
     &    ' largest cutpoint is set to 20.')
	endif
       sc1(0)=1.
c  Input the survival function.  The value should be between 0 and 1
       read(5,*)
       read(5,*)
       read(5,*) (sc1(i), i=1,lc1)
       if( sc1(lc1) .le. 0.00001 ) then
	  sc1(lc1)=0.00001
	  write(6, 901) tc1(lc1)
       endif
c  Input the censoring distribution for group 2
c Input the number of intervals to specify the censoring survival for the 
c treatment group, up to ',i2,' are allowed:
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*) lc2
       if( lc2 .le. 0) then
	 lc2=1
	 write(6,*) ' The number of intervals is set to 1.'
       endif
       tc2(0)=0.0
       read(5,*)
       read(5,*)
       read(5,*) (tc2(i),i=1,lc2)
       if( tc2(lc2) .le. 0.00000001 ) then
         do 426 i=1,lc2
426        tc2(i)=20.*float(i)/float(lc2)
         write(6,92) 
       endif
       sc2(0)=1.
       read(5,*)
       read(5,*)
       read(5,*) (sc2(i), i=1,lc2)
       if( sc2(lc2) .le. 0.00001 ) then
	  sc2(lc2)=0.00001
	  write(6, 901) tc2(lc2)
	endif

c  Input the drop-in (out) rate
c Input the number of intervals to specify the drop-in rate(s): Up to nl are allowed
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*)
       read(5,*) ld1
	if( ld1 .le. 0 ) then
	   ld1=1
	   write(6,*) ' The number of intervals is set to 1.'
        endif
c Input the time point(s) for specifying drop-in rate(s):             
       read(5,*)
       read(5,*) (td1(i), i=1,ld1)
c Input the drop-in rate(s):
       read(5,*)
       read(5,*) (rdin(i), i=1,ld1)
c Input the number of intervals to specify the drop-out rate(s): Up to nl are allowed:
       read(5,*)
       read(5,*)
       read(5,*) 
       read(5,*) ld2
	if( ld2 .le. 0 ) then
	   ld2=1
	   write(6,*) ' The number of intervals is set to 1.'
        endif
c Input the time point(s) for specifying drop-out rate(s):             
       read(5,*) 
       read(5,*) (td2(i), i=1,ld2)
c Input the ',i2,' drop-out rate(s):')
       read(5,*) 
       read(5,*) (rdout(i), i=1,ld2)
c	 read(5,*)
c      read(5,*)
c      read(5,*)
c If method 4 is specified, the following lines read the input info. on
c futility stopping
c       if( jj .eq. 4) then
       futility=1.
c          read(5,*) futility
c  futility = 2 specifies stopping for futility
c          if (futility .gt. 1.5) then
c             if( nside .eq. 2) then
c               nside = 1
c               write(6,*) ' The test is reset to one-sided.'
c             endif
c             read(5,*)
c             read(5,*) (cc(i),i=1,k)
c          endif
c       endif
c-------------
c  Start printout
c
       write(4,785) iseed
785    format('The seed = ', I10) 
       tmp=ran3(iseed)
       write(4,786) k
786    format('Number of interim analyses = ', I2)
       write(4,980) (tlook(i),i=1,k)
980     format('Times of the interim analysis = ', 12(x,f5.2))
       write(4,789) (rho(i,1),i=1,ih)
789    format('The rho(s) of the statistic(s)  = ',x, 5(f8.2,x) )
       write(4,788) (rho(i,2),i=1,ih)
788    format('The tau(s) of the statistic(s) = ',x, 5(f8.2,x) )
       write(4,790) accut, accrat
790    format('Accrual period and accrual rate = ', x, 2(f7.2,x) )
	write(4,997) pn, nside
997     format('Type-I error is ',f5.3,' and ',
     &  i1,' sided test has been specified.')
	write(4,*) '(x,S(x)) of the baseline survival function are:'
	write(4,998) (t(i),surt(i),i=0,l)
998     format(x,(x,'(',f8.3,x,f7.3,')'))
       if( haznot .gt. 1.5 ) then
        write(4,*) 'The alternative survival distribution is specified'
        write(4,*) 'by the hazard ratio to the baseline survival:'
         write(4,996) t(0), rt(1), rhaz(1) 
        if( ihaz .gt. 1 ) then
	 write(4,996) ( rt(i-1), rt(i), rhaz(i), i=2,ihaz)
996      format(x,'hazard ratio in (',f7.2,x,f7.2,') is ', f8.4)
        endif
       else
	write(4,*) '(x,S(x)) of the alternative survival function are:'
	write(4,998) (t2(i),surt2(i),i=0,l2)
       endif
	write(4,5500)
5500   format(' Censoring survival function for the control group',
     & ' are:')
	write(4,998) (tc1(i),sc1(i),i=0,lc1)
	write(4,5501) 
5501   format(' Censoring survival function for the treatment group',
     & ' are:')
	write(4,998) (tc2(i),sc2(i),i=0,lc2)
	write(4,*) 'The drop-in rate(s) are:'
	write(4,998) (td1(i), rdin(i),i=1,ld1)
	write(4,*) 'The drop-out rate(s) are:'
	write(4,998) (td2(i), rdout(i),i=1,ld2)
       if( jj .eq. 2) then
	 if( i_usefunction .eq. 1) then
       write(4,*) 'The OBrien & Fleming use function is specified.'
	 else
	    if( i_usefunction .eq. 2) then
	      write(4,*) 'The Pocock use function is specifieisd.'
	    else
	      if( i_usefunction .eq. 3 ) then
	    write(4,*) 'The linear use function is specified.'
	      else
            write(4,*) 'The Siegmund (Haybittle-type) use function is'
	    write(4,*) 'specified. The b used is:', b
	      endif
	    endif
	 endif
       endif
       do 888 i=1,k
	 va(i)=1.
	 bbv(i)=1.
888    continue
c
c   Calculating variances under the null distribution 
c
       mm=int(accut*accrat)
	 write(4,335) mm
335    format('--------------------------------------------------',
     &       /,'The total number of subjects simulated =', i5)
      write(4,999) n
999   format('The number of simulations = ',I8)
      do 66 i=1,ih
	write(4, 375)
375   format('-------------------------------------------------------')
	if( i .eq. 1) write(4,773) i
	if( i .eq. 2) write(4,573) i
	if( i .eq. 3) write(4,273) i
	if( i .gt. 3) write(4,373) i
773   format('Pre-simulation information for the',i2,'st statistic are:'
     &)
573   format('Pre-simulation information for the',i2,'nd statistic are:'
     &)
273   format('Pre-simulation information for the',i2,'rd statistic are:'
     &)
373   format('Pre-simulation information for the',i2,'th statistic are:'
     &)
      rho1=rho(i,1)
	rho2=rho(i,2)
	call vari(nk,accrat,accut,rho1,rho2,k,tlook,v)
        if( jj .eq. 3 ) then
	   e(k)=pn/3.
c  In here we use the theoretical variance to check if the b is too small
105        do 106 j=1,k
106        bb(j)=b
           nbtoosmall=0
	   call probbdy3(nk,k,v,e,nside,nbtoosmall,bb)
           if( nbtoosmall .eq. 1) then
	      write(6,2023)   
2023  format('b is too small.  I will increase b by 0.1')
	      b=b+.1
	      goto 105
	     else
	       write(4,*) 'The b used in the Haybittle boundary is:', b 
           endif
          write(4,2021) bb(k)
2021  format('The c based on the asymptotic null distribution:', F6.2)

        endif
	  write(4,374)  (v(j),j=1,k)
374	format('The variance under the null are:'/x, 15(f7.2,2x))

       if(v0(i) .le. 0.00000001) v0(i)=v(1)*.5
       if(v1(i) .le. 0.00000001) v1(i)=v(k)
c      if ( k .eq. 1 ) goto 66
c       write(4,774) v0(i), v1(i)
c774   format('The smallest variance needed and the maximum allowable',
c     & ' variance are:'
c     &/x, 2(f7.2,2x))
66     continue
	
      if ( jj .eq. 2 .and. i_usefunction .eq. 4 ) then
	      if ( b .lt. 0.001 ) then
		 b=3
		 write(6,*) ' b is set to 3.'
              endif
	      tmp=0.
	      do 33 i=1,ih
		tmpp=alog(v1(i)/v0(i))
		if( tmpp .gt. tmp ) then
	       	   tmp=tmpp
		endif
33            continue
34            if( (b-1/b)*phi(b)*tmp+4.0*phi(b)/b .gt. pn ) then
		 write(6,*) ' The present b is too small.'
                 write(6,*) ' I will increase b by 0.1.'
                 b = b+.1
		 goto 34
	      endif
      endif		      
      if( jj .eq. 4 ) then
c  In this case, different statistics will use the same boundary
	    do 104 i=1,ih
	    do 104 ik=1,k
104         be(i,ik)=bb(ik)
       endif
c
	 if( jj .eq. 3 ) then
           write(6,334) real(n*mm)*(alog(real(k))*.2 + 1.)/1600000.
	 else
	     if( jj .eq. 2) then
	          write(6,334) real(n*mm*k)/1600000.
	     else
	          if ( jj .eq. 1 ) then
	        write(6,334) real(n*mm*(k+1))/3200000.
	          else
	       write(6,334) real(n*mm)*(alog(real(k))*.2 + 1.)/3200000.
	          endif
	     endif
	 endif
       write(6,*) 'The output is stored in the file named below:'
       write(6,*) outfile
334   format('Start simulation, it will take approximately',f7.2,' min'
     & /'on a PC equiped with a P-II 300 processor.') 

	write(4, 375)

c
c----------------------------------------------------
c     Start simulation
c     initialize the icount
100     do 20 i=1,ih
	   sumtime(i)=0.
20      icount(i)=0
        nbtoosmall=0
c
c  start simulations here, ii is the number of simulation already done
c
	 do 2000 ii=1, n
c
c    generate data
c
	call gd(nl,m,mm,x,c,y,z,accut,accrat,iseed,ld1,
     &   ld2,td1,td2,rdin,rdout, haznot)
c
c   if the i-th stat has crossed boundary, then idone(i)=1
	  do 30 i=1, ih
	   idone(i)=0
30      continue
c
c   start calculations. ilook is the number of looks currently.

       do 1000 ilook=1,k 
c   test to see if all the statistics have crossed boundary
	  isum=0

	 do 40 i=1,ih
40       isum=isum+idone(i)
	  if( isum .eq. ih ) goto 2001
c
c     if isum not equal ih, some statistics havn't crossed boundary, continue...
c
c     make censoring, depending on the looking time. 
c
	tlk=tlook(ilook)
	 call handle(m,x,c,y,z,mm,xx,zz,ms,idelta,tlk)
c       xx----contains the current sample.
c       zz----the corresponding covariates.
c       delta----censoring indicator.
c       mm-----is the old sample saize
c       ms is the new sample size 
c
	 if ( ms .eq. 0 ) then
	     write(4,*) 'no life enters before', tlook(ilook)
	     write(4,*) 'no. of simulation=', ii
	     goto 1000
	  endif
	 call sort(m,ms,xx,zz,idelta)
c
c    MCALSTAT calculates the ih standardized statistics
	CALL  MCALSTAT(m,ms,xx,zz,idelta,ih,rho,stat,var,md)
	if( jj .eq. 4 ) goto 47
	do 45 i=1,ih
	   ve(i,ilook)=var(i)
45      continue
c
c    getbound will determine the current boundary based on the methods,
c    usefunction, spending error, previous boundary, v0, v1 and variances.
c    The boundary points will be 5 if var 
c    is smaller than v0 or is larger than v1
c
       call getbound(nk, nih, ilook, jj, i_usefunction, e, be, v0, v1, 
     &  ve, ih, pn, k, b, nside, nbtoosmall)
c
47     do 50 i=1,ih
       if( idone(i) .eq. 0 ) then
	 if( ilook .lt. k ) then
	     if (stat(i) .ge. be(i,ilook) ) then
		   icount(i)=icount(i)+1
		   sumtime(i)=sumtime(i)+tlook(ilook)
		   idone(i)=1
	     else
		   if( nside .eq. 2 .and. stat(i) .le. -be(i,ilook)) then
		      icount(i)=icount(i)+1
		      sumtime(i)=sumtime(i)+tlook(ilook)
		      idone(i)=1
             else
		      if( futility .gt. 1.5 .and. stat(i) .le. cc(ilook)) 
     &then
		          sumtime(i)=sumtime(i)+tlook(ilook)
		          idone(i)=1
                endif
		   endif
	     endif
	 else
	     if (stat(i) .ge. be(i,k) ) then
		  icount(i)=icount(i)+1
		  idone(i)=1
	     else
		   if( nside .eq. 2 .and. stat(i) .le. -be(i,k)) then
		          icount(i)=icount(i)+1
		          idone(i)=1
              else
		      if( futility .gt. 1.5 .and. stat(i) .le. cc(k)) 
     &              idone(i)=1
		      endif 
	        endif
	        sumtime(i)=sumtime(i)+tlook(k)
	     endif
       endif      
50     continue
c
1000     continue
2001     continue
2000     continue
c 
c
c prints out results from simulations
c

c-------------
	  write(4,1011)
	  write(4,1008) jj
1008       format('The results of simulation for method:', i2)
        if( futility .gt. 1.5 )
     &  write(4,*) 'Stop for futility is specified.' 
	write(4,1011)
	rn=real(n)
	do 300 i=1,ih
	   pp = real(icount(i))/rn
	  if( i .eq. 1) write(4,1009) 
	  if( i .eq. 2) write(4,2009) 
	  if( i .eq. 3) write(4,3009) 
	  if( i .gt. 3) write(4,4009) i 
	  write(4,1010) rho(i,1), rho(i,2), pp, sqrt(pp*(1.-pp)/rn)
	  if ( jj .eq. 4 ) then
c     write(4,1012) (be(i,j),j=1,k)
c          else
	     write(4,2012) (bb(j),j=1,k)
	     if (futility .gt. 1.5) write(4,2013) (cc(j),j=1,k)
          endif
	  if( jj .le. 1) then
	  write(4,661) k, (e(j)-e(j-1),j=1,k)
661   format('The type-I error is splitted as e(1),..., e(',i2,').',
     &  ' They are:',/x, 12(f6.4,x) ) 
	  endif
	  write(4,1013) sumtime(i)/rn
	  write(4,1011)
1009    format('---------Result of the 1st statistic:-----------')
2009    format('---------Result of the 2nd statistic:-----------')
3009    format('---------Result of the 3rd statistic:-----------')
4009    format('---------Result of the ',i2,'th statistic:-----------')
1011    format(80X)
1010    format('rho =',F4.1,3x,'tau =',F4.1,3x,' power =',F6.4,3x,
     &  's.d. = ',F7.5)
1012    format('The last boundary in the test is:',/,12(f6.2,x))
2012    format('The boundary values in the test are:',/,12(f6.2,x))
2013    format('The futility boundary values in the test are:',
     &   /,12(f6.2,x))
1013    format('The average stopping time = ',F6.3)
300     continue
	  if( jj .eq. 3 .and. binput .lt. b ) then
	     write(4,4022) binput, b
4022	format(/,'********************** Warning ***********************',
     &    /,' The b in the input file is', F6.2,', which is too small.',
     &     /,'The b used in the simulation is', F6.2, '. Please keep ',
     &     /,'this in mind in your considerations.',
     &     /, '******************************************************')
	  endif
      end



c--------------
       subroutine getbound(nk, nih, ilook, jj, i_use, e, be, v0, v1,
     &  ve, ih, pn, k, b, nside, nbtoosmall)
c--------------
       implicit none
       integer nk, nih, ilook, jj, i_use, nbtoosmall, nside
       integer ifault, j, i, ih, k
       real v0(nih), v1(nih), e(0:nk), be(nih,nk), ve(nih,nk)
       real bb(12), b, va(12), phi, alnorm, gauinv
       real zp, zpn, v(12), pn
       logical upper


       do 1000 i=1, ih

       if( ve(i,ilook) .lt. v0(i) ) then
	  if( ilook .eq. k .and. k .gt. 1) then
	    write(6,*) 'The', i, '   th v0 is too large.'
	    bb(k)=gauinv(1.-pn/real(nside),ifault)
	    if( ifault .eq. 0 ) then
	       bb(k)=5.0
	    endif
	  else
	    bb(ilook)=5.0
	  endif
       else
	  v(1)=ve(i,1)
	  if( ilook .gt. 1) then
	    do 20 j=1, ilook-1
	      bb(j)=be(i,j)
	      if( ve(i,j+1) .ge. v(j)+.0001 ) then
		v(j+1)=ve(i,j+1)
	      else
		v(j+1)=v(j)+.0001
	      endif
20          continue
	  endif

	  do 622 j=1,ilook
	      if( v(j) .le. v1(i) ) then
		va(j)=v(j)/v1(i)
	      else
		va(j)=1.
	      endif
622       continue
	  
	 if( jj .eq. 2) then
	   if( va(ilook) .lt. .99 .and. ilook .lt. k ) then
	   if( i_use .eq. 1 ) then
c  write(4,*) 'The OBrien & Fleming type use function is specified.'
	      upper=.true. 
	      zpn=gauinv(1.-pn*.5,ifault)
	      zp=zpn/sqrt(va(ilook))
	      e(ilook)=2.*alnorm(zp,upper)
	   else
	      if ( i_use .eq. 2) then
c The Pocock type use function is specifieisd.
		e(ilook)=pn*alog(1.+1.71828*va(ilook))
	      else
		 if( i_use .eq. 3) then
c  The linear use function is specified.
		     e(ilook)=pn*va(ilook)
		 else
c  The Haybittle use function is specified
		   e(ilook)=(b-1/b)*phi(b)*alog(v(ilook)/v0(i))
     &                         +4.0*phi(b)/b
		 endif
	      endif
	   endif
	   call probbdy3(nk,ilook,v,e,nside,nbtoosmall,bb)
	   else
	     if( va(ilook) .ge. .99 .and. ilook .lt. k) then
		   bb(ilook)=5.0
		   if( ilook .gt. 1 ) then
		      ve(i,ilook)=(ve(i,ilook)+ve(i,ilook-1))*.5 
                   else
		      ve(i,ilook)=ve(i,ilook-1)*.5
                   endif
	     else
c   In here, ilook=k 
		   e(k)=pn
		   call probbdy3(nk,k,v,e,nside,nbtoosmall,bb)
	     endif
	   endif
	 else
	      if( jj .eq. 3) then
		if( ilook .eq. k) then
		  e(ilook)=pn
		  call probbdy3(nk,ilook,v,e,nside,nbtoosmall,bb)
		else
		  bb(ilook)=b
		endif
	      else
		if( jj .eq. 1 ) then
		  call probbdy3(nk,ilook,v,e,nside,nbtoosmall,bb)
		endif
	      endif
	 endif
       endif
       if( nbtoosmall .eq. 1) then
	   be(i,ilook)=5. 
       else
	   be(i,ilook)=bb(ilook)
       endif
1000   continue 
       end


c-------------
c      this prog. generates staggered entry data.
c      no censoring.  Uniform entry.
c      baseline sample is specified by function surv.
c-------------
	subroutine  gd(nl,m,mm,x,c,y,z,accut,accrat,iseed,
     & ld1, ld2, td1, td2, rdin, rdout, haznot)
c-------------
       implicit none
       integer m, i, mm, iseed, nl, ld1,ld2, j
       real x(m), y(m), z(m), c(m), td1(nl), rdin(nl), rdout(nl)
       real accut, accrat, tmp, td2(nl), haznot
       real ran3, rinvsu, rinvsu2, hinvsu, rinvc1su, rinvc2su
       
c           mm is the sample size.
c         x is the survival time, y is the entry time and z is the covariates,
c         where first half of the  covariate is 0 and the
c         second half of the covariate is 1.  c is the censoring variable.
c
       do 100 i=1,mm/2
	   y(i)=ran3(iseed)*accut
	   tmp=ran3(iseed)
	   x(i)=rinvsu(tmp)
c The crossover are supposed to happen at td1(j), j=1,ld1
	   do 50 j=1,ld1
	      if( x(i) .gt. td1(j) ) then
	         if( ran3(iseed) .le. rdin(j) ) then
40                if( haznot .gt. 1.5 ) then
		           x(i)=hinvsu(tmp)
                  else
		           x(i)=rinvsu2(tmp)
                  endif
		        if ( x(i) .le. td1(j)) then
		           tmp=ran3(iseed)
		           goto 40
		        endif
		        goto 60
		      endif
	      endif
50       continue
60       c(i)=rinvc1su(ran3(iseed))
	   z(i)=0.
100    continue
       if ( haznot .gt. 1.5 ) then
       do 101 i=mm/2+1, mm
	   y(i)=ran3(iseed)*accut
	   tmp=ran3(iseed)
	   x(i)=hinvsu(tmp)
c The crossover are supposed to happen at td2(j), j=1,ld2
	   do 70 j=1,ld2
	     if( x(i) .gt. td2(j) ) then
	       if( ran3(iseed) .le. rdout(j) ) then
80                x(i)=rinvsu(tmp)
		  if ( x(i) .le. td2(j)) then
		     tmp=ran3(iseed)
		     goto 80
		  endif
		  goto 90
		endif
	      endif
70         continue
90         c(i)=rinvc2su(ran3(iseed))
	     z(i)=1.
101    continue
       else
       do 201 i=mm/2+1, mm
	   y(i)=ran3(iseed)*accut
	   tmp=ran3(iseed)
	   x(i)=rinvsu2(tmp)
c The crossover are supposed to happen at td2(j), j=1,ld2
	   do 270 j=1,ld2
	     if( x(i) .gt. td2(j) ) then
	       if( ran3(iseed) .le. rdout(j) ) then
280                x(i)=rinvsu(tmp)
       	     if ( x(i) .le. td2(j)) then
		        tmp=ran3(iseed)
		        goto 280
		     endif
		     goto 290
		   endif
	     endif
270      continue
290         c(i)=rinvc2su(ran3(iseed))
	   z(i)=1.
201    continue
       endif
       end

	real function surv(x)
c       This routine takes the value of a survival function
c       at 0<t(1)<...<t(l) and x and then returns the value
c       of Surv(x)
c 
	implicit none
	integer nl, i
	parameter(nl=20)
	real x, tp

	integer l, ihaz, lc1, lc2, l2
	real t(0:nl), surt(0:nl), rt(nl), rhaz(nl)
	real t2(0:nl), surt2(0:nl)
	real tc1(0:nl), tc2(0:nl), sc1(0:nl), sc2(0:nl)

        common l, t, surt, ihaz, rt, rhaz, l2, t2, surt2
	common lc1, lc2, tc1, tc2, sc1, sc2
     
	if( x .le. t(0) ) then
	   surv=1.
	else
	   i=1
10         if( x .le. t(i) ) then
              tp=(alog(surt(i-1))-alog(surt(i)))/(t(i)-t(i-1))
              surv=surt(i-1)*exp(-tp*(x-t(i-1)))
	   else
	      i=i+1
	      if(i .le. l) goto 10
		tp=(alog(surt(l-1))-alog(surt(l)))/(t(l)-t(l-1))
              surv=surt(l)*exp(-tp*(x-t(l)))
	   endif
	 endif
	 end


	real function c1surv(x)
c       This routine takes the value of group1 censoring survival function
c       at 0<tc1(1)<...<tc1(lc1) and x and then returns the value
c       of c1surv(x)
c 
	implicit none
	integer nl
	parameter(nl=20)
	real x, t(0:nl), surt(0:nl), rt(nl), rhaz(nl), tp
	real t2(0:nl), surt2(0:nl), haznot
	integer l, i, ihaz
       integer lc1, lc2, l2
       real tc1(0:nl), tc2(0:nl), sc1(0:nl), sc2(0:nl)


        common l, t, surt, ihaz, rt, rhaz, l2, t2, surt2
	common lc1, lc2, tc1, tc2, sc1, sc2
     
	if( x .le. tc1(0) ) then
	   c1surv=1.
	else
	   i=1
10         if( x .le. tc1(i) ) then

              tp=(alog(sc1(i-1))-alog(sc1(i)))/(tc1(i)-tc1(i-1)) 
              c1surv=sc1(i-1)*exp(-tp*(x-tc1(i-1)))
	   else
	      i=i+1
	      if(i .le. lc1) goto 10
              tp=(alog(sc1(lc1-1))-alog(sc1(lc1)))/(tc1(lc1)-tc1(lc1-1))
              c1surv=sc1(lc1)*exp(-tp*(x-tc1(lc1)))
	   endif
	 endif
	 end


	real function c2surv(x)
c       This routine takes the value of group2 censoring survival function
c       at 0<tc1(1)<...<tc1(lc1) and x and then returns the value
c       of c2surv(x)
c 
	implicit none
	integer nl
	parameter(nl=20)
	real x, t(0:nl), surt(0:nl), rt(nl), rhaz(nl), tp
	real t2(0:nl), surt2(0:nl), haznot
	integer l, i, ihaz, l2
       integer lc1, lc2
       real tc1(0:nl), tc2(0:nl), sc1(0:nl), sc2(0:nl)


	common l, t, surt, ihaz, rt, rhaz, l2, t2, surt2
	common lc1, lc2, tc1, tc2, sc1, sc2
     
	if( x .le. tc2(0) ) then
	   c2surv=1.
	else
	   i=1
10         if( x .le. tc2(i) ) then
              tp=(alog(sc2(i-1))-alog(sc2(i)))/(tc2(i)-tc2(i-1)) 
              c2surv=sc2(i-1)*exp(-tp*(x-tc2(i-1)))
	   else
	      i=i+1
	      if(i .le. lc2) goto 10

              tp=(alog(sc2(lc2-1))-alog(sc2(lc2)))/(tc2(lc2)-tc2(lc2-1))
              c2surv=sc2(lc2)*exp(-tp*(x-tc2(lc2)))
	   endif
	 endif
	 end
c
c
	real function psi(x, rho1, rho2)
c       This routine returns the weight function in the
c       Bete family of statistics
c 
	implicit none
	real x, t1, t2, tp, rho1, rho2, surv, c1surv, c2surv 

	t1=c1surv(x)
	t2=c2surv(x)
	tp=t1+t2
	if( tp .ge. .00001 ) then
	  psi=surv(x)**(2*rho1)*(1-surv(x))**(2*rho2)*t1*t2/tp
	else
	  psi=0.0
	endif
	end
c
c
	real function rinvsu(u)
c       This routine takes the value of a survival function
c       at 0<t(1)<...<t(l) and u and then returns the value
c       of inverse-surv(u), where 0 <= u <= 1.
c
	implicit none
	integer nl
	parameter(nl=20)
	real u, t(0:nl), surt(0:nl), rt(nl), rhaz(nl), dsurt
	real t2(0:nl), surt2(0:nl), haznot
	integer l, i, ihaz, l2
        integer lc1, lc2
        real tc1(0:nl), tc2(0:nl), sc1(0:nl), sc2(0:nl)
        real tp

        common l, t, surt, ihaz, rt, rhaz, l2, t2, surt2
        common lc1, lc2, tc1, tc2, sc1, sc2

        if( u .ge. 1. ) then
           rinvsu=t(0)
        else
           i=1
10         dsurt=surt(i-1)-surt(i)
           if( u .ge. surt(i) .and. dsurt .ge. .000001 ) then
             tp=(alog(surt(i-1))-alog(surt(i)))/(t(i)-t(i-1))
             rinvsu=t(i-1)-(alog(u/surt(i-1)))/tp
           else
              i=i+1
              if(i .le. l) goto 10
              tp=(alog(surt(l-1))-alog(surt(l)))/(t(l)-t(l-1))
             rinvsu=t(l)-(alog(u/surt(l)))/tp
           endif
         endif
         end
c
	real function rinvsu2(u)
c       This routine takes the value of a survival function
c       at 0<t(1)<...<t(l) and u and then returns the value
c       of inverse-surv(u), where 0 <= u <= 1.
c
	implicit none
	integer nl
	parameter(nl=20)
	real u, t(0:nl), surt(0:nl), rt(nl), rhaz(nl), dsurt
	real t2(0:nl), surt2(0:nl), haznot
	integer l, i, ihaz, l2
       integer lc1, lc2
       real tc1(0:nl), tc2(0:nl), sc1(0:nl), sc2(0:nl)
       real tp
	
	common l, t, surt, ihaz, rt, rhaz, l2, t2, surt2
       common lc1, lc2, tc1, tc2, sc1, sc2

	if( u .ge. 1. ) then
	   rinvsu2=t2(0)
	else
	   i=1
10         dsurt=surt2(i-1)-surt2(i)
	   if( u .ge. surt2(i) .and. dsurt .ge. .000001 ) then

             tp=(alog(surt2(i-1))-alog(surt2(i)))/(t2(i)-t2(i-1))
             rinvsu2=t2(i-1)-(alog(u/surt2(i-1)))/tp
	   else
	      i=i+1
	      if(i .le. l2) goto 10
              tp=(alog(surt2(l2-1))-alog(surt2(l2)))/(t2(l2)-t2(l2-1))
              rinvsu2=t2(l2)-(alog(u/surt2(l2)))/tp
	   endif
	 endif
	 end

c
	real function rinvc1su(u)
c       This routine takes the value of a survival function
c       at 0<tc1(1)<...<tc1(lc1) and u and then returns the value
c       of inverse-c1surv(u), where 0 <= u <= 1.
c
	implicit none
	integer nl
	parameter(nl=20)
	real u, t(0:nl), surt(0:nl), rt(nl), rhaz(nl), dsurt
	real t2(0:nl), surt2(0:nl), haznot
	integer l, i, ihaz, l2
       integer lc1, lc2
       real tc1(0:nl), tc2(0:nl), sc1(0:nl), sc2(0:nl)
       real tp
	
	common l, t, surt, ihaz, rt, rhaz, l2, t2, surt2
      common lc1, lc2, tc1, tc2, sc1, sc2

	if( u .ge. 1. ) then
	   rinvc1su=tc1(0)
	else
	   i=1
10       dsurt=sc1(i-1)-sc1(i)
	   if( u .ge. sc1(i) .and. dsurt .ge. .000001 ) then
             tp=(alog(sc1(i-1))-alog(sc1(i)))/(tc1(i)-tc1(i-1))
             rinvc1su=tc1(i-1)-(alog(u/sc1(i-1)))/tp
	   else
	      i=i+1
	      if(i .le. lc1) goto 10
            tp=(alog(sc1(lc1-1))-alog(sc1(lc1)))/(tc1(lc1)-tc1(lc1-1))
            rinvc1su=tc1(lc1)-(alog(u/sc1(lc1)))/tp
	   endif
	endif
	 end
c
	real function rinvc2su(u)
c       This routine takes the value of a survival function
c       at 0<tc2(1)<...<tc2(lc2) and u and then returns the value
c       of inverse-c1surv(u), where 0 <= u <= 1.
c
	implicit none
	integer nl
	parameter(nl=20)
	real u, t(0:nl), surt(0:nl), rt(nl), rhaz(nl), dsurt
	real t2(0:nl), surt2(0:nl), haznot
	integer l, i, ihaz, l2
       integer lc1, lc2
       real tc1(0:nl), tc2(0:nl), sc1(0:nl), sc2(0:nl)
       real tp
	
	common l, t, surt, ihaz, rt, rhaz, l2, t2, surt2
       common lc1, lc2, tc1, tc2, sc1, sc2

	if( u .ge. 1. ) then
	   rinvc2su=tc2(0)
	else
	   i=1
10       dsurt=sc2(i-1)-sc2(i)
	   if( u .ge. sc2(i) .and. dsurt .ge. .000001 ) then

             tp=(alog(sc2(i-1))-alog(sc2(i)))/(tc2(i)-tc2(i-1))
             rinvc2su=tc2(i-1)-(alog(u/sc2(i-1)))/tp
	   else
	      i=i+1
	      if(i .le. lc2) goto 10
            tp=(alog(sc2(lc2-1))-alog(sc2(lc2)))/(tc2(lc2)-tc2(lc2-1))
            rinvc2su=tc2(lc2)-(alog(u/sc2(lc2)))/tp
	   endif
	 endif
	 end
c
c
	real function hinvsu(u)
c       This routine takes the discrete baseline survival function
c       and the hazard ratio to get the inverse of the alternative
c       survival function
c
c       u, l, t surt are as in rinvsur
c       rt = the end points of intervals of constant hazard ratio
c       rhaz = the hazard ratios, rhaz(m+1)=0, rhaz(i)>0 if i>0.
c
	implicit none
	integer nl
	parameter(nl=20)
	real u, t(0:nl), surt(0:nl), rt(nl), rhaz(nl)
	real t2(0:nl), surt2(0:nl), haznot
	real rmlp, tmp, tmpp, rinvsu, surv
	integer l, i, ihaz, l2
       integer lc1, lc2
       real tc1(0:nl), tc2(0:nl), sc1(0:nl), sc2(0:nl)

	common l, t, surt, ihaz, rt, rhaz, l2, t2, surt2
       common lc1, lc2, tc1, tc2, sc1, sc2

	if( u .ge. 1. ) then
	   hinvsu=t(0)
	else
	   i=1
	   rmlp=1.
10         tmp=(u/rmlp)**(1./rhaz(i))
	   tmpp=surv(rt(i))
	   if( tmp .ge. tmpp ) then
	     hinvsu=rinvsu(tmp) 
	   else
	     i=i+1
	     if( tmpp .lt. .000001) write(4,*) 'tmpp too small'
	     if( rmlp .lt. .000001) write(4,*) 'rmlp too small'
	     if( i .le. ihaz ) then
	       rmlp=rmlp*tmpp**(rhaz(i-1)-rhaz(i))
	       goto 10 
	     else
	       rmlp=rmlp*tmpp**(rhaz(ihaz)-1.)
	       hinvsu=rinvsu(u/rmlp)
	     endif
	   endif
	endif
	end

	subroutine vari(nk,accrat,accut,rho1,rho2,k,tlook,var)
c
c       This routine takes the length of accrual time --accut--
c       and the rho in self's (G -rho) statistics (assumes that the accrual
c       is uniform) and calculates the variance at k different
c       looking time. At most nk looks are permitted.
c
	implicit none
	integer nk, nl
	parameter (nl=20)
	real accrat, accut, rho1, rho2, tlook(nk), var(nk)
	real sd, vinterval
	integer k, i, j 

	integer l, ihaz, lc1, lc2, l2
	real t(0:nl), surt(0:nl), rt(nl), rhaz(nl)
	real t2(0:nl), surt2(0:nl), haznot
	real tc1(0:nl), tc2(0:nl), sc1(0:nl), sc2(0:nl)

	common l, t, surt, ihaz, rt, rhaz, l2, t2, surt2
	common lc1, lc2, tc1, tc2, sc1, sc2

	do 100 i=1, k
	var(i)=0.0
	j=1
10      sd=(surt(j-1)-surt(j))/(t(j)-t(j-1))
	if( t(j) .ge. tlook(i) ) then
	    var(i)=var(i) + vinterval( t(j-1),tlook(i), tlook(i),
     &                accut, rho1, rho2)*sd
	else
	    var(i)=var(i) + vinterval( t(j-1), t(j), tlook(i),
     &                accut, rho1, rho2)*sd
	    if( j .lt. l ) then
c             this is to prevent to increase j but surt(j)=0
	      j=j+1
	      goto 10
	    endif
	endif
90      continue
	var(i)=var(i)/2.*accrat
100     continue
	end


	real function vinterval(v, w, tl, accut, rho1, rho2)
c  This calculates the variance contributed from the interval (v,w)
	implicit none
	real v, w, tl, rho1, rho2, tmp, psi, hh, tmpp, accut
	real tlmaccut

	tlmaccut = tl - accut
	vinterval= 0.0
	hh =.01
	tmpp= v + hh
	tmp = v + hh*.5
10      if( tmpp .le. w ) then  
	  if( tmp .le. tlmaccut) then
	     vinterval=vinterval+psi(tmp, rho1, rho2)*accut*hh
	  else
	     vinterval=vinterval+psi(tmp, rho1, rho2)*(tl-tmp)*hh
	  endif
	  tmpp=tmpp + hh
	  tmp=tmp + hh
	  goto 10
	else
	  tmpp=tmpp - hh
	  tmp=(tmpp + w)*.5
	  if( tmp .le. tlmaccut) then
	     vinterval=vinterval+psi(tmp, rho1, rho2)*accut*hh
	  else
	    vinterval=vinterval+psi(tmp, rho1, rho2)*(tl-tmp)*(w-tmpp)
	  endif
	endif
	end

c
c
c
	subroutine probbdy3(nk,kk,v,e,nside,nbtoosmall,bb)
c  prob computes the boundary crossing probability corresponding to a 
c  square root boundary b and returns the last value c (Haybittle,1977)
c  such that the total
c  probability equals to e(kk). The variances satisfy 
c        0<v(1)<v(2)<...<v(k)
c
c       bb  = the value of previous bdy
c       pn  = overall significance lever
c       v   = variances
c       k   = number of data review points (2 <= k <= 12)
c       bb(k)   = last bdy value, output
c       
      implicit none
      integer  nk, nside
      dimension    fn(0:1,12), ybn(0:1), var(0:1)
      dimension    v(nk),bb(nk),bbv(12),va(12)
      dimension    w(12),abss(12),e(0:nk)
      real    del, h, sump, bb, ybn, fn, v, e
      real    phi, abss,w,var, bbv, va, tmpp
      real    dsq, den, den1, a1, a3, p1, p2, varsqt
      real    x, c, h2, gauinv
      integer nbtoosmall, kk, i, ifault
      integer i1, i0, j, l
c
      data abss/.064056892862605626085, .191118867473616309159, 
     &          .315042679696163374387, .433793507626045138487, 
     &          .545421471388839535658, .648093651936975569252,
     &          .740124191578554364244, .820001985973902921954, 
     &          .886415527004401034213, .938274552002732758524, 
     &          .974728555971309498198, .995187219997021360180/
      data w/.127938195346752156974, .125837456346828296121, 
     &       .121670472927803391204, .115505668053725601353, 
     &       .107444270115965634783, .097618652104113888270,
     &       .086190161531953275917, .073346481411080305734, 
     &       .059298584915436780746, .044277438817419806169, 
     &       .028531388628933663181, .012341229799987199547/

c
c       fn(ik,j) = a matrix defining the density function at
c                   t(ik) at the j-th point in the interval
c                   between yan(ik) and ybn(ik).
c       ybn(ik)  = the upper limit of integration
c       h        = grid size used in the last round of numerical integration
c
        nbtoosmall=0
	if(kk .eq. 1) then
	   bb(1)=gauinv(1.-e(1)/real(nside), ifault)
	   if( ifault .eq. 1 ) then
	     bb(1)=5.0
	   endif
	   return
	endif

c  The boundary will be calculated as if 2 sided test has been specified
        if ( nside .eq. 1) then
	   tmpp=.5-e(kk)
	else
	   tmpp=(1.-e(kk))*.5
	endif

        if( v(kk) .gt. 0.0001 ) then
   	  do 60 i=1,kk
	     va(i)=v(i)/v(kk)
60        continue
        else
   	  do 61 i=1,kk
	     va(i)=real(i)/real(kk)
61        continue
	endif

c   limit on the number of grid points at each lever of density function
c   initialization
	i=2
	i1=1
	i0=0
	var(0)=va(1)
19      if( var(0) .lt. .00000001) then
c  In this case, the first increment is small, no integration is done
	  if( i .lt. kk ) then
	     var(0)=va(i)
	     i=i+1
	     goto 19
          else
c In this case, it is like you have only one look
	    bb(kk)=gauinv(1.-e(kk)/real(nside), ifault)
	    if( ifault .eq. 1 ) then
	      bb(kk)=5.0
	    endif
	    return
	  endif  
        endif
	varsqt=sqrt(var(0))
c   initialize the vector fn
	bbv(i-1)=bb(i-1)*varsqt
	do 10 j=1,12
	fn(0,j)=phi(abss(j)*bb(i-1))/varsqt
10      continue
	ybn(0)=bbv(i-1)
c
c      if k=2 then the next iteration is the last one
c
15     if (i .eq. kk) goto 100                  
c
c          fn(i1,j) is the prob. density  
c
20      var(i1)=va(i)
	bbv(i)=bb(i)*sqrt(va(i))
	if(var(i1)-var(i0) .lt. .00001 ) then
	  i=i+1
          goto 15
	endif  
c
c         prepare to integrate
	del=var(i1)-var(i0)
	dsq=sqrt(del)
c
       do 30 j=1,12
	fn(i1,j)=0.
	do 31 l=1,12
       fn(i1,j)=fn(i1,j)+fn(i0,l)*(phi((ybn(i0)*abss(l)-bbv(i)*abss(j))/
     & dsq) + phi((-ybn(i0)*abss(l)-bbv(i)*abss(j))/dsq))*w(l)
31      continue
	fn(i1,j)=fn(i1,j)*ybn(i0)/dsq
30     continue
       ybn(i1)=bbv(i)
c
c       prepare to do the next integration
	i=i+1
	i1=1-i1
	i0=1-i0
	goto 15
c
100     continue
c
c     begin to get the density at the last point.  
c
	h=.005
	h2=.0025
	sump=0.
	var(i1)=va(i)
c  var(i1) should be 1
	if(var(i1)-var(i0) .lt. 0.00001) then
c          write(4,*) ' the last increment is too small.'
           c=5
	endif
	del=var(i1)-var(i0)
	dsq=sqrt(del)
	x=0.
105     den=0.
	do 110 l=1,12
	den=den+fn(i0,l)*(phi((ybn(i0)*abss(l)-x)/dsq)
     &  + phi((-ybn(i0)*abss(l)-x)/dsq))*w(l)
110     continue
	den=den*ybn(i0)/dsq

	x=x+h
c
106     den1=0.
	do 111 l=1,12
	den1=den1+fn(i0,l)*(phi((ybn(i0)*abss(l)-x)/dsq)
     &  + phi((-ybn(i0)*abss(l)-x)/dsq))*w(l)
111     continue
	den1=den1*ybn(i0)/dsq

c
c         Pick up the next grid
	sump=sump + (den+den1)*h2
c
c         Test to see if the probability has been exceed
c
	if( sump .gt. tmpp ) goto 140
	x=x+h
	if( x .gt. 5. ) then
	  nbtoosmall=1
	  bb(kk)=5.
	  return
	endif
	den=den1
	goto 106
c
140     continue
c   Begin interpolation process for the last point
c   which produce the estimate for c
	p2=sump
	a1=x-h 
	p1=sump-(den+den1)*h2
c   a3 is the interpolated point
	a3=a1+h*(tmpp-p1)/(p2-p1)
c
c   do one more interpolation
	den1=0.
	do 112 l=1,12
	den1=den1+fn(i0,l)*(phi((ybn(i0)*abss(l)-a3)/dsq)
     &  + phi((-ybn(i0)*abss(l)-a3)/dsq))*w(l)
112     continue
	den1=den1*ybn(i0)/dsq

	a3=a1+2.*(tmpp-p1)/(den+den1)
c
c   get the output value
	bb(kk)=a3
       return
       end 


      real FUNCTION GAUINV(P, IFAULT)
C
C        ALGORITHM AS 70 APPL. STATST. (1974) VOL.23, NO.1
C
C        GAUINV FINDS PERCENTAGE POINTS  OF THE NORMAL DISTRIBUTION
C
      DATA ZERO, ONE, HALF, ALIMIT /0.0, 1.0, 0.5, 1.0E-20/
C
      DATA             P0,   P1,             P2,                P3
     *   / -.322232431088, -1.0, -.342242088547, -.204231210245E-1/
C
      DATA                P4,               Q0,            Q1
     *   / -.453642210148E-4, .993484626060E-1, .588581570495/
C
      DATA            Q2,            Q3,              Q4
     *   / .531103462366, .103537752850, .38560700634E-2 /
C
      IFAULT = 1
      GAUINV = ZERO
      PS = P
      IF (PS .GT. HALF) PS =ONE - PS
      IF (PS .LT. ALIMIT) RETURN
      IFAULT = 0
      IF (PS .EQ. HALF) RETURN
      YI = SQRT( ALOG( ONE / (PS*PS)))
      GAUINV = YI + ((((YI *P4 + P3)* YI +P2) * YI +P1) * YI +P0)
     *            / ((((YI *Q4 + Q3)* YI +Q2) * YI +Q1) * YI +Q0)
      IF (P .LT. HALF) GAUINV = -GAUINV
      RETURN
      END
c
c
       real FUNCTION ALNORM(X, UPPER)
C      1-p=phi(x)
C      ALGORITHM AS 66 APPL, STATST, (1979) VOL.22, NO.3
C
C      EVALUATESN THE TAIL AREA OF THE STANDARDISED NORMAL CURVE
C      FROM X TO INFINITY IF UPPER IF .TRUE. OR
C      FROM MINUS INFINITY TO X IF  UPPER IS .FALSE.
C
      REAL LTONE, UTZERO, ZERO, HALF, ONE, CON, Z, Y, X
      LOGICAL UPPER, UP
C
      DATA LTONE, UTZERO /7.0, 18.66/
      DATA ZERO, HALF, ONE, CON / 0.0, 0.5, 1.0,1.28/
      UP = UPPER
      Z = X
      IF (Z .GT. ZERO) GOTO 10
      UP = .NOT.UP
      Z = -Z
10    IF (Z.LE.LTONE.OR.UP.AND. Z.LE.UTZERO) GOTO 20
      ALNORM = ZERO
      GOTO 40
20    Y = HALF * Z * Z
      IF (Z.GT.CON) GOTO 30
C
      ALNORM = HALF - Z * (0.398942280444 - 0.399903438504 * Y /
     1  (Y + 5.75885480458 - 29.8213557808 /
     2  (Y + 2.62433121679 + 48.6959930692 /
     3  (Y + 5.92885724438))))
      GOTO 40
C
30    ALNORM = 0.398942280385 * EXP(-Y) /
     1  (Z - 3.8052E-8 + 1.00000615302 /
     2  (Z + 3.98064794E-4 + 1.98615381364 /
     3  (Z + 0.151679116635 + 5.29330324926 /
     4  (Z + 4.8385912808 - 15.1508972451 /
     5  (Z + 0.742380924027 + 30.789933034 / (Z + 3.99019417011))))))
C
40    IF (.NOT.UP) ALNORM = ONE - ALNORM
      END
c
c
      REAL FUNCTION PHI(X)

      REAL ARG, HALF, SQ2P, X, XLOW, ZERO
      REAL EXP

      DATA XLOW /-87.0/, SQ2P /0.91893853320467274/, ZERO /0.0/,
     *  HALF /0.5/
      PHI = ZERO
      ARG = -HALF * X * X - SQ2P 
      IF (ARG .GT. XLOW) PHI = EXP(ARG)
      RETURN
      END

c-------------
c-------------
	  subroutine handle(m,x,c,y,z,mm,xx,zz,ms,idelta,tk)
c-------------
c          this rountine put the minimal of x and t-y into xx
c          and if t-y < x then idelta=0
c          mm is the original sample size and ms is the new sample size.
c-------------
       implicit none
       integer m
       real x(m), y(m), z(m), xx(m), zz(m), c(m)
       real tk, tmp
       integer idelta(m), ms, mm, i

       ms=0
       do 10 i=1, mm
       tmp=min(tk-y(i), c(i))
       if ( tmp .gt. 0. ) then
	   ms=ms+1
	   zz(ms)=z(i)
	   if ( tmp .ge. x(i) ) then 
	      xx(ms)=x(i)
	      idelta(ms)=1
	   else
	      xx(ms)=tmp
	      idelta(ms)=0
	   endif
       endif
10     continue
       end

c-------------
c          this is the shell sort program
c          sort the dble array element xx(1) to xx(ms) and point1, point2 sorted
c          accordingly.
c-------------
       subroutine sort(m,ms,xx,point1,point2)
c-------------
       implicit none
       logical inordr
       integer m
       real xx(m), point1(m), t
       integer point2(m), i, ms, ndelta, k
       
       ndelta=ms
10      if(ndelta .gt. 1) then
	  ndelta=ndelta/2
20          inordr=.true.
	       do 30 i=1,ms-ndelta
		  if(xx(i).gt.xx(i+ndelta)) then
		     if (point2(i) .ne. point2(i+ndelta)) then
			k=point2(i)
			point2(i)=point2(i+ndelta)
			point2(i+ndelta)=k
		      endif
		     if (point1(i) .ne. point1(i+ndelta)) then
		       t=point1(i)
		       point1(i)=point1(i+ndelta)
		       point1(i+ndelta)=t
		      endif
		 t=xx(i)
		 xx(i)=xx(i+ndelta)
		 xx(i+ndelta)=t
		 inordr=.false.
		  endif
30       continue
	     if ( .not. inordr) goto 20
	  goto 10
       endif
       return
       end
c-------------
	 subroutine MCALSTAT(m,ms,xx,zz,idelta,ih,rho,stat,var,md)
c-------------
c      this subroutine calculates the G-rho statistics
c-------------
	implicit none
	integer nih, m
	parameter(nih=5)
	integer idelta(m), i, ms, ih, md, j
	real xx(m), zz(m), rho(nih,2), stat(nih), var(nih)
	real a, a1, h, term, vterm
  
c          ih is the number of statistics need to be cal. in here.

	  a=float(ms)
	  a1=0.
	  md=0
	  h=1.
	   do 14 i=1, ms
	   a1=a1 + zz(i)
14       continue
	 do 15 j=1,ih
	    var(j)=0.
15       stat(j)=0.
c          h denotes the value of Kaplan-Meier est.
	   do 10 i=1, ms
	   if ( idelta(i) .eq. 1) then
	       h=h*a/(a+1.)
	       md=md+1
	       term=zz(i)-a1/a
	       vterm=a1*(a-a1)/a/a
	      do 16 j=1,ih
		 stat(j)=stat(j)+(h**rho(j,1))*((1-h)**rho(j,2))*term
		 var(j)=var(j)+(h**(2.*rho(j,1)))
     &   *((1-h)**(2.*rho(j,2)))*vterm
16       continue
	   endif
	   a=a-1.
	   a1=a1-zz(i)
10       continue
c
c
c      md is the number of death
	if (md .eq. 0 ) goto 40
	do 20 j=1,ih
c  the negative sign here makes the statistic has positive drift when treatment
c  group has smaller hazard, or treatment is better than the control.   
	stat(j)=-stat(j)/sqrt(var(j))
20      continue
40      return
	end
c
c
c
c
c  ran3 for iris: use static vars.  Seed idum must be a integer.
c  Initial seed: put idum to some large negative odd number like -1999
c
      real function ran3(idum)
c         implicit real*4(m)
c         parameter (mbig=4000000.,mseed=1618033.,mz=0.,fac=2.5e-7)
      parameter (mbig=1000000000,mseed=161803398,mz=0,fac=1.e-9)
      static ma(55)
      static iff,inext,inextp
      data iff /0/
      if(idum.lt.0.or.iff.eq.0)then
	iff=1
	mj=mseed-iabs(idum)
	mj=mod(mj,mbig)
	ma(55)=mj
	mk=1
	do 11 i=1,54
	  ii=mod(21*i,55)
	  ma(ii)=mk
	  mk=mj-mk
	  if(mk.lt.mz)mk=mk+mbig
	  mj=ma(ii)
11      continue
	do 13 k=1,4
	  do 12 i=1,55
	    ma(i)=ma(i)-ma(1+mod(i+30,55))
	    if(ma(i).lt.mz)ma(i)=ma(i)+mbig
12        continue
13      continue
	inext=0
	inextp=31
	idum=1
      endif
      inext=inext+1
      if(inext.eq.56)inext=1
      inextp=inextp+1
      if(inextp.eq.56)inextp=1
      mj=ma(inext)-ma(inextp)
      if(mj.lt.mz)mj=mj+mbig
      ma(inext)=mj
      ran3=mj*fac
c
      return
      end
