!     Last change:  DOS   2 Aug 2000    9:36 pm
!                                       *** copyright 2000 ***
! *** filename spkbl5.f95               *** John F. Monahan **
!                                       **********************
      program pspkbl5
!  test of spkbl5 "spikeball" -- generates random points
!   on the surface of a p-dimensional sphere
!    -- fifth-order integration rule -- see Mysovskikh
!
!  test using integral of Langevin distribution
!   see Watson, Statistics on Spheres p.101 and compare with
!    Abramowitz & Stegun, Handbook of Math Fns, Tables 9.8 and 10.8
!
      integer, parameter :: n = 1000   ! number of replications
      integer, parameter :: nsc = 5    ! number of scale params
      real, dimension(6) :: mu,v
      real, dimension(nsc) :: kap,s,ss
      real, dimension(6,56) :: a      ! points are columns
      real t,w,w1,w2
      integer i,j,k,p,np,n1
!
      data kap/ .3, .5, .8, 1., 1.2 /  ! five scale params
!
  20  format(/' Integration of Langevin distributions with p=',i3,    &
     & ' and',i6,' replications'/9x,'kappa',9x,'mean',9x,'std err')
  21  format(i6,f8.2,2f15.8)
!                   interface block
      interface
      subroutine spkbl5(p,a,np,n1,w1,w2)
      integer, intent(in) :: p
      real, dimension(:,:), intent(out) :: a
      integer, intent(out) :: np,n1
      real, intent(out) :: w1,w2
      end subroutine spkbl5
      end interface
!     
      open( unit=6, file='spkbl5.out')
!
!                   initialize random number generator
      t = ran(5151917)
!                   vary dimension from 2 to 6
      do p = 2,6
      write(6,20) p,n
!                   parameter vector
!            mu is center of distribution (?), on surface
!              so || mu || = 1
      t = sqrt( real(p*(p+1)*(2*p+1)) / 6. )
      do k = 1,p
      mu(k) = real(k) / t
      end do  ! loop on k
!                   initialize sums
      s = 0.
      ss = 0.
!                   n replications (rotations of abscissas)
      do i = 1,n
!                   initialize for this rotation
      v = 0.
!
      call spkbl5(p,a,np,n1,w1,w2)
!                   generate this spikeball on surface of sphere
      do j = 1,np
!                   compute inner product of mu and vector
      t = 0.
      do k = 1,p
      t = t + mu(k)*a(k,j)
      end do  ! loop on k
!                   now compute integral for 5 values of kappa
      if( j .le. n1 ) then
        w = w1
                      else
        w = w2
                      end if ! ( j .le. n1 )
      do k = 1,nsc
      v(k) = v(k) + w*exp( - kap(k)*t )
      end do  ! loop on k
!                   end of rotations
      end do  ! loop on j
      do k = 1,nsc
      s(k) = s(k) + v(k)
      if( i .gt. 1 ) ss(k) = ss(k) +                                  &
     &        (real(i)*v(k)-s(k))*(real(i)*v(k)-s(k))/real(i*(i-1))
      end do  ! loop on k
!                   end of replication loop 
      end do  ! loop on i
!                   gets means, std errs and rescale as in tables
      do k = 1,nsc
      s(k) =  s(k) / real(n)
      ss(k) = ss(k) / real((n-1)*n)
!                   note that various functions are tabled for
!                    different values of p with different scalings
      t = 1.
      if( p .eq. 2 ) t = exp( - kap(k) )
      if( p .eq. 4 ) t = kap(k) * exp( - kap(k) ) / 2.
      if( p .eq. 5 ) t = kap(k) / 3.
      if( p .eq. 6 ) t = 1./8.
      s(k) = t * s(k) 
      ss(k) = t * sqrt( ss(k) )
      write(6,21) k,kap(k),s(k),ss(k)
      end do  ! loop on k
!                   end of dimension loop
      end do  ! loop on p
      stop
      end program pspkbl5
      subroutine spkbl5(p,a,np,n1,w1,w2)
