Skip to content
Snippets Groups Projects
Commit 56c388ef authored by Tom Fischer's avatar Tom Fischer
Browse files

my initial contribution to OGS-6

parent 328cd770
No related branches found
No related tags found
No related merge requests found
Showing
with 733 additions and 0 deletions
# 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)
#include "CPUTimeTimer.h"
void CPUTimeTimer::start()
{
_start = clock();
}
void CPUTimeTimer::stop()
{
_stop = clock();
}
double CPUTimeTimer::elapsed()
{
return (_stop-_start)/(double)(CLOCKS_PER_SEC);
}
#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
/**
* \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
#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));
}
#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
#ifndef TIMEMEASUREMENT_H
#define TIMEMEASUREMENT_H
class TimeMeasurementBase
{
public:
virtual void start () = 0;
virtual void stop () = 0;
virtual double elapsed () = 0;
virtual ~TimeMeasurementBase () {};
};
#endif
# 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 )
# 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)
#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
#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
/*
* 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_ */
/*
* 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
/*
* 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_ */
#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
#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;
}
}
#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
#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
File added
## 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
)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment