Commit a43c7fbe authored by Thomas Jahns's avatar Thomas Jahns 🤸
Browse files

Add function to generate random integers.

parent ffb16f6e
......@@ -52,6 +52,7 @@ libscalesppm_a_SOURCES=ppm/scales_ppm.f90 \
solver/solver_all.f90 \
\
core/core.c core/core.h \
core/yarandom.c core/yarandom.h \
core/ppm_base.f90 core/ppm_set_default_comm.f90 \
core/ppm_extents.f90 core/ppm_extents_c.c \
core/ppm_strided_extents.f90 core/ppm_strided_extents_c.c \
......
......@@ -53,9 +53,13 @@
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#ifdef _OPENMP
#include <omp.h>
#endif
#include <stdio.h>
#include <cfortran.h>
#include "core.h"
#include "yarandom.h"
MPI_Comm PPM_default_comm =
#ifdef USE_MPI
......@@ -126,6 +130,43 @@ PPM_abort_default(MPI_Comm *comm, const char *msg, const char *source, int line)
PPM_abort_func PPM_abort = PPM_abort_default;
int
IRand()
{
int i = (int)ppm_ya_random();
if (i == -2147483648) i = 0;
return i;
}
FCALLSCFUN0(INT,IRand,IRAND,irand)
void
aIRand(int *a, int n)
{
int i;
for (i = 0; i < n; ++i)
a[i] = IRand();
}
FCALLSCSUB2(aIRand,A_RAND_I,a_rand_i,PINT,INT)
void
initIRand_f(MPI_Fint *comm_f, int random_seed)
{
MPI_Comm comm_c;
#if defined(USE_MPI) && defined(HAVE_MPI_COMM_F2C)
int flag = 0;
if (MPI_Initialized(&flag) == MPI_SUCCESS && flag)
comm_c = MPI_Comm_f2c((MPI_Fint)*comm_f);
#else
comm_c = *comm_f;
#endif
ppm_ya_rand_init(&comm_c, random_seed);
}
FCALLSCSUB2(initIRand_f,INITIALIZE_IRAND,initialize_irand,PVOID,INT)
/*
* Local Variables:
* license-project-url: "https://www.dkrz.de/redmine/projects/show/scales-ppm"
......
......@@ -50,6 +50,8 @@
*
* Code:
*/
#ifndef PPM_CORE_H
#define PPM_CORE_H
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
......@@ -69,6 +71,7 @@ extern MPI_Comm PPM_default_comm;
extern void
PPM_set_default_comm(MPI_Comm comm);
#endif
/*
* Local Variables:
* license-project-url: "https://www.dkrz.de/redmine/projects/show/scales-ppm"
......
......@@ -60,8 +60,34 @@ MODULE ppm_base
CHARACTER(*), INTENT(in) :: msg, source
END SUBROUTINE ppm_abort
END INTERFACE
!> irand returns integers in the range irand_min..irand_max
INTEGER, PARAMETER :: irand_min=-2147483647, irand_max=2147483647
!> this function shall be implemented in an OpenMP-thread-safe means,
!! currently it's not
INTERFACE
FUNCTION irand()
INTEGER :: irand
END FUNCTION irand
END INTERFACE
INTERFACE
SUBROUTINE initialize_irand(comm, random_seed)
INTEGER, INTENT(inout) :: comm
INTEGER, INTENT(in) :: random_seed
END SUBROUTINE initialize_irand
END INTERFACE
!> unfortunately, Fortrans random number generator is only prepared
!! to produce REAL-type results, this add similar capabilities for
!! INTEGER results
INTERFACE a_rand
SUBROUTINE a_rand_i(a, n)
INTEGER, INTENT(out) :: a(*)
INTEGER, INTENT(in) :: n
END SUBROUTINE a_rand_i
MODULE PROCEDURE a_rand_i_1d
END INTERFACE a_rand
PUBLIC :: ppm_default_comm, set_default_comm
PUBLIC :: abort_ppm, assertion
PUBLIC :: irand, initialize_irand, a_rand
CONTAINS
SUBROUTINE abort_ppm(msg, source, line, comm)
CHARACTER(len=*), INTENT(in) :: source, msg
......@@ -108,6 +134,11 @@ CONTAINS
END IF
END SUBROUTINE assertion
SUBROUTINE a_rand_i_1d(a)
INTEGER, INTENT(out) :: a(:)
CALL a_rand(a(1:SIZE(a)), SIZE(a))
END SUBROUTINE a_rand_i_1d
END MODULE ppm_base
!
! Local Variables:
......
/* Some changes for interoperability with OpenMP
* Copyright (c) 2011 by Thomas Jahns <jahns@dkrz.de>
*/
/* yarandom.c -- Yet Another Random Number Generator.
* Copyright (c) 1997, 1998, 2003 by Jamie Zawinski <jwz@jwz.org>
*
* Permission to use, copy, modify, distribute, and sell this software and its
* documentation for any purpose is hereby granted without fee, provided that
* the above copyright notice appear in all copies and that both that
* copyright notice and this permission notice appear in supporting
* documentation. No representations are made about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*/
/* The unportable mess that is rand(), random(), drand48() and friends led me
to ask Phil Karlton <karlton@netscape.com> what the Right Thing to Do was.
He responded with this. It is non-cryptographically secure, reasonably
random (more so than anything that is in any C library), and very fast.
I don't understand how it works at all, but he says "look at Knuth,
Vol. 2 (original edition), page 26, Algorithm A. In this case n=55,
k=20 and m=2^32."
So there you have it.
---------------------------
Note: xlockmore 4.03a10 uses this very simple RNG:
if ((seed = seed % 44488 * 48271 - seed / 44488 * 3399) < 0)
seed += 2147483647;
return seed-1;
of which it says
``Dr. Park's algorithm published in the Oct. '88 ACM "Random Number
Generators: Good Ones Are Hard To Find" His version available at
ftp://cs.wm.edu/pub/rngs.tar Present form by many authors.''
Karlton says: ``the usual problem with that kind of RNG turns out to
be unexepected short cycles for some word lengths.''
Karlton's RNG is faster, since it does three adds and two stores, while the
xlockmore RNG does two multiplies, two divides, three adds, and one store.
Compiler optimizations make a big difference here:
gcc -O: difference is 1.2x.
gcc -O2: difference is 1.4x.
gcc -O3: difference is 1.5x.
SGI cc -O: difference is 2.4x.
SGI cc -O2: difference is 2.4x.
SGI cc -O3: difference is 5.1x.
Irix 6.2; Indy r5k; SGI cc version 6; gcc version 2.7.2.1.
*/
#include <stdlib.h>
#include <unistd.h> /* for getpid() */
#include <sys/time.h> /* for gettimeofday() */
#include <inttypes.h>
#ifdef _OPENMP
#include <omp.h>
#endif
#include "core/yarandom.h"
# undef gt_ya_rand_init
/* The following 'random' numbers are taken from CRC, 18th Edition, page 622.
Each array element was taken from the corresponding line in the table,
except that a[0] was from line 100. 8s and 9s in the table were simply
skipped. The high order digit was taken mod 4.
*/
enum {
VectorSize = 55
};
static const uint32_t a_source[VectorSize] = {
035340171546, 010401501101, 022364657325, 024130436022, 002167303062, /* 5 */
037570375137, 037210607110, 016272055420, 023011770546, 017143426366, /* 10 */
014753657433, 021657231332, 023553406142, 004236526362, 010365611275, /* 14 */
007117336710, 011051276551, 002362132524, 001011540233, 012162531646, /* 20 */
007056762337, 006631245521, 014164542224, 032633236305, 023342700176, /* 25 */
002433062234, 015257225043, 026762051606, 000742573230, 005366042132, /* 30 */
012126416411, 000520471171, 000725646277, 020116577576, 025765742604, /* 35 */
007633473735, 015674255275, 017555634041, 006503154145, 021576344247, /* 40 */
014577627653, 002707523333, 034146376720, 030060227734, 013765414060, /* 45 */
036072251540, 007255221037, 024364674123, 006200353166, 010126373326, /* 50 */
015664104320, 016401041535, 016215305520, 033115351014, 017411670323 /* 55 */
};
#ifdef _OPENMP
static uint32_t (*a)[VectorSize];
static int (*vidx)[2];
#else
static uint32_t a[1][VectorSize];
static int vidx[1][2];
#endif
uint32_t
ppm_ya_random (void)
{
int ret, tid
#ifdef _OPENMP
= omp_get_thread_num();
#else
= 0;
#endif
int i1 = vidx[tid][0], i2 = vidx[tid][1];
ret = a[tid][i1] + a[tid][i2];
a[tid][i1] = ret;
if (++i1 >= VectorSize) i1 = 0;
if (++i2 >= VectorSize) i2 = 0;
vidx[tid][0] = i1;
vidx[tid][1] = i2;
return ret;
}
unsigned int
ppm_ya_rand_init(MPI_Comm *comm, int seed_arg)
{
unsigned seed = (unsigned)seed_arg;
int i, tid;
#ifdef _OPENMP
tid = omp_get_thread_num();
#else
tid = 0;
#endif
if (seed == 0)
{
struct timeval tp;
struct timezone tzp;
gettimeofday(&tp, &tzp);
/* ignore overflow */
seed = (999*tp.tv_sec) + (1001*tp.tv_usec) + (1003 * getpid()) + tid;
}
#ifdef _OPENMP
#pragma omp single
{
int n = omp_get_num_threads();
if (vidx) free(vidx);
if (a) free(a);
vidx = malloc(n * sizeof(*vidx));
a = malloc(n * sizeof(*a));
if (!vidx || !a)
PPM_abort(comm, "error in ppm_ya_rand_init", __FILE__, __LINE__);
}
#endif
/* allow to call gt_ya_rand_init() multiple times */
for (i = 0; i < VectorSize; i++)
a[tid][i] = a_source[i];
a[tid][0] += seed;
for (i = 1; i < VectorSize; i++)
{
seed = a[tid][i-1]*1001 + seed*999;
a[tid][i] += seed;
}
vidx[tid][0] = a[tid][0] % VectorSize;
vidx[tid][1] = (vidx[tid][0] + 024) % VectorSize;
return seed;
}
/* xscreensaver, Copyright (c) 1997, 1998, 2003 by Jamie Zawinski <jwz@jwz.org>
*
* Permission to use, copy, modify, distribute, and sell this software and its
* documentation for any purpose is hereby granted without fee, provided that
* the above copyright notice appear in all copies and that both that
* copyright notice and this permission notice appear in supporting
* documentation. No representations are made about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*/
#ifndef YARANDOM_H
#define YARANDOM_H
#include <inttypes.h>
#include "core.h"
#undef random
#undef rand
#undef drand48
#undef srandom
#undef srand
#undef srand48
#undef frand
#undef RAND_MAX
uint32_t ppm_ya_random (void);
unsigned int ppm_ya_rand_init(MPI_Comm *comm, int);
#define RAND_MAX 0x7FFFFFFF
#define random() ((long) (ppm_ya_random() & RAND_MAX))
/*#define srandom(i) ya_rand_init(0)*/
/* Define these away to keep people from using the wrong APIs in GenomeTools.
*/
#define rand __ERROR_use_random_not_rand_in_GenomeTools__
#define drand48 __ERROR_use_frand_not_drand48_in_GenomeTools__
#define srandom __ERROR_do_not_call_srandom_in_GenomeTools__
#define srand __ERROR_do_not_call_srand_in_GenomeTools__
#define srand48 __ERROR_do_not_call_srand48_in_GenomeTools__
#if defined (__GNUC__) && (__GNUC__ >= 2)
/* Implement frand using GCC's statement-expression extension. */
# define frand(f) \
__extension__ \
({ double tmp = ((((double) random()) * ((double) (f))) / \
((double) ((unsigned int)~0))); \
tmp < 0 ? (-tmp) : tmp; })
#else /* not GCC2 - implement frand using a global variable.*/
static double _frand_tmp_;
# define frand(f) \
(_frand_tmp_ = ((((double) random()) * ((double) (f))) / \
((double) ((unsigned int)~0))), \
_frand_tmp_ < 0 ? (-_frand_tmp_) : _frand_tmp_)
#endif /* not GCC2 */
#endif
......@@ -43,7 +43,11 @@
! Code:
!
MODULE scales_ppm
USE ppm_base, ONLY: abort_ppm, ppm_default_comm, set_default_comm
#ifdef _OPENMP
USE omp_lib
#endif
USE ppm_base, ONLY: abort_ppm, ppm_default_comm, set_default_comm, &
initialize_irand, irand, ppm_default_comm
USE ppm_extents, ONLY: extent, extent_size, extent_end, &
extent_start, rebased_extent, &
extent_set_iinterval, extent_from_iinterval, char, iinterval, &
......@@ -74,8 +78,9 @@ MODULE scales_ppm
PUBLIC :: abort_ppm, initialize_scales_ppm
PUBLIC :: char
CONTAINS
SUBROUTINE initialize_scales_ppm(default_comm)
INTEGER, OPTIONAL, INTENT(in) :: default_comm
SUBROUTINE initialize_scales_ppm(default_comm, random_seed)
INTEGER, OPTIONAL, INTENT(in) :: default_comm, random_seed
INTEGER :: random_seed_arg
IF (PRESENT(default_comm)) THEN
......@@ -83,7 +88,16 @@ CONTAINS
CALL set_default_comm(default_comm)
END IF
IF (PRESENT(random_seed)) THEN; random_seed_arg=random_seed; ELSE
#ifdef _OPENMP
random_seed_arg = 9 + omp_get_thread_num()
#else
random_seed_arg = 9
#endif
END IF
CALL initialize_irand(ppm_default_comm, random_seed_arg)
END SUBROUTINE initialize_scales_ppm
END MODULE scales_ppm
!
! Local Variables:
......
......@@ -32,7 +32,7 @@ AM_FCFLAGS = $(FPP_INCOPT)$(top_srcdir)/include/f90 \
$(FPP_INCOPT)../include/f77 \
$(FC_MOD_FLAG)../src $(MPI_FC_INCLUDE)
noinst_PROGRAMS = test_qsort test_uniform_partition
noinst_PROGRAMS = test_qsort test_uniform_partition test_irand
if USE_MPI
noinst_PROGRAMS+= test_strided_extents
endif
......@@ -49,6 +49,10 @@ test_strided_extents_SOURCES = test_strided_extents.f90
test_strided_extents_LDADD = ../src/libscalesppm.a $(MPI_FC_LIB)
test_irand_SOURCES = test_irand.f90
test_irand_LDADD = ../src/libscalesppm.a $(MPI_FC_LIB)
EXTRA_DIST=test_strided_extents.f90
./$(DEPDIR)/FC.deps: $(SOURCES) Makefile
......
PROGRAM test_irand
#ifdef _OPENMP
USE omp_lib
#endif
USE ppm_base, ONLY: irand, initialize_irand, ppm_default_comm, a_rand
IMPLICIT NONE
INTEGER, ALLOCATABLE :: irand_res(:, :)
INTEGER :: num_rands, num_threads, tid, i, ref
LOGICAL :: inequality_found
inequality_found = .FALSE.
!$omp parallel private(tid, num_threads, i, ref) &
!$omp reduction(.OR.:inequality_found)
#ifdef _OPENMP
num_threads = omp_get_num_threads()
tid = omp_get_thread_num()
#else
num_threads = 1
tid = 0
#endif
CALL initialize_irand(ppm_default_comm, 9)
!$omp master
num_rands = MOD(IAND(irand(), 2147483647), 10000)
ALLOCATE(irand_res(num_rands, 0:num_threads-1))
!$omp end master
CALL initialize_irand(ppm_default_comm, 9)
DO i = 1, num_rands
irand_res(i, tid) = irand()
END DO
!$omp barrier
ref = MOD(tid + 1, num_threads)
inequality_found = inequality_found &
.OR. ANY(irand_res(:, tid) /= irand_res(:, ref))
CALL a_rand(irand_res(:, tid))
!$omp barrier
inequality_found = inequality_found &
.OR. ANY(irand_res(:, tid) /= irand_res(:, ref))
!$omp end parallel
IF (inequality_found) THEN
WRITE (0, *) "Inequality found in random numbers!"
STOP 1
END IF
END PROGRAM test_irand
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment