diff --git a/AssemblerLib/ComponentGlobalIndexDict.h b/AssemblerLib/ComponentGlobalIndexDict.h
index 59aba9e82166a69f027f80a1732d7b7187bb93d2..e5fc2de61c9d05c1af181785ba0cbf408dd4f862 100644
--- a/AssemblerLib/ComponentGlobalIndexDict.h
+++ b/AssemblerLib/ComponentGlobalIndexDict.h
@@ -22,6 +22,7 @@
#include <boost/multi_index_container.hpp>
#include "MeshLib/Location.h"
+#include "ProcessLib/NumericsConfig.h"
namespace AssemblerLib
{
@@ -38,21 +39,21 @@ struct Line
std::size_t comp_id;
// Position in global matrix or vector
- std::size_t global_index;
+ GlobalIndexType global_index;
- Line(MeshLib::Location const& l, std::size_t c, std::size_t i)
+ Line(MeshLib::Location const& l, std::size_t c, GlobalIndexType i)
: location(l), comp_id(c), global_index(i)
{}
Line(MeshLib::Location const& l, std::size_t c)
: location(l), comp_id(c),
- global_index(std::numeric_limits<std::size_t>::max())
+ global_index(std::numeric_limits<GlobalIndexType>::max())
{}
explicit Line(MeshLib::Location const& l)
: location(l),
comp_id(std::numeric_limits<std::size_t>::max()),
- global_index(std::numeric_limits<std::size_t>::max())
+ global_index(std::numeric_limits<GlobalIndexType>::max())
{}
friend std::ostream& operator<<(std::ostream& os, Line const& l)
@@ -112,7 +113,7 @@ typedef boost::multi_index::multi_index_container<
boost::multi_index::ordered_non_unique
<
boost::multi_index::tag<ByGlobalIndex>,
- boost::multi_index::member<Line, std::size_t, &Line::global_index>
+ boost::multi_index::member<Line, GlobalIndexType, &Line::global_index>
>
>
> ComponentGlobalIndexDict;
diff --git a/AssemblerLib/LocalToGlobalIndexMap.h b/AssemblerLib/LocalToGlobalIndexMap.h
index 0a60576cb75f95aeb03d2e227ce780e7e99514a6..aeb52880bdcfe5c9836cfc26989400ab2a9d0154 100644
--- a/AssemblerLib/LocalToGlobalIndexMap.h
+++ b/AssemblerLib/LocalToGlobalIndexMap.h
@@ -35,7 +35,7 @@ namespace AssemblerLib
class LocalToGlobalIndexMap
{
public:
- typedef MathLib::RowColumnIndices<std::size_t> RowColumnIndices;
+ typedef MathLib::RowColumnIndices<GlobalIndexType> RowColumnIndices;
typedef RowColumnIndices::LineIndex LineIndex;
public:
@@ -78,7 +78,7 @@ public:
LineIndex rowIndices(std::size_t const mesh_item_id) const;
LineIndex columnIndices(std::size_t const mesh_item_id) const;
- std::size_t getGlobalIndex(MeshLib::Location const& l,
+ GlobalIndexType getGlobalIndex(MeshLib::Location const& l,
std::size_t const c) const
{
return _mesh_component_map.getGlobalIndex(l, c);
diff --git a/AssemblerLib/MeshComponentMap.cpp b/AssemblerLib/MeshComponentMap.cpp
index 93a867638d2be02f3a77e5236ce1f4b7d4b99ad5..fe6cd48305141c7eebf825fc59a375f5850a46fe 100644
--- a/AssemblerLib/MeshComponentMap.cpp
+++ b/AssemblerLib/MeshComponentMap.cpp
@@ -21,13 +21,14 @@ namespace AssemblerLib
using namespace detail;
-std::size_t const MeshComponentMap::nop = std::numeric_limits<std::size_t>::max();
+GlobalIndexType const MeshComponentMap::nop =
+ std::numeric_limits<GlobalIndexType>::max();
MeshComponentMap::MeshComponentMap(
const std::vector<MeshLib::MeshSubsets*> &components, ComponentOrder order)
{
// construct dict (and here we number global_index by component type)
- std::size_t global_index = 0;
+ GlobalIndexType global_index = 0;
std::size_t comp_id = 0;
for (auto c = components.cbegin(); c != components.cend(); ++c)
{
@@ -81,9 +82,9 @@ MeshComponentMap::getSubset(std::vector<MeshLib::MeshSubsets*> const& components
return MeshComponentMap(subset_dict);
}
-void MeshComponentMap::renumberByLocation(std::size_t offset)
+void MeshComponentMap::renumberByLocation(GlobalIndexType offset)
{
- std::size_t global_index = offset;
+ GlobalIndexType global_index = offset;
auto &m = _dict.get<ByLocation>(); // view as sorted by mesh item
for (auto itr_mesh_item=m.begin(); itr_mesh_item!=m.end(); ++itr_mesh_item)
@@ -113,7 +114,7 @@ Line MeshComponentMap::getLine(Location const& l,
return *itr;
}
-std::size_t MeshComponentMap::getGlobalIndex(Location const& l,
+GlobalIndexType MeshComponentMap::getGlobalIndex(Location const& l,
std::size_t const comp_id) const
{
auto const &m = _dict.get<ByLocationAndComponent>();
@@ -121,25 +122,25 @@ std::size_t MeshComponentMap::getGlobalIndex(Location const& l,
return itr!=m.end() ? itr->global_index : nop;
}
-std::vector<std::size_t> MeshComponentMap::getGlobalIndices(const Location &l) const
+std::vector<GlobalIndexType> MeshComponentMap::getGlobalIndices(const Location &l) const
{
auto const &m = _dict.get<ByLocation>();
auto const p = m.equal_range(Line(l));
- std::vector<std::size_t> global_indices;
+ std::vector<GlobalIndexType> global_indices;
for (auto itr=p.first; itr!=p.second; ++itr)
global_indices.push_back(itr->global_index);
return global_indices;
}
template <>
-std::vector<std::size_t>
+std::vector<GlobalIndexType>
MeshComponentMap::getGlobalIndices<ComponentOrder::BY_LOCATION>(
std::vector<Location> const &ls) const
{
// Create vector of global indices sorted by location containing all
// locations given in ls parameter.
- std::vector<std::size_t> global_indices;
+ std::vector<GlobalIndexType> global_indices;
global_indices.reserve(ls.size());
auto const &m = _dict.get<ByLocation>();
@@ -154,12 +155,12 @@ MeshComponentMap::getGlobalIndices<ComponentOrder::BY_LOCATION>(
}
template <>
-std::vector<std::size_t>
+std::vector<GlobalIndexType>
MeshComponentMap::getGlobalIndices<ComponentOrder::BY_COMPONENT>(
std::vector<Location> const &ls) const
{
// vector of (Component, global Index) pairs.
- typedef std::pair<std::size_t, std::size_t> CIPair;
+ typedef std::pair<std::size_t, GlobalIndexType> CIPair;
std::vector<CIPair> pairs;
pairs.reserve(ls.size());
@@ -181,7 +182,7 @@ MeshComponentMap::getGlobalIndices<ComponentOrder::BY_COMPONENT>(
if (!std::is_sorted(pairs.begin(), pairs.end(), CIPairLess))
std::stable_sort(pairs.begin(), pairs.end(), CIPairLess);
- std::vector<std::size_t> global_indices;
+ std::vector<GlobalIndexType> global_indices;
global_indices.reserve(pairs.size());
for (auto p = pairs.cbegin(); p != pairs.cend(); ++p)
global_indices.push_back(p->second);
diff --git a/AssemblerLib/MeshComponentMap.h b/AssemblerLib/MeshComponentMap.h
index a4a3dbfb81ba14e55fb7700fbd7a9c52c45b757e..1ed5b0c13740b17569338a438c2aae333a166cd1 100644
--- a/AssemblerLib/MeshComponentMap.h
+++ b/AssemblerLib/MeshComponentMap.h
@@ -70,7 +70,7 @@ public:
/// | Location | ComponentID | GlobalIndex |
/// | -------- | ----------- | ----------- |
/// | l | comp_id | gi |
- std::size_t getGlobalIndex(Location const &l, std::size_t const comp_id) const;
+ GlobalIndexType getGlobalIndex(Location const &l, std::size_t const comp_id) const;
/// Global indices for all components at the given location \c l.
///
@@ -82,7 +82,7 @@ public:
/// | l | comp_id_1 | gi23 |
/// | ... | ... | ... |
/// | l | comp_id_k | gi45 |
- std::vector<std::size_t> getGlobalIndices(const Location &l) const;
+ std::vector<GlobalIndexType> getGlobalIndices(const Location &l) const;
/// Global indices for all components at all given locations \c ls ordered
/// as required by the template parameter ORDER.
@@ -113,11 +113,11 @@ public:
/// | ... | ... | ... |
/// | l_n | comp_id_n | gi89 |
template <ComponentOrder ORDER>
- std::vector<std::size_t> getGlobalIndices(const std::vector<Location> &ls) const;
+ std::vector<GlobalIndexType> getGlobalIndices(const std::vector<Location> &ls) const;
/// A value returned if no global index was found for the requested
/// location/component. The value is implementation dependent.
- static std::size_t const nop;
+ static GlobalIndexType const nop;
#ifndef NDEBUG
const detail::ComponentGlobalIndexDict& getDictionary() const
@@ -146,7 +146,7 @@ private:
/// \return a copy of the line.
detail::Line getLine(Location const& l, std::size_t const component_id) const;
- void renumberByLocation(std::size_t offset=0);
+ void renumberByLocation(GlobalIndexType offset=0);
private:
detail::ComponentGlobalIndexDict _dict;
diff --git a/AssemblerLib/VectorMatrixBuilder.h b/AssemblerLib/VectorMatrixBuilder.h
index dc0c747225dcda9a135dc328af7397710061e8a8..4ead7cac2866dbe02ec8eb10d7120fc2d8a914fd 100644
--- a/AssemblerLib/VectorMatrixBuilder.h
+++ b/AssemblerLib/VectorMatrixBuilder.h
@@ -13,8 +13,6 @@
#ifndef ASSEMBLERLIB_VECTORMATRIXBUILDER_H_
#define ASSEMBLERLIB_VECTORMATRIXBUILDER_H_
-#include "AssemblerLib/MeshComponentMap.h"
-
namespace AssemblerLib
{
diff --git a/MathLib/LinAlg/Dense/DenseVector.h b/MathLib/LinAlg/Dense/DenseVector.h
index e924b15564b0a9323188f0545d90b317e3c8a912..e630ddf98073b316bbc2b623e8543f19325b461b 100644
--- a/MathLib/LinAlg/Dense/DenseVector.h
+++ b/MathLib/LinAlg/Dense/DenseVector.h
@@ -31,6 +31,7 @@ class DenseVector : public std::valarray<T>
{
public:
typedef T FP_T;
+ using IndexType = std::size_t; // The type of valarray indices.
public:
using std::valarray<T>::operator=;
diff --git a/MathLib/LinAlg/Dense/GlobalDenseMatrix.h b/MathLib/LinAlg/Dense/GlobalDenseMatrix.h
index 46f13c7e50d30786d41315d848ead9b2c88e4d94..60e97121c7dee47eea71fb5ba7950f1b70dafc74 100644
--- a/MathLib/LinAlg/Dense/GlobalDenseMatrix.h
+++ b/MathLib/LinAlg/Dense/GlobalDenseMatrix.h
@@ -30,6 +30,7 @@ class GlobalDenseMatrix: public DenseMatrix<FP_TYPE, IDX_TYPE>
{
public:
typedef FP_TYPE FP_T;
+ using IndexType = IDX_TYPE;
public:
/// Dense square matrix constructor.
diff --git a/MathLib/LinAlg/Eigen/EigenMatrix.h b/MathLib/LinAlg/Eigen/EigenMatrix.h
index ad5758afc24ba7d038225ba9203c69b3510d89fa..4cbd970c03d1b4800ffd2ed887dd4e1fddc0d20c 100644
--- a/MathLib/LinAlg/Eigen/EigenMatrix.h
+++ b/MathLib/LinAlg/Eigen/EigenMatrix.h
@@ -10,7 +10,6 @@
#ifndef EIGENMATRIX_H_
#define EIGENMATRIX_H_
-#include <cassert>
#ifndef NDEBUG
#include <fstream>
#include <string>
@@ -33,6 +32,7 @@ class EigenMatrix final
{
public:
using RawMatrixType = Eigen::SparseMatrix<double, Eigen::RowMajor>;
+ using IndexType = RawMatrixType::Index;
/**
* constructor
@@ -61,14 +61,14 @@ public:
void setZero()
{
auto const N = _mat.nonZeros();
- for (std::remove_const<decltype(N)>::type i=0; i<N; i++)
+ for (auto i = decltype(N){0}; i<N; i++)
_mat.valuePtr()[i] = 0;
// don't use _mat.setZero(). it makes a matrix uncompressed
}
/// set a value to the given entry. If the entry doesn't exist, this class
/// dynamically allocates it.
- int setValue(std::size_t row, std::size_t col, double val)
+ int setValue(IndexType row, IndexType col, double val)
{
_mat.coeffRef(row, col) = val;
return 0;
@@ -76,9 +76,8 @@ public:
/// add a value to the given entry. If the entry doesn't exist, the value is
/// inserted.
- int add(std::size_t row, std::size_t col, double val)
+ int add(IndexType row, IndexType col, double val)
{
- assert(row < getNRows() && col < getNCols());
_mat.coeffRef(row, col) += val;
return 0;
}
@@ -86,7 +85,7 @@ public:
/// Add sub-matrix at positions \c row_pos and same column positions as the
/// given row positions. If the entry doesn't exist, the value is inserted.
template<class T_DENSE_MATRIX>
- void add(std::vector<std::size_t> const& row_pos,
+ void add(std::vector<IndexType> const& row_pos,
const T_DENSE_MATRIX &sub_matrix,
double fkt = 1.0)
{
@@ -96,7 +95,7 @@ public:
/// Add sub-matrix at positions given by \c indices. If the entry doesn't exist,
/// this class inserts the value.
template<class T_DENSE_MATRIX>
- void add(RowColumnIndices<std::size_t> const& indices,
+ void add(RowColumnIndices<IndexType> const& indices,
const T_DENSE_MATRIX &sub_matrix,
double fkt = 1.0)
{
@@ -112,19 +111,18 @@ public:
/// @param sub_matrix a sub-matrix to be added
/// @param fkt a scaling factor applied to all entries in the sub-matrix
template <class T_DENSE_MATRIX>
- void add(std::vector<std::size_t> const& row_pos,
- std::vector<std::size_t> const& col_pos, const T_DENSE_MATRIX &sub_matrix,
+ void add(std::vector<IndexType> const& row_pos,
+ std::vector<IndexType> const& col_pos, const T_DENSE_MATRIX &sub_matrix,
double fkt = 1.0);
/// get value. This function returns zero if the element doesn't exist.
- double get(std::size_t row, std::size_t col) const
+ double get(IndexType row, IndexType col) const
{
- assert(row < getNRows() && col < getNCols());
return _mat.coeff(row, col);
}
/// get value. This function returns zero if the element doesn't exist.
- double operator() (std::size_t row, std::size_t col) const
+ double operator() (IndexType row, IndexType col) const
{
return get(row, col);
}
@@ -170,18 +168,17 @@ protected:
RawMatrixType _mat;
};
-
-template<class T_DENSE_MATRIX>
-void
-EigenMatrix::add(std::vector<std::size_t> const& row_pos, std::vector<std::size_t> const& col_pos,
- const T_DENSE_MATRIX &sub_matrix, double fkt)
+template <class T_DENSE_MATRIX>
+void EigenMatrix::add(std::vector<IndexType> const& row_pos,
+ std::vector<IndexType> const& col_pos,
+ const T_DENSE_MATRIX& sub_matrix, double fkt)
{
- const std::size_t n_rows = row_pos.size();
- const std::size_t n_cols = col_pos.size();
- for (std::size_t i = 0; i < n_rows; i++) {
- const std::size_t row = row_pos[i];
- for (std::size_t j = 0; j < n_cols; j++) {
- const std::size_t col = col_pos[j];
+ auto const n_rows = row_pos.size();
+ auto const n_cols = col_pos.size();
+ for (auto i = decltype(n_rows){0}; i < n_rows; i++) {
+ auto const row = row_pos[i];
+ for (auto j = decltype(n_cols){0}; j < n_cols; j++) {
+ auto const col = col_pos[j];
add(row, col, fkt * sub_matrix(i, j));
}
}
@@ -191,4 +188,3 @@ EigenMatrix::add(std::vector<std::size_t> const& row_pos, std::vector<std::size_
} // end namespace MathLib
#endif
-
diff --git a/MathLib/LinAlg/Eigen/EigenVector.h b/MathLib/LinAlg/Eigen/EigenVector.h
index 0c27be964077d94202afead0c832cbe0d472bfac..b723b15beae0d06e28ea4ee3d0213f2e25ae3cae 100644
--- a/MathLib/LinAlg/Eigen/EigenVector.h
+++ b/MathLib/LinAlg/Eigen/EigenVector.h
@@ -17,6 +17,7 @@
#endif
#include <Eigen/Eigen>
+#include <Eigen/Sparse>
namespace MathLib
{
@@ -27,6 +28,11 @@ class EigenVector final
public:
using RawVectorType = Eigen::VectorXd;
+ // The Index type of the Eigen::VectorXd class differs from the
+ // Eigen::SparseMatrix<double> index type. Maybe an Eigen::SparseVector is a
+ // more appropriate RawVectorType for the global vectors.
+ using IndexType = Eigen::SparseMatrix<double>::Index;
+
/// Constructor for initialization of the number of rows
/// @param length number of rows
explicit EigenVector(std::size_t length) : _vec(length) {}
@@ -50,30 +56,30 @@ public:
EigenVector& operator*= (double v) { _vec *= v; return *this; }
/// access entry
- double const & operator[] (std::size_t rowId) const { return _vec[rowId]; }
- double& operator[] (std::size_t rowId) { return _vec[rowId]; }
+ double const & operator[] (IndexType rowId) const { return _vec[rowId]; }
+ double& operator[] (IndexType rowId) { return _vec[rowId]; }
/// get entry
- double get(std::size_t rowId) const
+ double get(IndexType rowId) const
{
return _vec[rowId];
}
/// set entry
- void set(std::size_t rowId, double v)
+ void set(IndexType rowId, double v)
{
_vec[rowId] = v;
}
/// add entry
- void add(std::size_t rowId, double v)
+ void add(IndexType rowId, double v)
{
_vec[rowId] += v;
}
/// add entries
template<class T_SUBVEC>
- void add(const std::vector<std::size_t> &pos, const T_SUBVEC &sub_vec)
+ void add(const std::vector<IndexType> &pos, const T_SUBVEC &sub_vec)
{
for (std::size_t i=0; i<pos.size(); ++i) {
this->add(pos[i], sub_vec[i]);
diff --git a/MathLib/LinAlg/Lis/LisMatrix.cpp b/MathLib/LinAlg/Lis/LisMatrix.cpp
index 6b3e420fd2eefd8c7671d91b31542a94d65cc3f8..db5a238a81902d1322617a282a6ab15babcc13e8 100644
--- a/MathLib/LinAlg/Lis/LisMatrix.cpp
+++ b/MathLib/LinAlg/Lis/LisMatrix.cpp
@@ -37,8 +37,12 @@ LisMatrix::LisMatrix(std::size_t n_rows, LisOption::MatrixType mat_type)
checkLisError(ierr);
}
-LisMatrix::LisMatrix(std::size_t n_rows, int nnz, int* row_ptr, int* col_idx, double* data)
-: _n_rows(n_rows), _mat_type(LisOption::MatrixType::CRS), _is_assembled(false), _use_external_arrays(true)
+LisMatrix::LisMatrix(std::size_t n_rows, int nnz, IndexType *row_ptr,
+ IndexType *col_idx, double *data)
+ : _n_rows(n_rows),
+ _mat_type(LisOption::MatrixType::CRS),
+ _is_assembled(false),
+ _use_external_arrays(true)
{
int ierr = lis_matrix_create(0, &_AA);
checkLisError(ierr);
@@ -81,7 +85,7 @@ void LisMatrix::setZero()
_is_assembled = false;
}
-int LisMatrix::setValue(std::size_t rowId, std::size_t colId, double v)
+int LisMatrix::setValue(IndexType rowId, IndexType colId, double v)
{
lis_matrix_set_value(LIS_INS_VALUE, rowId, colId, v, _AA);
if (rowId==colId)
@@ -90,7 +94,7 @@ int LisMatrix::setValue(std::size_t rowId, std::size_t colId, double v)
return 0;
}
-int LisMatrix::add(std::size_t rowId, std::size_t colId, double v)
+int LisMatrix::add(IndexType rowId, IndexType colId, double v)
{
lis_matrix_set_value(LIS_ADD_VALUE, rowId, colId, v, _AA);
if (rowId==colId)
diff --git a/MathLib/LinAlg/Lis/LisMatrix.h b/MathLib/LinAlg/Lis/LisMatrix.h
index 4d9d725c6943a30a2e6768880e88a576e82fa4a9..1756accbde367ad0ce5d07c490cf7a5b934f5c2e 100644
--- a/MathLib/LinAlg/Lis/LisMatrix.h
+++ b/MathLib/LinAlg/Lis/LisMatrix.h
@@ -44,6 +44,8 @@ struct SetMatrixSparsity<LisMatrix, SPARSITY_PATTERN>;
*/
class LisMatrix
{
+public:
+ using IndexType = LIS_INT;
public:
/**
* constructor
@@ -62,7 +64,8 @@ public:
* @param col_idx array of column indexes
* @param data the non-zero entry values
*/
- LisMatrix(std::size_t n_rows, int nonzero, int* row_ptr, int* col_idx, double* data);
+ LisMatrix(std::size_t n_rows, int nonzero, IndexType* row_ptr, IndexType* col_idx,
+ double* data);
/**
*
@@ -85,10 +88,10 @@ public:
void setZero();
/// set entry
- int setValue(std::size_t rowId, std::size_t colId, double v);
+ int setValue(IndexType rowId, IndexType colId, double v);
/// add value
- int add(std::size_t rowId, std::size_t colId, double v);
+ int add(IndexType rowId, IndexType colId, double v);
/// printout this equation for debugging
void write(const std::string &filename) const;
@@ -105,7 +108,7 @@ public:
/// Add sub-matrix at positions \c row_pos and same column positions as the
/// given row positions.
template<class T_DENSE_MATRIX>
- void add(std::vector<std::size_t> const& row_pos,
+ void add(std::vector<IndexType> const& row_pos,
const T_DENSE_MATRIX &sub_matrix,
double fkt = 1.0)
{
@@ -114,18 +117,17 @@ public:
/// Add sub-matrix at positions given by \c indices.
template<class T_DENSE_MATRIX>
- void add(RowColumnIndices<std::size_t> const& indices,
+ void add(RowColumnIndices<IndexType> const& indices,
const T_DENSE_MATRIX &sub_matrix,
double fkt = 1.0)
{
this->add(indices.rows, indices.columns, sub_matrix, fkt);
}
- ///
template <class T_DENSE_MATRIX>
- void add(std::vector<std::size_t> const& row_pos,
- std::vector<std::size_t> const& col_pos, const T_DENSE_MATRIX &sub_matrix,
- double fkt = 1.0);
+ void add(std::vector<IndexType> const& row_pos,
+ std::vector<IndexType> const& col_pos,
+ const T_DENSE_MATRIX& sub_matrix, double fkt = 1.0);
/// get this matrix type
LisOption::MatrixType getMatrixType() const { return _mat_type; }
@@ -139,8 +141,8 @@ private:
LIS_MATRIX _AA;
LIS_VECTOR _diag;
bool _is_assembled;
- LIS_INT _is; ///< location where the partial matrix _AA starts in global matrix.
- LIS_INT _ie; ///< location where the partial matrix _AA ends in global matrix.
+ IndexType _is; ///< location where the partial matrix _AA starts in global matrix.
+ IndexType _ie; ///< location where the partial matrix _AA ends in global matrix.
bool _use_external_arrays;
// friend function
@@ -150,17 +152,19 @@ private:
friend struct SetMatrixSparsity;
};
-template<class T_DENSE_MATRIX>
-void
-LisMatrix::add(std::vector<std::size_t> const& row_pos, std::vector<std::size_t> const& col_pos,
- const T_DENSE_MATRIX &sub_matrix, double fkt)
+template <class T_DENSE_MATRIX>
+void LisMatrix::add(std::vector<IndexType> const& row_pos,
+ std::vector<IndexType> const& col_pos,
+ const T_DENSE_MATRIX& sub_matrix, double fkt)
{
- const std::size_t n_rows = row_pos.size();
- const std::size_t n_cols = col_pos.size();
- for (std::size_t i = 0; i < n_rows; i++) {
- const std::size_t row = row_pos[i];
- for (std::size_t j = 0; j < n_cols; j++) {
- const std::size_t col = col_pos[j];
+ auto const n_rows = row_pos.size();
+ auto const n_cols = col_pos.size();
+ for (auto i = decltype(n_rows){0}; i < n_rows; i++)
+ {
+ auto const row = row_pos[i];
+ for (auto j = decltype(n_cols){0}; j < n_cols; j++)
+ {
+ auto const col = col_pos[j];
add(row, col, fkt * sub_matrix(i, j));
}
}
@@ -176,12 +180,12 @@ struct SetMatrixSparsity<LisMatrix, SPARSITY_PATTERN>
void operator()(LisMatrix &matrix, SPARSITY_PATTERN const& sparsity_pattern)
{
- std::size_t n_rows = matrix.getNRows();
- std::vector<int> row_sizes;
+ auto const n_rows = matrix.getNRows();
+ std::vector<LisMatrix::IndexType> row_sizes;
row_sizes.reserve(n_rows);
// LIS needs 1 more entry, otherewise it starts reallocating arrays.
- for (std::size_t i = 0; i < n_rows; i++)
+ for (auto i = decltype(n_rows){0}; i < n_rows; i++)
row_sizes.push_back(sparsity_pattern.getNodeDegree(i) + 1);
int ierr = lis_matrix_malloc(matrix._AA, 0, row_sizes.data());
@@ -193,4 +197,3 @@ void operator()(LisMatrix &matrix, SPARSITY_PATTERN const& sparsity_pattern)
} // MathLib
#endif //LISMATRIX_H_
-
diff --git a/MathLib/LinAlg/Lis/LisVector.cpp b/MathLib/LinAlg/Lis/LisVector.cpp
index 72af5d6447e2b0d95d91e43c4f0aa055b5e337eb..924471175d56f25e8ac73bc4971bad711d46d876 100644
--- a/MathLib/LinAlg/Lis/LisVector.cpp
+++ b/MathLib/LinAlg/Lis/LisVector.cpp
@@ -12,27 +12,28 @@
*
*/
+#include <cassert>
+
#include "LisVector.h"
#include "LisCheck.h"
namespace MathLib
{
-
LisVector::LisVector(std::size_t length)
{
- lis_vector_create(0, &_vec);
- lis_vector_set_size(_vec, 0, length);
+ lis_vector_create(0, &_vec);
+ lis_vector_set_size(_vec, 0, length);
}
LisVector::LisVector(std::size_t length, double* data)
{
- lis_vector_create(0, &_vec);
- lis_vector_set_size(_vec, 0, length);
- for (std::size_t i=0; i<length; i++)
- lis_vector_set_value(LIS_INS_VALUE, i, data[i], _vec);
+ lis_vector_create(0, &_vec);
+ lis_vector_set_size(_vec, 0, length);
+ for (std::size_t i = 0; i < length; i++)
+ lis_vector_set_value(LIS_INS_VALUE, i, data[i], _vec);
}
-LisVector::LisVector(LisVector const &src)
+LisVector::LisVector(LisVector const& src)
{
lis_vector_duplicate(src._vec, &_vec);
lis_vector_copy(src._vec, _vec);
@@ -40,45 +41,44 @@ LisVector::LisVector(LisVector const &src)
LisVector::~LisVector()
{
- lis_vector_destroy(_vec);
+ lis_vector_destroy(_vec);
}
-LisVector& LisVector::operator= (const LisVector &src)
+LisVector& LisVector::operator=(const LisVector& src)
{
lis_vector_copy(src._vec, _vec);
- return *this;
+ return *this;
}
-void LisVector::operator+= (const LisVector& v)
+void LisVector::operator+=(const LisVector& v)
{
lis_vector_axpy(1.0, v._vec, _vec);
}
-void LisVector::operator-= (const LisVector& v)
+void LisVector::operator-=(const LisVector& v)
{
lis_vector_axpy(-1.0, v._vec, _vec);
}
-LisVector& LisVector::operator= (double v)
+LisVector& LisVector::operator=(double v)
{
- lis_vector_set_all(v, _vec);
- return *this;
+ lis_vector_set_all(v, _vec);
+ return *this;
}
std::size_t LisVector::size() const
{
- LIS_INT dummy;
- LIS_INT size;
+ IndexType dummy;
+ IndexType size;
int const ierr = lis_vector_get_size(_vec, &dummy, &size);
checkLisError(ierr);
+ assert(size >= 0); // For safe implicit conversion to std::size_t.
return size;
}
-void LisVector::write (const std::string &filename) const
+void LisVector::write(const std::string& filename) const
{
lis_output_vector(_vec, LIS_FMT_PLAIN, const_cast<char*>(filename.c_str()));
}
-
-} // MathLib
-
+} // MathLib
diff --git a/MathLib/LinAlg/Lis/LisVector.h b/MathLib/LinAlg/Lis/LisVector.h
index 68ef9637e6b609f768a13f7eaae45b79c36d3050..e2b3f82ef2cfd2a235ea8b67821f763b55f4a4a4 100644
--- a/MathLib/LinAlg/Lis/LisVector.h
+++ b/MathLib/LinAlg/Lis/LisVector.h
@@ -22,97 +22,95 @@
namespace MathLib
{
-
/**
* \brief Lis vector wrapper class
*/
class LisVector
{
+public:
+ using IndexType = LIS_INT;
public:
/**
* Constructor for initialization of the number of rows
* @param length number of rows
*/
- explicit LisVector(std::size_t length);
-
- /**
- * Constructor using the given raw data
- * @param length the length of the vector
- * @param data the raw data
- */
- LisVector(std::size_t length, double* data);
-
- /// copy constructor
- LisVector(LisVector const &src);
-
- /**
- *
- */
- virtual ~LisVector();
-
- /// return a vector length
- std::size_t size() const;
-
- /// return a start index of the active data range
- std::size_t getRangeBegin() const { return 0;}
-
- /// return an end index of the active data range
- std::size_t getRangeEnd() const { return this->size(); }
-
- /// set all values in this vector
- LisVector& operator= (double v);
-
- /// access entry
- double operator[] (std::size_t rowId) const { return get(rowId); }
-
- /// get entry
- double get(std::size_t rowId) const
- {
- double v = .0;
- lis_vector_get_value(_vec, rowId, &v);
- return v;
- }
-
- /// set entry
- void set(std::size_t rowId, double v)
- {
- lis_vector_set_value(LIS_INS_VALUE, rowId, v, _vec);
- }
-
- /// add entry
- void add(std::size_t rowId, double v)
- {
- lis_vector_set_value(LIS_ADD_VALUE, rowId, v, _vec);
- }
-
- /// printout this equation for debugging
- void write (const std::string &filename) const;
-
- /// return a raw Lis vector object
- LIS_VECTOR& getRawVector() {return _vec; }
-
- /// vector operation: set data
- LisVector& operator= (const LisVector &src);
-
- /// vector operation: add
- void operator+= (const LisVector& v);
-
- /// vector operation: subtract
- void operator-= (const LisVector& v);
-
- ///
- template<class T_SUBVEC>
- void add(const std::vector<std::size_t> &pos, const T_SUBVEC &sub_vec)
- {
- for (std::size_t i=0; i<pos.size(); ++i) {
- this->add(pos[i], sub_vec[i]);
- }
- }
+ explicit LisVector(std::size_t length);
+
+ /**
+ * Constructor using the given raw data
+ * @param length the length of the vector
+ * @param data the raw data
+ */
+ LisVector(std::size_t length, double* data);
+
+ /// copy constructor
+ LisVector(LisVector const& src);
+
+ /**
+ *
+ */
+ virtual ~LisVector();
+
+ /// return a vector length
+ std::size_t size() const;
+
+ /// return a start index of the active data range
+ std::size_t getRangeBegin() const { return 0; }
+ /// return an end index of the active data range
+ std::size_t getRangeEnd() const { return this->size(); }
+ /// set all values in this vector
+ LisVector& operator=(double v);
+
+ /// access entry
+ double operator[](IndexType rowId) const { return get(rowId); }
+ /// get entry
+ double get(IndexType rowId) const
+ {
+ double v = .0;
+ lis_vector_get_value(_vec, rowId, &v);
+ return v;
+ }
+
+ /// set entry
+ void set(IndexType rowId, double v)
+ {
+ lis_vector_set_value(LIS_INS_VALUE, rowId, v, _vec);
+ }
+
+ /// add entry
+ void add(IndexType rowId, double v)
+ {
+ lis_vector_set_value(LIS_ADD_VALUE, rowId, v, _vec);
+ }
+
+ /// printout this equation for debugging
+ void write(const std::string& filename) const;
+
+ /// return a raw Lis vector object
+ LIS_VECTOR& getRawVector() { return _vec; }
+ /// vector operation: set data
+ LisVector& operator=(const LisVector& src);
+
+ /// vector operation: add
+ void operator+=(const LisVector& v);
+
+ /// vector operation: subtract
+ void operator-=(const LisVector& v);
+
+ ///
+ template <class T_SUBVEC>
+ void add(const std::vector<IndexType>& pos, const T_SUBVEC& sub_vec)
+ {
+ for (std::size_t i = 0; i < pos.size(); ++i)
+ {
+ this->add(pos[i], sub_vec[i]);
+ }
+ }
+
private:
- LIS_VECTOR _vec;
+ LIS_VECTOR _vec;
};
-} // MathLib
-
-#endif //LISVECTOR_H_
+} // MathLib
+#endif // LISVECTOR_H_
diff --git a/MathLib/LinAlg/PETSc/PETScMatrix.h b/MathLib/LinAlg/PETSc/PETScMatrix.h
index 43706b873052adec1d3ed7b26dc27ff0e26c582e..4f830e865bb86f78d8d3472689a2c68010ca03e2 100644
--- a/MathLib/LinAlg/PETSc/PETScMatrix.h
+++ b/MathLib/LinAlg/PETSc/PETScMatrix.h
@@ -33,6 +33,9 @@ namespace MathLib
*/
class PETScMatrix
{
+ public:
+ using IndexType = PetscInt;
+
public:
/*!
\brief Constructor for a square matrix partitioning with more options
diff --git a/MathLib/LinAlg/PETSc/PETScVector.h b/MathLib/LinAlg/PETSc/PETScVector.h
index b940d0d861b8df3c7ad67c7f29f502c8dec629f4..ee7d4482829bb5304194aa2d89f53b2e59387542 100644
--- a/MathLib/LinAlg/PETSc/PETScVector.h
+++ b/MathLib/LinAlg/PETSc/PETScVector.h
@@ -34,6 +34,9 @@ namespace MathLib
*/
class PETScVector
{
+ public:
+ using IndexType = PetscInt;
+
public:
/*!
diff --git a/ProcessLib/GroundwaterFlowProcess.h b/ProcessLib/GroundwaterFlowProcess.h
index 649ab89c5e8e418f178b75d8c9af967d4903d749..42d6fed50d6fb283f5dd1a473642ebf9c999c5ad 100644
--- a/ProcessLib/GroundwaterFlowProcess.h
+++ b/ProcessLib/GroundwaterFlowProcess.h
@@ -237,7 +237,7 @@ public:
{
MeshLib::Location const l(_mesh.getID(),
MeshLib::MeshItemType::Node, i);
- std::size_t const global_index =
+ auto const global_index =
_local_to_global_index_map->getGlobalIndex(
l, 0); // 0 is the component id.
_x->set(global_index,
diff --git a/ProcessLib/NumericsConfig.h b/ProcessLib/NumericsConfig.h
index 3bc349a638077cfa004b3864ce113c1c56ff97e5..2b4fcaf86eeb23baa27c8d5beb2428d8a7013498 100644
--- a/ProcessLib/NumericsConfig.h
+++ b/ProcessLib/NumericsConfig.h
@@ -103,5 +103,18 @@ using GlobalSetupType =
// Check the configuration
//
static_assert(std::is_class<GlobalSetupType>::value,
- "GlobalSetupType was not defined.");
+ "GlobalSetupType was not defined.");
+static_assert(std::is_integral<detail::GlobalMatrixType::IndexType>::value,
+ "The index type for global matrices is not an integral type.");
+static_assert(std::is_integral<detail::GlobalVectorType::IndexType>::value,
+ "The index type for global vectors is not an integral type.");
+static_assert(std::is_same<detail::GlobalMatrixType::IndexType,
+ detail::GlobalVectorType::IndexType>::value,
+ "The global matrix and vector index types do not match.");
+// Both types are integral types and equal, define a single GlobalIndexType.
+
+/// A type used for indexing of global vectors and matrices. It is equal to the
+/// GlobalMatrixType::IndexType and the GlobalVectorType::IndexType.
+using GlobalIndexType = detail::GlobalMatrixType::IndexType;
+
#endif // APPLICATIONS_NUMERICSCONFIG_H_
diff --git a/Tests/MathLib/TestGlobalMatrixInterface.cpp b/Tests/MathLib/TestGlobalMatrixInterface.cpp
index 430704887e4f52983060d20daaf5ca1ccf9445e3..a1741f33a10643da6495a42e35b4aed3a99bd89c 100644
--- a/Tests/MathLib/TestGlobalMatrixInterface.cpp
+++ b/Tests/MathLib/TestGlobalMatrixInterface.cpp
@@ -15,21 +15,20 @@
#include <gtest/gtest.h>
-#include "MathLib/LinAlg/Dense/DenseMatrix.h"
-#include "MathLib/LinAlg/Dense/GlobalDenseMatrix.h"
-#include "MathLib/LinAlg/FinalizeMatrixAssembly.h"
-
-#ifdef OGS_USE_EIGEN
+#if defined(USE_LIS)
+#include "MathLib/LinAlg/Lis/LisMatrix.h"
+#elif defined(USE_PETSC)
+#include "MathLib/LinAlg/PETSc/PETScMatrix.h"
+#elif defined(OGS_USE_EIGEN)
#include "MathLib/LinAlg/Eigen/EigenMatrix.h"
+#else
+#include "MathLib/LinAlg/Dense/GlobalDenseMatrix.h"
#endif
-#ifdef USE_LIS
-#include "MathLib/LinAlg/Lis/LisMatrix.h"
-#endif
+#include "MathLib/LinAlg/Dense/DenseMatrix.h"
+#include "MathLib/LinAlg/FinalizeMatrixAssembly.h"
-#ifdef USE_PETSC
-#include "MathLib/LinAlg/PETSc/PETScMatrix.h"
-#endif
+#include "ProcessLib/NumericsConfig.h"
namespace
{
@@ -47,7 +46,7 @@ void checkGlobalMatrixInterface(T_MATRIX &m)
m.setZero();
MathLib::DenseMatrix<double> local_m(2, 2, 1.0);
- std::vector<std::size_t> vec_pos(2);
+ std::vector<GlobalIndexType> vec_pos(2);
vec_pos[0] = 1;
vec_pos[1] = 3;
m.add(vec_pos, vec_pos, local_m);
@@ -84,8 +83,8 @@ void checkGlobalMatrixInterfaceMPI(T_MATRIX &m, T_VECTOR &v)
loc_m(1, 0) = 3.;
loc_m(1, 1) = 4.;
- std::vector<int> row_pos(2);
- std::vector<int> col_pos(2);
+ std::vector<GlobalIndexType> row_pos(2);
+ std::vector<GlobalIndexType> col_pos(2);
row_pos[0] = 2 * mrank;
row_pos[1] = 2 * mrank + 1;
col_pos[0] = row_pos[0];
@@ -143,8 +142,8 @@ void checkGlobalRectangularMatrixInterfaceMPI(T_MATRIX &m, T_VECTOR &v)
loc_m(1, 1) = 2.;
loc_m(1, 2) = 3.;
- std::vector<int> row_pos(2);
- std::vector<int> col_pos(3);
+ std::vector<GlobalIndexType> row_pos(2);
+ std::vector<GlobalIndexType> col_pos(3);
row_pos[0] = 2 * mrank;
row_pos[1] = 2 * mrank + 1;
col_pos[0] = 3 * mrank;
@@ -167,29 +166,13 @@ void checkGlobalRectangularMatrixInterfaceMPI(T_MATRIX &m, T_VECTOR &v)
} // end namespace
-TEST(Math, CheckInterface_GlobalDenseMatrix)
-{
- MathLib::GlobalDenseMatrix<double> m(10, 10);
- checkGlobalMatrixInterface(m);
-}
-
-#ifdef OGS_USE_EIGEN
-TEST(Math, CheckInterface_EigenMatrix)
-{
- MathLib::EigenMatrix m(10);
- checkGlobalMatrixInterface(m);
-}
-#endif
-
-#ifdef USE_LIS
+#if defined(USE_LIS)
TEST(Math, CheckInterface_LisMatrix)
{
MathLib::LisMatrix m(10);
checkGlobalMatrixInterface(m);
}
-#endif
-
-#ifdef USE_PETSC // or MPI
+#elif defined(USE_PETSC)
TEST(MPITest_Math, CheckInterface_PETScMatrix_Local_Size)
{
MathLib::PETScMatrixOption opt;
@@ -217,7 +200,6 @@ TEST(MPITest_Math, CheckInterface_PETScMatrix_Global_Size)
checkGlobalMatrixInterfaceMPI(A, x);
}
-// Test rectangular matrix
TEST(MPITest_Math, CheckInterface_PETSc_Rectangular_Matrix_Local_Size)
{
MathLib::PETScMatrixOption opt;
@@ -244,6 +226,16 @@ TEST(MPITest_Math, CheckInterface_PETSc_Rectangular_Matrix_Global_Size)
checkGlobalRectangularMatrixInterfaceMPI(A, x);
}
-
-#endif // end of: ifdef USE_PETSC // or MPI
-
+#elif defined(OGS_USE_EIGEN)
+TEST(Math, CheckInterface_EigenMatrix)
+{
+ MathLib::EigenMatrix m(10);
+ checkGlobalMatrixInterface(m);
+}
+#else
+TEST(Math, CheckInterface_GlobalDenseMatrix)
+{
+ MathLib::GlobalDenseMatrix<double> m(10, 10);
+ checkGlobalMatrixInterface(m);
+}
+#endif
diff --git a/Tests/MathLib/TestGlobalVectorInterface.cpp b/Tests/MathLib/TestGlobalVectorInterface.cpp
index caa06479069ac8c2065978600a09f4204fbd73fb..59473457786b2cb79891e9b4169bebfd26e6c61c 100644
--- a/Tests/MathLib/TestGlobalVectorInterface.cpp
+++ b/Tests/MathLib/TestGlobalVectorInterface.cpp
@@ -16,21 +16,19 @@
#include <gtest/gtest.h>
#include "../TestTools.h"
-#include "MathLib/LinAlg/Dense/DenseVector.h"
-#include "MathLib/LinAlg/FinalizeVectorAssembly.h"
-
-#ifdef OGS_USE_EIGEN
-#include "MathLib/LinAlg/Eigen/EigenVector.h"
-#endif
-
-#ifdef USE_LIS
+#if defined(USE_LIS)
#include "MathLib/LinAlg/Lis/LisVector.h"
-#endif
-
-#ifdef USE_PETSC
+#elif defined(USE_PETSC)
#include "MathLib/LinAlg/PETSc/PETScVector.h"
+#elif defined(OGS_USE_EIGEN)
+#include "MathLib/LinAlg/Eigen/EigenVector.h"
+#else
+#include "MathLib/LinAlg/Dense/DenseVector.h"
#endif
+#include "MathLib/LinAlg/FinalizeVectorAssembly.h"
+#include "ProcessLib/NumericsConfig.h"
+
namespace
{
@@ -66,7 +64,7 @@ void checkGlobalVectorInterface()
ASSERT_EQ(2.0, y.get(0));
std::vector<double> local_vec(2, 1.0);
- std::vector<std::size_t> vec_pos(2);
+ std::vector<GlobalIndexType> vec_pos(2);
vec_pos[0] = 0;
vec_pos[1] = 3;
y.add(vec_pos, local_vec);
@@ -114,7 +112,7 @@ void checkGlobalVectorInterfaceMPI()
ASSERT_EQ(40., y.getNorm());
std::vector<double> local_vec(2, 10.0);
- std::vector<int> vec_pos(2);
+ std::vector<GlobalIndexType> vec_pos(2);
vec_pos[0] = r0; // any index in [0,15]
vec_pos[1] = r0+1; // any index in [0,15]
@@ -194,29 +192,25 @@ void checkGlobalVectorInterfaceMPI()
} // end namespace
-TEST(Math, CheckInterface_DenseVector)
-{
- checkGlobalVectorInterface<MathLib::DenseVector<double> >();
-}
-
-#ifdef OGS_USE_EIGEN
-TEST(Math, CheckInterface_EigenVector)
-{
- checkGlobalVectorInterface<MathLib::EigenVector >();
-}
-#endif
-
-#ifdef USE_LIS
+//--------------------------------------------
+#if defined(USE_LIS)
TEST(Math, CheckInterface_LisVector)
{
checkGlobalVectorInterface<MathLib::LisVector >();
}
-#endif
-
-//--------------------------------------------
-#ifdef USE_PETSC
+#elif defined(USE_PETSC)
TEST(MPITest_Math, CheckInterface_PETScVector)
{
checkGlobalVectorInterfaceMPI<MathLib::PETScVector >();
}
+#elif defined(OGS_USE_EIGEN)
+TEST(Math, CheckInterface_EigenVector)
+{
+ checkGlobalVectorInterface<MathLib::EigenVector >();
+}
+#else
+TEST(Math, CheckInterface_DenseVector)
+{
+ checkGlobalVectorInterface<MathLib::DenseVector<double> >();
+}
#endif