!  computes random rotation of fifth degree quadrature points on surface
!  of sphere in p dimensions -- gives np points stored as columns in a
!
!  ** note different results for p = 2 or 3 because of duplications
!
!   p    integer    dimension of sphere
!   a    real       points on surface: point k is in a( ,k)
!   np   integer    number of points generated
!   n1   integer    number of points with first weighting w1
!   w1   real       weight of first n1 = 2*(p+1) points
!   w2   real       weight of last np-n1 = p*(p+1) points
!
!  for fifth degree quadrature, weights of points are as follows
!   first 2*p+2 points, weight is (7-p)*p/(2*(p+1)*(p+1)*(p+2))
!   last  p*(p+1) points, weight is 2*(p-1)*(p-1)/(p*(p+1)*(p+1)*(p+2))
!
! J F Monahan (June 1995) Dept of Statistics, N C State U, Raleigh 27695
!     Recoded February 2000 for Fortran 95
      implicit none
      integer, intent(in) :: p
      real, dimension(:,:), intent(out) :: a      ! dummy dimensions
      integer, intent(out) :: np,n1
      real, intent(out) :: w1,w2
      real, dimension(p) :: x
      integer i,ii,j,k,l,pp1
      real s,d,vi,vj
      real gnroul
!                   compute some constants
      pp1 = p + 1
!                   initialize to zero first
      a = 0.
!

!                   first put in the p+1 simplex vertices
      do i = 1,p
!                   i is indexing component
      vi = sqrt( real( (p+1)*(p-i+1) ) / real( p*(p-i+2) ) )
      vj = - sqrt( real(p+1) / real( (p-i+1)*p*(p-i+2) ) )
      do j = 1,pp1
!                   j is indexing vertices
      if( j .lt. i ) a(i,j) = 0.d0
      if( j .eq. i ) a(i,j) = vi
      if( j .gt. i ) a(i,j) = vj
!                   next p+1 points are negatives
      a(i,pp1+j) = - a(i,j)
      end do  ! loop on j
      end do  ! loop on i
!                   have simplex vertices, now midpoints and project
!                   set up for special case of p = 2
      if( p .eq. 2 ) then
      n1 = 6
      np = 6
      w1 = 1./6.
      w2 = w1
                     else     !   if( p .eq. 2 )
!                   continue with midpoints except for p = 2
      l = 2*(p+1)
!                   increment to reflect last points
      n1 = (p*(p+1))/2
      do i = 1,pp1
      do j = 1,i
!                   do all pairs (i,j) but not i=j
      if( i .ne. j ) then
!                   count points
      s = 0.d0
      l = l + 1
      do k = 1,p
      a(k,l) = a(k,i) + a(k,j)
      s = s + a(k,l)*a(k,l)
      end do  ! loop on k
!                   normalize - or - project to sphere
      s = sqrt(s)
      do k = 1,p
      a(k,l) = a(k,l) / s
      if( p .ne. 3 ) a(k,l+n1) = - a(k,l)         ! note anomaly with 3
      end do  ! loop on k
!                   done with this point
                     end if   ! ( i .ne. j )
      end do  ! loop on j
      end do  ! loop on i
!                   fix special case of p = 3
      if( p .eq. 3 ) then
      n1 = 8
      np = 14
      w1 = 3./40.
      w2 = 1./15.
                     else     ! if( p .eq. 3 )
!                   have created all np = (p+1)*(p+2) points
      n1 = 2*(p+1)
      np = (p+1)*(p+2)
      w1 = real( (7-p)*p ) / real( 2*(p+1)*(p+1)*(p+2) )
      w2 = real( 2*(p-1)*(p-1) ) / real( p*(p+1)*(p+1)*(p+2) )
                     end if   !   if( p .eq. 3 )
!                   have n1, np and weights, do the rotations
                     end if   !   if( p .eq. 2 )
!
!                   got points, start rotations
!
!                   multiply each vector v(k) = a( ,k), k = 1, ..., np
!                    by Householder tfms in the right order
!                    need H(1) * H(2) * ... * H(p-1) * v(k)
!                    so start with *** p-1 ***
      do ii = 2,p
      i = p - ii + 1
!                   generate random rotation by Householder transforms
!                   defined by iid normal elements
      s = 0.
      do j = i,p
      x(j) = gnroul(i+j)
      s = s + x(j)*x(j)
      end do  ! loop on j
      d = s
      s = sqrt(s)
      d = d - s*x(i)          ! use negative to insure R
      x(i) = x(i) - s         !  has positive diagonals
!                   x now holds vector u, d holds constant
!
!                   apply rotation to axial points
      do k = 1,np
!                   Householder requires one inner product
      s = 0.
      do j = i,p
      s = s + x(j)*a(j,k)
      end do  ! loop on j
      s = s / d
      do j = i,p
      a(j,k) = a(j,k) - s*x(j)
      end do  ! loop on j
      end do  ! loop on k
!                   done with rotation H(i)
      end do  ! loop on ii or i
      return
      end subroutine spkbl5
      REAL FUNCTION GNROUL(IX)
!  RATIO OF UNIFORMS ALGORITHM FOR GENERATING NORMAL(0,1) RV'S
!   USING LEVA'S QUADRATIC INNER AND OUTER BOUNDS
!
!  A J KINDERMAN AND J F MONAHAN (1977) "COMPUTER GENERATION OF RANDOM
!   VARIABLES USING THE RATIO OF UNIFORM DEVIATES," ACM TRANSACTIONS ON
!   MATHEMATICAL SOFTWARE, VOLUME 3, PP.257-260
!
!  J L LEVA (1992) "A FAST NORMAL RANDOM NUMBER GENERATOR," ACM
!   TRANSACTIONS ON MATHEMATICAL SOFTWARE, VOLUME 18, PP. 449-453.
!
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: IX         ! DUMMY ARGUMENT
      REAL, PARAMETER :: R = 1.7155277  ! FIRST CONSTANT R = SQRT(2/E)
!                   CONSTANTS S,T CENTER OF ELLIPSES
      REAL, PARAMETER :: S = .449871
      REAL, PARAMETER :: T = -.386595
!                   SHAPE PARAMETERS OF ELLIPSES
      REAL, PARAMETER :: A = .19600
      REAL, PARAMETER :: B = .25472
!                   RADII OF ELLIPSES
      REAL, PARAMETER :: R1 = .27597
      REAL, PARAMETER :: R2 = .27846
      REAL RAN
      REAL U,V,X,Y,Q
!                   START / RESTART HERE
      DO            ! NOTE UNRESTRICTED DO
!                   GENERATE (U,V) UNIFORMLY IN BOX
      U = RAN(1)
      V = R*(RAN(2) - 0.5)
      X = U - S
      Y = ABS(V) - T
!                   COMPUTE QUADRATIC PIECE
      Q = X*X + Y*(A*Y - B*X)
!                   COMPUTE DEVIATE BEFORE TESTS
      GNROUL = V / U
!                   INNER BOUND -- QUICK ACCEPT CHECK
      IF( Q .LE. R1 ) RETURN
!                   OUTER BOUND -- QUICK REJECT CHECK
      IF( A .GT. R2 ) CYCLE
!                   FINAL RATIO OF UNIFORMS CHECK
      IF( GNROUL*GNROUL .LE. -4.*LOG(U) ) RETURN
!                   REJECT -- START OVER
      END DO        ! END OF UNRESTRICTED DO
      END FUNCTION GNROUL
      REAL FUNCTION RAN(IDUM)
!  PORTABLE IMPLEMENTATION OF UNIFORM PSEUDORANDOM NUMBER GENERATOR
!    LEWIS, GOODMAN, & MILLER MULTIPLICATIVE GENERATOR
!     X(N+1) = MOD( 16807*X(N), 2**31-1 )
!
!  P. BRANTLEY, B.L. FOX, L. SCHRAGE (1983) A GUIDE TO SIMULATION
!   SPRINGER-VERLAG, NEW YORK.  PP. 200-202.
!  UPDATED VERSION OF
!  LINUS SCHRAGE (1979)'A MORE PORTABLE FORTRAN RANDOM NUMBER GENERATOR' 
!   ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE, VOLUME 5, PP. 132-138.
!
!  ARGUMENT 
!   IDUM    INTEGER  FIRST CALL SETS SEED, IGNORED IN SUBSEQUENT CALLS
!
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: IDUM
      REAL, PARAMETER :: RP = 4.656612875E-10     ! 1/P
      INTEGER, PARAMETER :: A = 16807             ! MULTIPLIER
      INTEGER, PARAMETER :: B = 127773            ! B = P / A
      INTEGER, PARAMETER :: C = 2836              ! C = P MOD A
      INTEGER, PARAMETER :: P = 2147483647        ! MODULUS 2**31 - 1
      INTEGER, SAVE :: IX = 0
      INTEGER K1
!
!                   IF NOT FIRST CALL, THEN SKIP SETTING SEED
      IF( IX .EQ. 0) IX = IDUM
!                   WRITE NUMBER AS ALPHA*2**16 + BETA
      K1 = IX / B
      IX = A*( IX - K1*B) - K1*C
!                   ABOVE DOES A*IX MOD B -K1*C
      IF( IX .LT. 0 ) IX = IX + P 
!                   RP BELOW IS RECIPROCAL OF P
      RAN = REAL(IX)*RP
      RETURN
      END FUNCTION RAN
! *** end of filename spkbl5.f95        *********************