From c00289ba543121a78c5ab07a8e45385cc12fb9a8 Mon Sep 17 00:00:00 2001
From: TiborGY <gyori.tibor@stud.u-szeged.hu>
Date: Sat, 1 Jun 2019 21:30:06 +0200
Subject: [PATCH 1/3] upload thread safety test folder

---
 cpp_thread_test/Makefile                   |  14 +++
 cpp_thread_test/cpp_thread_safety_common.h |  55 +++++++++++
 cpp_thread_test/dgemm_thread_safety.cpp    |  92 +++++++++++++++++++
 cpp_thread_test/dgemv_thread_safety.cpp    | 101 +++++++++++++++++++++
 4 files changed, 262 insertions(+)
 create mode 100644 cpp_thread_test/Makefile
 create mode 100644 cpp_thread_test/cpp_thread_safety_common.h
 create mode 100644 cpp_thread_test/dgemm_thread_safety.cpp
 create mode 100644 cpp_thread_test/dgemv_thread_safety.cpp

diff --git a/cpp_thread_test/Makefile b/cpp_thread_test/Makefile
new file mode 100644
index 000000000..81e3470ef
--- /dev/null
+++ b/cpp_thread_test/Makefile
@@ -0,0 +1,14 @@
+include ../Makefile.rule
+
+all :: dgemv_tester dgemm_tester
+
+dgemv_tester :
+	$(CXX) $(COMMON_OPT) -Wall -Wextra -Wshadow -fopenmp -std=c++11 dgemv_thread_safety.cpp ../libopenblas.a -lpthread -o dgemv_tester
+	./dgemv_tester
+
+dgemm_tester : dgemv_tester
+	$(CXX) $(COMMON_OPT) -Wall -Wextra -Wshadow -fopenmp -std=c++11 dgemm_thread_safety.cpp ../libopenblas.a -lpthread -o dgemm_tester
+	./dgemm_tester
+
+clean ::
+	rm -f dgemv_tester dgemm_tester
diff --git a/cpp_thread_test/cpp_thread_safety_common.h b/cpp_thread_test/cpp_thread_safety_common.h
new file mode 100644
index 000000000..60ab5bb2f
--- /dev/null
+++ b/cpp_thread_test/cpp_thread_safety_common.h
@@ -0,0 +1,55 @@
+inline void pauser(){
+    /// a portable way to pause a program
+    std::string dummy;
+    std::cout << "Press enter to continue...";
+    std::getline(std::cin, dummy);
+}
+
+void FillMatrices(std::vector<std::vector<double>>& matBlock, std::mt19937_64& PRNG, std::uniform_real_distribution<double>& rngdist, const blasint randomMatSize, const uint32_t numConcurrentThreads, const uint32_t numMat){
+	for(uint32_t i=0; i<numMat; i++){
+		for(uint32_t j = 0; j < static_cast<uint32_t>(randomMatSize*randomMatSize); j++){
+			matBlock[i][j] = rngdist(PRNG);
+		}
+	}
+	for(uint32_t i=numMat; i<(numConcurrentThreads*numMat); i+=numMat){
+		for(uint32_t j=0; j<numMat; j++){
+			matBlock[i+j] = matBlock[j];
+		}
+	}
+}
+
+void FillVectors(std::vector<std::vector<double>>& vecBlock, std::mt19937_64& PRNG, std::uniform_real_distribution<double>& rngdist, const blasint randomMatSize, const uint32_t numConcurrentThreads, const uint32_t numVec){
+	for(uint32_t i=0; i<numVec; i++){
+		for(uint32_t j = 0; j < static_cast<uint32_t>(randomMatSize); j++){
+			vecBlock[i][j] = rngdist(PRNG);
+		}
+	}
+	for(uint32_t i=numVec; i<(numConcurrentThreads*numVec); i+=numVec){
+		for(uint32_t j=0; j<numVec; j++){
+			vecBlock[i+j] = vecBlock[j];
+		}
+	}
+}
+
+std::mt19937_64 InitPRNG(){
+	std::random_device rd;
+	std::mt19937_64 PRNG(rd()); //seed PRNG using /dev/urandom or similar OS provided RNG
+	std::uniform_real_distribution<double> rngdist{-1.0, 1.0};
+	//make sure the internal state of the PRNG is properly mixed by generating 10M random numbers
+	//PRNGs often have unreliable distribution uniformity and other statistical properties before their internal state is sufficiently mixed
+	for (uint32_t i=0;i<10000000;i++) rngdist(PRNG);
+	return PRNG;
+}
+
+void PrintMatrices(const std::vector<std::vector<double>>& matBlock, const blasint randomMatSize, const uint32_t numConcurrentThreads, const uint32_t numMat){
+	for (uint32_t i=0;i<numConcurrentThreads*numMat;i++){
+		std::cout<<i<<std::endl;
+		for (uint32_t j = 0; j < static_cast<uint32_t>(randomMatSize); j++){
+			for (uint32_t k = 0; k < static_cast<uint32_t>(randomMatSize); k++){
+				std::cout<<matBlock[i][j*randomMatSize + k]<<"  ";
+			}
+			std::cout<<std::endl;
+		}
+		std::cout<<std::endl;
+	}
+}
diff --git a/cpp_thread_test/dgemm_thread_safety.cpp b/cpp_thread_test/dgemm_thread_safety.cpp
new file mode 100644
index 000000000..cecf794fa
--- /dev/null
+++ b/cpp_thread_test/dgemm_thread_safety.cpp
@@ -0,0 +1,92 @@
+#include <iostream>
+#include <vector>
+#include <random>
+#include <future>
+#include <omp.h>
+#include "../cblas.h"
+#include "cpp_thread_safety_common.h"
+
+void launch_cblas_dgemm(double* A, double* B, double* C, const blasint randomMatSize){
+	cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, randomMatSize, randomMatSize, randomMatSize, 1.0, A, randomMatSize, B, randomMatSize, 0.1, C, randomMatSize);
+}
+
+int main(int argc, char* argv[]){
+	blasint randomMatSize = 1024; //dimension of the random square matrices used
+	uint32_t numConcurrentThreads = 52; //number of concurrent calls of the functions being tested
+	uint32_t numTestRounds = 16; //number of testing rounds before success exit
+	
+	if (argc > 4){
+		std::cout<<"ERROR: too many arguments for thread safety tester"<<std::endl;
+		abort();
+	}
+	
+	if(argc == 4){
+		std::vector<std::string> cliArgs;
+		for (int i = 1; i < argc; i++){
+			cliArgs.push_back(argv[i]);
+			std::cout<<argv[i]<<std::endl;
+		}
+		randomMatSize = std::stoul(cliArgs[0]);
+		numConcurrentThreads = std::stoul(cliArgs[1]);
+		numTestRounds = std::stoul(cliArgs[2]);
+	}
+	
+	std::uniform_real_distribution<double> rngdist{-1.0, 1.0};
+	std::vector<std::vector<double>> matBlock(numConcurrentThreads*3);
+	std::vector<std::future<void>> futureBlock(numConcurrentThreads);
+	
+	std::cout<<"*----------------------------*\n";
+	std::cout<<"| DGEMM thread safety tester |\n";
+	std::cout<<"*----------------------------*\n";
+	std::cout<<"Size of random matrices(N=M=K): "<<randomMatSize<<'\n';
+	std::cout<<"Number of concurrent calls into OpenBLAS : "<<numConcurrentThreads<<'\n';
+	std::cout<<"Number of testing rounds : "<<numTestRounds<<'\n';
+	std::cout<<"This test will need "<<(static_cast<uint64_t>(randomMatSize*randomMatSize)*numConcurrentThreads*3*8)/static_cast<double>(1024*1024)<<" MiB of RAM\n"<<std::endl;
+	
+	std::cout<<"Initializing random number generator..."<<std::flush;
+	std::mt19937_64 PRNG = InitPRNG();
+	std::cout<<"done\n";
+	
+	std::cout<<"Preparing to test CBLAS DGEMM thread safety\n";
+	std::cout<<"Allocating matrices..."<<std::flush;
+	for(uint32_t i=0; i<(numConcurrentThreads*3); i++){
+		matBlock[i].resize(randomMatSize*randomMatSize);
+	}
+	std::cout<<"done\n";
+	//pauser();
+	std::cout<<"Filling matrices with random numbers..."<<std::flush;
+	FillMatrices(matBlock, PRNG, rngdist, randomMatSize, numConcurrentThreads, 3);
+	//PrintMatrices(matBlock, randomMatSize, numConcurrentThreads, 3);
+	std::cout<<"done\n";
+	std::cout<<"Testing CBLAS DGEMM thread safety\n";
+	omp_set_num_threads(numConcurrentThreads);
+	for(uint32_t R=0; R<numTestRounds; R++){
+		std::cout<<"DGEMM round #"<<R<<std::endl;
+		std::cout<<"Launching "<<numConcurrentThreads<<" threads simultaneously using OpenMP..."<<std::flush;
+		#pragma omp parallel for default(none) shared(futureBlock, matBlock, randomMatSize, numConcurrentThreads)
+		for(uint32_t i=0; i<numConcurrentThreads; i++){
+			futureBlock[i] = std::async(std::launch::async, launch_cblas_dgemm, &matBlock[i*3][0], &matBlock[i*3+1][0], &matBlock[i*3+2][0], randomMatSize);
+			//launch_cblas_dgemm( &matBlock[i][0], &matBlock[i+1][0], &matBlock[i+2][0]);
+		}
+		std::cout<<"done\n";
+		std::cout<<"Waiting for threads to finish..."<<std::flush;
+		for(uint32_t i=0; i<numConcurrentThreads; i++){
+			futureBlock[i].get();
+		}
+		std::cout<<"done\n";
+		//PrintMatrices(matBlock, randomMatSize, numConcurrentThreads, 3);
+		std::cout<<"Comparing results from different threads..."<<std::flush;
+		for(uint32_t i=3; i<(numConcurrentThreads*3); i+=3){ //i is the index of matrix A, for a given thread
+			for(uint32_t j = 0; j < static_cast<uint32_t>(randomMatSize*randomMatSize); j++){
+				if (std::abs(matBlock[i+2][j] - matBlock[2][j]) > 1.0E-13){ //i+2 is the index of matrix C, for a given thread
+					std::cout<<"ERROR: one of the threads returned a different result! Index : "<<i+2<<std::endl;
+					std::cout<<"CBLAS DGEMM thread safety test FAILED!"<<std::endl;
+					return -1;
+				}
+			}
+		}
+		std::cout<<"OK!\n"<<std::endl;
+	}
+	std::cout<<"CBLAS DGEMM thread safety test PASSED!\n"<<std::endl;
+	return 0;
+}
diff --git a/cpp_thread_test/dgemv_thread_safety.cpp b/cpp_thread_test/dgemv_thread_safety.cpp
new file mode 100644
index 000000000..22505d03f
--- /dev/null
+++ b/cpp_thread_test/dgemv_thread_safety.cpp
@@ -0,0 +1,101 @@
+#include <iostream>
+#include <vector>
+#include <random>
+#include <future>
+#include <omp.h>
+#include "../cblas.h"
+#include "cpp_thread_safety_common.h"
+
+void launch_cblas_dgemv(double* A, double* x, double* y, const blasint randomMatSize){
+	const blasint inc = 1;
+	cblas_dgemv(CblasColMajor, CblasNoTrans, randomMatSize, randomMatSize, 1.0, A, randomMatSize, x, inc, 0.1, y, inc);
+}
+
+int main(int argc, char* argv[]){
+	blasint randomMatSize = 1024; //dimension of the random square matrices and vectors being used
+	uint32_t numConcurrentThreads = 52; //number of concurrent calls of the functions being tested
+	uint32_t numTestRounds = 16; //number of testing rounds before success exit
+	
+	if (argc > 4){
+		std::cout<<"ERROR: too many arguments for thread safety tester"<<std::endl;
+		abort();
+	}
+	if(argc == 4){
+		std::vector<std::string> cliArgs;
+		for (int i = 1; i < argc; i++){
+			cliArgs.push_back(argv[i]);
+			std::cout<<argv[i]<<std::endl;
+		}
+		randomMatSize = std::stoul(cliArgs.at(0));
+		numConcurrentThreads = std::stoul(cliArgs.at(1));
+		numTestRounds = std::stoul(cliArgs.at(2));
+	}
+	
+	std::uniform_real_distribution<double> rngdist{-1.0, 1.0};
+	std::vector<std::vector<double>> matBlock(numConcurrentThreads);
+	std::vector<std::vector<double>> vecBlock(numConcurrentThreads*2);
+	std::vector<std::future<void>> futureBlock(numConcurrentThreads);
+	
+	std::cout<<"*----------------------------*\n";
+	std::cout<<"| DGEMV thread safety tester |\n";
+	std::cout<<"*----------------------------*\n";
+	std::cout<<"Size of random matrices and vectors(N=M): "<<randomMatSize<<'\n';
+	std::cout<<"Number of concurrent calls into OpenBLAS : "<<numConcurrentThreads<<'\n';
+	std::cout<<"Number of testing rounds : "<<numTestRounds<<'\n';
+	std::cout<<"This test will need "<<((static_cast<uint64_t>(randomMatSize*randomMatSize)*numConcurrentThreads*8)+(static_cast<uint64_t>(randomMatSize)*numConcurrentThreads*8*2))/static_cast<double>(1024*1024)<<" MiB of RAM\n"<<std::endl;
+	
+	std::cout<<"Initializing random number generator..."<<std::flush;
+	std::mt19937_64 PRNG = InitPRNG();
+	std::cout<<"done\n";
+	
+	std::cout<<"Preparing to test CBLAS DGEMV thread safety\n";
+	std::cout<<"Allocating matrices..."<<std::flush;
+	for(uint32_t i=0; i<numConcurrentThreads; i++){
+		matBlock.at(i).resize(randomMatSize*randomMatSize);
+	}
+	std::cout<<"done\n";
+	std::cout<<"Allocating vectors..."<<std::flush;
+	for(uint32_t i=0; i<(numConcurrentThreads*2); i++){
+		vecBlock.at(i).resize(randomMatSize);
+	}
+	std::cout<<"done\n";
+	//pauser();
+	
+	std::cout<<"Filling matrices with random numbers..."<<std::flush;
+	FillMatrices(matBlock, PRNG, rngdist, randomMatSize, numConcurrentThreads, 1);
+	//PrintMatrices(matBlock, randomMatSize, numConcurrentThreads);
+	std::cout<<"done\n";
+	std::cout<<"Filling vectors with random numbers..."<<std::flush;
+	FillVectors(vecBlock, PRNG, rngdist, randomMatSize, numConcurrentThreads, 2);
+	std::cout<<"done\n";
+	
+	std::cout<<"Testing CBLAS DGEMV thread safety"<<std::endl;
+	omp_set_num_threads(numConcurrentThreads);
+	for(uint32_t R=0; R<numTestRounds; R++){
+		std::cout<<"DGEMV round #"<<R<<std::endl;
+		std::cout<<"Launching "<<numConcurrentThreads<<" threads simultaneously using OpenMP..."<<std::flush;
+		#pragma omp parallel for default(none) shared(futureBlock, matBlock, vecBlock, randomMatSize, numConcurrentThreads)
+		for(uint32_t i=0; i<numConcurrentThreads; i++){
+			futureBlock[i] = std::async(std::launch::async, launch_cblas_dgemv, &matBlock[i][0], &vecBlock[i*2][0], &vecBlock[i*2+1][0], randomMatSize);
+		}
+		std::cout<<"done\n";
+		std::cout<<"Waiting for threads to finish..."<<std::flush;
+		for(uint32_t i=0; i<numConcurrentThreads; i++){
+			futureBlock[i].get();
+		}
+		std::cout<<"done\n";
+		std::cout<<"Comparing results from different threads..."<<std::flush;
+		for(uint32_t i=2; i<(numConcurrentThreads*2); i+=2){ //i is the index of vector x, for a given thread
+			for(uint32_t j = 0; j < static_cast<uint32_t>(randomMatSize); j++){
+				if (std::abs(vecBlock[i+1][j] - vecBlock[1][j]) > 1.0E-13){ //i+1 is the index of vector y, for a given thread
+					std::cout<<"ERROR: one of the threads returned a different result! Index : "<<i+1<<std::endl;
+					std::cout<<"CBLAS DGEMV thread safety test FAILED!"<<std::endl;
+					return -1;
+				}
+			}
+		}
+		std::cout<<"OK!\n"<<std::endl;
+	}
+	std::cout<<"CBLAS DGEMV thread safety test PASSED!\n"<<std::endl;
+	return 0;
+}

From 1aded698212cca0e6fd447306ae346460cf88616 Mon Sep 17 00:00:00 2001
From: TiborGY <gyori.tibor@stud.u-szeged.hu>
Date: Sat, 1 Jun 2019 21:32:52 +0200
Subject: [PATCH 2/3] hook up c++ thread safety test (main Makefile)

---
 Makefile | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/Makefile b/Makefile
index 21096f893..20ef1e868 100644
--- a/Makefile
+++ b/Makefile
@@ -34,7 +34,7 @@ endif
 
 LAPACK_NOOPT := $(filter-out -O0 -O1 -O2 -O3 -Ofast,$(LAPACK_FFLAGS))
 
-SUBDIRS_ALL = $(SUBDIRS) test ctest utest exports benchmark ../laswp ../bench
+SUBDIRS_ALL = $(SUBDIRS) test ctest utest exports benchmark ../laswp ../bench cpp_thread_test
 
 .PHONY : all libs netlib $(RELA) test ctest shared install
 .NOTPARALLEL : all libs $(RELA) prof lapack-test install blas-test
@@ -127,6 +127,9 @@ ifndef NO_FBLAS
 endif
 ifndef NO_CBLAS
 	$(MAKE) -C ctest all
+ifeq ($(CPP_THREAD_SAFETY_TEST), 1)
+	$(MAKE) -C cpp_thread_test all
+endif
 endif
 endif
 

From 16f3df5d3551ff705d5d23dcdf26853114fb6956 Mon Sep 17 00:00:00 2001
From: TiborGY <gyori.tibor@stud.u-szeged.hu>
Date: Sat, 1 Jun 2019 21:36:41 +0200
Subject: [PATCH 3/3]  add c++ thread test option to Makefile.rule

---
 Makefile.rule | 15 +++++++++++++++
 1 file changed, 15 insertions(+)

diff --git a/Makefile.rule b/Makefile.rule
index 7c128fb49..209934991 100644
--- a/Makefile.rule
+++ b/Makefile.rule
@@ -220,6 +220,21 @@ COMMON_PROF = -pg
 # SYMBOLPREFIX=
 # SYMBOLSUFFIX=
 
+# Run a C++ based thread safety tester after the build is done.
+# This is mostly intended as a developer feature to spot regressions, but users and
+# package maintainers can enable this if they have doubts about the thread safety of
+# the library, given the configuration in this file.
+# By default, the thread safety tester launches 52 concurrent calculations at the same
+# time.
+#
+# Please note that the test uses ~1300 MiB of RAM for the DGEMM test.
+#
+# The test requires CBLAS to be built, a C++11 capable compiler and the presence of
+# an OpenMP implementation. If you are cross-compiling this test will probably not
+# work at all.
+#
+# CPP_THREAD_SAFETY_TEST = 1
+
 #
 #  End of user configuration
 #