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

Merge pull request #1199 from norihiro-w/remove-GlobalDenseMatVec

Remove unused GlobalDenseMatrix and DenseVector
parents 9e9cf967 f909ff77
No related branches found
No related tags found
No related merge requests found
Showing
with 15 additions and 631 deletions
...@@ -11,8 +11,6 @@ ...@@ -11,8 +11,6 @@
#include <vector> #include <vector>
#include "MathLib/LinAlg/Dense/DenseTools.h"
#ifdef OGS_USE_EIGEN #ifdef OGS_USE_EIGEN
#include "MathLib/LinAlg/Eigen/EigenTools.h" #include "MathLib/LinAlg/Eigen/EigenTools.h"
#endif // OGS_USE_EIGEN #endif // OGS_USE_EIGEN
......
/**
* \file
* \author Norihiro Watanabe
* \date 2013-04-16
* \brief Implementation of dense matrix and vector utility functions.
*
* \copyright
* Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*
*/
#include "DenseTools.h"
namespace MathLib
{
void applyKnownSolution(DenseMatrix<double>& eqsA, DenseVector<double>& eqsRHS,
DenseMatrix<double>::IDX_T k, double val)
{
const std::size_t n_rows = eqsA.getNRows();
const std::size_t n_cols = eqsA.getNCols();
// set all entries of the k-th row of the matrix to zero
// except the diagonal entry that is set to one
for (std::size_t j = 0; j < n_cols; j++)
eqsA(k, j) = .0;
eqsA(k, k) = 1.0;
// b_i -= A(i,k)*val, i!=k and
// set the entries of the k-th column of the matrix to zero
// except the diagonal entry A(k,k)
for (std::size_t i = 0; i < k; ++i)
{
eqsRHS[i] -= eqsA(i, k) * val;
eqsA(i, k) = 0.0;
}
for (std::size_t i = k + 1; i < n_rows; ++i)
{
eqsRHS[i] -= eqsA(i, k) * val;
eqsA(i, k) = 0.0;
}
// b_k = val
eqsRHS[k] = val;
}
void applyKnownSolution(
DenseMatrix<double>& A, DenseVector<double>& b, DenseVector<double>& /*x*/,
const std::vector<DenseMatrix<double>::IDX_T>& vec_knownX_id,
const std::vector<double>& vec_knownX_x)
{
const std::size_t n_bc = vec_knownX_id.size();
for (std::size_t i_bc = 0; i_bc < n_bc; i_bc++)
{
applyKnownSolution(A, b, vec_knownX_id[i_bc], vec_knownX_x[i_bc]);
}
}
} // MathLib
/**
* \file
* \author Norihiro Watanabe
* \date 2013-04-16
* \brief Declaration of dense matrix and vector utility functions.
*
* \copyright
* Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*
*/
#ifndef DENSETOOLS_H_
#define DENSETOOLS_H_
#include <vector>
#include "DenseMatrix.h"
#include "DenseVector.h"
namespace MathLib
{
/**
* Apply known solutions to a system of linear equations.
*
* This function introduces the given constrains by diagonalizing a coefficient matrix.
* Symmetry of the matrix is preserved.
*
* @param A Coefficient matrix
* @param b RHS vector
* @param vec_knownX_id a vector of known solution entry IDs
* @param vec_knownX_x a vector of known solutions
*/
void applyKnownSolution(
DenseMatrix<double>& A, DenseVector<double>& b, DenseVector<double>& /*x*/,
const std::vector<DenseMatrix<double>::IDX_T>& vec_knownX_id,
const std::vector<double>& vec_knownX_x);
/**
* Apply known solutions to a system of linear equations \f$A x = b\f$.
*
* This function introduces the given constrain into the system of linear
* equations. For this purpose it modifies the entries of \f$A\f$ in the
* \f$k\f$-th row and column, i.e. all those entries are set to zero except
* the diagonal entry that is set to one. The right
* hand side \f$b\f$ is modified, too. The symmetry of \f$A\f$ is preserved.
*
* @param A Coefficient matrix
* @param b RHS vector
* @param row_id a known solution entry ID
* @param val a known solution
*/
void applyKnownSolution(DenseMatrix<double>& A, DenseVector<double>& b,
DenseMatrix<double>::IDX_T row_id, double val);
} // MathLib
#endif //DENSETOOLS_H_
/**
* \file
* \author Norihiro Watanabe
* \date 2013-04-16
* \brief Definition of the DenseVector class.
*
* \copyright
* Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*
*/
#ifndef DENSEVECTOR_H_
#define DENSEVECTOR_H_
#include <vector>
#include <valarray>
#include <fstream>
#include <iterator>
namespace MathLib
{
/**
* Dense vector class
*/
template <typename T>
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=;
using std::valarray<T>::operator+=;
using std::valarray<T>::operator-=;
using std::valarray<T>::operator[];
/**
* Constructor for initialization of the number of rows
* @param nrows number of rows
*/
explicit DenseVector(std::size_t nrows=0)
: std::valarray<T>(nrows)
{}
/// 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(); }
/// get entry
T get(std::size_t i) const { return (*this)[i]; }
/// set a value to entry
void set(std::size_t i, T v) { (*this)[i] = v; }
/// add a value to entry
void add(std::size_t i, T v) { (*this)[i] += v; }
/**
* add a sub vector
* @param pos positions of each sub-vector entry in this vector
* @param sub_vec sub-vector
*/
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]);
}
}
/// Copy vector values.
void copyValues(std::vector<double>& u)
{
assert(u.size() == this->size());
std::copy(this->cbegin(), this->cend(), u.begin());
}
/**
* writes the matrix entries into a file
* @param filename output file name
*/
void write (const std::string &filename) const
{
std::ofstream os(filename);
os << *this;
os.close();
}
/// vector operation: add
void operator+= (const DenseVector<T>& v)
{
*this += static_cast<std::valarray<T> >(v);
}
/// vector operation: subtract
void operator-= (const DenseVector<T>& v)
{
*this -= static_cast<std::valarray<T> >(v);
}
};
/**
* writes a vector content into the output stream
* @param os the output stream
*/
template <typename T>
std::ostream& operator<<(std::ostream& os, DenseVector<T> const & v)
{
std::copy(std::begin(v), std::end(v), std::ostream_iterator<T>(os, "\n"));
return os;
}
}
#endif /* DENSEVECTOR_H_ */
/**
* @file GlobalDenseMatrix-impl.h
* @author Thomas Fischer
* @date Jun 10, 2013
* @brief
*
* @copyright
* Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/LICENSE.txt
*/
namespace MathLib
{
template<typename FP_TYPE, typename IDX_TYPE>
GlobalDenseMatrix<FP_TYPE, IDX_TYPE>::GlobalDenseMatrix(IDX_TYPE rows) :
DenseMatrix<FP_TYPE, IDX_TYPE>(rows, rows)
{}
template<typename FP_TYPE, typename IDX_TYPE>
GlobalDenseMatrix<FP_TYPE, IDX_TYPE>::GlobalDenseMatrix(IDX_TYPE rows, IDX_TYPE cols) :
DenseMatrix<FP_TYPE, IDX_TYPE>(rows, cols)
{}
template<typename FP_TYPE, typename IDX_TYPE>
GlobalDenseMatrix<FP_TYPE, IDX_TYPE>::GlobalDenseMatrix (IDX_TYPE rows, IDX_TYPE cols, const FP_TYPE& val) :
DenseMatrix<FP_TYPE, IDX_TYPE>(rows, cols, val)
{}
template<typename FP_TYPE, typename IDX_TYPE>
GlobalDenseMatrix<FP_TYPE, IDX_TYPE>::GlobalDenseMatrix(const GlobalDenseMatrix &src) :
DenseMatrix<FP_TYPE, IDX_TYPE>(src.getNRows(), src.getNCols())
{
std::copy(src._data, src._data + this->_n_rows * this->_n_cols, this->_data);
}
template<typename FP_TYPE, typename IDX_TYPE>
void
GlobalDenseMatrix<FP_TYPE, IDX_TYPE>::setZero()
{
std::fill(this->_data, this->_data + this->_n_rows * this->_n_cols, static_cast<FP_TYPE>(0));
}
template<typename FP_TYPE, typename IDX_TYPE>
bool
GlobalDenseMatrix<FP_TYPE, IDX_TYPE>::setValue(IDX_TYPE row, IDX_TYPE col, FP_TYPE val)
{
if (row >= this->_n_rows || col >= this->_n_cols)
return false;
this->operator()(row, col) = val;
return true;
}
template<typename FP_TYPE, typename IDX_TYPE>
bool
GlobalDenseMatrix<FP_TYPE, IDX_TYPE>::add(IDX_TYPE row, IDX_TYPE col, FP_TYPE val)
{
if (row >= this->_n_rows || col >= this->_n_cols)
return false;
this->operator()(row, col) += val;
return true;
}
template<typename FP_TYPE, typename IDX_TYPE>
template<class T_DENSE_MATRIX>
void
GlobalDenseMatrix<FP_TYPE, IDX_TYPE>::add(std::vector<IDX_TYPE> const& row_pos, std::vector<IDX_TYPE> const& col_pos,
const T_DENSE_MATRIX &sub_matrix, FP_TYPE 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 IDX_TYPE row = row_pos[i];
for (std::size_t j = 0; j < n_cols; j++) {
const IDX_TYPE col = col_pos[j];
add(row, col, fkt * sub_matrix(i, j));
}
}
}
template<typename FP_TYPE, typename IDX_TYPE>
void
GlobalDenseMatrix<FP_TYPE, IDX_TYPE>::multiply( const DenseVector<FP_TYPE> &x, DenseVector<FP_TYPE> &y) const
{
this->axpy (1.0, &x[0], 0.0, &y[0]);
}
} // end namespace MathLib
/**
* @file GlobalDenseMatrix.h
* @author Thomas Fischer
* @date Jun 4, 2013
* @brief
*
* @copyright
* Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/LICENSE.txt
*/
#ifndef GLOBALDENSEMATRIX_H_
#define GLOBALDENSEMATRIX_H_
#include <algorithm>
#include <vector>
#include "DenseMatrix.h"
#include "DenseVector.h"
#include "MathLib/LinAlg/RowColumnIndices.h"
namespace MathLib
{
template<typename FP_TYPE, typename IDX_TYPE = std::size_t>
class GlobalDenseMatrix: public DenseMatrix<FP_TYPE, IDX_TYPE>
{
public:
typedef FP_TYPE FP_T;
using IndexType = IDX_TYPE;
public:
/// Dense square matrix constructor.
GlobalDenseMatrix (IDX_TYPE rows);
/// Dense rectangular matrix constructor.
GlobalDenseMatrix (IDX_TYPE rows, IDX_TYPE cols);
GlobalDenseMatrix (IDX_TYPE rows, IDX_TYPE cols, const FP_TYPE& val);
GlobalDenseMatrix (const GlobalDenseMatrix &src);
virtual ~GlobalDenseMatrix() {}
/// return a start index of the active data range
IDX_TYPE getRangeBegin() const { return static_cast<IDX_TYPE>(0); }
/// return an end index of the active data range
IDX_TYPE getRangeEnd() const { return this->_n_rows; }
/**
* Method setZero() set all matrix entries to zero.
*/
virtual void setZero();
/**
* Set the value of the matrix entry (row,col) to val.
* @param row The index of the row of the matrix.
* @param col The index of the column of the matrix.
* @param val The value that shoud be set.
* @return False if row index or column index are to large, else true.
*/
virtual bool setValue(IDX_TYPE row, IDX_TYPE col, FP_TYPE val);
/**
* Method adds a value to the entry at position (row,col).
* @param row The index of the row of the matrix.
* @param col The index of the column of the matrix.
* @param val The value that shoud be added.
* @return False if row index or column index are to large, else true.
*/
virtual bool add(IDX_TYPE row, IDX_TYPE col, FP_TYPE val);
/// 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<IDX_TYPE> const& row_pos,
const T_DENSE_MATRIX &sub_matrix,
FP_TYPE fkt = static_cast<FP_TYPE>(1.0))
{
this->add(row_pos, row_pos, sub_matrix, fkt);
}
template<class T_DENSE_MATRIX>
void add(RowColumnIndices<IDX_TYPE> const& indices,
const T_DENSE_MATRIX &sub_matrix,
FP_TYPE fkt = static_cast<FP_TYPE>(1.0))
{
this->add(indices.rows, indices.columns, sub_matrix, fkt);
}
template<class T_DENSE_MATRIX>
void add(std::vector<IDX_TYPE> const& row_pos,
std::vector<IDX_TYPE> const& col_pos, const T_DENSE_MATRIX &sub_matrix,
FP_TYPE fkt = static_cast<FP_TYPE>(1.0));
/// y = mat * x
void multiply( const DenseVector<FP_TYPE> &x, DenseVector<FP_TYPE> &y) const;
};
} // end namespace MathLib
#include "GlobalDenseMatrix-impl.h"
#endif /* GLOBALDENSEMATRIX_H_ */
...@@ -13,7 +13,6 @@ ...@@ -13,7 +13,6 @@
#include "CRSTools.h" #include "CRSTools.h"
#include "CRSMatrix.h" #include "CRSMatrix.h"
#include "../Dense/DenseVector.h"
namespace MathLib namespace MathLib
{ {
......
...@@ -16,8 +16,7 @@ ...@@ -16,8 +16,7 @@
#include <cmath> #include <cmath>
#include <vector> #include <vector>
#include "MathLib/LinAlg/Dense/DenseMatrix.h" #include <Eigen/Dense>
#include "MathLib/LinAlg/Dense/DenseVector.h"
#include "MeshLib/Elements/Element.h" #include "MeshLib/Elements/Element.h"
#include "MeshLib/Mesh.h" #include "MeshLib/Mesh.h"
...@@ -26,8 +25,9 @@ ...@@ -26,8 +25,9 @@
template<typename IndexType>struct SteadyDiffusion2DExample1 template<typename IndexType>struct SteadyDiffusion2DExample1
{ {
using LocalMatrixType = MathLib::DenseMatrix<double>; using LocalMatrixType =
using LocalVectorType = MathLib::DenseVector<double>; Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
using LocalVectorType = Eigen::VectorXd;
template <typename GlobalMatrix, typename GlobalVector> template <typename GlobalMatrix, typename GlobalVector>
class LocalAssemblerData class LocalAssemblerData
......
...@@ -101,9 +101,6 @@ REGISTER_TYPED_TEST_CASE_P(AssemblerLibVectorMatrixBuilder, ...@@ -101,9 +101,6 @@ REGISTER_TYPED_TEST_CASE_P(AssemblerLibVectorMatrixBuilder,
DISABLED_createVector, DISABLED_createMatrix); DISABLED_createVector, DISABLED_createMatrix);
#endif #endif
#include "MathLib/LinAlg/Dense/DenseVector.h"
#include "MathLib/LinAlg/Dense/GlobalDenseMatrix.h"
#ifdef USE_LIS #ifdef USE_LIS
#include "MathLib/LinAlg/Lis/LisVector.h" #include "MathLib/LinAlg/Lis/LisVector.h"
#include "MathLib/LinAlg/Lis/LisMatrix.h" #include "MathLib/LinAlg/Lis/LisMatrix.h"
...@@ -120,8 +117,9 @@ REGISTER_TYPED_TEST_CASE_P(AssemblerLibVectorMatrixBuilder, ...@@ -120,8 +117,9 @@ REGISTER_TYPED_TEST_CASE_P(AssemblerLibVectorMatrixBuilder,
#endif // OGS_USE_EIGEN #endif // OGS_USE_EIGEN
typedef ::testing::Types typedef ::testing::Types
< AssemblerLib::VectorMatrixBuilder< <
MathLib::GlobalDenseMatrix<double>, MathLib::DenseVector<double>> AssemblerLib::VectorMatrixBuilder<
MathLib::EigenMatrix, MathLib::EigenVector>
#ifdef USE_LIS #ifdef USE_LIS
, AssemblerLib::VectorMatrixBuilder< , AssemblerLib::VectorMatrixBuilder<
MathLib::LisMatrix, MathLib::LisVector> MathLib::LisMatrix, MathLib::LisVector>
...@@ -130,10 +128,6 @@ typedef ::testing::Types ...@@ -130,10 +128,6 @@ typedef ::testing::Types
, AssemblerLib::VectorMatrixBuilder< , AssemblerLib::VectorMatrixBuilder<
MathLib::PETScMatrix, MathLib::PETScVector> MathLib::PETScMatrix, MathLib::PETScVector>
#endif // USE_PETSC #endif // USE_PETSC
#ifdef OGS_USE_EIGEN
, AssemblerLib::VectorMatrixBuilder<
MathLib::EigenMatrix, MathLib::EigenVector>
#endif // OGS_USE_EIGEN
> TestTypes; > TestTypes;
INSTANTIATE_TYPED_TEST_CASE_P(templated, AssemblerLibVectorMatrixBuilder, INSTANTIATE_TYPED_TEST_CASE_P(templated, AssemblerLibVectorMatrixBuilder,
......
...@@ -18,8 +18,6 @@ ...@@ -18,8 +18,6 @@
#include <gtest/gtest.h> #include <gtest/gtest.h>
#include "MathLib/LinAlg/Dense/DenseMatrix.h"
#include "MathLib/LinAlg/Dense/DenseVector.h"
#include "MathLib/LinAlg/Solvers/GaussAlgorithm.h" #include "MathLib/LinAlg/Solvers/GaussAlgorithm.h"
TEST(MathLib, DenseGaussAlgorithm) TEST(MathLib, DenseGaussAlgorithm)
...@@ -97,69 +95,3 @@ TEST(MathLib, DenseGaussAlgorithm) ...@@ -97,69 +95,3 @@ TEST(MathLib, DenseGaussAlgorithm)
delete [] b3_copy; delete [] b3_copy;
} }
TEST(MathLib, DenseGaussAlgorithmDenseVector)
{
std::size_t n_rows(50);
std::size_t n_cols(n_rows);
MathLib::DenseMatrix<double,std::size_t> mat(n_rows, n_cols);
// *** fill matrix with arbitrary values
// ** initialize random seed
srand ( static_cast<unsigned>(time(nullptr)) );
// ** loop over rows and columns
for (std::size_t i(0); i<n_rows; i++) {
for (std::size_t j(0); j<n_cols; j++) {
mat(i,j) = rand()/static_cast<double>(RAND_MAX);
}
}
// *** create solution vector, set all entries to 0.0
MathLib::DenseVector<double> x(n_cols);
std::fill(std::begin(x), std::end(x), 0.0);
MathLib::DenseVector<double> b0(mat * x);
// *** create other right hand sides,
// set all entries of the solution vector to 1.0
std::fill(std::begin(x), std::end(x), 1.0);
MathLib::DenseVector<double> b1(mat * x);
std::generate(std::begin(x), std::end(x), std::rand);
MathLib::DenseVector<double> b2(mat * x);
// right hand side and solution vector with random entries
MathLib::DenseVector<double> b3(mat * x);
MathLib::DenseVector<double> b3_copy(mat * x);
MathLib::DenseVector<double> x3 (n_cols);
std::generate(std::begin(x3),std::end(x3), std::rand);
MathLib::GaussAlgorithm<MathLib::DenseMatrix<double, std::size_t>,
MathLib::DenseVector<double>> gauss;
// solve with b0 as right hand side
gauss.solve(mat, b0, true);
for (std::size_t i(0); i<n_rows; i++) {
ASSERT_NEAR(b0[i], 0.0, 1e5 * std::numeric_limits<float>::epsilon());
}
// solve with b1 as right hand side
gauss.solve(mat, b1, false);
for (std::size_t i(0); i<n_rows; i++) {
ASSERT_NEAR(b1[i], 1.0, std::numeric_limits<float>::epsilon());
}
// solve with b2 as right hand side
gauss.solve(mat, b2, false);
for (std::size_t i(0); i<n_rows; i++) {
ASSERT_NEAR(fabs(b2[i]-x[i])/fabs(x[i]), 0.0, std::numeric_limits<float>::epsilon());
}
gauss.solve(mat, b3, x3, false);
for (std::size_t i(0); i<n_rows; i++) {
ASSERT_NEAR(fabs(x3[i]-x[i])/fabs(x[i]), 0.0, std::numeric_limits<float>::epsilon());
}
// assure entries of vector b3 are not changed
for (std::size_t i(0); i<n_rows; i++) {
ASSERT_NEAR(fabs(b3[i]-b3_copy[i])/fabs(b3[i]), 0.0, std::numeric_limits<float>::epsilon());
}
}
/**
* @file TestGlobalDenseMatrix.cpp
* @author Thomas Fischer
* @date Jun 5, 2013
* @brief
*
* @copyright
* Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/LICENSE.txt
*/
#include <limits>
#include <gtest/gtest.h>
#include "MathLib/LinAlg/Dense/GlobalDenseMatrix.h"
TEST(MathLib, GlobalDenseMatrix)
{
const std::size_t n_rows(5);
const std::size_t n_cols(5);
MathLib::GlobalDenseMatrix<double, std::size_t> mat0(n_rows,n_cols);
MathLib::GlobalDenseMatrix<double, std::size_t> mat1(n_rows,n_cols);
MathLib::GlobalDenseMatrix<double, std::size_t> mat2(n_rows,n_cols-1);
for (std::size_t i(0); i<n_rows; i++) {
for (std::size_t j(0); j<n_cols; j++) {
mat0(i,j) = 1.0 / (i+1.0+j+1.0);
}
}
mat1.setZero();
mat1 = mat0;
for (std::size_t i(0); i<n_rows; i++) {
for (std::size_t j(0); j<n_cols; j++) {
ASSERT_NEAR(1.0/(i+j+2.0), mat1(i,j), std::numeric_limits<double>::epsilon());
}
}
ASSERT_THROW(mat2 = mat1, std::range_error);
}
...@@ -15,9 +15,7 @@ ...@@ -15,9 +15,7 @@
#include <gtest/gtest.h> #include <gtest/gtest.h>
#include "MathLib/LinAlg/Dense/DenseVector.h"
#include "MathLib/LinAlg/Dense/DenseMatrix.h" #include "MathLib/LinAlg/Dense/DenseMatrix.h"
#include "MathLib/LinAlg/Dense/GlobalDenseMatrix.h"
#include "MathLib/LinAlg/FinalizeMatrixAssembly.h" #include "MathLib/LinAlg/FinalizeMatrixAssembly.h"
#include "MathLib/LinAlg/ApplyKnownSolution.h" #include "MathLib/LinAlg/ApplyKnownSolution.h"
#include "MathLib/LinAlg/Solvers/GaussAlgorithm.h" #include "MathLib/LinAlg/Solvers/GaussAlgorithm.h"
...@@ -70,7 +68,7 @@ void setMatrix9x9(T_Mat &mat) ...@@ -70,7 +68,7 @@ void setMatrix9x9(T_Mat &mat)
template<typename IntType> struct Example1 template<typename IntType> struct Example1
{ {
MathLib::GlobalDenseMatrix<double> mat; MathLib::EigenMatrix mat;
std::vector<IntType> vec_dirichlet_bc_id; std::vector<IntType> vec_dirichlet_bc_id;
std::vector<double> vec_dirichlet_bc_value; std::vector<double> vec_dirichlet_bc_value;
static const std::size_t dim_eqs = 9; static const std::size_t dim_eqs = 9;
...@@ -110,7 +108,7 @@ void checkLinearSolverInterface(T_MATRIX &A, BaseLib::ConfigTree const& ls_optio ...@@ -110,7 +108,7 @@ void checkLinearSolverInterface(T_MATRIX &A, BaseLib::ConfigTree const& ls_optio
{ {
for (std::size_t j=0; j<ex1.dim_eqs; j++) for (std::size_t j=0; j<ex1.dim_eqs; j++)
{ {
double v = ex1.mat(i, j); double v = ex1.mat.get(i, j);
if (v!=.0) if (v!=.0)
A.add(i, j, v); A.add(i, j, v);
} }
...@@ -202,21 +200,6 @@ void checkLinearSolverInterface(T_MATRIX& A, T_VECTOR& b, ...@@ -202,21 +200,6 @@ void checkLinearSolverInterface(T_MATRIX& A, T_VECTOR& b,
} // end namespace } // end namespace
TEST(MathLib, CheckInterface_GaussAlgorithm)
{
boost::property_tree::ptree t_root;
BaseLib::ConfigTree conf(t_root, "",
BaseLib::ConfigTree::onerror, BaseLib::ConfigTree::onwarning);
using Example = Example1<std::size_t>;
typedef MathLib::GaussAlgorithm<MathLib::GlobalDenseMatrix<double>, MathLib::DenseVector<double> > LinearSolverType;
MathLib::GlobalDenseMatrix<double> A(Example::dim_eqs, Example::dim_eqs);
checkLinearSolverInterface<MathLib::GlobalDenseMatrix<double>,
MathLib::DenseVector<double>, LinearSolverType, std::size_t>(
A, conf);
}
#ifdef OGS_USE_EIGEN #ifdef OGS_USE_EIGEN
TEST(Math, CheckInterface_Eigen) TEST(Math, CheckInterface_Eigen)
{ {
......
/**
* @file TestMatrixAssemblyTemplate.cpp
* @author Thomas Fischer
* @date Jun 13, 2013
* @brief
*
* @copyright
* Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/LICENSE.txt
*/
#include <gtest/gtest.h>
#include "MathLib/LinAlg/FinalizeMatrixAssembly.h"
#include "MathLib/LinAlg/Dense/GlobalDenseMatrix.h"
TEST(MathLib, GlobalDenseMatrixAssembly)
{
const std::size_t n_rows(5);
const std::size_t n_cols(5);
MathLib::GlobalDenseMatrix<double, std::size_t> mat0(n_rows,n_cols);
// assembly entries
for (std::size_t i(0); i<n_rows; i++) {
for (std::size_t j(0); j<n_cols; j++) {
mat0(i,j) = 1.0 / (i+1.0+j+1.0);
}
}
ASSERT_TRUE(finalizeMatrixAssembly(mat0));
}
...@@ -11,8 +11,9 @@ ...@@ -11,8 +11,9 @@
#include <gtest/gtest.h> #include <gtest/gtest.h>
#include <valarray>
#include "MathLib/LinAlg/Dense/DenseMatrix.h" #include "MathLib/LinAlg/Dense/DenseMatrix.h"
#include "MathLib/LinAlg/Dense/DenseVector.h"
#include "MathLib/LinAlg/Solvers/GaussAlgorithm.h" #include "MathLib/LinAlg/Solvers/GaussAlgorithm.h"
#include "MathLib/Nonlinear/NewtonRaphson.h" #include "MathLib/Nonlinear/NewtonRaphson.h"
#include "Tests/TestTools.h" #include "Tests/TestTools.h"
...@@ -21,7 +22,7 @@ namespace ...@@ -21,7 +22,7 @@ namespace
{ {
typedef MathLib::DenseMatrix<double> MatrixType; typedef MathLib::DenseMatrix<double> MatrixType;
typedef MathLib::DenseVector<double> VectorType; typedef std::valarray<double> VectorType;
typedef MathLib::GaussAlgorithm<MatrixType,VectorType> DenseSolverType; typedef MathLib::GaussAlgorithm<MatrixType,VectorType> DenseSolverType;
template<class F_JACOBIAN> template<class F_JACOBIAN>
......
...@@ -13,8 +13,9 @@ ...@@ -13,8 +13,9 @@
#include <gtest/gtest.h> #include <gtest/gtest.h>
#include <valarray>
#include "MathLib/LinAlg/Dense/DenseMatrix.h" #include "MathLib/LinAlg/Dense/DenseMatrix.h"
#include "MathLib/LinAlg/Dense/DenseVector.h"
#include "MathLib/LinAlg/Solvers/GaussAlgorithm.h" #include "MathLib/LinAlg/Solvers/GaussAlgorithm.h"
#include "MathLib/Nonlinear/Picard.h" #include "MathLib/Nonlinear/Picard.h"
#include "Tests/TestTools.h" #include "Tests/TestTools.h"
...@@ -23,7 +24,7 @@ namespace ...@@ -23,7 +24,7 @@ namespace
{ {
typedef MathLib::DenseMatrix<double> MatrixType; typedef MathLib::DenseMatrix<double> MatrixType;
typedef MathLib::DenseVector<double> VectorType; typedef std::valarray<double> VectorType;
typedef MathLib::GaussAlgorithm<MatrixType,VectorType> DenseSolverType; typedef MathLib::GaussAlgorithm<MatrixType,VectorType> DenseSolverType;
......
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