diff --git a/MathLib/LinAlg/Solvers/BiCGStab.cpp b/MathLib/LinAlg/Solvers/BiCGStab.cpp
index a63da21968505448aa71a3cecc8389e37f8d6f0e..7f6a4ef3deb685f84ecb47141ad857585016b1de 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 1f7821b7f29494955e229bfa7d0d9923013ce725..cdf699d336a27bffcf621a140e109ae03880cce2 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 05c2f1d977f0943cf644b04b4473a172c4e55c39..0c37460a6603aaddb6cfe46c608ebb4b3d8958d3 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 36cd6917402f478be69b95e0bc1f90878894dc9e..4b07debdf6a01ac5293e0e8ceb9d8810245bde76 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 45aeaf5a7348c4f412dab9a24e22993fc520065d..871431f40d8023210ac1eca11b6138d12ff3a3d8 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 f7acd6e3e010d1c9a919f359ae732c379183c38a..732c560b425b78483a16c2179fc91540d1658831 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 a2c889cf3926acd1967fbdc0fadd95e299438283..f2995a353fc7907b834c0c3bc067bbed5d5595ed 100644
--- a/MathLib/LinAlg/Sparse/CRSMatrix.h
+++ b/MathLib/LinAlg/Sparse/CRSMatrix.h
@@ -11,6 +11,12 @@
#include <string>
#include <fstream>
#include <iostream>
+#include <cassert>
+
+// Base
+#include "swap.h"
+
+// MathLib
#include "SparseMatrixBase.h"
#include "sparse.h"
#include "amuxCRS.h"
@@ -18,29 +24,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 +57,292 @@ 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
{}
+ /**
+ * get the number of non-zero entries
+ * @return number of non-zero entries
+ */
+ IDX_TYPE getNNZ() const { return _row_ptr[MatrixBase::_n_rows]; }
+
+ /**
+ * This is the constant access operator to a non-zero matrix entry.
+ * Precondition: the entries have to be in the sparsity pattern!
+ * @param row the row number
+ * @param col the column number
+ * @return The corresponding matrix entry.
+ */
+ FP_TYPE operator() (IDX_TYPE row, IDX_TYPE col) const
+ {
+ assert(0 <= row && row < MatrixBase::_n_rows);
+
+ // linear search - for matrices with many entries per row binary search is much faster
+ const IDX_TYPE idx_end (_row_ptr[row+1]);
+ IDX_TYPE j(_row_ptr[row]), k;
+
+ while (j<idx_end && (k=_col_idx[j]) <= col) {
+ if (k == col) {
+ return _data[j];
+ }
+ j++;
+ }
+ return 0.0;
+ }
+
+ /**
+ * This is the non constant access operator to a non-zero matrix entry.
+ * Precondition: the entries have to be in the sparsity pattern!
+ * @param row the row number
+ * @param col the column number
+ * @return The corresponding matrix entry that can be modified by the user.
+ */
+ FP_TYPE & operator() (IDX_TYPE row, IDX_TYPE col)
+ {
+ assert(0 <= row && row < MatrixBase::_n_rows);
+
+ // linear search - for matrices with many entries per row binary search is much faster
+ const IDX_TYPE idx_end (_row_ptr[row+1]);
+ IDX_TYPE j(_row_ptr[row]), k;
+
+ while (j<idx_end && (k=_col_idx[j]) <= col) {
+ if (k == col) {
+ return _data[j];
+ }
+ j++;
+ }
+ }
+
+ /**
+ * get const access to the row pointer array of CRS matrix
+ * @return the index array _row_ptr
+ */
+ IDX_TYPE const* getRowPtrArray() const { return _row_ptr; }
+
+ /**
+ * get const access to the column index array of CRS matrix
+ * @return the index array _col_idx
+ */
+ IDX_TYPE const* getColIdxArray ()const { return _col_idx; }
+
+ /**
+ * get the matrix entries within an array of CRS matrix
+ * @return
+ */
+ FP_TYPE const* getEntryArray() const { return _data; }
+
+ /**
+ * erase rows and columns from sparse matrix
+ * @param n_rows_cols number of rows / columns to remove
+ * @param rows sorted list of rows that should be removed
+ * @param cols sorted list of columns that should be removed
+ */
+ void eraseEntries (IDX_TYPE n_rows_cols,
+ IDX_TYPE const*const rows, IDX_TYPE const*const cols)
+ {
+// std::cout << "original: " << std::endl;
+// printMat();
+ IDX_TYPE n_cols(MatrixBase::_n_rows);
+ //*** remove the rows
+ removeRows(n_rows_cols, rows);
+// std::cout << "removed rows: " << std::endl;
+// printMat();
+ //*** transpose
+ transpose(n_cols);
+// std::cout << "transposed: " << std::endl;
+// printMat();
+ //*** remove columns in original means removing rows in the transposed
+ removeRows(n_rows_cols, cols);
+// std::cout << "removed rows: " << std::endl;
+// printMat();
+ //*** transpose again
+ transpose(MatrixBase::_n_rows);
+// std::cout << "transposed: " << std::endl;
+// printMat();
+// //*** apply changes to the column indices
+// correctColumnIndices(n_rows_cols, cols);
+ }
+
+ /**
+ * get the j-th column of the sparse matrix
+ * @param j the column number that should be returned
+ * @param column_entries the column entries (have to be allocated
+ */
+ void getColumn(IDX_TYPE j, FP_TYPE* column_entries) const
+ {
+ for (IDX_TYPE k(0); k<MatrixBase::_n_rows; k++) {
+ const IDX_TYPE end_row(_row_ptr[k+1]);
+ IDX_TYPE i(_row_ptr[k+1]);
+ while (i<end_row && _col_idx[i] != j) {
+ i++;
+ }
+ if (i==end_row) {
+ column_entries[k] = 0.0;
+ } else {
+ column_entries[k] = _data[i];
+ }
+ }
+ }
+
protected:
- unsigned *_row_ptr;
- unsigned *_col_idx;
- T* _data;
+ void removeRows (IDX_TYPE n_rows_cols, IDX_TYPE const*const rows)
+ {
+ //*** determine the number of new rows and the number of entries without the rows
+ const IDX_TYPE n_new_rows(MatrixBase::_n_rows - n_rows_cols);
+ IDX_TYPE *row_ptr_new(new IDX_TYPE[n_new_rows+1]);
+ row_ptr_new[0] = 0;
+ IDX_TYPE row_cnt (1), erase_row_cnt(0);
+ for (unsigned k(0); k<MatrixBase::_n_rows; k++) {
+ if (k != rows[erase_row_cnt]) {
+ row_ptr_new[row_cnt] = _row_ptr[k+1] - _row_ptr[k];
+ row_cnt++;
+ } else {
+ erase_row_cnt++;
+ }
+ }
+
+ //*** sum up the entries
+ for (IDX_TYPE k(0); k<n_new_rows; k++) {
+ row_ptr_new[k+1] = row_ptr_new[k+1] + row_ptr_new[k];
+ }
+
+ //*** create new memory for col_idx and data
+ IDX_TYPE nnz_new(row_ptr_new[n_new_rows]);
+ IDX_TYPE *col_idx_new (new IDX_TYPE[nnz_new]);
+ FP_TYPE *data_new (new FP_TYPE[nnz_new]);
+
+ //*** copy the entries
+ // initialization
+ IDX_TYPE *row_ptr_new_tmp(new IDX_TYPE[n_new_rows+1]);
+ for (unsigned k(0); k<=n_new_rows; k++) {
+ row_ptr_new_tmp[k] = row_ptr_new[k];
+ }
+ erase_row_cnt = 0;
+ row_cnt = 0;
+ // copy column index and data entries
+ for (IDX_TYPE k(0); k<MatrixBase::_n_rows; k++) {
+ if (k != rows[erase_row_cnt]) {
+ const IDX_TYPE end (_row_ptr[k+1]);
+ // walk through row
+ for (IDX_TYPE j(_row_ptr[k]); j<end; j++) {
+ col_idx_new[row_ptr_new_tmp[row_cnt]] = _col_idx[j];
+ data_new[row_ptr_new_tmp[row_cnt]] = _data[j];
+ row_ptr_new_tmp[row_cnt]++;
+ }
+ row_cnt++;
+ } else {
+ erase_row_cnt++;
+ }
+ }
+
+ MatrixBase::_n_rows -= n_rows_cols;
+ BASELIB::swap (row_ptr_new, _row_ptr);
+ BASELIB::swap (col_idx_new, _col_idx);
+ BASELIB::swap (data_new, _data);
+
+ delete [] row_ptr_new;
+ delete [] col_idx_new;
+ delete [] data_new;
+ }
+
+ void transpose (IDX_TYPE n_cols)
+ {
+ // create a helper array row_ptr_nnz
+ IDX_TYPE *row_ptr_nnz(new IDX_TYPE[n_cols+1]);
+ for (IDX_TYPE k(0); k <= n_cols; k++) {
+ row_ptr_nnz[k] = 0;
+ }
+
+ // count entries per row in the transposed matrix
+ IDX_TYPE nnz(_row_ptr[MatrixBase::_n_rows]);
+ for (IDX_TYPE k(0); k < nnz; k++) {
+ row_ptr_nnz[_col_idx[k]]++;
+ }
+
+ // create row_ptr_trans
+ IDX_TYPE *row_ptr_trans(new IDX_TYPE[n_cols + 1]);
+ row_ptr_trans[0] = 0;
+ for (IDX_TYPE k(0); k < n_cols; k++) {
+ row_ptr_trans[k+1] = row_ptr_trans[k] + row_ptr_nnz[k];
+ }
+
+ // make a copy of row_ptr_trans
+ for (IDX_TYPE k(0); k <= n_cols; k++) {
+ row_ptr_nnz[k] = row_ptr_trans[k];
+ }
+
+ // create arrays col_idx_trans and data_trans
+ assert(nnz == row_ptr_trans[n_cols]);
+ IDX_TYPE *col_idx_trans(new IDX_TYPE[nnz]);
+ FP_TYPE *data_trans(new FP_TYPE[nnz]);
+
+ // fill arrays col_idx_trans and data_trans
+ for (IDX_TYPE i(0); i < MatrixBase::_n_rows; i++) {
+ const IDX_TYPE row_end(_row_ptr[i + 1]);
+ for (IDX_TYPE j(_row_ptr[i]); j < row_end; j++) {
+ const IDX_TYPE k(_col_idx[j]);
+ col_idx_trans[row_ptr_nnz[k]] = i;
+ data_trans[row_ptr_nnz[k]] = _data[j];
+ row_ptr_nnz[k]++;
+ }
+ }
+
+ MatrixBase::_n_rows = n_cols;
+ BASELIB::swap(row_ptr_trans, _row_ptr);
+ BASELIB::swap(col_idx_trans, _col_idx);
+ BASELIB::swap(data_trans, _data);
+
+ delete[] row_ptr_nnz;
+ delete[] row_ptr_trans;
+ delete[] col_idx_trans;
+ delete[] data_trans;
+ }
+
+ void correctColumnIndices(IDX_TYPE n_rows_cols, IDX_TYPE const*const cols)
+ {
+ // create the mapping
+ const IDX_TYPE size(MatrixBase::_n_rows + n_rows_cols);
+ IDX_TYPE *mapping(new IDX_TYPE[size]);
+ IDX_TYPE k(0);
+ for (IDX_TYPE i(0); i < n_rows_cols; i++) {
+ while (k < cols[i]) {
+ mapping[k] = k - i;
+ k++;
+ }
+ }
+ for (IDX_TYPE k(cols[n_rows_cols - 1]); k < size; k++) {
+ mapping[k] = k - n_rows_cols;
+ }
+
+ // apply mapping to _col_idx
+ const IDX_TYPE nnz(_row_ptr[MatrixBase::_n_rows]);
+ for (IDX_TYPE k(0); k < nnz; k++) {
+ _col_idx[k] = mapping[_col_idx[k]];
+ }
+
+ delete[] mapping;
+ }
+
+ void printMat() const
+ {
+ for (IDX_TYPE k(0); k<MatrixBase::_n_rows; k++) {
+ std::cout << k << ": " << std::flush;
+ const IDX_TYPE row_end(_row_ptr[k+1]);
+ for (IDX_TYPE j(_row_ptr[k]); j<row_end; j++) {
+ std::cout << _col_idx[j] << " " << std::flush;
+ }
+ std::cout << std::endl;
+ }
+ }
+
+ 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 7bee253808a0dff9e7112d3d094d8d77a8f04f19..49fad7602837f2e4947c49b3ee8db9a4806cb734 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 408304ae8b69f0eb1b44b397567c0ea770854c05..b2467bb86e2c28729e66c8f2e2e48898357ead40 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 eeabaac76173e1d5c676c626a8cb1e7694c19816..32ca857f505bd713bb3c771623da92dedc8e196c 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/CRSTranspose.h b/MathLib/LinAlg/Sparse/CRSTranspose.h
index 6b3eba986045d4327dd39bb15fc5dabd5f81e6ea..010a49e5a31126c57d14232f3185573857892aff 100644
--- a/MathLib/LinAlg/Sparse/CRSTranspose.h
+++ b/MathLib/LinAlg/Sparse/CRSTranspose.h
@@ -2,17 +2,17 @@
* CRSTranspose.h
*
* Created on: Sep 27, 2011
- * Author: fischeth
+ * Author: TF
*/
#ifndef CRSTRANSPOSE_H_
#define CRSTRANSPOSE_H_
-void CS_transp(unsigned n, double* A, unsigned* jA, unsigned *iA,
- double* B, unsigned *jB, unsigned *iB)
+void CS_transp(unsigned n, unsigned *iA, unsigned* jA, double* A,
+ unsigned *iB, unsigned *jB, double* B)
{
unsigned nnz = iA[n];
- unsigned *inz = new unsigned[n];
+ unsigned *inz(new unsigned[n]);
for (unsigned i=0; i<n; ++i) inz[i] = 0;
diff --git a/MathLib/LinAlg/Sparse/SparseMatrixBase.h b/MathLib/LinAlg/Sparse/SparseMatrixBase.h
index d6c65de2747ab2728aaf953f9ef4a7454b94e9ac..51d612da5c5a6171040a14df0a8e8c73d8f97555 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/MathLib/LinAlg/Sparse/amuxCRS.cpp b/MathLib/LinAlg/Sparse/amuxCRS.cpp
index 79b847d2c135c6f60086705c27f753fe6b147741..88c112409771098affce9b7d4559b3341dc09cb6 100644
--- a/MathLib/LinAlg/Sparse/amuxCRS.cpp
+++ b/MathLib/LinAlg/Sparse/amuxCRS.cpp
@@ -14,20 +14,6 @@
namespace MathLib {
-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,
@@ -49,14 +35,14 @@ struct MatMultThreadParam {
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);
+ 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++) {
@@ -110,24 +96,23 @@ void amuxCRSParallelPThreads (double a,
#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)
+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;
- }
- }
+ 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,
@@ -139,17 +124,17 @@ void amuxCRSSym (double a,
}
for (unsigned i(0); i<n; i++) {
- unsigned j (iA[i]);
- // handle diagonal
- if (jA[j] == 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]];
- 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];
- }
+ y[jA[j]] += A[j] * x[i];
+ }
}
for (unsigned i(0); i<n; i++) {
diff --git a/MathLib/LinAlg/Sparse/amuxCRS.h b/MathLib/LinAlg/Sparse/amuxCRS.h
index b1845218f8c5e8f3ad6f9a48c448b7e1eada5a67..fbe9d2ad97ec3f026b4c0f0f153a4eaa34ffa9cc 100644
--- a/MathLib/LinAlg/Sparse/amuxCRS.h
+++ b/MathLib/LinAlg/Sparse/amuxCRS.h
@@ -3,9 +3,19 @@
namespace MathLib {
-void amuxCRS (double a,
- unsigned n, unsigned const * const iA, unsigned const * const jA,
- double const * const A, double const * const x, double* y);
+template<typename FP_TYPE, typename IDX_TYPE>
+void amuxCRS(FP_TYPE a, IDX_TYPE n, IDX_TYPE const * const iA, IDX_TYPE const * const jA,
+ FP_TYPE const * const A, FP_TYPE const * const x, FP_TYPE* y)
+{
+ for (IDX_TYPE i(0); i < n; i++) {
+ y[i] = 0.0;
+ const IDX_TYPE end(iA[i + 1]);
+ for (IDX_TYPE j(iA[i]); j < end; j++) {
+ y[i] += A[j] * x[jA[j]];
+ }
+ y[i] *= a;
+ }
+}
void amuxCRSParallelPThreads (double a,
unsigned n, unsigned const * const iA, unsigned const * const jA,
diff --git a/MathLib/sparse.h b/MathLib/sparse.h
index fa3f6623659f7ec59afca7531d7b9028af31975e..3278ea93c699bee6803959e543d1dc22371ff09d 100644
--- a/MathLib/sparse.h
+++ b/MathLib/sparse.h
@@ -4,51 +4,48 @@
#include <iostream>
#include <cassert>
-extern void CS_write(char*, unsigned, unsigned*, unsigned*, double*);
+extern void CS_write(char*, unsigned, unsigned const*, unsigned const*, double const*);
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)
+template<class T> void CS_write(std::ostream &os, unsigned n, unsigned const* iA, unsigned const* jA, T const* 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));
+ 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)
+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));
+ 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;
+ // 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
}
diff --git a/SimpleTests/MatrixTests/CMakeLists.txt b/SimpleTests/MatrixTests/CMakeLists.txt
index bda3c0ca41a4d64b35dab5b17a6418d155cb8afb..50238c4a25c2bcd8a3082b9bdc6127ba148e0bca 100644
--- a/SimpleTests/MatrixTests/CMakeLists.txt
+++ b/SimpleTests/MatrixTests/CMakeLists.txt
@@ -19,6 +19,13 @@ ADD_EXECUTABLE( MatMult
${HEADERS}
)
+# Create the executable
+ADD_EXECUTABLE( MatTestRemoveRowsCols
+ MatTestRemoveRowsCols.cpp
+ ${SOURCES}
+ ${HEADERS}
+)
+
IF (WIN32)
TARGET_LINK_LIBRARIES(MatMult Winmm.lib)
ENDIF (WIN32)
@@ -28,3 +35,12 @@ TARGET_LINK_LIBRARIES ( MatMult
MathLib
)
+IF (WIN32)
+ TARGET_LINK_LIBRARIES(MatTestRemoveRowsCols Winmm.lib)
+ENDIF (WIN32)
+
+TARGET_LINK_LIBRARIES ( MatTestRemoveRowsCols
+ Base
+ MathLib
+)
+
diff --git a/SimpleTests/MatrixTests/MatTestRemoveRowsCols.cpp b/SimpleTests/MatrixTests/MatTestRemoveRowsCols.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1d27d44c9956fdc6ec144f1459d76baaaa382f1f
--- /dev/null
+++ b/SimpleTests/MatrixTests/MatTestRemoveRowsCols.cpp
@@ -0,0 +1,71 @@
+/*
+ * MatTestRemoveRowsCols.cpp
+ *
+ * Created on: Nov 8, 2011
+ * Author: TF
+ */
+
+#include <fstream>
+#include <iostream>
+
+// Base
+#include "RunTimeTimer.h"
+#include "CPUTimeTimer.h"
+
+// MathLib
+#include "LinAlg/Sparse/CRSMatrix.h"
+
+int main(int argc, char *argv[])
+{
+ if (argc < 3) {
+ std::cout << "Usage: " << argv[0] << " input-matrix output-matrix" << std::endl;
+ return 1;
+ }
+
+ std::string fname_mat (argv[1]);
+ 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);
+ in.close();
+ 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;
+ }
+
+ MathLib::CRSMatrix<double, unsigned> mat (n, iA, jA, A);
+
+ const unsigned n_rows_cols_to_erase(3);
+ unsigned *rows_to_erase(new unsigned[n_rows_cols_to_erase]);
+ unsigned *cols_to_erase(new unsigned[n_rows_cols_to_erase]);
+
+ for (unsigned k(0); k<n_rows_cols_to_erase; k++) {
+ rows_to_erase[k] = (k+1)*2;
+ cols_to_erase[k] = (k+1)*2;
+ }
+
+ mat.eraseEntries(n_rows_cols_to_erase, rows_to_erase, cols_to_erase);
+ delete[] rows_to_erase;
+ delete[] cols_to_erase;
+
+ fname_mat = argv[2];
+ std::ofstream out (fname_mat.c_str(), std::ios::binary);
+ CS_write (out, mat.getNRows(), mat.getRowPtrArray(), mat.getColIdxArray(), mat.getEntryArray());
+ out.close();
+}
diff --git a/SimpleTests/SolverTests/ConjugateGradientUnpreconditioned.cpp b/SimpleTests/SolverTests/ConjugateGradientUnpreconditioned.cpp
index 2413a1502def9a6f29460fd0d8bf7edd2f33532c..52f1ba239f7a459c1901ee9382cdf36ee3571b28 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) {