libgfortran/intrinsics/random.c

   1 /* Implementation of the RANDOM intrinsics
   2    Copyright 2002, 2004, 2005 Free Software Foundation, Inc.
   3    Contributed by Lars Segerlund <seger@linuxmail.org>
   4    and Steve Kargl.
   5
   6 This file is part of the GNU Fortran 95 runtime library (libgfortran).
   7
   8 Libgfortran is free software; you can redistribute it and/or
   9 modify it under the terms of the GNU General Public
  10 License as published by the Free Software Foundation; either
  11 version 2 of the License, or (at your option) any later version.
  12
  13 In addition to the permissions in the GNU General Public License, the
  14 Free Software Foundation gives you unlimited permission to link the
  15 compiled version of this file into combinations with other programs,
  16 and to distribute those combinations without any restriction coming
  17 from the use of this file.  (The General Public License restrictions
  18 do apply in other respects; for example, they cover modification of
  19 the file, and distribution when not linked into a combine
  20 executable.)
  21
  22 Ligbfortran is distributed in the hope that it will be useful,
  23 but WITHOUT ANY WARRANTY; without even the implied warranty of
  24 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  25 GNU General Public License for more details.
  26
  27 You should have received a copy of the GNU General Public
  28 License along with libgfortran; see the file COPYING.  If not,
  29 write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
  30 Boston, MA 02110-1301, USA.  */
  31
  32 #include "config.h"
  33 #include "libgfortran.h"
  34 #include "../io/io.h"
  35
  36 extern void random_r4 (GFC_REAL_4 *);
  37 iexport_proto(random_r4);
  38
  39 extern void random_r8 (GFC_REAL_8 *);
  40 iexport_proto(random_r8);
  41
  42 extern void arandom_r4 (gfc_array_r4 *);
  43 export_proto(arandom_r4);
  44
  45 extern void arandom_r8 (gfc_array_r8 *);
  46 export_proto(arandom_r8);
  47
  48 #ifdef __GTHREAD_MUTEX_INIT
  49 static __gthread_mutex_t random_lock = __GTHREAD_MUTEX_INIT;
  50 #else
  51 static __gthread_mutex_t random_lock;
  52 #endif
  53
  54 #if 0
  55
  56 /*  The Mersenne Twister code is currently commented out due to
  57
  58     (1) Simple user specified seeds lead to really bad sequences for
  59         nearly 100000 random numbers.
  60     (2) open(), read(), and close() are not properly declared via header
  61         files.
  62     (3) The global index i is abused and causes unexpected behavior with
  63         GET and PUT.
  64     (4) See PR 15619.
  65
  66   The algorithm was taken from the paper :
  67
  68         Mersenne Twister:       623-dimensionally equidistributed
  69                                 uniform pseudorandom generator.
  70
  71         by:     Makoto Matsumoto
  72                 Takuji Nishimura
  73
  74         Which appeared in the: ACM Transactions on Modelling and Computer
  75         Simulations: Special Issue on Uniform Random Number
  76         Generation. ( Early in 1998 ).  */
  77
  78
  79 #include <stdio.h>
  80 #include <stdlib.h>
  81 #include <sys/types.h>
  82 #include <sys/stat.h>
  83 #include <fcntl.h>
  84
  85 #ifdef HAVE_UNISTD_H
  86 #include <unistd.h>
  87 #endif
  88
  89 /*Use the 'big' generator by default ( period -> 2**19937 ).  */
  90
  91 #define MT19937
  92
  93 /* Define the necessary constants for the algorithm.  */
  94
  95 #ifdef  MT19937
  96 enum constants
  97 {
  98   N = 624, M = 397, R = 19, TU = 11, TS = 7, TT = 15, TL = 17
  99 };
 100 #define M_A     0x9908B0DF
 101 #define T_B     0x9D2C5680
 102 #define T_C     0xEFC60000
 103 #else
 104 enum constants
 105 {
 106   N = 351, M = 175, R = 19, TU = 11, TS = 7, TT = 15, TL = 17
 107 };
 108 #define M_A     0xE4BD75F5
 109 #define T_B     0x655E5280
 110 #define T_C     0xFFD58000
 111 #endif
 112
 113 static int i = N;
 114 static unsigned int seed[N];
 115
 116 /* This is the routine which handles the seeding of the generator,
 117    and also reading and writing of the seed.  */
 118
 119 void
 120 random_seed (GFC_INTEGER_4 *size, gfc_array_i4 *put, gfc_array_i4 *get)
 121 {
 122   __gthread_mutex_lock (&random_lock);
 123
 124   /* Initialize the seed in system dependent manner.  */
 125   if (get == NULL && put == NULL && size == NULL)
 126     {
 127       int fd;
 128       fd = open ("/dev/urandom", O_RDONLY);
 129       if (fd < 0)
 130         {
 131           /* We dont have urandom.  */
 132           GFC_UINTEGER_4 s = (GFC_UINTEGER_4) seed;
 133           for (i = 0; i < N; i++)
 134             {
 135               s = s * 29943829 - 1;
 136               seed[i] = s;
 137             }
 138         }
 139       else
 140         {
 141           /* Using urandom, might have a length issue.  */
 142           read (fd, &seed[0], sizeof (GFC_UINTEGER_4) * N);
 143           close (fd);
 144           i = N;
 145         }
 146       goto return_unlock;
 147     }
 148
 149   /* Return the size of the seed */
 150   if (size != NULL)
 151     {
 152       *size = N;
 153       goto return_unlock;
 154     }
 155
 156   /* if we have gotten to this pount we have a get or put
 157    * now we check it the array fulfills the demands in the standard .
 158    */
 159
 160   /* Set the seed to PUT data */
 161   if (put != NULL)
 162     {
 163       /* if the rank of the array is not 1 abort */
 164       if (GFC_DESCRIPTOR_RANK (put) != 1)
 165         abort ();
 166
 167       /* if the array is too small abort */
 168       if (((put->dim[0].ubound + 1 - put->dim[0].lbound)) < N)
 169         abort ();
 170
 171       /* If this is the case the array is a temporary */
 172       if (put->dim[0].stride == 0)
 173         goto return_unlock;
 174
 175       /*  This code now should do correct strides. */
 176       for (i = 0; i < N; i++)
 177         seed[i] = put->data[i * put->dim[0].stride];
 178     }
 179
 180   /* Return the seed to GET data */
 181   if (get != NULL)
 182     {
 183       /* if the rank of the array is not 1 abort */
 184       if (GFC_DESCRIPTOR_RANK (get) != 1)
 185         abort ();
 186
 187       /* if the array is too small abort */
 188       if (((get->dim[0].ubound + 1 - get->dim[0].lbound)) < N)
 189         abort ();
 190
 191       /* If this is the case the array is a temporary */
 192       if (get->dim[0].stride == 0)
 193         goto return_unlock;
 194
 195       /*  This code now should do correct strides. */
 196       for (i = 0; i < N; i++)
 197         get->data[i * get->dim[0].stride] = seed[i];
 198     }
 199
 200  random_unlock:
 201   __gthread_mutex_unlock (&random_lock);
 202 }
 203 iexport(random_seed);
 204
 205 /* Here is the internal routine which generates the random numbers
 206    in 'batches' based upon the need for a new batch.
 207    It's an integer based routine known as 'Mersenne Twister'.
 208    This implementation still lacks 'tempering' and a good verification,
 209    but gives very good metrics.  */
 210
 211 static void
 212 random_generate (void)
 213 {
 214   /* 32 bits.  */
 215   GFC_UINTEGER_4 y;
 216
 217   /* Generate batch of N.  */
 218   int k, m;
 219   for (k = 0, m = M; k < N - 1; k++)
 220     {
 221       y = (seed[k] & (-1 << R)) | (seed[k + 1] & ((1u << R) - 1));
 222       seed[k] = seed[m] ^ (y >> 1) ^ (-(GFC_INTEGER_4) (y & 1) & M_A);
 223       if (++m >= N)
 224         m = 0;
 225     }
 226
 227   y = (seed[N - 1] & (-1 << R)) | (seed[0] & ((1u << R) - 1));
 228   seed[N - 1] = seed[M - 1] ^ (y >> 1) ^ (-(GFC_INTEGER_4) (y & 1) & M_A);
 229   i = 0;
 230 }
 231
 232 /* A routine to return a REAL(KIND=4).  */
 233
 234 void
 235 random_r4 (GFC_REAL_4 * harv)
 236 {
 237   __gthread_mutex_lock (&random_lock);
 238
 239   /* Regenerate if we need to.  */
 240   if (i >= N)
 241     random_generate ();
 242
 243   /* Convert uint32 to REAL(KIND=4).  */
 244   *harv = (GFC_REAL_4) ((GFC_REAL_4) (GFC_UINTEGER_4) seed[i++] /
 245                         (GFC_REAL_4) (~(GFC_UINTEGER_4) 0));
 246   __gthread_mutex_unlock (&random_lock);
 247 }
 248 iexport(random_r4);
 249
 250 /* A routine to return a REAL(KIND=8).  */
 251
 252 void
 253 random_r8 (GFC_REAL_8 * harv)
 254 {
 255   __gthread_mutex_lock (&random_lock);
 256
 257   /* Regenerate if we need to, may waste one 32-bit value.  */
 258   if ((i + 1) >= N)
 259     random_generate ();
 260
 261   /* Convert two uint32 to a REAL(KIND=8).  */
 262   *harv = ((GFC_REAL_8) ((((GFC_UINTEGER_8) seed[i+1]) << 32) + seed[i])) /
 263           (GFC_REAL_8) (~(GFC_UINTEGER_8) 0);
 264   i += 2;
 265   __gthread_mutex_unlock (&random_lock);
 266 }
 267 iexport(random_r8);
 268
 269 /* Code to handle arrays will follow here.  */
 270
 271 /* REAL(KIND=4) REAL array.  */
 272
 273 void
 274 arandom_r4 (gfc_array_r4 * harv)
 275 {
 276   index_type count[GFC_MAX_DIMENSIONS];
 277   index_type extent[GFC_MAX_DIMENSIONS];
 278   index_type stride[GFC_MAX_DIMENSIONS];
 279   index_type stride0;
 280   index_type dim;
 281   GFC_REAL_4 *dest;
 282   int n;
 283
 284   dest = harv->data;
 285
 286   if (harv->dim[0].stride == 0)
 287     harv->dim[0].stride = 1;
 288
 289   dim = GFC_DESCRIPTOR_RANK (harv);
 290
 291   for (n = 0; n < dim; n++)
 292     {
 293       count[n] = 0;
 294       stride[n] = harv->dim[n].stride;
 295       extent[n] = harv->dim[n].ubound + 1 - harv->dim[n].lbound;
 296       if (extent[n] <= 0)
 297         return;
 298     }
 299
 300   stride0 = stride[0];
 301
 302   __gthread_mutex_lock (&random_lock);
 303
 304   while (dest)
 305     {
 306       /* Set the elements.  */
 307
 308       /* regenerate if we need to */
 309       if (i >= N)
 310         random_generate ();
 311
 312       /* Convert uint32 to float in a hopefully g95 compiant manner */
 313       *dest = (GFC_REAL_4) ((GFC_REAL_4) (GFC_UINTEGER_4) seed[i++] /
 314                             (GFC_REAL_4) (~(GFC_UINTEGER_4) 0));
 315
 316       /* Advance to the next element.  */
 317       dest += stride0;
 318       count[0]++;
 319       /* Advance to the next source element.  */
 320       n = 0;
 321       while (count[n] == extent[n])
 322         {
 323           /* When we get to the end of a dimension,
 324              reset it and increment
 325              the next dimension.  */
 326           count[n] = 0;
 327           /* We could precalculate these products,
 328              but this is a less
 329              frequently used path so proabably not worth it.  */
 330           dest -= stride[n] * extent[n];
 331           n++;
 332           if (n == dim)
 333             {
 334               dest = NULL;
 335               break;
 336             }
 337           else
 338             {
 339               count[n]++;
 340               dest += stride[n];
 341             }
 342         }
 343     }
 344
 345   __gthread_mutex_unlock (&random_lock);
 346 }
 347
 348 /* REAL(KIND=8) array.  */
 349
 350 void
 351 arandom_r8 (gfc_array_r8 * harv)
 352 {
 353   index_type count[GFC_MAX_DIMENSIONS];
 354   index_type extent[GFC_MAX_DIMENSIONS];
 355   index_type stride[GFC_MAX_DIMENSIONS];
 356   index_type stride0;
 357   index_type dim;
 358   GFC_REAL_8 *dest;
 359   int n;
 360
 361   dest = harv->data;
 362
 363   if (harv->dim[0].stride == 0)
 364     harv->dim[0].stride = 1;
 365
 366   dim = GFC_DESCRIPTOR_RANK (harv);
 367
 368   for (n = 0; n < dim; n++)
 369     {
 370       count[n] = 0;
 371       stride[n] = harv->dim[n].stride;
 372       extent[n] = harv->dim[n].ubound + 1 - harv->dim[n].lbound;
 373       if (extent[n] <= 0)
 374         return;
 375     }
 376
 377   stride0 = stride[0];
 378
 379   __gthread_mutex_lock (&random_lock);
 380
 381   while (dest)
 382     {
 383       /* Set the elements.  */
 384
 385       /* regenerate if we need to, may waste one 32-bit value */
 386       if ((i + 1) >= N)
 387         random_generate ();
 388
 389       /* Convert two uint32 to a REAL(KIND=8).  */
 390       *dest = ((GFC_REAL_8) ((((GFC_UINTEGER_8) seed[i+1]) << 32) + seed[i])) /
 391               (GFC_REAL_8) (~(GFC_UINTEGER_8) 0);
 392       i += 2;
 393
 394       /* Advance to the next element.  */
 395       dest += stride0;
 396       count[0]++;
 397       /* Advance to the next source element.  */
 398       n = 0;
 399       while (count[n] == extent[n])
 400         {
 401           /* When we get to the end of a dimension,
 402              reset it and increment
 403              the next dimension.  */
 404           count[n] = 0;
 405           /* We could precalculate these products,
 406              but this is a less
 407              frequently used path so proabably not worth it.  */
 408           dest -= stride[n] * extent[n];
 409           n++;
 410           if (n == dim)
 411             {
 412               dest = NULL;
 413               break;
 414             }
 415           else
 416             {
 417               count[n]++;
 418               dest += stride[n];
 419             }
 420         }
 421     }
 422
 423   __gthread_mutex_unlock (&random_lock);
 424 }
 425
 426 #else
 427
 428 /* George Marsaglia's KISS (Keep It Simple Stupid) random number generator.
 429
 430    This PRNG combines:
 431
 432    (1) The congruential generator x(n)=69069*x(n-1)+1327217885 with a period
 433        of 2^32,
 434    (2) A 3-shift shift-register generator with a period of 2^32-1,
 435    (3) Two 16-bit multiply-with-carry generators with a period of
 436        597273182964842497 > 2^59.
 437
 438    The overall period exceeds 2^123.
 439
 440    http://www.ciphersbyritter.com/NEWS4/RANDC.HTM#369F6FCA.74C7C041@stat.fsu.edu
 441
 442    The above web site has an archive of a newsgroup posting from George
 443    Marsaglia with the statement:
 444
 445    Subject: Random numbers for C: Improvements.
 446    Date: Fri, 15 Jan 1999 11:41:47 -0500
 447    From: George Marsaglia <geo@stat.fsu.edu>
 448    Message-ID: <369F6FCA.74C7C041@stat.fsu.edu>
 449    References: <369B5E30.65A55FD1@stat.fsu.edu>
 450    Newsgroups: sci.stat.math,sci.math,sci.math.numer-analysis
 451    Lines: 93
 452
 453    As I hoped, several suggestions have led to
 454    improvements in the code for RNG's I proposed for
 455    use in C. (See the thread "Random numbers for C: Some
 456    suggestions" in previous postings.) The improved code
 457    is listed below.
 458
 459    A question of copyright has also been raised.  Unlike
 460    DIEHARD, there is no copyright on the code below. You
 461    are free to use it in any way you want, but you may
 462    wish to acknowledge the source, as a courtesy.
 463
 464 "There is no copyright on the code below." included the original
 465 KISS algorithm.  */
 466
 467 #define GFC_SL(k, n)    ((k)^((k)<<(n)))
 468 #define GFC_SR(k, n)    ((k)^((k)>>(n)))
 469
 470 static const GFC_INTEGER_4 kiss_size = 4;
 471 #define KISS_DEFAULT_SEED {123456789, 362436069, 521288629, 916191069}
 472 static const GFC_UINTEGER_4 kiss_default_seed[4] = KISS_DEFAULT_SEED;
 473 static GFC_UINTEGER_4 kiss_seed[4] = KISS_DEFAULT_SEED;
 474
 475 /* kiss_random_kernel() returns an integer value in the range of
 476    (0, GFC_UINTEGER_4_HUGE].  The distribution of pseudorandom numbers
 477    should be uniform.  */
 478
 479 static GFC_UINTEGER_4
 480 kiss_random_kernel(void)
 481 {
 482   GFC_UINTEGER_4 kiss;
 483
 484   kiss_seed[0] = 69069 * kiss_seed[0] + 1327217885;
 485   kiss_seed[1] = GFC_SL(GFC_SR(GFC_SL(kiss_seed[1],13),17),5);
 486   kiss_seed[2] = 18000 * (kiss_seed[2] & 65535) + (kiss_seed[2] >> 16);
 487   kiss_seed[3] = 30903 * (kiss_seed[3] & 65535) + (kiss_seed[3] >> 16);
 488   kiss = kiss_seed[0] + kiss_seed[1] + (kiss_seed[2] << 16) + kiss_seed[3];
 489
 490   return kiss;
 491 }
 492
 493 /*  This function produces a REAL(4) value from the uniform distribution
 494     with range [0,1).  */
 495
 496 void
 497 random_r4 (GFC_REAL_4 *x)
 498 {
 499   GFC_UINTEGER_4 kiss;
 500
 501   __gthread_mutex_lock (&random_lock);
 502   kiss = kiss_random_kernel ();
 503   /* Burn a random number, so the REAL*4 and REAL*8 functions
 504      produce similar sequences of random numbers.  */
 505   kiss_random_kernel ();
 506   *x = normalize_r4_i4 (kiss, ~(GFC_UINTEGER_4) 0);
 507   __gthread_mutex_unlock (&random_lock);
 508 }
 509 iexport(random_r4);
 510
 511 /*  This function produces a REAL(8) value from the uniform distribution
 512     with range [0,1).  */
 513
 514 void
 515 random_r8 (GFC_REAL_8 *x)
 516 {
 517   GFC_UINTEGER_8 kiss;
 518
 519   __gthread_mutex_lock (&random_lock);
 520   kiss = ((GFC_UINTEGER_8)kiss_random_kernel ()) << 32;
 521   kiss += kiss_random_kernel ();
 522   *x = normalize_r8_i8 (kiss, ~(GFC_UINTEGER_8) 0);
 523   __gthread_mutex_unlock (&random_lock);
 524 }
 525 iexport(random_r8);
 526
 527 /*  This function fills a REAL(4) array with values from the uniform
 528     distribution with range [0,1).  */
 529
 530 void
 531 arandom_r4 (gfc_array_r4 *x)
 532 {
 533   index_type count[GFC_MAX_DIMENSIONS];
 534   index_type extent[GFC_MAX_DIMENSIONS];
 535   index_type stride[GFC_MAX_DIMENSIONS];
 536   index_type stride0;
 537   index_type dim;
 538   GFC_REAL_4 *dest;
 539   GFC_UINTEGER_4 kiss;
 540   int n;
 541
 542   dest = x->data;
 543
 544   if (x->dim[0].stride == 0)
 545     x->dim[0].stride = 1;
 546
 547   dim = GFC_DESCRIPTOR_RANK (x);
 548
 549   for (n = 0; n < dim; n++)
 550     {
 551       count[n] = 0;
 552       stride[n] = x->dim[n].stride;
 553       extent[n] = x->dim[n].ubound + 1 - x->dim[n].lbound;
 554       if (extent[n] <= 0)
 555         return;
 556     }
 557
 558   stride0 = stride[0];
 559
 560   __gthread_mutex_lock (&random_lock);
 561
 562   while (dest)
 563     {
 564       /* random_r4 (dest); */
 565       kiss = kiss_random_kernel ();
 566       /* Burn a random number, so the REAL*4 and REAL*8 functions
 567          produce similar sequences of random numbers.  */
 568       kiss_random_kernel ();
 569       *dest = normalize_r4_i4 (kiss, ~(GFC_UINTEGER_4) 0);
 570
 571       /* Advance to the next element.  */
 572       dest += stride0;
 573       count[0]++;
 574       /* Advance to the next source element.  */
 575       n = 0;
 576       while (count[n] == extent[n])
 577         {
 578           /* When we get to the end of a dimension, reset it and increment
 579              the next dimension.  */
 580           count[n] = 0;
 581           /* We could precalculate these products, but this is a less
 582              frequently used path so probably not worth it.  */
 583           dest -= stride[n] * extent[n];
 584           n++;
 585           if (n == dim)
 586             {
 587               dest = NULL;
 588               break;
 589             }
 590           else
 591             {
 592               count[n]++;
 593               dest += stride[n];
 594             }
 595         }
 596     }
 597   __gthread_mutex_unlock (&random_lock);
 598 }
 599
 600 /*  This function fills a REAL(8) array with values from the uniform
 601     distribution with range [0,1).  */
 602
 603 void
 604 arandom_r8 (gfc_array_r8 *x)
 605 {
 606   index_type count[GFC_MAX_DIMENSIONS];
 607   index_type extent[GFC_MAX_DIMENSIONS];
 608   index_type stride[GFC_MAX_DIMENSIONS];
 609   index_type stride0;
 610   index_type dim;
 611   GFC_REAL_8 *dest;
 612   GFC_UINTEGER_8 kiss;
 613   int n;
 614
 615   dest = x->data;
 616
 617   if (x->dim[0].stride == 0)
 618     x->dim[0].stride = 1;
 619
 620   dim = GFC_DESCRIPTOR_RANK (x);
 621
 622   for (n = 0; n < dim; n++)
 623     {
 624       count[n] = 0;
 625       stride[n] = x->dim[n].stride;
 626       extent[n] = x->dim[n].ubound + 1 - x->dim[n].lbound;
 627       if (extent[n] <= 0)
 628         return;
 629     }
 630
 631   stride0 = stride[0];
 632
 633   __gthread_mutex_lock (&random_lock);
 634
 635   while (dest)
 636     {
 637       /* random_r8 (dest); */
 638       kiss = ((GFC_UINTEGER_8)kiss_random_kernel ()) << 32;
 639       kiss += kiss_random_kernel ();
 640       *dest = normalize_r8_i8 (kiss, ~(GFC_UINTEGER_8) 0);
 641
 642       /* Advance to the next element.  */
 643       dest += stride0;
 644       count[0]++;
 645       /* Advance to the next source element.  */
 646       n = 0;
 647       while (count[n] == extent[n])
 648         {
 649           /* When we get to the end of a dimension, reset it and increment
 650              the next dimension.  */
 651           count[n] = 0;
 652           /* We could precalculate these products, but this is a less
 653              frequently used path so probably not worth it.  */
 654           dest -= stride[n] * extent[n];
 655           n++;
 656           if (n == dim)
 657             {
 658               dest = NULL;
 659               break;
 660             }
 661           else
 662             {
 663               count[n]++;
 664               dest += stride[n];
 665             }
 666         }
 667     }
 668   __gthread_mutex_unlock (&random_lock);
 669 }
 670
 671 /* random_seed is used to seed the PRNG with either a default
 672    set of seeds or user specified set of seeds.  random_seed
 673    must be called with no argument or exactly one argument.  */
 674
 675 void
 676 random_seed (GFC_INTEGER_4 *size, gfc_array_i4 *put, gfc_array_i4 *get)
 677 {
 678   int i;
 679
 680   __gthread_mutex_lock (&random_lock);
 681
 682   if (size == NULL && put == NULL && get == NULL)
 683     {
 684       /* From the standard: "If no argument is present, the processor assigns
 685          a processor-dependent value to the seed."  */
 686       kiss_seed[0] = kiss_default_seed[0];
 687       kiss_seed[1] = kiss_default_seed[1];
 688       kiss_seed[2] = kiss_default_seed[2];
 689       kiss_seed[3] = kiss_default_seed[3];
 690     }
 691
 692   if (size != NULL)
 693     *size = kiss_size;
 694
 695   if (put != NULL)
 696     {
 697       /* If the rank of the array is not 1, abort.  */
 698       if (GFC_DESCRIPTOR_RANK (put) != 1)
 699         runtime_error ("Array rank of PUT is not 1.");
 700
 701       /* If the array is too small, abort.  */
 702       if (((put->dim[0].ubound + 1 - put->dim[0].lbound)) < kiss_size)
 703         runtime_error ("Array size of PUT is too small.");
 704
 705       if (put->dim[0].stride == 0)
 706         put->dim[0].stride = 1;
 707
 708       /*  This code now should do correct strides.  */
 709       for (i = 0; i < kiss_size; i++)
 710         kiss_seed[i] =(GFC_UINTEGER_4) put->data[i * put->dim[0].stride];
 711     }
 712
 713   /* Return the seed to GET data.  */
 714   if (get != NULL)
 715     {
 716       /* If the rank of the array is not 1, abort.  */
 717       if (GFC_DESCRIPTOR_RANK (get) != 1)
 718         runtime_error ("Array rank of GET is not 1.");
 719
 720       /* If the array is too small, abort.  */
 721       if (((get->dim[0].ubound + 1 - get->dim[0].lbound)) < kiss_size)
 722         runtime_error ("Array size of GET is too small.");
 723
 724       if (get->dim[0].stride == 0)
 725         get->dim[0].stride = 1;
 726
 727       /*  This code now should do correct strides.  */
 728       for (i = 0; i < kiss_size; i++)
 729         get->data[i * get->dim[0].stride] = (GFC_INTEGER_4) kiss_seed[i];
 730     }
 731
 732   __gthread_mutex_unlock (&random_lock);
 733 }
 734 iexport(random_seed);
 735
 736 #endif /* mersenne twister */
 737
 738 #ifndef __GTHREAD_MUTEX_INIT
 739 static void __attribute__((constructor))
 740 init (void)
 741 {
 742   __GTHREAD_MUTEX_INIT_FUNCTION (&random_lock);
 743 }
 744 #endif