From 66cd7b7ee3d99731664c984b2425d2241feb6de5 Mon Sep 17 00:00:00 2001
From: Thomas Fischer <thomas.fischer@ufz.de>
Date: Thu, 3 Nov 2011 15:39:10 +0100
Subject: [PATCH] a) added to class templates SparseMatrixBase and CRSMatrix
another template parameter for index sets b) added method getNNZ() to class
template CRSMatrix
---
MathLib/LinAlg/Solvers/BiCGStab.cpp | 2 +-
MathLib/LinAlg/Solvers/BiCGStab.h | 2 +-
MathLib/LinAlg/Solvers/CG.cpp | 2 +-
MathLib/LinAlg/Solvers/CG.h | 4 +--
MathLib/LinAlg/Solvers/GMRes.cpp | 4 +--
MathLib/LinAlg/Solvers/GMRes.h | 2 +-
MathLib/LinAlg/Sparse/CRSMatrix.h | 29 ++++++++++---------
MathLib/LinAlg/Sparse/CRSMatrixDiagPrecond.h | 4 +--
MathLib/LinAlg/Sparse/CRSMatrixOpenMP.h | 10 +++----
MathLib/LinAlg/Sparse/CRSMatrixPThreads.h | 12 ++++----
MathLib/LinAlg/Sparse/SparseMatrixBase.h | 8 ++---
.../ConjugateGradientUnpreconditioned.cpp | 2 +-
12 files changed, 43 insertions(+), 38 deletions(-)
diff --git a/MathLib/LinAlg/Solvers/BiCGStab.cpp b/MathLib/LinAlg/Solvers/BiCGStab.cpp
index a63da219685..7f6a4ef3deb 100644
--- a/MathLib/LinAlg/Solvers/BiCGStab.cpp
+++ b/MathLib/LinAlg/Solvers/BiCGStab.cpp
@@ -12,7 +12,7 @@
namespace MathLib {
-unsigned BiCGStab(CRSMatrix<double> const& A, double* const b, double* const x,
+unsigned BiCGStab(CRSMatrix<double, unsigned> const& A, double* const b, double* const x,
double& eps, unsigned& nsteps)
{
const unsigned N(A.getNRows());
diff --git a/MathLib/LinAlg/Solvers/BiCGStab.h b/MathLib/LinAlg/Solvers/BiCGStab.h
index 1f7821b7f29..cdf699d336a 100644
--- a/MathLib/LinAlg/Solvers/BiCGStab.h
+++ b/MathLib/LinAlg/Solvers/BiCGStab.h
@@ -14,7 +14,7 @@
namespace MathLib {
-unsigned BiCGStab(CRSMatrix<double> const& A, double* const b, double* const x,
+unsigned BiCGStab(CRSMatrix<double, unsigned> const& A, double* const b, double* const x,
double& eps, unsigned& nsteps);
} // end namespace MathLib
diff --git a/MathLib/LinAlg/Solvers/CG.cpp b/MathLib/LinAlg/Solvers/CG.cpp
index 05c2f1d977f..0c37460a660 100644
--- a/MathLib/LinAlg/Solvers/CG.cpp
+++ b/MathLib/LinAlg/Solvers/CG.cpp
@@ -27,7 +27,7 @@
namespace MathLib {
-unsigned CG(CRSMatrix<double> const * mat, double const * const b,
+unsigned CG(CRSMatrix<double,unsigned> const * mat, double const * const b,
double* const x, double& eps, unsigned& nsteps, unsigned num_threads)
{
unsigned N = mat->getNRows();
diff --git a/MathLib/LinAlg/Solvers/CG.h b/MathLib/LinAlg/Solvers/CG.h
index 36cd6917402..4b07debdf6a 100644
--- a/MathLib/LinAlg/Solvers/CG.h
+++ b/MathLib/LinAlg/Solvers/CG.h
@@ -11,9 +11,9 @@
namespace MathLib {
// forward declaration
-template <typename T> class CRSMatrix;
+template <typename PF_TYPE, typename IDX_TYPE> class CRSMatrix;
-unsigned CG(CRSMatrix<double> const * mat, double const * const b,
+unsigned CG(CRSMatrix<double,unsigned> const * mat, double const * const b,
double* const x, double& eps, unsigned& nsteps, unsigned num_threads = 1);
} // end namespace MathLib
diff --git a/MathLib/LinAlg/Solvers/GMRes.cpp b/MathLib/LinAlg/Solvers/GMRes.cpp
index 45aeaf5a734..871431f40d8 100644
--- a/MathLib/LinAlg/Solvers/GMRes.cpp
+++ b/MathLib/LinAlg/Solvers/GMRes.cpp
@@ -37,7 +37,7 @@ inline void applPlRot(double& dx, double& dy, double cs, double sn)
}
// solve H y = s and update x += MVy
-static void update(const CRSMatrix<double>& A, unsigned k, double* H,
+static void update(const CRSMatrix<double,unsigned>& A, unsigned k, double* H,
unsigned ldH, double* s, double* V, double* x)
{
const size_t n(A.getNRows());
@@ -59,7 +59,7 @@ static void update(const CRSMatrix<double>& A, unsigned k, double* H,
delete[] y;
}
-unsigned GMRes(const CRSMatrix<double>& A, double* const b, double* const x,
+unsigned GMRes(const CRSMatrix<double,unsigned>& A, double* const b, double* const x,
double& eps, unsigned m, unsigned& nsteps)
{
double resid;
diff --git a/MathLib/LinAlg/Solvers/GMRes.h b/MathLib/LinAlg/Solvers/GMRes.h
index f7acd6e3e01..732c560b425 100644
--- a/MathLib/LinAlg/Solvers/GMRes.h
+++ b/MathLib/LinAlg/Solvers/GMRes.h
@@ -12,7 +12,7 @@
namespace MathLib {
-unsigned GMRes(const CRSMatrix<double>& mat, double* const b, double* const x,
+unsigned GMRes(const CRSMatrix<double,unsigned>& mat, double* const b, double* const x,
double& eps, unsigned m, unsigned& steps);
} // end namespace MathLib
diff --git a/MathLib/LinAlg/Sparse/CRSMatrix.h b/MathLib/LinAlg/Sparse/CRSMatrix.h
index a2c889cf392..493a00af4a3 100644
--- a/MathLib/LinAlg/Sparse/CRSMatrix.h
+++ b/MathLib/LinAlg/Sparse/CRSMatrix.h
@@ -18,29 +18,30 @@
namespace MathLib {
-template<class T> class CRSMatrix : public SparseMatrixBase<T>
+template<typename FP_TYPE, typename IDX_TYPE>
+class CRSMatrix: public SparseMatrixBase<FP_TYPE, IDX_TYPE>
{
public:
CRSMatrix(std::string const &fname) :
- SparseMatrixBase<T>(),
+ SparseMatrixBase<FP_TYPE, IDX_TYPE>(),
_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;
+ CS_read(in, SparseMatrixBase<FP_TYPE, IDX_TYPE>::_n_rows, _row_ptr, _col_idx, _data);
+ SparseMatrixBase<FP_TYPE, IDX_TYPE>::_n_cols = SparseMatrixBase<FP_TYPE, IDX_TYPE>::_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(IDX_TYPE n, IDX_TYPE *iA, IDX_TYPE *jA, FP_TYPE* A) :
+ SparseMatrixBase<FP_TYPE, IDX_TYPE>(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)
+ CRSMatrix(IDX_TYPE n1) :
+ SparseMatrixBase<FP_TYPE, IDX_TYPE>(n1, n1), _row_ptr(NULL), _col_idx(NULL), _data(NULL)
{}
virtual ~CRSMatrix()
@@ -50,18 +51,20 @@ public:
delete [] _data;
}
- virtual void amux(T d, T const * const x, T *y) const
+ virtual void amux(FP_TYPE d, FP_TYPE const * const x, FP_TYPE *y) const
{
amuxCRS(d, MatrixBase::_n_rows, _row_ptr, _col_idx, _data, x, y);
}
- virtual void precondApply(T* x) const
+ virtual void precondApply(FP_TYPE* x) const
{}
+ IDX_TYPE getNNZ() const { return _row_ptr[MatrixBase::_n_rows]; }
+
protected:
- unsigned *_row_ptr;
- unsigned *_col_idx;
- T* _data;
+ IDX_TYPE *_row_ptr;
+ IDX_TYPE *_col_idx;
+ FP_TYPE* _data;
};
} // end namespace MathLib
diff --git a/MathLib/LinAlg/Sparse/CRSMatrixDiagPrecond.h b/MathLib/LinAlg/Sparse/CRSMatrixDiagPrecond.h
index 7bee253808a..49fad760283 100644
--- a/MathLib/LinAlg/Sparse/CRSMatrixDiagPrecond.h
+++ b/MathLib/LinAlg/Sparse/CRSMatrixDiagPrecond.h
@@ -9,11 +9,11 @@
namespace MathLib {
-class CRSMatrixDiagPrecond : public CRSMatrix<double>
+class CRSMatrixDiagPrecond : public CRSMatrix<double, unsigned>
{
public:
CRSMatrixDiagPrecond (std::string const &fname)
- : CRSMatrix<double>(fname), _inv_diag(new double[_n_rows])
+ : CRSMatrix<double, unsigned>(fname), _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;
diff --git a/MathLib/LinAlg/Sparse/CRSMatrixOpenMP.h b/MathLib/LinAlg/Sparse/CRSMatrixOpenMP.h
index 408304ae8b6..b2467bb86e2 100644
--- a/MathLib/LinAlg/Sparse/CRSMatrixOpenMP.h
+++ b/MathLib/LinAlg/Sparse/CRSMatrixOpenMP.h
@@ -15,18 +15,18 @@
namespace MathLib {
-template<class T> class CRSMatrixOpenMP : public CRSMatrix<T> {
+template<class T> class CRSMatrixOpenMP : public CRSMatrix<T, unsigned> {
public:
CRSMatrixOpenMP(std::string const &fname, unsigned num_of_threads) :
- CRSMatrix<T>(fname), _num_of_threads (num_of_threads)
+ CRSMatrix<T, unsigned>(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)
+ CRSMatrix<T, unsigned>(n, iA, jA, A), _num_of_threads (num_of_threads)
{}
CRSMatrixOpenMP(unsigned n1) :
- CRSMatrix<T>(n1), _num_of_threads (1)
+ CRSMatrix<T, unsigned>(n1), _num_of_threads (1)
{}
virtual ~CRSMatrixOpenMP()
@@ -34,7 +34,7 @@ public:
virtual void amux(T d, T const * const x, T *y) const
{
- amuxCRSParallelOpenMP(d, MatrixBase::_n_rows, CRSMatrix<T>::_row_ptr, CRSMatrix<T>::_col_idx, CRSMatrix<T>::_data, x, y, _num_of_threads);
+ amuxCRSParallelOpenMP(d, MatrixBase::_n_rows, CRSMatrix<T,unsigned>::_row_ptr, CRSMatrix<T,unsigned>::_col_idx, CRSMatrix<T,unsigned>::_data, x, y, _num_of_threads);
}
private:
diff --git a/MathLib/LinAlg/Sparse/CRSMatrixPThreads.h b/MathLib/LinAlg/Sparse/CRSMatrixPThreads.h
index eeabaac7617..32ca857f505 100644
--- a/MathLib/LinAlg/Sparse/CRSMatrixPThreads.h
+++ b/MathLib/LinAlg/Sparse/CRSMatrixPThreads.h
@@ -17,19 +17,19 @@
namespace MathLib {
-template<class T> class CRSMatrixPThreads : public CRSMatrix<T>
+template<class T> class CRSMatrixPThreads : public CRSMatrix<T,unsigned>
{
public:
CRSMatrixPThreads(std::string const &fname, unsigned num_of_threads) :
- CRSMatrix<T>(fname), _num_of_threads (num_of_threads)
+ CRSMatrix<T,unsigned>(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)
+ CRSMatrix<T,unsigned>(n, iA, jA, A), _num_of_threads (num_of_threads)
{}
CRSMatrixPThreads(unsigned n1) :
- CRSMatrix<T>(n1), _num_of_threads (1)
+ CRSMatrix<T,unsigned>(n1), _num_of_threads (1)
{}
virtual ~CRSMatrixPThreads()
@@ -37,7 +37,9 @@ public:
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);
+ amuxCRSParallelPThreads(d, SparseMatrixBase<T, unsigned>::_n_rows,
+ CRSMatrix<T, unsigned>::_row_ptr, CRSMatrix<T, unsigned>::_col_idx,
+ CRSMatrix<T, unsigned>::_data, x, y, _num_of_threads);
}
protected:
diff --git a/MathLib/LinAlg/Sparse/SparseMatrixBase.h b/MathLib/LinAlg/Sparse/SparseMatrixBase.h
index d6c65de2747..51d612da5c5 100644
--- a/MathLib/LinAlg/Sparse/SparseMatrixBase.h
+++ b/MathLib/LinAlg/Sparse/SparseMatrixBase.h
@@ -5,13 +5,13 @@
namespace MathLib {
-template<class T> class SparseMatrixBase : public MatrixBase
+template<typename FP_TYPE, typename IDX_TYPE> class SparseMatrixBase : public MatrixBase
{
public:
- SparseMatrixBase(unsigned n1, unsigned n2) : MatrixBase (n1,n2) {}
+ SparseMatrixBase(IDX_TYPE n1, IDX_TYPE n2) : MatrixBase (n1,n2) {}
SparseMatrixBase() : MatrixBase () {}
- virtual void amux(T d, T const * const x, T *y) const = 0; // y +=d*Ax
- virtual ~SparseMatrixBase() { }
+ virtual void amux(FP_TYPE d, FP_TYPE const * const x, FP_TYPE *y) const = 0; // y +=d*Ax
+ virtual ~SparseMatrixBase() {};
};
} // end namespace MathLib
diff --git a/SimpleTests/SolverTests/ConjugateGradientUnpreconditioned.cpp b/SimpleTests/SolverTests/ConjugateGradientUnpreconditioned.cpp
index 2413a1502de..52f1ba239f7 100644
--- a/SimpleTests/SolverTests/ConjugateGradientUnpreconditioned.cpp
+++ b/SimpleTests/SolverTests/ConjugateGradientUnpreconditioned.cpp
@@ -12,7 +12,7 @@ int main(int argc, char *argv[])
// *** reading matrix in crs format from file
std::string fname("/work/fischeth/data/testmat.bin");
// std::ifstream in(fname.c_str(), std::ios::binary);
- MathLib::CRSMatrix<double> *mat (new MathLib::CRSMatrix<double>(fname));
+ MathLib::CRSMatrix<double, unsigned> *mat (new MathLib::CRSMatrix<double, unsigned>(fname));
/* double *A(NULL);
unsigned *iA(NULL), *jA(NULL), n;
if (in) {
--
GitLab