From b774669268341acc22f352f845b6bc34f0f6d418 Mon Sep 17 00:00:00 2001
From: Alberto Ruiz <aruiz@um.es>
Date: Mon, 1 Dec 2014 13:35:54 +0100
Subject: change rand() to drand48_r

(I have tried also random_r but got segfaults on exit) (?)
---
 packages/base/THANKS.md          |  6 +++++-
 packages/base/src/C/vector-aux.c | 22 +++++++++++++++-------
 2 files changed, 20 insertions(+), 8 deletions(-)

(limited to 'packages')

diff --git a/packages/base/THANKS.md b/packages/base/THANKS.md
index fcbdb8e..90e6199 100644
--- a/packages/base/THANKS.md
+++ b/packages/base/THANKS.md
@@ -171,5 +171,9 @@ module reorganization, monadic mapVectorM, and many other improvements.
 
 - Niklas Hambüchen improved the documentation.
 
-- "erdeszt" optimized "conv" using a direct vector reverse
+- "erdeszt" optimized "conv" using a direct vector reverse.
+
+- John Shahbazian added support for openBLAS.
+
+- "yongqli" reported the bug in randomVector (rand() is not thread-safe).
 
diff --git a/packages/base/src/C/vector-aux.c b/packages/base/src/C/vector-aux.c
index a7eaa08..d5b6d4d 100644
--- a/packages/base/src/C/vector-aux.c
+++ b/packages/base/src/C/vector-aux.c
@@ -700,8 +700,15 @@ int saveMatrix(char * file, char * format, KDMAT(a)){
 
 ////////////////////////////////////////////////////////////////////////////////
 
+inline double urandom(struct drand48_data * buffer) {
+    double res;
+    drand48_r(buffer,&res);
+    return res;
+}
+
+
 // http://c-faq.com/lib/gaussian.html
-double gaussrand()
+double gaussrand(struct drand48_data * buffer)
 {
 	static double V1, V2, S;
 	static int phase = 0;
@@ -709,8 +716,8 @@ double gaussrand()
 
 	if(phase == 0) {
 		do {
-			double U1 = (double)rand() / RAND_MAX;
-			double U2 = (double)rand() / RAND_MAX;
+			double U1 = urandom(buffer);
+			double U2 = urandom(buffer);
 
 			V1 = 2 * U1 - 1;
 			V2 = 2 * U2 - 1;
@@ -726,19 +733,20 @@ double gaussrand()
 	return X;
 }
 
-int random_vector(int seed, int code, DVEC(r)) {
-    srand(seed);
+int random_vector(unsigned int seed, int code, DVEC(r)) {
+    struct drand48_data buffer;
+    srand48_r(seed,&buffer);
     int k;
     switch (code) {
       case 0: { // uniform
         for (k=0; k<rn; k++) {
-            rp[k] = (double)rand()/RAND_MAX;
+            rp[k] = urandom(&buffer);
         }
         OK
       }
       case 1: { // gaussian
         for (k=0; k<rn; k++) {
-            rp[k] = gaussrand();
+            rp[k] = gaussrand(&buffer);
         }
         OK
       }
-- 
cgit v1.2.3