MODULE rand_utility
  USE precisioninfo
  USE iso_c_binding
  USE error_codes
  IMPLICIT NONE
#include "macro.h"
  
  INTERFACE
     SUBROUTINE init_mt(d_status, d_texture) BIND(C, name="init_mt_fermi")
       USE iso_c_binding
       TYPE(c_ptr) :: d_status
       TYPE(c_ptr), DIMENSION(3) :: d_texture
     END SUBROUTINE init_mt
     TYPE(c_ptr) FUNCTION init_mt_status() BIND(C, name="init_mt_status_fermi")
       USE iso_c_binding
     END FUNCTION init_mt_status
     FUNCTION init_mt_texture(d_status) RESULT(res) BIND(C, name="init_mt_texture_fermi") 
       USE iso_c_binding
       TYPE(c_ptr), VALUE :: d_status
       TYPE(c_ptr) :: res
     END FUNCTION  init_mt_texture
     SUBROUTINE finalize_mt(d_status, d_texture) BIND(C, name="finalize_mt")
       USE iso_c_binding
       TYPE(c_ptr) :: d_status
       TYPE(c_ptr) :: d_texture
     END SUBROUTINE finalize_mt
     SUBROUTINE generate_mt_data(d_status, d_texture, num_data, d_data) BIND(C, name="generate_mt_data_fermi")
       USE iso_c_binding
       TYPE(c_ptr), VALUE :: d_status, d_texture
       INTEGER(c_int), VALUE :: num_data
       REAL(c_double), DIMENSION(*), device :: d_data
     END SUBROUTINE generate_mt_data
  END INTERFACE

CONTAINS
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! initialize random number generator (RNG)
  ! using Mersenne-Twister
  ! return: STATE, INFO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  SUBROUTINE init_rndgen(STATE, INFO)
    INTEGER, DIMENSION(:), ALLOCATABLE :: STATE
    INTEGER :: INFO
    !local
    INTEGER :: ii
    INTEGER :: newseed, nrun
    REAL :: tempsec, tempsec2
    INTEGER, PARAMETER :: LSTATE=633,LSEED=624
    INTEGER :: SEED(LSEED)
    !! Define random number generator parameters
    INFO = 0
    nrun = 1

    IF (.NOT. ALLOCATED(STATE)) ALLOCATE(STATE(LSTATE))

    !! Seed the random number generator
    tempsec = 372.d0*SECNDS(0.0)-1.d0
    tempsec2 = 2.d0*SECNDS(0.0)-1.d0
    newseed = int(tempsec) + mod(int(tempsec2),127)*975121
    seed(1) = (2*newseed+1)+nrun*10
    DO ii=2, LSEED
       SEED(ii) = seed(ii-1)+3*ii
    enddo
    PRINT *, "Use clock random seed"
    !!     Initialize the base generator. 
    CALL DRANDINITIALIZE(3,0,SEED,LSEED,STATE,LSTATE,INFO)

    RETURN
  END SUBROUTINE init_rndgen

  !============================================
  ! This function always use the same seed value
  ! which is designed to test the program
  !============================================
  SUBROUTINE init_rndgen_testonly(STATE, INFO, rand_seed)
    INTEGER, DIMENSION(:), ALLOCATABLE :: STATE
    INTEGER :: INFO
    INTEGER, OPTIONAL :: rand_seed
    !local
    INTEGER :: ii, nrun
    INTEGER, PARAMETER :: LSTATE=633,LSEED=624
    INTEGER :: SEED(LSEED)
    !! Define random number generator parameters
    INFO = 0
    nrun = 1

    IF (.NOT. ALLOCATED(STATE)) ALLOCATE(STATE(LSTATE))

    !! Seed the random number generator
    IF (PRESENT(rand_seed)) THEN
       seed(1) = rand_seed
    ELSE
       seed(1) = 1
    END IF
    DO ii=2, LSEED
       SEED(ii) = seed(ii-1)+3*ii
    enddo
    PRINT *, "Use fixed-random number"

    !!     Initialize the base generator. 
    CALL DRANDINITIALIZE(3,0,SEED,LSEED,STATE,LSTATE,INFO)

    RETURN
  END SUBROUTINE init_rndgen_testonly

END MODULE rand_utility

