!     Last change:  DOS   2 Aug 2000    9:38 pm
!                                       *** copyright 2000 ***
! *** filename spkblh.f95               *** John F. Monahan **
!                                       **********************
      program pspkblh
!  test of spkbll "spikeball" -- generates random points antipodally
!   on the surface of a p-dimensional sphere
!
!  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(10) :: mu,v
      real, dimension(nsc) :: kap,s,ss
      real, dimension(10,20) :: a      ! points are columns
      real t
      integer i,j,k,p,p2
!
      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 spkblh(p,a)
      integer, intent(in) :: p
      real, dimension(:,:), intent(out) :: a
      end subroutine spkblh
      end interface
!     
      open( unit=6, file='spkblh.out')
!
!                   initialize random number generator
      t = ran(5151917)
!                   vary dimension from 2 to 6
      do p = 2,6
      write(6,20) p,n
      p2 = 2 * p
!                   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 spkblh(p,a)
!                   generate this spikeball on surface of sphere
      do j = 1,p2
!                   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
      do k = 1,nsc
      v(k) = v(k) + exp( - kap(k)*t )
      end do  ! loop on k
!                   end of rotations
      end do  ! loop on j
      do k = 1,nsc
      v(k) = v(k) / (2.*p)
      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 pspkblh
      subroutine spkblh(p,a)
!  compute random rotation of antipodal points on surface of a sphere
!  in p dimensions  -- gives  2*p points stored as columns in a
!
!   p    integer    dimension of sphere
!   a    real       points on surface: point k is in a( ,k)
!
! J F Monahan (June 1994) Dept of Statistics, N C State U, Raleigh 27695
!     recoded  June 1998
!     Recoded February, April 2000 for Fortran 95
      implicit none
      integer, intent(in) :: p
      real, dimension(:,:), intent(out) :: a      ! dummy dimensions
      real, dimension(p) :: x
      integer i,ii,j,k,p2
      real s,d
      real gnroul
!                   compute some constants
      p2 = p + p
!                   initialize to zero first
      a = 0.
!                   now put 1 and -1 for poles
      do j = 1,p
      a(j,j) = 1.
      a(j,j+p) = -1.
      end do  ! loop on j
!                   matrix a is now I(p), -I(p) 
      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,p2
!                   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
!                   done with this rotation
      end do  ! loop on j
      end do  ! loop on ii
      return
      end subroutine spkblh
      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 spkblh.f95        *********************