diff --git a/Base/CMakeLists.txt b/Base/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..2acb5e21be9845f6427db8933eb2a281c50370cf
--- /dev/null
+++ b/Base/CMakeLists.txt
@@ -0,0 +1,18 @@
+# Source files
+SET( HEADERS
+ Configure.h.in
+ TimeMeasurementBase.h
+ RunTimeTimer.h
+ CPUTimeTimer.h
+)
+
+SET( SOURCES
+ RunTimeTimer.cpp
+ CPUTimeTimer.cpp
+)
+
+# Create the library
+ADD_LIBRARY( Base STATIC ${HEADERS} ${SOURCES} )
+
+SET_TARGET_PROPERTIES(Base PROPERTIES LINKER_LANGUAGE CXX)
+
diff --git a/Base/CPUTimeTimer.cpp b/Base/CPUTimeTimer.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ad8377cd9348611c414574d9cf914b08ece7a7ff
--- /dev/null
+++ b/Base/CPUTimeTimer.cpp
@@ -0,0 +1,17 @@
+#include "CPUTimeTimer.h"
+
+void CPUTimeTimer::start()
+{
+ _start = clock();
+}
+
+void CPUTimeTimer::stop()
+{
+ _stop = clock();
+}
+
+double CPUTimeTimer::elapsed()
+{
+ return (_stop-_start)/(double)(CLOCKS_PER_SEC);
+}
+
diff --git a/Base/CPUTimeTimer.h b/Base/CPUTimeTimer.h
new file mode 100644
index 0000000000000000000000000000000000000000..8b5eb18b6022e395f7b5beda3dadd7472218f734
--- /dev/null
+++ b/Base/CPUTimeTimer.h
@@ -0,0 +1,21 @@
+#ifndef CPUTIMETIMER_H
+#define CPUTIMETIMER_H
+
+#include <ctime>
+
+#include "TimeMeasurementBase.h"
+
+class CPUTimeTimer
+{
+public:
+ virtual void start();
+ virtual void stop();
+ virtual double elapsed();
+ ~CPUTimeTimer() {};
+private:
+ clock_t _start;
+ clock_t _stop;
+};
+
+#endif
+
diff --git a/Base/Configure.h.in b/Base/Configure.h.in
new file mode 100644
index 0000000000000000000000000000000000000000..aead27c31e567cdc1e141fa27ebce45f7e77f5a5
--- /dev/null
+++ b/Base/Configure.h.in
@@ -0,0 +1,15 @@
+/**
+ * \file Configure.h.in
+ *
+ * #defines which gets set through CMake
+ */
+#ifndef CONFIGURE_H
+#define CONFIGURE_H
+
+#define SOURCEPATH "${CMAKE_SOURCE_DIR}"
+
+#cmakedefine HAVE_PTHREADS
+#cmakedefine PROCESSOR_COUNT "${PROCESSOR_COUNT}"
+
+#endif // CONFIGURE_H
+
diff --git a/Base/RunTimeTimer.cpp b/Base/RunTimeTimer.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b89563245665c5178af08395f02c339cf2890ed4
--- /dev/null
+++ b/Base/RunTimeTimer.cpp
@@ -0,0 +1,17 @@
+#include "RunTimeTimer.h"
+
+void RunTimeTimer::start()
+{
+ gettimeofday(&_start, 0);
+}
+
+void RunTimeTimer::stop()
+{
+ gettimeofday(&_stop, 0);
+}
+
+double RunTimeTimer::elapsed()
+{
+ return (_stop.tv_sec + _stop.tv_usec/1000000.0 - (_start.tv_sec + _start.tv_usec/1000000.0));
+}
+
diff --git a/Base/RunTimeTimer.h b/Base/RunTimeTimer.h
new file mode 100644
index 0000000000000000000000000000000000000000..871cfb5a601a0212ab67e77dc7846dd03c4672a1
--- /dev/null
+++ b/Base/RunTimeTimer.h
@@ -0,0 +1,21 @@
+#ifndef RUNTIMETIMER_H
+#define RUNTIMETIMER_H
+
+#include <sys/time.h>
+
+#include "TimeMeasurementBase.h"
+
+class RunTimeTimer
+{
+public:
+ virtual void start();
+ virtual void stop();
+ virtual double elapsed();
+ ~RunTimeTimer() {};
+private:
+ timeval _start;
+ timeval _stop;
+};
+
+#endif
+
diff --git a/Base/TimeMeasurementBase.h b/Base/TimeMeasurementBase.h
new file mode 100644
index 0000000000000000000000000000000000000000..36c081fe3db1afe184da77f95482b40ba6210d73
--- /dev/null
+++ b/Base/TimeMeasurementBase.h
@@ -0,0 +1,14 @@
+#ifndef TIMEMEASUREMENT_H
+#define TIMEMEASUREMENT_H
+
+class TimeMeasurementBase
+{
+public:
+ virtual void start () = 0;
+ virtual void stop () = 0;
+ virtual double elapsed () = 0;
+ virtual ~TimeMeasurementBase () {};
+};
+
+#endif
+
diff --git a/CMakeLists.txt b/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..56d6bf2f70ae4cc4994bf0e6d9276105abfbfdd4
--- /dev/null
+++ b/CMakeLists.txt
@@ -0,0 +1,33 @@
+# Specify minimum CMake version
+cmake_minimum_required(VERSION 2.6)
+
+# Project name
+project( OGS-6 )
+
+# Set cmake module path
+#SET(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/CMakeConfiguration")
+
+### For GNU C/CXX
+IF(CMAKE_COMPILER_IS_GNUCXX OR CMAKE_COMPILER_IS_GNUCC)
+ IF( NOT CMAKE_BUILD_TYPE STREQUAL "Debug" )
+ MESSAGE(STATUS "Set GCC release flags")
+ SET(CMAKE_CXX_FLAGS "-O3 -march=native -mtune=native -msse4.2 -DNDEBUG")
+ ENDIF()
+ # -g
+ SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-deprecated -Wall -Wextra -fno-nonansi-builtins")
+ SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fopenmp")
+ ADD_DEFINITIONS(
+ -DGCC
+ )
+ENDIF(CMAKE_COMPILER_IS_GNUCXX OR CMAKE_COMPILER_IS_GNUCC)
+
+# Set build directories
+SET( EXECUTABLE_OUTPUT_PATH ${PROJECT_BINARY_DIR}/bin )
+SET( LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/lib )
+
+# Add subdirectories with the projects
+ADD_SUBDIRECTORY( Base )
+ADD_SUBDIRECTORY( MathLib )
+ADD_SUBDIRECTORY( SimpleTests/MatrixTests )
+
+
diff --git a/MathLib/CMakeLists.txt b/MathLib/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..07597457dc76a27b94258c6492ece9bcf88e750a
--- /dev/null
+++ b/MathLib/CMakeLists.txt
@@ -0,0 +1,25 @@
+# Source files
+SET( HEADERS
+ amuxCRS.h
+ CRSMatrix.h
+ CRSMatrixPThreads.h
+ CRSMatrixOpenMP.h
+ CRSSymMatrix.h
+ sparse.h
+ SparseMatrixBase.h
+)
+
+SET( SOURCES
+ amuxCRS.cpp
+)
+
+INCLUDE_DIRECTORIES (
+ .
+ ${PROJECT_BINARY_DIR}/Base
+)
+
+# Create the library
+ADD_LIBRARY( MathLib STATIC ${HEADERS} ${SOURCES} )
+
+SET_TARGET_PROPERTIES(MathLib PROPERTIES LINKER_LANGUAGE CXX)
+
diff --git a/MathLib/CRSMatrix.h b/MathLib/CRSMatrix.h
new file mode 100644
index 0000000000000000000000000000000000000000..f040e654276b375e4fb2c61e2ada60845275e180
--- /dev/null
+++ b/MathLib/CRSMatrix.h
@@ -0,0 +1,55 @@
+#ifndef CRSMATRIX_H
+#define CRSMATRIX_H
+
+#include <string>
+#include <fstream>
+#include <iostream>
+#include "SparseMatrixBase.h"
+#include "sparse.h"
+#include "amuxCRS.h"
+
+template<class T> class CRSMatrix : public SparseMatrixBase<T>
+{
+public:
+ CRSMatrix(std::string const &fname) :
+ SparseMatrixBase<T>(),
+ _row_ptr(NULL), _col_idx(NULL), _data(NULL)
+ {
+ std::ifstream in(fname.c_str(), std::ios::in | std::ios::binary);
+ if (in) {
+ CS_read(in, SparseMatrixBase<T>::_n_rows, _row_ptr, _col_idx, _data);
+ SparseMatrixBase<T>::_n_cols = SparseMatrixBase<T>::_n_rows;
+ in.close();
+ } else {
+ std::cout << "cannot open " << fname << std::endl;
+ }
+ }
+
+ CRSMatrix(unsigned n, unsigned *iA, unsigned *jA, T* A) :
+ SparseMatrixBase<T>(n,n), _row_ptr(iA), _col_idx(jA), _data(A)
+ {}
+
+ CRSMatrix(unsigned n1) :
+ SparseMatrixBase<T>(n1, n1), _row_ptr(NULL), _col_idx(NULL), _data(NULL)
+ {}
+
+ virtual ~CRSMatrix()
+ {
+ delete [] _row_ptr;
+ delete [] _col_idx;
+ delete [] _data;
+ }
+
+ virtual void amux(T d, T const * const x, T *y) const
+ {
+ amuxCRS(d, SparseMatrixBase<T>::_n_rows, _row_ptr, _col_idx, _data, x, y);
+ }
+
+protected:
+ unsigned *_row_ptr;
+ unsigned *_col_idx;
+ T* _data;
+};
+
+#endif
+
diff --git a/MathLib/CRSMatrixDiagPrecond.h b/MathLib/CRSMatrixDiagPrecond.h
new file mode 100644
index 0000000000000000000000000000000000000000..2335eedf001f7cf4b6af353402fc9a9fd2fc2bf8
--- /dev/null
+++ b/MathLib/CRSMatrixDiagPrecond.h
@@ -0,0 +1,40 @@
+#ifndef CRSMATRIXDIAGPRECOND_H
+#define CRSMATRIXDIAGPRECOND_H
+
+#include <omp.h>
+
+#include "CRSMatrix.h"
+#include "sparse.h"
+
+class CRSMatrixDiagPrecond : public CRSMatrix<double>
+{
+public:
+ CRSMatrixDiagPrecond (std::string const &fname, unsigned num_of_threads=1)
+ : CRSMatrix<double>(fname, num_of_threads), _inv_diag(NULL)
+ {
+ _inv_diag = new double[_n_rows];
+ if (!generateDiagPrecond (_n_rows, _row_ptr, _col_idx, _data, _inv_diag)) {
+ std::cout << "Could not create diagonal preconditioner" << std::endl;
+ }
+ }
+
+ void precondApply(double* x) const {
+ omp_set_num_threads( _num_of_threads );
+ {
+ unsigned k;
+ #pragma omp parallel for
+ for (k=0; k<_n_rows; ++k) {
+ x[k] = _inv_diag[k]*x[k];
+ }
+ }
+ }
+
+ ~CRSMatrixDiagPrecond() {
+ delete [] _inv_diag;
+ }
+private:
+ double *_inv_diag;
+};
+
+#endif
+
diff --git a/MathLib/CRSMatrixOpenMP.h b/MathLib/CRSMatrixOpenMP.h
new file mode 100644
index 0000000000000000000000000000000000000000..009be74ee2991377c9cc3e9cf70158d906160c51
--- /dev/null
+++ b/MathLib/CRSMatrixOpenMP.h
@@ -0,0 +1,42 @@
+/*
+ * CRSMatrixOpenMP.h
+ *
+ * Created on: Aug 8, 2011
+ * Author: TF
+ */
+
+#ifndef CRSMATRIXOPENMP_H_
+#define CRSMATRIXOPENMP_H_
+
+#include <string>
+
+#include "CRSMatrix.h"
+#include "amuxCRS.h"
+
+template<class T> class CRSMatrixOpenMP : public CRSMatrix<T> {
+public:
+ CRSMatrixOpenMP(std::string const &fname, unsigned num_of_threads) :
+ CRSMatrix<T>(fname), _num_of_threads (num_of_threads)
+ {}
+
+ CRSMatrixOpenMP(unsigned n, unsigned *iA, unsigned *jA, T* A, unsigned num_of_threads) :
+ CRSMatrix<T>(n, iA, jA, A), _num_of_threads (num_of_threads)
+ {}
+
+ CRSMatrixOpenMP(unsigned n1) :
+ CRSMatrix<T>(n1), _num_of_threads (1)
+ {}
+
+ virtual ~CRSMatrixOpenMP()
+ {}
+
+ virtual void amux(T d, T const * const x, T *y) const
+ {
+ amuxCRSParallelOpenMP(d, SparseMatrixBase<T>::_n_rows, CRSMatrix<T>::_row_ptr, CRSMatrix<T>::_col_idx, CRSMatrix<T>::_data, x, y, _num_of_threads);
+ }
+
+private:
+ unsigned _num_of_threads;
+};
+
+#endif /* CRSMATRIXOPENMP_H_ */
diff --git a/MathLib/CRSMatrixPThreads.h b/MathLib/CRSMatrixPThreads.h
new file mode 100644
index 0000000000000000000000000000000000000000..7c56d9b89488bd1ca1516acc08889cbcf2f67f16
--- /dev/null
+++ b/MathLib/CRSMatrixPThreads.h
@@ -0,0 +1,46 @@
+/*
+ * CRSMatrixPThreads.h
+ *
+ * Created on: Aug 2, 2011
+ * Author: TF
+ */
+
+#ifndef CRSMATRIXPTHREADS_H
+#define CRSMATRIXPTHREADS_H
+
+#include <string>
+
+#include "SparseMatrixBase.h"
+#include "sparse.h"
+#include "CRSMatrix.h"
+#include "amuxCRS.h"
+
+template<class T> class CRSMatrixPThreads : public CRSMatrix<T>
+{
+public:
+ CRSMatrixPThreads(std::string const &fname, unsigned num_of_threads) :
+ CRSMatrix<T>(fname), _num_of_threads (num_of_threads)
+ {}
+
+ CRSMatrixPThreads(unsigned n, unsigned *iA, unsigned *jA, T* A, unsigned num_of_threads) :
+ CRSMatrix<T>(n, iA, jA, A), _num_of_threads (num_of_threads)
+ {}
+
+ CRSMatrixPThreads(unsigned n1) :
+ CRSMatrix<T>(n1), _num_of_threads (1)
+ {}
+
+ virtual ~CRSMatrixPThreads()
+ {}
+
+ virtual void amux(T d, T const * const x, T *y) const
+ {
+ amuxCRSParallelPThreads(d, SparseMatrixBase<T>::_n_rows, CRSMatrix<T>::_row_ptr, CRSMatrix<T>::_col_idx, CRSMatrix<T>::_data, x, y, _num_of_threads);
+ }
+
+protected:
+ unsigned _num_of_threads;
+};
+
+#endif
+
diff --git a/MathLib/CRSSymMatrix.h b/MathLib/CRSSymMatrix.h
new file mode 100644
index 0000000000000000000000000000000000000000..37a75ed86f7966a1d3d543a81bf0ba6ddafeb6c0
--- /dev/null
+++ b/MathLib/CRSSymMatrix.h
@@ -0,0 +1,64 @@
+/*
+ * CRSSymMatrix.h
+ *
+ * Created on: Jul 22, 2011
+ * Author: fischeth
+ */
+
+#ifndef CRSSYMMATRIX_H_
+#define CRSSYMMATRIX_H_
+
+#include "CRSMatrix.h"
+
+template<class T> class CRSSymMatrix : public CRSMatrix<T>
+{
+public:
+ CRSSymMatrix(std::string const &fname)
+ : CRSMatrix<T> (fname)
+ {
+ unsigned nnz (0);
+
+ // count number of non-zeros in the upper triangular part
+ for (unsigned i = 0; i < SparseMatrixBase<T>::_n_rows; i++) {
+ unsigned idx = CRSMatrix<T>::_row_ptr[i+1];
+ for (unsigned j = CRSMatrix<T>::_row_ptr[i]; j < idx; j++)
+ if (CRSMatrix<T>::_col_idx[j] >= i)
+ ++nnz;
+ }
+
+ double *A_new (new double[nnz]);
+ unsigned *jA_new (new unsigned[nnz]);
+ unsigned *iA_new (new unsigned[SparseMatrixBase<T>::_n_rows+1]);
+
+ iA_new[0] = nnz = 0;
+
+ for (unsigned i = 0; i < SparseMatrixBase<T>::_n_rows; i++) {
+ const unsigned idx (CRSMatrix<T>::_row_ptr[i+1]);
+ for (unsigned j = CRSMatrix<T>::_row_ptr[i]; j < idx; j++) {
+ if (CRSMatrix<T>::_col_idx[j] >= i) {
+ A_new[nnz] = CRSMatrix<T>::_data[j];
+ jA_new[nnz++] = CRSMatrix<T>::_col_idx[j];
+ }
+ }
+ iA_new[i+1] = nnz;
+ }
+
+ std::swap(CRSMatrix<T>::_row_ptr, iA_new);
+ std::swap(CRSMatrix<T>::_col_idx, jA_new);
+ std::swap(CRSMatrix<T>::_data, A_new);
+
+ delete[] iA_new;
+ delete[] jA_new;
+ delete[] A_new;
+ }
+
+ virtual ~CRSSymMatrix() {}
+
+ void amux(T d, T const * const x, T *y) const
+ {
+ amuxCRSSym (d, SparseMatrixBase<T>::_n_rows, CRSMatrix<T>::_row_ptr, CRSMatrix<T>::_col_idx, CRSMatrix<T>::_data, x, y);
+ }
+
+};
+
+#endif /* CRSSYMMATRIX_H_ */
diff --git a/MathLib/SparseMatrixBase.h b/MathLib/SparseMatrixBase.h
new file mode 100644
index 0000000000000000000000000000000000000000..831818a1b13e7231d32183101612e0e7912ee4cf
--- /dev/null
+++ b/MathLib/SparseMatrixBase.h
@@ -0,0 +1,21 @@
+#ifndef SPARSEMATRIXBASE_H
+#define SPARSEMATRIXBASE_H
+
+template<class T> class SparseMatrixBase
+{
+public:
+ SparseMatrixBase(unsigned n1, unsigned n2) : _n_rows(n1), _n_cols(n2) { }
+ SparseMatrixBase() : _n_rows(0), _n_cols(0) { }
+ virtual void amux(T d, T const * const x, T *y) const = 0; // y +=d*Ax
+ virtual ~SparseMatrixBase() { }
+ unsigned getNRows () const { return _n_rows; }
+ unsigned getNCols () const { return _n_cols; }
+
+protected:
+ unsigned _n_rows;
+ unsigned _n_cols;
+};
+
+
+#endif
+
diff --git a/MathLib/amuxCRS.cpp b/MathLib/amuxCRS.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..336527d0bd6d12a44ddac62f1c25f12ff6cee420
--- /dev/null
+++ b/MathLib/amuxCRS.cpp
@@ -0,0 +1,148 @@
+#include "amuxCRS.h"
+#include <omp.h>
+#include <omp.h>
+#include <pthread.h>
+
+void amuxCRS (double a,
+ unsigned n, unsigned const * const iA, unsigned const * const jA,
+ double const * const A, double const * const x, double* y)
+{
+ for (unsigned i(0); i<n; i++) {
+ y[i] = 0.0;
+ const unsigned end (iA[i+1]);
+ for (unsigned j(iA[i]); j<end; j++) {
+ y[i] += A[j] * x[jA[j]];
+ }
+ y[i] *= a;
+ }
+}
+
+struct MatMultThreadParam {
+ MatMultThreadParam (double scalar_factor, unsigned beg_row, unsigned end_row,
+ unsigned const * const iA, unsigned const * const jA,
+ double const * const A, double const * const x, double* y) :
+ _a (scalar_factor), _beg_row(beg_row), _end_row(end_row),
+ _row_ptr(iA), _col_idx(jA), _data(A), _x(x), _y(y)
+ {}
+
+ double _a;
+ unsigned _beg_row;
+ unsigned _end_row;
+ unsigned const * const _row_ptr;
+ unsigned const * const _col_idx;
+ double const * const _data;
+ double const * const _x;
+ double * _y;
+};
+
+extern "C" {
+void* amuxCRSpthread (void* ptr)
+{
+ MatMultThreadParam *thread_param ((MatMultThreadParam*)(ptr));
+ const double a (thread_param->_a);
+ const unsigned beg_row (thread_param->_beg_row);
+ const unsigned end_row (thread_param->_end_row);
+ unsigned const * const iA (thread_param->_row_ptr);
+ unsigned const * const jA (thread_param->_col_idx);
+ double const * const A (thread_param->_data);
+ double const * const x (thread_param->_x);
+ double* y(thread_param->_y);
+
+ for (unsigned i(beg_row); i<end_row; i++) {
+ y[i] = 0.0;
+ const unsigned end (iA[i+1]);
+ for (unsigned j(iA[i]); j<end; j++) {
+ y[i] += A[j] * x[jA[j]];
+ }
+ y[i] *= a;
+ }
+ return NULL;
+}
+} // end extern "C"
+
+void amuxCRSParallelPThreads (double a,
+ unsigned n, unsigned const * const iA, unsigned const * const jA,
+ double const * const A, double const * const x, double* y,
+ unsigned num_of_pthreads)
+{
+#ifdef HAVE_PTHREADS
+ // fill thread data objects
+ MatMultThreadParam** thread_param_array (new MatMultThreadParam*[num_of_pthreads]);
+ double step_size ((double)(n)/(double)(num_of_pthreads));
+ for (unsigned k(0); k<num_of_pthreads; k++) {
+ const unsigned beg (static_cast<unsigned>(k*step_size));
+ const unsigned end (static_cast<unsigned>((k+1)*step_size));
+ thread_param_array[k] = new MatMultThreadParam (a, beg, end, iA, jA, A, x, y);
+ }
+
+ // allocate thread_array and return value array
+ pthread_t *thread_array (new pthread_t[num_of_pthreads]);
+ int *ret_vals (new int[num_of_pthreads]);
+
+ // create threads
+ for (unsigned k(0); k<num_of_pthreads; k++) {
+ ret_vals[k] = pthread_create( &(thread_array[k]), NULL, amuxCRSpthread, thread_param_array[k]);
+ }
+
+ // join threads
+ for (unsigned k(0); k<num_of_pthreads; k++) {
+ pthread_join (thread_array[k], NULL);
+ }
+
+ delete [] ret_vals;
+ for (unsigned k(0); k<num_of_pthreads; k++)
+ delete thread_param_array[k];
+ delete [] thread_param_array;
+ delete [] thread_array;
+#else
+ amuxCRS (a, n, iA, jA, A, x, y);
+#endif
+}
+
+void amuxCRSParallelOpenMP (double a,
+ unsigned n, unsigned const * const iA, unsigned const * const jA,
+ double const * const A, double const * const x, double* y,
+ unsigned num_of_omp_threads)
+{
+ unsigned i;
+ omp_set_num_threads(num_of_omp_threads);
+ {
+ #pragma omp parallel for
+ for (i=0; i<n; i++) {
+ y[i] = 0.0;
+ const unsigned end (iA[i+1]);
+ for (unsigned j(iA[i]); j<end; j++) {
+ y[i] += A[j] * x[jA[j]];
+ }
+ y[i] *= a;
+ }
+ }
+}
+
+void amuxCRSSym (double a,
+ unsigned n, unsigned const * const iA, unsigned const * const jA,
+ double const * const A, double const * const x, double* y)
+{
+ for (unsigned i(0); i<n; i++) {
+ y[i] = 0.0;
+ }
+
+ for (unsigned i(0); i<n; i++) {
+ unsigned j (iA[i]);
+ // handle diagonal
+ if (jA[j] == i) {
+ y[i] += A[j] * x[jA[j]];
+ j++;
+ }
+ const unsigned end (iA[i+1]);
+ for (; j<end; j++) {
+ y[i] += A[j] * x[jA[j]];
+ y[jA[j]] += A[j] * x[i];
+ }
+ }
+
+ for (unsigned i(0); i<n; i++) {
+ y[i] *= a;
+ }
+}
+
diff --git a/MathLib/amuxCRS.h b/MathLib/amuxCRS.h
new file mode 100644
index 0000000000000000000000000000000000000000..24ec77d6fa41dac6530ab06821e5cf96a3213246
--- /dev/null
+++ b/MathLib/amuxCRS.h
@@ -0,0 +1,22 @@
+#ifndef AMUXCRS_H
+#define AMUXCRS_H
+
+void amuxCRS (double a,
+ unsigned n, unsigned const * const iA, unsigned const * const jA,
+ double const * const A, double const * const x, double* y);
+
+void amuxCRSParallelPThreads (double a,
+ unsigned n, unsigned const * const iA, unsigned const * const jA,
+ double const * const A, double const * const x, double* y,
+ unsigned num_of_pthreads);
+
+void amuxCRSParallelOpenMP (double a,
+ unsigned n, unsigned const * const iA, unsigned const * const jA,
+ double const * const A, double const * const x, double* y,
+ unsigned num_of_omp_threads);
+
+void amuxCRSSym (double a,
+ unsigned n, unsigned const * const iA, unsigned const * const jA,
+ double const * const A, double const * const x, double* y);
+
+#endif
diff --git a/MathLib/sparse.h b/MathLib/sparse.h
new file mode 100644
index 0000000000000000000000000000000000000000..44ecf09e79b58078d235b8bfccd847d8102b4cc4
--- /dev/null
+++ b/MathLib/sparse.h
@@ -0,0 +1,88 @@
+#ifndef SPARSE_H
+#define SPARSE_H
+
+#include <iostream>
+#include <cassert>
+
+extern void CS_write(char*, unsigned, unsigned*, unsigned*, double*);
+extern void CS_read(char*, unsigned&, unsigned*&, unsigned*&, double*&);
+
+template<class T> void
+CS_write(std::ostream &os, unsigned n, unsigned* iA, unsigned* jA, T* A)
+{
+ os.write((char*) &n, sizeof(unsigned));
+ os.write((char*) iA, (n+1)*sizeof(unsigned));
+ os.write((char*) jA, iA[n]*sizeof(unsigned));
+ os.write((char*) A, iA[n]*sizeof(T));
+}
+
+template<class T> void
+CS_read(std::istream &is, unsigned &n, unsigned* &iA, unsigned* &jA, T* &A)
+{
+ is.read((char*) &n, sizeof(unsigned));
+ if (iA != NULL) {
+ delete [] iA;
+ delete [] jA;
+ delete [] A;
+ }
+ iA = new unsigned[n+1];
+ assert(iA!=NULL);
+ is.read((char*) iA, (n+1)*sizeof(unsigned));
+
+ jA = new unsigned[iA[n]];
+ assert(jA!=NULL);
+ is.read((char*) jA, iA[n]*sizeof(unsigned));
+
+ A = new T[iA[n]];
+ assert(A!=NULL);
+ is.read((char*) A, iA[n]*sizeof(T));
+
+#ifndef NDEBUG
+ // do simple checks
+ if (iA[0]!=0)
+ std::cerr << std::endl
+ << "CRS matrix: array iA doesn't start with 0" << std::endl;
+
+ unsigned i = 0;
+ while (i<iA[n] && jA[i]<n) ++i;
+ if (i<iA[n])
+ std::cerr << std::endl
+ << "CRS matrix: the " << i << "th entry of jA has the value "
+ << jA[i] << ", which is out of bounds." << std::endl;
+#endif
+}
+
+void CS_transp(unsigned n, double* A, unsigned* jA, unsigned *iA,
+ double* B, unsigned *jB, unsigned *iB)
+{
+ unsigned nnz = iA[n];
+ unsigned *inz = new unsigned[n];
+
+ for (unsigned i=0; i<n; ++i) inz[i] = 0;
+
+ // compute number of entries of each column in A
+ for (unsigned l=0; l<nnz; ++l) inz[jA[l]]++;
+
+ // create iB
+ iB[0] = nnz = 0;
+ for (unsigned i=0; i<n; ++i) {
+ nnz += inz[i];
+ inz[i] = 0;
+ iB[i+1] = nnz;
+ }
+
+ // create arrays jB, B
+ unsigned l = iA[0];
+ for (unsigned i=0; i<n; ++i) {
+ while (l<iA[i+1]) {
+ unsigned j = jA[l];
+ unsigned k = iB[j] + inz[j]++;
+ B[k] = A[l++];
+ jB[k] = i;
+ }
+ }
+
+ delete [] inz;
+}
+
+#endif
diff --git a/SimpleTests/MatrixTests/.MatMult.cpp.swp b/SimpleTests/MatrixTests/.MatMult.cpp.swp
new file mode 100644
index 0000000000000000000000000000000000000000..808689bb3f0aefe083b5aa63d6b6f520e7d2869f
Binary files /dev/null and b/SimpleTests/MatrixTests/.MatMult.cpp.swp differ
diff --git a/SimpleTests/MatrixTests/CMakeLists.txt b/SimpleTests/MatrixTests/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..08a920431b8ac1f63950cc94cc0133be33015ec8
--- /dev/null
+++ b/SimpleTests/MatrixTests/CMakeLists.txt
@@ -0,0 +1,26 @@
+## pthread ##
+SET ( CMAKE_THREAD_PREFER_PTHREAD On )
+FIND_PACKAGE( Threads )
+IF ( CMAKE_USE_PTHREADS_INIT )
+ SET (HAVE_PTHREADS TRUE)
+ MESSAGE (STATUS "pthread library found." )
+ENDIF (CMAKE_USE_PTHREADS_INIT )
+
+INCLUDE_DIRECTORIES(
+ .
+ ../../Base/
+ ../../MathLib/
+)
+
+# Create the executable
+ADD_EXECUTABLE( MatMult
+ MatMult.cpp
+ ${SOURCES}
+ ${HEADERS}
+)
+
+TARGET_LINK_LIBRARIES ( MatMult
+ Base
+ MathLib
+)
+
diff --git a/SimpleTests/MatrixTests/MatMult.cpp b/SimpleTests/MatrixTests/MatMult.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..182e119d6deef48289d4e1fcf7037af6cef90236
--- /dev/null
+++ b/SimpleTests/MatrixTests/MatMult.cpp
@@ -0,0 +1,116 @@
+#include <fstream>
+#include <iostream>
+#include <cmath>
+#include <limits>
+#include <omp.h>
+#include <cstdlib>
+#include "sparse.h"
+#include "CRSMatrix.h"
+#include "CRSMatrixOpenMP.h"
+#include "CRSMatrixPThreads.h"
+#include "RunTimeTimer.h"
+#include "CPUTimeTimer.h"
+
+int main(int argc, char *argv[])
+{
+ if (argc < 4) {
+ std::cout << "Usage: " << argv[0] << " num_of_threads matrix number_of_multiplications resultfile" << std::endl;
+ exit (1);
+ }
+
+ // read number of threads
+ unsigned n_threads (1);
+ n_threads = atoi (argv[1]);
+
+ // read the number of multiplication to execute
+ unsigned n_mults (0);
+ n_mults = atoi (argv[3]);
+
+ std::string fname_mat (argv[2]);
+
+ bool verbose (true);
+
+ // *** reading matrix in crs format from file
+ std::ifstream in(fname_mat.c_str(), std::ios::in | std::ios::binary);
+ double *A(NULL);
+ unsigned *iA(NULL), *jA(NULL), n;
+ if (in) {
+ if (verbose) {
+ std::cout << "reading matrix from " << fname_mat << " ... " << std::flush;
+ }
+ RunTimeTimer timer;
+ timer.start();
+ CS_read(in, n, iA, jA, A);
+ timer.stop();
+ if (verbose) {
+ std::cout << "ok, " << timer.elapsed() << " s)" << std::endl;
+ }
+ } else {
+ std::cout << "error reading matrix from " << fname_mat << std::endl;
+ }
+ unsigned nnz(iA[n]);
+ if (verbose) {
+ std::cout << "Parameters read: n=" << n << ", nnz=" << nnz << std::endl;
+ }
+
+/*
+ // test: fill all nnz entries with 1
+ for (unsigned k(0); k<nnz; k++) A[k] = 2.0;
+*/
+
+ CRSMatrix<double> mat (n, iA, jA, A);
+// CRSMatrixOpenMP<double> mat (n, iA, jA, A, n_threads);
+// CRSMatrixPThreads<double> mat (n, iA, jA, A, n_threads);
+
+ double *x(new double[n]);
+ double *y(new double[n]);
+
+ for (unsigned k(0); k<n; ++k)
+ x[k] = 1.0;
+
+ if (verbose) {
+ std::cout << "matrix vector multiplication with Toms amuxCRS (" << n_threads << " threads) ... " << std::flush;
+ }
+ RunTimeTimer run_timer;
+ CPUTimeTimer cpu_timer;
+ run_timer.start();
+ cpu_timer.start();
+ for (size_t k(0); k<n_mults; k++) {
+// std::cout << "mult " << k << " ... " << std::flush;
+ mat.amux (1.0, x, y);
+// std::cout << "ok" << std::endl;
+ }
+ cpu_timer.stop();
+ run_timer.stop();
+
+/*
+ // test result of matrix vector mult
+ std::cout << "simple test of matrix vector mult ... " << std::flush;
+ for (unsigned k(0); k<n; k++) {
+ if (std::fabs(y[k] - 2*(iA[k+1]-iA[k])) > std::numeric_limits<double>::min()) {
+ std::cout << "value y[" << k << "] = " << y[k] << " not equal to number of entries in row " << k << ", that is " << iA[k+1]-iA[k] << std::endl;
+ }
+ }
+*/
+
+ if (verbose) {
+ std::cout << "done ["
+<< cpu_timer.elapsed() << " sec cpu time], [" << run_timer.elapsed() << " sec run time]" << std::endl;
+ } else {
+ if (argc == 5) {
+ std::ofstream result_os (argv[4], std::ios::app);
+ if (result_os) {
+ result_os << cpu_timer.elapsed() << "\t" << run_timer.elapsed() << std::endl;
+ }
+ result_os.close();
+ } else {
+ std::cout << cpu_timer.elapsed() << "\t" << run_timer.elapsed() << std::endl;
+ }
+ }
+
+ delete [] x;
+ delete [] y;
+
+ return 0;
+}
+