]> git.ipfire.org Git - thirdparty/gcc.git/blob - libgfortran/intrinsics/random.c
re PR fortran/30964 (optional arguments to random_seed)
[thirdparty/gcc.git] / libgfortran / intrinsics / random.c
1 /* Implementation of the RANDOM intrinsics
2 Copyright 2002, 2004, 2005, 2006, 2007 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 <gthr.h>
35 #include <string.h>
36
37 extern void random_r4 (GFC_REAL_4 *);
38 iexport_proto(random_r4);
39
40 extern void random_r8 (GFC_REAL_8 *);
41 iexport_proto(random_r8);
42
43 extern void arandom_r4 (gfc_array_r4 *);
44 export_proto(arandom_r4);
45
46 extern void arandom_r8 (gfc_array_r8 *);
47 export_proto(arandom_r8);
48
49 #ifdef HAVE_GFC_REAL_10
50
51 extern void random_r10 (GFC_REAL_10 *);
52 iexport_proto(random_r10);
53
54 extern void arandom_r10 (gfc_array_r10 *);
55 export_proto(arandom_r10);
56
57 #endif
58
59 #ifdef HAVE_GFC_REAL_16
60
61 extern void random_r16 (GFC_REAL_16 *);
62 iexport_proto(random_r16);
63
64 extern void arandom_r16 (gfc_array_r16 *);
65 export_proto(arandom_r16);
66
67 #endif
68
69 #ifdef __GTHREAD_MUTEX_INIT
70 static __gthread_mutex_t random_lock = __GTHREAD_MUTEX_INIT;
71 #else
72 static __gthread_mutex_t random_lock;
73 #endif
74
75 /* Helper routines to map a GFC_UINTEGER_* to the corresponding
76 GFC_REAL_* types in the range of [0,1). If GFC_REAL_*_RADIX are 2
77 or 16, respectively, we mask off the bits that don't fit into the
78 correct GFC_REAL_*, convert to the real type, then multiply by the
79 correct offset.
80 */
81
82
83 static inline void
84 rnumber_4 (GFC_REAL_4 *f, GFC_UINTEGER_4 v)
85 {
86 GFC_UINTEGER_4 mask;
87 #if GFC_REAL_4_RADIX == 2
88 mask = ~ (GFC_UINTEGER_4) 0u << (32 - GFC_REAL_4_DIGITS);
89 #elif GFC_REAL_4_RADIX == 16
90 mask = ~ (GFC_UINTEGER_4) 0u << ((8 - GFC_REAL_4_DIGITS) * 4);
91 #else
92 #error "GFC_REAL_4_RADIX has unknown value"
93 #endif
94 v = v & mask;
95 *f = (GFC_REAL_4) v * (GFC_REAL_4) 0x1.p-32f;
96 }
97
98 static inline void
99 rnumber_8 (GFC_REAL_8 *f, GFC_UINTEGER_8 v)
100 {
101 GFC_UINTEGER_8 mask;
102 #if GFC_REAL_8_RADIX == 2
103 mask = ~ (GFC_UINTEGER_8) 0u << (64 - GFC_REAL_8_DIGITS);
104 #elif GFC_REAL_8_RADIX == 16
105 mask = ~ (GFC_UINTEGER_8) 0u << (16 - GFC_REAL_8_DIGITS) * 4);
106 #else
107 #error "GFC_REAL_8_RADIX has unknown value"
108 #endif
109 v = v & mask;
110 *f = (GFC_REAL_8) v * (GFC_REAL_8) 0x1.p-64;
111 }
112
113 #ifdef HAVE_GFC_REAL_10
114
115 static inline void
116 rnumber_10 (GFC_REAL_10 *f, GFC_UINTEGER_8 v)
117 {
118 GFC_UINTEGER_8 mask;
119 #if GFC_REAL_10_RADIX == 2
120 mask = ~ (GFC_UINTEGER_8) 0u << (64 - GFC_REAL_10_DIGITS);
121 #elif GFC_REAL_10_RADIX == 16
122 mask = ~ (GFC_UINTEGER_10) 0u << ((16 - GFC_REAL_10_DIGITS) * 4);
123 #else
124 #error "GFC_REAL_10_RADIX has unknown value"
125 #endif
126 v = v & mask;
127 *f = (GFC_REAL_10) v * (GFC_REAL_10) 0x1.p-64;
128 }
129 #endif
130
131 #ifdef HAVE_GFC_REAL_16
132
133 /* For REAL(KIND=16), we only need to mask off the lower bits. */
134
135 static inline void
136 rnumber_16 (GFC_REAL_16 *f, GFC_UINTEGER_8 v1, GFC_UINTEGER_8 v2)
137 {
138 GFC_UINTEGER_8 mask;
139 #if GFC_REAL_16_RADIX == 2
140 mask = ~ (GFC_UINTEGER_8) 0u << (128 - GFC_REAL_16_DIGITS);
141 #elif GFC_REAL_16_RADIX == 16
142 mask = ~ (GFC_UINTEGER_8) 0u << ((32 - GFC_REAL_16_DIGITS) * 4);
143 #else
144 #error "GFC_REAL_16_RADIX has unknown value"
145 #endif
146 v2 = v2 & mask;
147 *f = (GFC_REAL_16) v1 * (GFC_REAL_16) 0x1.p-64
148 + (GFC_REAL_16) v2 * (GFC_REAL_16) 0x1.p-128;
149 }
150 #endif
151 /* libgfortran previously had a Mersenne Twister, taken from the paper:
152
153 Mersenne Twister: 623-dimensionally equidistributed
154 uniform pseudorandom generator.
155
156 by Makoto Matsumoto & Takuji Nishimura
157 which appeared in the: ACM Transactions on Modelling and Computer
158 Simulations: Special Issue on Uniform Random Number
159 Generation. ( Early in 1998 ).
160
161 The Mersenne Twister code was replaced due to
162
163 (1) Simple user specified seeds lead to really bad sequences for
164 nearly 100000 random numbers.
165 (2) open(), read(), and close() were not properly declared via header
166 files.
167 (3) The global index i was abused and caused unexpected behavior with
168 GET and PUT.
169 (4) See PR 15619.
170
171
172 libgfortran currently uses George Marsaglia's KISS (Keep It Simple Stupid)
173 random number generator. This PRNG combines:
174
175 (1) The congruential generator x(n)=69069*x(n-1)+1327217885 with a period
176 of 2^32,
177 (2) A 3-shift shift-register generator with a period of 2^32-1,
178 (3) Two 16-bit multiply-with-carry generators with a period of
179 597273182964842497 > 2^59.
180
181 The overall period exceeds 2^123.
182
183 http://www.ciphersbyritter.com/NEWS4/RANDC.HTM#369F6FCA.74C7C041@stat.fsu.edu
184
185 The above web site has an archive of a newsgroup posting from George
186 Marsaglia with the statement:
187
188 Subject: Random numbers for C: Improvements.
189 Date: Fri, 15 Jan 1999 11:41:47 -0500
190 From: George Marsaglia <geo@stat.fsu.edu>
191 Message-ID: <369F6FCA.74C7C041@stat.fsu.edu>
192 References: <369B5E30.65A55FD1@stat.fsu.edu>
193 Newsgroups: sci.stat.math,sci.math,sci.math.numer-analysis
194 Lines: 93
195
196 As I hoped, several suggestions have led to
197 improvements in the code for RNG's I proposed for
198 use in C. (See the thread "Random numbers for C: Some
199 suggestions" in previous postings.) The improved code
200 is listed below.
201
202 A question of copyright has also been raised. Unlike
203 DIEHARD, there is no copyright on the code below. You
204 are free to use it in any way you want, but you may
205 wish to acknowledge the source, as a courtesy.
206
207 "There is no copyright on the code below." included the original
208 KISS algorithm. */
209
210 /* We use three KISS random number generators, with different
211 seeds.
212 As a matter of Quality of Implementation, the random numbers
213 we generate for different REAL kinds, starting from the same
214 seed, are always the same up to the precision of these types.
215 We do this by using three generators with different seeds, the
216 first one always for the most significant bits, the second one
217 for bits 33..64 (if present in the REAL kind), and the third one
218 (called twice) for REAL(16).
219 */
220
221 #define GFC_SL(k, n) ((k)^((k)<<(n)))
222 #define GFC_SR(k, n) ((k)^((k)>>(n)))
223
224 /* Reference for the seed:
225 From: "George Marsaglia" <g...@stat.fsu.edu>
226 Newsgroups: sci.math
227 Message-ID: <e7CcnWxczriWssCjXTWc3A@comcast.com>
228
229 The KISS RNG uses four seeds, x, y, z, c,
230 with 0<=x<2^32, 0<y<2^32, 0<=z<2^32, 0<=c<698769069
231 except that the two pairs
232 z=0,c=0 and z=2^32-1,c=698769068
233 should be avoided.
234 */
235
236 #define KISS_DEFAULT_SEED_1 123456789, 362436069, 521288629, 316191069
237 #define KISS_DEFAULT_SEED_2 987654321, 458629013, 582859209, 438195021
238 #ifdef HAVE_GFC_REAL_16
239 #define KISS_DEFAULT_SEED_3 573658661, 185639104, 582619469, 296736107
240 #endif
241
242 static GFC_UINTEGER_4 kiss_seed[] = {
243 KISS_DEFAULT_SEED_1,
244 KISS_DEFAULT_SEED_2,
245 #ifdef HAVE_GFC_REAL_16
246 KISS_DEFAULT_SEED_3
247 #endif
248 };
249
250 static GFC_UINTEGER_4 kiss_default_seed[] = {
251 KISS_DEFAULT_SEED_1,
252 KISS_DEFAULT_SEED_2,
253 #ifdef HAVE_GFC_REAL_16
254 KISS_DEFAULT_SEED_3
255 #endif
256 };
257
258 static const GFC_INTEGER_4 kiss_size = sizeof(kiss_seed)/sizeof(kiss_seed[0]);
259
260 static GFC_UINTEGER_4 * const kiss_seed_1 = kiss_seed;
261 static GFC_UINTEGER_4 * const kiss_seed_2 = kiss_seed + 4;
262
263 #ifdef HAVE_GFC_REAL_16
264 static GFC_UINTEGER_4 * const kiss_seed_3 = kiss_seed + 8;
265 #endif
266
267 /* kiss_random_kernel() returns an integer value in the range of
268 (0, GFC_UINTEGER_4_HUGE]. The distribution of pseudorandom numbers
269 should be uniform. */
270
271 static GFC_UINTEGER_4
272 kiss_random_kernel(GFC_UINTEGER_4 * seed)
273 {
274 GFC_UINTEGER_4 kiss;
275
276 seed[0] = 69069 * seed[0] + 1327217885;
277 seed[1] = GFC_SL(GFC_SR(GFC_SL(seed[1],13),17),5);
278 seed[2] = 18000 * (seed[2] & 65535) + (seed[2] >> 16);
279 seed[3] = 30903 * (seed[3] & 65535) + (seed[3] >> 16);
280 kiss = seed[0] + seed[1] + (seed[2] << 16) + seed[3];
281
282 return kiss;
283 }
284
285 /* This function produces a REAL(4) value from the uniform distribution
286 with range [0,1). */
287
288 void
289 random_r4 (GFC_REAL_4 *x)
290 {
291 GFC_UINTEGER_4 kiss;
292
293 __gthread_mutex_lock (&random_lock);
294 kiss = kiss_random_kernel (kiss_seed_1);
295 rnumber_4 (x, kiss);
296 __gthread_mutex_unlock (&random_lock);
297 }
298 iexport(random_r4);
299
300 /* This function produces a REAL(8) value from the uniform distribution
301 with range [0,1). */
302
303 void
304 random_r8 (GFC_REAL_8 *x)
305 {
306 GFC_UINTEGER_8 kiss;
307
308 __gthread_mutex_lock (&random_lock);
309 kiss = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32;
310 kiss += kiss_random_kernel (kiss_seed_2);
311 rnumber_8 (x, kiss);
312 __gthread_mutex_unlock (&random_lock);
313 }
314 iexport(random_r8);
315
316 #ifdef HAVE_GFC_REAL_10
317
318 /* This function produces a REAL(10) value from the uniform distribution
319 with range [0,1). */
320
321 void
322 random_r10 (GFC_REAL_10 *x)
323 {
324 GFC_UINTEGER_8 kiss;
325
326 __gthread_mutex_lock (&random_lock);
327 kiss = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32;
328 kiss += kiss_random_kernel (kiss_seed_2);
329 rnumber_10 (x, kiss);
330 __gthread_mutex_unlock (&random_lock);
331 }
332 iexport(random_r10);
333
334 #endif
335
336 /* This function produces a REAL(16) value from the uniform distribution
337 with range [0,1). */
338
339 #ifdef HAVE_GFC_REAL_16
340
341 void
342 random_r16 (GFC_REAL_16 *x)
343 {
344 GFC_UINTEGER_8 kiss1, kiss2;
345
346 __gthread_mutex_lock (&random_lock);
347 kiss1 = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32;
348 kiss1 += kiss_random_kernel (kiss_seed_2);
349
350 kiss2 = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_3)) << 32;
351 kiss2 += kiss_random_kernel (kiss_seed_3);
352
353 rnumber_16 (x, kiss1, kiss2);
354 __gthread_mutex_unlock (&random_lock);
355 }
356 iexport(random_r16);
357
358
359 #endif
360 /* This function fills a REAL(4) array with values from the uniform
361 distribution with range [0,1). */
362
363 void
364 arandom_r4 (gfc_array_r4 *x)
365 {
366 index_type count[GFC_MAX_DIMENSIONS];
367 index_type extent[GFC_MAX_DIMENSIONS];
368 index_type stride[GFC_MAX_DIMENSIONS];
369 index_type stride0;
370 index_type dim;
371 GFC_REAL_4 *dest;
372 GFC_UINTEGER_4 kiss;
373 int n;
374
375 dest = x->data;
376
377 dim = GFC_DESCRIPTOR_RANK (x);
378
379 for (n = 0; n < dim; n++)
380 {
381 count[n] = 0;
382 stride[n] = x->dim[n].stride;
383 extent[n] = x->dim[n].ubound + 1 - x->dim[n].lbound;
384 if (extent[n] <= 0)
385 return;
386 }
387
388 stride0 = stride[0];
389
390 __gthread_mutex_lock (&random_lock);
391
392 while (dest)
393 {
394 /* random_r4 (dest); */
395 kiss = kiss_random_kernel (kiss_seed_1);
396 rnumber_4 (dest, kiss);
397
398 /* Advance to the next element. */
399 dest += stride0;
400 count[0]++;
401 /* Advance to the next source element. */
402 n = 0;
403 while (count[n] == extent[n])
404 {
405 /* When we get to the end of a dimension, reset it and increment
406 the next dimension. */
407 count[n] = 0;
408 /* We could precalculate these products, but this is a less
409 frequently used path so probably not worth it. */
410 dest -= stride[n] * extent[n];
411 n++;
412 if (n == dim)
413 {
414 dest = NULL;
415 break;
416 }
417 else
418 {
419 count[n]++;
420 dest += stride[n];
421 }
422 }
423 }
424 __gthread_mutex_unlock (&random_lock);
425 }
426
427 /* This function fills a REAL(8) array with values from the uniform
428 distribution with range [0,1). */
429
430 void
431 arandom_r8 (gfc_array_r8 *x)
432 {
433 index_type count[GFC_MAX_DIMENSIONS];
434 index_type extent[GFC_MAX_DIMENSIONS];
435 index_type stride[GFC_MAX_DIMENSIONS];
436 index_type stride0;
437 index_type dim;
438 GFC_REAL_8 *dest;
439 GFC_UINTEGER_8 kiss;
440 int n;
441
442 dest = x->data;
443
444 dim = GFC_DESCRIPTOR_RANK (x);
445
446 for (n = 0; n < dim; n++)
447 {
448 count[n] = 0;
449 stride[n] = x->dim[n].stride;
450 extent[n] = x->dim[n].ubound + 1 - x->dim[n].lbound;
451 if (extent[n] <= 0)
452 return;
453 }
454
455 stride0 = stride[0];
456
457 __gthread_mutex_lock (&random_lock);
458
459 while (dest)
460 {
461 /* random_r8 (dest); */
462 kiss = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32;
463 kiss += kiss_random_kernel (kiss_seed_2);
464 rnumber_8 (dest, kiss);
465
466 /* Advance to the next element. */
467 dest += stride0;
468 count[0]++;
469 /* Advance to the next source element. */
470 n = 0;
471 while (count[n] == extent[n])
472 {
473 /* When we get to the end of a dimension, reset it and increment
474 the next dimension. */
475 count[n] = 0;
476 /* We could precalculate these products, but this is a less
477 frequently used path so probably not worth it. */
478 dest -= stride[n] * extent[n];
479 n++;
480 if (n == dim)
481 {
482 dest = NULL;
483 break;
484 }
485 else
486 {
487 count[n]++;
488 dest += stride[n];
489 }
490 }
491 }
492 __gthread_mutex_unlock (&random_lock);
493 }
494
495 #ifdef HAVE_GFC_REAL_10
496
497 /* This function fills a REAL(10) array with values from the uniform
498 distribution with range [0,1). */
499
500 void
501 arandom_r10 (gfc_array_r10 *x)
502 {
503 index_type count[GFC_MAX_DIMENSIONS];
504 index_type extent[GFC_MAX_DIMENSIONS];
505 index_type stride[GFC_MAX_DIMENSIONS];
506 index_type stride0;
507 index_type dim;
508 GFC_REAL_10 *dest;
509 GFC_UINTEGER_8 kiss;
510 int n;
511
512 dest = x->data;
513
514 dim = GFC_DESCRIPTOR_RANK (x);
515
516 for (n = 0; n < dim; n++)
517 {
518 count[n] = 0;
519 stride[n] = x->dim[n].stride;
520 extent[n] = x->dim[n].ubound + 1 - x->dim[n].lbound;
521 if (extent[n] <= 0)
522 return;
523 }
524
525 stride0 = stride[0];
526
527 __gthread_mutex_lock (&random_lock);
528
529 while (dest)
530 {
531 /* random_r10 (dest); */
532 kiss = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32;
533 kiss += kiss_random_kernel (kiss_seed_2);
534 rnumber_10 (dest, kiss);
535
536 /* Advance to the next element. */
537 dest += stride0;
538 count[0]++;
539 /* Advance to the next source element. */
540 n = 0;
541 while (count[n] == extent[n])
542 {
543 /* When we get to the end of a dimension, reset it and increment
544 the next dimension. */
545 count[n] = 0;
546 /* We could precalculate these products, but this is a less
547 frequently used path so probably not worth it. */
548 dest -= stride[n] * extent[n];
549 n++;
550 if (n == dim)
551 {
552 dest = NULL;
553 break;
554 }
555 else
556 {
557 count[n]++;
558 dest += stride[n];
559 }
560 }
561 }
562 __gthread_mutex_unlock (&random_lock);
563 }
564
565 #endif
566
567 #ifdef HAVE_GFC_REAL_16
568
569 /* This function fills a REAL(16) array with values from the uniform
570 distribution with range [0,1). */
571
572 void
573 arandom_r16 (gfc_array_r16 *x)
574 {
575 index_type count[GFC_MAX_DIMENSIONS];
576 index_type extent[GFC_MAX_DIMENSIONS];
577 index_type stride[GFC_MAX_DIMENSIONS];
578 index_type stride0;
579 index_type dim;
580 GFC_REAL_16 *dest;
581 GFC_UINTEGER_8 kiss1, kiss2;
582 int n;
583
584 dest = x->data;
585
586 dim = GFC_DESCRIPTOR_RANK (x);
587
588 for (n = 0; n < dim; n++)
589 {
590 count[n] = 0;
591 stride[n] = x->dim[n].stride;
592 extent[n] = x->dim[n].ubound + 1 - x->dim[n].lbound;
593 if (extent[n] <= 0)
594 return;
595 }
596
597 stride0 = stride[0];
598
599 __gthread_mutex_lock (&random_lock);
600
601 while (dest)
602 {
603 /* random_r16 (dest); */
604 kiss1 = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32;
605 kiss1 += kiss_random_kernel (kiss_seed_2);
606
607 kiss2 = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_3)) << 32;
608 kiss2 += kiss_random_kernel (kiss_seed_3);
609
610 rnumber_16 (dest, kiss1, kiss2);
611
612 /* Advance to the next element. */
613 dest += stride0;
614 count[0]++;
615 /* Advance to the next source element. */
616 n = 0;
617 while (count[n] == extent[n])
618 {
619 /* When we get to the end of a dimension, reset it and increment
620 the next dimension. */
621 count[n] = 0;
622 /* We could precalculate these products, but this is a less
623 frequently used path so probably not worth it. */
624 dest -= stride[n] * extent[n];
625 n++;
626 if (n == dim)
627 {
628 dest = NULL;
629 break;
630 }
631 else
632 {
633 count[n]++;
634 dest += stride[n];
635 }
636 }
637 }
638 __gthread_mutex_unlock (&random_lock);
639 }
640
641 #endif
642
643 /* random_seed is used to seed the PRNG with either a default
644 set of seeds or user specified set of seeds. random_seed
645 must be called with no argument or exactly one argument. */
646
647 void
648 random_seed_i4 (GFC_INTEGER_4 *size, gfc_array_i4 *put, gfc_array_i4 *get)
649 {
650 int i;
651
652 __gthread_mutex_lock (&random_lock);
653
654 /* Check that we only have one argument present. */
655 if ((size ? 1 : 0) + (put ? 1 : 0) + (get ? 1 : 0) > 1)
656 runtime_error ("RANDOM_SEED should have at most one argument present.");
657
658 /* From the standard: "If no argument is present, the processor assigns
659 a processor-dependent value to the seed." */
660 if (size == NULL && put == NULL && get == NULL)
661 for (i = 0; i < kiss_size; i++)
662 kiss_seed[i] = kiss_default_seed[i];
663
664 if (size != NULL)
665 *size = kiss_size;
666
667 if (put != NULL)
668 {
669 /* If the rank of the array is not 1, abort. */
670 if (GFC_DESCRIPTOR_RANK (put) != 1)
671 runtime_error ("Array rank of PUT is not 1.");
672
673 /* If the array is too small, abort. */
674 if (((put->dim[0].ubound + 1 - put->dim[0].lbound)) < kiss_size)
675 runtime_error ("Array size of PUT is too small.");
676
677 /* This code now should do correct strides. */
678 for (i = 0; i < kiss_size; i++)
679 kiss_seed[i] = (GFC_UINTEGER_4) put->data[i * put->dim[0].stride];
680 }
681
682 /* Return the seed to GET data. */
683 if (get != NULL)
684 {
685 /* If the rank of the array is not 1, abort. */
686 if (GFC_DESCRIPTOR_RANK (get) != 1)
687 runtime_error ("Array rank of GET is not 1.");
688
689 /* If the array is too small, abort. */
690 if (((get->dim[0].ubound + 1 - get->dim[0].lbound)) < kiss_size)
691 runtime_error ("Array size of GET is too small.");
692
693 /* This code now should do correct strides. */
694 for (i = 0; i < kiss_size; i++)
695 get->data[i * get->dim[0].stride] = (GFC_INTEGER_4) kiss_seed[i];
696 }
697
698 __gthread_mutex_unlock (&random_lock);
699 }
700 iexport(random_seed_i4);
701
702
703 void
704 random_seed_i8 (GFC_INTEGER_8 *size, gfc_array_i8 *put, gfc_array_i8 *get)
705 {
706 int i;
707
708 __gthread_mutex_lock (&random_lock);
709
710 /* Check that we only have one argument present. */
711 if ((size ? 1 : 0) + (put ? 1 : 0) + (get ? 1 : 0) > 1)
712 runtime_error ("RANDOM_SEED should have at most one argument present.");
713
714 /* From the standard: "If no argument is present, the processor assigns
715 a processor-dependent value to the seed." */
716 if (size == NULL && put == NULL && get == NULL)
717 for (i = 0; i < kiss_size; i++)
718 kiss_seed[i] = kiss_default_seed[i];
719
720 if (size != NULL)
721 *size = kiss_size / 2;
722
723 if (put != NULL)
724 {
725 /* If the rank of the array is not 1, abort. */
726 if (GFC_DESCRIPTOR_RANK (put) != 1)
727 runtime_error ("Array rank of PUT is not 1.");
728
729 /* If the array is too small, abort. */
730 if (((put->dim[0].ubound + 1 - put->dim[0].lbound)) < kiss_size / 2)
731 runtime_error ("Array size of PUT is too small.");
732
733 /* This code now should do correct strides. */
734 for (i = 0; i < kiss_size; i += 2)
735 memcpy (&kiss_seed[i], &(put->data[i * put->dim[0].stride]),
736 sizeof (GFC_UINTEGER_8));
737 }
738
739 /* Return the seed to GET data. */
740 if (get != NULL)
741 {
742 /* If the rank of the array is not 1, abort. */
743 if (GFC_DESCRIPTOR_RANK (get) != 1)
744 runtime_error ("Array rank of GET is not 1.");
745
746 /* If the array is too small, abort. */
747 if (((get->dim[0].ubound + 1 - get->dim[0].lbound)) < kiss_size / 2)
748 runtime_error ("Array size of GET is too small.");
749
750 /* This code now should do correct strides. */
751 for (i = 0; i < kiss_size; i += 2)
752 memcpy (&(get->data[i * get->dim[0].stride]), &kiss_seed[i],
753 sizeof (GFC_UINTEGER_8));
754 }
755
756 __gthread_mutex_unlock (&random_lock);
757 }
758 iexport(random_seed_i8);
759
760
761 #ifndef __GTHREAD_MUTEX_INIT
762 static void __attribute__((constructor))
763 init (void)
764 {
765 __GTHREAD_MUTEX_INIT_FUNCTION (&random_lock);
766 }
767 #endif