axpby as concrete method in DataContainer and BlockDataContainer (#489)

* axpby as concrete method in DataContainer and BlockDataContainer

* fixed axpby and added unittest

* PDHG to use axpby

* pass num_threads to axpby

* void commit

* add seed to random in test

* NUM_THREADS can be imported from ccpi.utilities

* added test to axpby with num_threads

7 files changed, 30 insertions, 103 deletions

diff --git a/src/Core/CMakeLists.txt b/src/Core/CMakeLists.txt
index e828fe5..9c9a89d 100644
--- a/src/Core/CMakeLists.txt
+++ b/src/Core/CMakeLists.txt
@@ -93,8 +93,9 @@ message("CMAKE_STATIC_LINKER_FLAGS ${CMAKE_STATIC_LINKER_FLAGS}")
 
 
 
-add_library(cilacc SHARED ${CMAKE_CURRENT_SOURCE_DIR}/axpby.c
-                          ${CMAKE_CURRENT_SOURCE_DIR}/FiniteDifferenceLibrary.c )
+add_library(cilacc SHARED ${CMAKE_CURRENT_SOURCE_DIR}/utilities.c
+                          ${CMAKE_CURRENT_SOURCE_DIR}/FiniteDifferenceLibrary.c
+                          ${CMAKE_CURRENT_SOURCE_DIR}/axpby.c )
 
 target_link_libraries(cilacc ${OpenMP_C_LIB_NAMES} )
 include_directories(cilacc PUBLIC 
diff --git a/src/Core/FiniteDifferenceLibrary.c b/src/Core/FiniteDifferenceLibrary.c
index fbf2646..244e170 100644
--- a/src/Core/FiniteDifferenceLibrary.c
+++ b/src/Core/FiniteDifferenceLibrary.c
@@ -16,17 +16,6 @@ DLL_EXPORT int openMPtest(int nThreads)
 	}
 	return nThreads_running;
 }
-void threads_setup(int nThreads_requested, int *nThreads_current)
-{
-#pragma omp parallel
-	{
-		if (omp_get_thread_num() == 0)
-		{
-			*nThreads_current = omp_get_num_threads();
-		}
-	}
-	omp_set_num_threads(nThreads_requested);
-}
 
 int fdiff_direct_neumann(const float *inimagefull, float *outimageXfull, float *outimageYfull, float *outimageZfull, float *outimageCfull, long nx, long ny, long nz, long nc)
 {
diff --git a/src/Core/axpby.c b/src/Core/axpby.c
index c4d162d..54a597f 100755
--- a/src/Core/axpby.c
+++ b/src/Core/axpby.c
@@ -1,87 +1,12 @@
 #include "axpby.h"
 
 
-DLL_EXPORT int padd(float * x, float * y, float * out, long size){
+DLL_EXPORT int saxpby(float * x, float * y, float * out, float a, float b, long size, int nThreads){
     long i = 0;
-#pragma omp parallel for
-    for (i=0; i < size; i++)
-    {
-        *(out + i ) = *(x + i) + *(y+i); 
-    }
-    return 0;
-}
-
-DLL_EXPORT int psubtract(float * x, float * y, float * out, long size){
-    long i = 0;
-#pragma omp parallel
-{
-//#pragma omp single
-//{
-//		printf("current number of threads %d\n", omp_get_num_threads());
-//}
-#pragma omp for
-    for (i=0; i < size; i++)
-    {
-        *(out + i ) = *(x + i) - *(y+i); 
-    }
-}
-    return 0;
-    
-}
-
-DLL_EXPORT int pmultiply(float * x, float * y, float * out, long size){
-    long i = 0;
-#pragma omp parallel for
-    for (i=0; i < size; i++)
-    {
-        *(out + i ) = *(x + i) * *(y+i); 
-    }
-    return 0;
-}
-
-DLL_EXPORT int pdivide(float * x, float * y, float * out, long size, float default_value)
-{
-    long i = 0;
-#pragma omp parallel for
-    for (i=0; i < size; i++)
-    {
-        *(out + i ) = *(y+i) ? *(x + i) / *(y+i) : default_value;
-    }
-    return 0;
-}
-DLL_EXPORT int ppower(float * x, float * y, float * out, long size){
-    long i = 0;
-#pragma omp parallel for
-    for (i=0; i < size; i++)
-    {
-        *(out + i ) = (float)pow(*(x + i) , *(y+i)) ; 
-    }
-    return 0;
-}
-
-DLL_EXPORT int pminimum(float * x, float * y, float * out, long size){
-    long i = 0;
-#pragma omp parallel for
-    for (i=0; i < size; i++)
-    {
-        *(out + i ) = *(y+i) > (*x+i) ? *(x + i) :  *(y+i); 
-    }
-    return 0;
-}
-
-DLL_EXPORT int pmaximum(float * x, float * y, float * out, long size) {
-	long i = 0;
-#pragma omp parallel for
-	for (i = 0; i < size; i++)
-	{
-		*(out + i) = *(y + i) < (*x + i) ? *(x + i) : *(y + i);
-	}
-	return 0;
-}
 
+    int nThreads_initial;
+	threads_setup(nThreads, &nThreads_initial);
 
-DLL_EXPORT int saxpby(float * x, float * y, float * out, float a, float b, long size){
-    long i = 0;
 #pragma omp parallel
 {
 #pragma omp for
@@ -90,11 +15,12 @@ DLL_EXPORT int saxpby(float * x, float * y, float * out, float a, float b, long
         *(out + i ) = a * ( *(x + i) ) + b * ( *(y + i) ); 
     }
 }
+    omp_set_num_threads(nThreads_initial);
     return 0;
     
 }
 
-DLL_EXPORT int daxpby(double * x, double * y, double * out, double a, double b, long size) {
+DLL_EXPORT int daxpby(double * x, double * y, double * out, double a, double b, long size, int nThreads) {
 	long i = 0;
 #pragma omp parallel
 	{
diff --git a/src/Core/include/FiniteDifferenceLibrary.h b/src/Core/include/FiniteDifferenceLibrary.h
index 6e426af..b8e6c4f 100644
--- a/src/Core/include/FiniteDifferenceLibrary.h
+++ b/src/Core/include/FiniteDifferenceLibrary.h
@@ -3,4 +3,5 @@
 #include <stdio.h>
 #include "omp.h"
 //#include "ipp.h"
-#include "dll_export.h"
\ No newline at end of file
+#include "dll_export.h"
+#include "utilities.h"
\ No newline at end of file
diff --git a/src/Core/include/axpby.h b/src/Core/include/axpby.h
index 2849547..e13d6e1 100644
--- a/src/Core/include/axpby.h
+++ b/src/Core/include/axpby.h
@@ -3,15 +3,8 @@
 #include <stdio.h>
 #include "omp.h"
 #include "dll_export.h"
+#include "utilities.h"
 
 
-DLL_EXPORT int padd(float * x, float * y, float * out, long size);
-DLL_EXPORT int psubtract(float * x, float * y, float * out, long size);
-DLL_EXPORT int pmultiply(float * x, float * y, float * out, long size);
-DLL_EXPORT int pdivide(float * x, float * y, float * out, long size, float default_value);
-DLL_EXPORT int ppower(float * x, float * y, float * out, long size);
-DLL_EXPORT int pminimum(float * x, float * y, float * out, long size);
-DLL_EXPORT int pmaximum(float * x, float * y, float * out, long size);
-
-DLL_EXPORT int saxpby(float * x, float * y, float * out, float a, float b, long size);
-DLL_EXPORT int daxpby(double * x, double * y, double * out, double a, double b, long size);
+DLL_EXPORT int saxpby(float * x, float * y, float * out, float a, float b, long size, int nThreads);
+DLL_EXPORT int daxpby(double * x, double * y, double * out, double a, double b, long size, int nThreads);
diff --git a/src/Core/include/utilities.h b/src/Core/include/utilities.h
new file mode 100644
index 0000000..c3003d6
--- /dev/null
+++ b/src/Core/include/utilities.h
@@ -0,0 +1,3 @@
+#include "omp.h"
+
+void threads_setup(int nThreads_requested, int *nThreads_current);
\ No newline at end of file
diff --git a/src/Core/utilities.c b/src/Core/utilities.c
new file mode 100644
index 0000000..86b23e8
--- /dev/null
+++ b/src/Core/utilities.c
@@ -0,0 +1,14 @@
+#include "utilities.h"
+
+
+void threads_setup(int nThreads_requested, int *nThreads_current)
+{
+#pragma omp parallel
+	{
+		if (omp_get_thread_num() == 0)
+		{
+			*nThreads_current = omp_get_num_threads();
+		}
+	}
+	omp_set_num_threads(nThreads_requested);
+}