diff --git a/MathLib/LinAlg/BLAS.cpp b/MathLib/LinAlg/BLAS.cpp
index c72cb1684840ee6444aec83bb7dc260883ab583a..b104ad290431a7265b293981ba6e9f1456d548b2 100644
--- a/MathLib/LinAlg/BLAS.cpp
+++ b/MathLib/LinAlg/BLAS.cpp
@@ -16,7 +16,9 @@
#include <Eigen/Core>
-namespace MathLib { namespace BLAS
+namespace MathLib
+{
+namespace BLAS
{
// Explicit specialization
@@ -70,7 +72,7 @@ namespace MathLib { namespace BLAS
void copy(PETScVector const& x, PETScVector& y)
{
- y = x;
+ VecCopy(x.getRawVector(), y.getRawVector());
}
void scale(PETScVector& x, double const a)
@@ -113,7 +115,9 @@ void componentwiseDivide(PETScVector& w,
template<>
double norm1(PETScVector const& x)
{
- return x.getNorm(MathLib::VecNormType::NORM1);
+ PetscScalar norm = 0.;
+ VecNorm(x.getRawVector(), NORM_1, &norm);
+ return norm;
}
// Explicit specialization
@@ -121,7 +125,9 @@ double norm1(PETScVector const& x)
template<>
double norm2(PETScVector const& x)
{
- return x.getNorm(MathLib::VecNormType::NORM2);
+ PetscScalar norm = 0.;
+ VecNorm(x.getRawVector(), NORM_2, &norm);
+ return norm;
}
// Explicit specialization
@@ -129,7 +135,9 @@ double norm2(PETScVector const& x)
template<>
double normMax(PETScVector const& x)
{
- return x.getNorm(MathLib::VecNormType::INFINITY_N);
+ PetscScalar norm = 0.;
+ VecNorm(x.getRawVector(), NORM_INFINITY, &norm);
+ return norm;
}
@@ -218,7 +226,7 @@ void copy(EigenVector const& x, EigenVector& y)
void scale(EigenVector& x, double const a)
{
- x *= a;
+ x.getRawVector() *= a;
}
// y = a*y + X
@@ -312,7 +320,7 @@ void axpy(EigenMatrix& Y, double const a, EigenMatrix const& X)
void matMult(EigenMatrix const& A, EigenVector const& x, EigenVector& y)
{
assert(&x != &y);
- A.multiply(x, y);
+ y.getRawVector() = A.getRawMatrix() * x.getRawVector();
}
// v3 = A*v1 + v2
diff --git a/MathLib/LinAlg/Eigen/EigenMatrix.h b/MathLib/LinAlg/Eigen/EigenMatrix.h
index 0fdbc3c47b5076141140d9423f35005004981d88..cfd7d645ee7cee1c32a634c9d81fe0487f5df410 100644
--- a/MathLib/LinAlg/Eigen/EigenMatrix.h
+++ b/MathLib/LinAlg/Eigen/EigenMatrix.h
@@ -142,12 +142,6 @@ public:
return _mat.diagonal().maxCoeff();
}
- /// y = mat * x
- void multiply(const EigenVector &x, EigenVector &y) const
- {
- y.getRawVector() = _mat * x.getRawVector();
- }
-
/// return always true, i.e. the matrix is always ready for use
bool isAssembled() const { return true; }
diff --git a/MathLib/LinAlg/Eigen/EigenVector.h b/MathLib/LinAlg/Eigen/EigenVector.h
index 83a35c6d4902d9a86f3ef492f0a4f1a3c39b09d3..667ff6ebfc1b3e1a46cbdea48da7ae661361515e 100644
--- a/MathLib/LinAlg/Eigen/EigenVector.h
+++ b/MathLib/LinAlg/Eigen/EigenVector.h
@@ -52,15 +52,9 @@ public:
/// return an end index of the active data range
std::size_t getRangeEnd() const { return size(); }
- /// set all values in this vector
- EigenVector& operator= (double v) { _vec.setConstant(v); return *this; }
-
// TODO preliminary
void setZero() { _vec.setZero(); }
- /// set all values in this vector
- EigenVector& operator*= (double v) { _vec *= v; return *this; }
-
/// access entry
double const & operator[] (IndexType rowId) const { return _vec[rowId]; }
double& operator[] (IndexType rowId) { return _vec[rowId]; }
@@ -110,15 +104,6 @@ public:
/// return a raw Eigen vector object
const RawVectorType& getRawVector() const {return _vec; }
- /// vector operation: set data
- EigenVector& operator= (const EigenVector &src) { _vec = src._vec; return *this; }
-
- /// vector operation: add
- void operator+= (const EigenVector& v) { _vec += v._vec; }
-
- /// vector operation: subtract
- void operator-= (const EigenVector& v) { _vec -= v._vec; }
-
private:
RawVectorType _vec;
};
diff --git a/MathLib/LinAlg/Lis/LisMatrix.cpp b/MathLib/LinAlg/Lis/LisMatrix.cpp
index 9b324f503aa76b875e60d8aa3d94dd0f35fdc370..307d3b62ab9868facac938867cc74e92faf689d0 100644
--- a/MathLib/LinAlg/Lis/LisMatrix.cpp
+++ b/MathLib/LinAlg/Lis/LisMatrix.cpp
@@ -133,16 +133,6 @@ double LisMatrix::getMaxDiagCoeff()
return abs_max_entry;
}
-void LisMatrix::multiply(const LisVector &x, LisVector &y) const
-{
- if (!_is_assembled)
- {
- OGS_FATAL("LisMatrix::multiply(): matrix not assembled.");
- }
- int ierr = lis_matvec(_AA, const_cast<LisVector*>(&x)->getRawVector(), y.getRawVector());
- checkLisError(ierr);
-}
-
bool finalizeMatrixAssembly(LisMatrix &mat)
{
LIS_MATRIX &A = mat.getRawMatrix();
diff --git a/MathLib/LinAlg/Lis/LisMatrix.h b/MathLib/LinAlg/Lis/LisMatrix.h
index 44546837234a87143248890d8d4a6d6251ca747e..389adfa73a48104ef2292281749cf4c235fc30e0 100644
--- a/MathLib/LinAlg/Lis/LisMatrix.h
+++ b/MathLib/LinAlg/Lis/LisMatrix.h
@@ -118,9 +118,6 @@ public:
/// return a raw Lis matrix object
LIS_MATRIX& getRawMatrix() { return _AA; }
- /// y = mat * x
- void multiply(const LisVector &x, LisVector &y) const;
-
/// Add sub-matrix at positions \c row_pos and same column positions as the
/// given row positions.
template<class T_DENSE_MATRIX>
diff --git a/MathLib/LinAlg/Lis/LisVector.cpp b/MathLib/LinAlg/Lis/LisVector.cpp
index b841046e29f0d0996b6f9e88206d2fe9ce34a303..71263edc79cf9075c3372c3291d18ea8125ca3b2 100644
--- a/MathLib/LinAlg/Lis/LisVector.cpp
+++ b/MathLib/LinAlg/Lis/LisVector.cpp
@@ -44,28 +44,6 @@ LisVector::~LisVector()
lis_vector_destroy(_vec);
}
-LisVector& LisVector::operator=(const LisVector& src)
-{
- lis_vector_copy(src._vec, _vec);
- return *this;
-}
-
-void LisVector::operator+=(const LisVector& v)
-{
- lis_vector_axpy(1.0, v._vec, _vec);
-}
-
-void LisVector::operator-=(const LisVector& v)
-{
- lis_vector_axpy(-1.0, v._vec, _vec);
-}
-
-LisVector& LisVector::operator=(double v)
-{
- lis_vector_set_all(v, _vec);
- return *this;
-}
-
std::size_t LisVector::size() const
{
IndexType dummy;
diff --git a/MathLib/LinAlg/Lis/LisVector.h b/MathLib/LinAlg/Lis/LisVector.h
index 202674ee0bbc396e5c3f93613715c7bc0a12bae3..35b1033e89cb4836d970550b1bc2e18a5599078c 100644
--- a/MathLib/LinAlg/Lis/LisVector.h
+++ b/MathLib/LinAlg/Lis/LisVector.h
@@ -58,11 +58,9 @@ public:
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);
// TODO preliminary
- void setZero() { *this = 0.0; }
+ void setZero() { lis_vector_set_all(0.0, _vec); }
/// access entry
double operator[](IndexType rowId) const { return get(rowId); }
@@ -91,14 +89,6 @@ public:
/// 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>
diff --git a/MathLib/LinAlg/PETSc/PETScMatrix.h b/MathLib/LinAlg/PETSc/PETScMatrix.h
index 48a6fc7992b35f4f1222af65b7133a4ad48b255b..ad8a3ac4fb0ce0921dffc3badb42359cb52ac4bb 100644
--- a/MathLib/LinAlg/PETSc/PETScMatrix.h
+++ b/MathLib/LinAlg/PETSc/PETScMatrix.h
@@ -146,17 +146,6 @@ class PETScMatrix
*/
void setRowsColumnsZero(std::vector<PetscInt> const& row_pos);
- /*
- \brief Perform operation \f$ y = A x \f$
- \param vec The given vector, e.g. \f$ x \f$
- \param vec_r The result vector, e.g. \f$ y \f$
- Both of the two arguments must be created prior to be used.
- */
- void multiply(const PETScVector &vec, PETScVector &vec_r)
- {
- MatMult(_A, vec.getData(), vec_r.getData() );
- }
-
/*!
\brief Set a single entry with a value.
\param i The row index.
diff --git a/MathLib/LinAlg/PETSc/PETScVector.cpp b/MathLib/LinAlg/PETSc/PETScVector.cpp
index b56c4031a0a2d6283683d42e849e0244b9baa6bc..4f8b3cf05ed48b1c1549e4911ed049d2e80576e1 100644
--- a/MathLib/LinAlg/PETSc/PETScVector.cpp
+++ b/MathLib/LinAlg/PETSc/PETScVector.cpp
@@ -191,29 +191,6 @@ void PETScVector::restoreArray(PetscScalar* array) const
VecRestoreArray(_v, &array);
}
-PetscScalar PETScVector::getNorm(MathLib::VecNormType nmtype) const
-{
- NormType petsc_norm = NORM_1;
- switch(nmtype)
- {
- case MathLib::VecNormType::NORM1:
- petsc_norm = NORM_1;
- break;
- case MathLib::VecNormType::NORM2:
- petsc_norm = NORM_2;
- break;
- case MathLib::VecNormType::INFINITY_N:
- petsc_norm = NORM_INFINITY;
- break;
- default:
- break;
- }
-
- PetscScalar norm = 0.;
- VecNorm(_v, petsc_norm, &norm);
- return norm;
-}
-
void PETScVector::viewer(const std::string &file_name, const PetscViewerFormat vw_format) const
{
PetscViewer viewer;
diff --git a/MathLib/LinAlg/PETSc/PETScVector.h b/MathLib/LinAlg/PETSc/PETScVector.h
index f2abd9aa187fb70898c55a410c8f5ca1da7e52e5..783e3f35593b7f7e1ca7923a4a9033e778701d34 100644
--- a/MathLib/LinAlg/PETSc/PETScVector.h
+++ b/MathLib/LinAlg/PETSc/PETScVector.h
@@ -113,12 +113,6 @@ class PETScVector
return _end_rank;
}
- /*!
- Get norm of vector
- \param nmtype Norm type, default Euclidean norm
- */
- PetscScalar getNorm(const MathLib::VecNormType nmtype = MathLib::VecNormType::NORM2) const;
-
/*!
Insert a single entry with value.
\param i Entry index
@@ -212,43 +206,14 @@ class PETScVector
return _v;
}
- /// Initialize the vector with a constant value
- void operator = (const PetscScalar val)
- {
- VecSet(_v, val);
- }
-
// TODO preliminary
- void setZero() { *this = 0.0; }
-
- /// Overloaded operator: assign
- PETScVector& operator = (const PETScVector &v_in)
- {
- if (!_v) shallowCopy(v_in);
- VecCopy(v_in._v, _v);
-
- return *this;
- }
+ void setZero() { VecSet(_v, 0.0); }
/// Disallow moving.
/// \todo This operator should be implemented properly when doing a
/// general cleanup of all matrix and vector classes.
PETScVector& operator = (PETScVector &&) = delete;
- /// Overloaded operator: add
- void operator += (const PETScVector& v_in)
- {
- if (!_v) shallowCopy(v_in);
- VecAXPY(_v, 1.0, v_in._v);
- }
-
- /// Overloaded operator: subtract
- void operator -= (const PETScVector& v_in)
- {
- if (!_v) shallowCopy(v_in);
- VecAXPY(_v, -1.0, v_in._v);
- }
-
//! Exposes the underlying PETSc vector.
PETSc_Vec getRawVector() { return _v; }
diff --git a/MathLib/LinAlg/UnifiedMatrixSetters.cpp b/MathLib/LinAlg/UnifiedMatrixSetters.cpp
index da657029005c6089f9605df3113948214b91b270..447d27b04d9e16e1990944c4d6bbf1fb6fe7412e 100644
--- a/MathLib/LinAlg/UnifiedMatrixSetters.cpp
+++ b/MathLib/LinAlg/UnifiedMatrixSetters.cpp
@@ -89,9 +89,27 @@ void setVector(Eigen::VectorXd& v, MatrixVectorTraits<Eigen::VectorXd>::Index co
namespace MathLib
{
-double norm(PETScVector const& x)
+double norm(PETScVector const& x, MathLib::VecNormType nmtype = MathLib::VecNormType::NORM2)
{
- return x.getNorm();
+ NormType petsc_norm = NORM_1;
+ switch(nmtype)
+ {
+ case MathLib::VecNormType::NORM1:
+ petsc_norm = NORM_1;
+ break;
+ case MathLib::VecNormType::NORM2:
+ petsc_norm = NORM_2;
+ break;
+ case MathLib::VecNormType::INFINITY_N:
+ petsc_norm = NORM_INFINITY;
+ break;
+ default:
+ break;
+ }
+
+ PetscScalar norm = 0.;
+ VecNorm(x.getRawVector(), petsc_norm, &norm);
+ return norm;
}
void setVector(PETScVector& v,
diff --git a/NumLib/ODESolver/TimeDiscretization.h b/NumLib/ODESolver/TimeDiscretization.h
index e60165cb2604c7312e6ec498756ca015c866961e..01e806b8f608b5e93e21f7ae3d223eab2ca02334 100644
--- a/NumLib/ODESolver/TimeDiscretization.h
+++ b/NumLib/ODESolver/TimeDiscretization.h
@@ -219,12 +219,12 @@ public:
void setInitialState(const double t0, Vector const& x0) override {
_t = t0;
- _x_old = x0;
+ MathLib::BLAS::copy(x0, _x_old);
}
void pushState(const double /*t*/, Vector const& x, InternalMatrixStorage const&) override
{
- _x_old = x;
+ MathLib::BLAS::copy(x, _x_old);
}
void nextTimestep(const double t, const double delta_t) override {
@@ -272,12 +272,12 @@ public:
void setInitialState(const double t0, Vector const& x0) override {
_t = t0;
_t_old = t0;
- _x_old = x0;
+ MathLib::BLAS::copy(x0, _x_old);
}
void pushState(const double /*t*/, Vector const& x, InternalMatrixStorage const&) override
{
- _x_old = x;
+ MathLib::BLAS::copy(x, _x_old);
}
void nextTimestep(const double t, const double delta_t) override {
@@ -350,12 +350,12 @@ public:
void setInitialState(const double t0, Vector const& x0) override {
_t = t0;
- _x_old = x0;
+ MathLib::BLAS::copy(x0, _x_old);
}
void pushState(const double, Vector const& x, InternalMatrixStorage const& strg) override
{
- _x_old = x;
+ MathLib::BLAS::copy(x, _x_old);
strg.pushMatrices();
}
@@ -455,6 +455,7 @@ public:
void pushState(const double, Vector const& x, InternalMatrixStorage const&) override
{
+ namespace BLAS = MathLib::BLAS;
// TODO use boost cirular buffer?
// until _xs_old is filled, lower-order BDF formulas are used.
@@ -462,7 +463,7 @@ public:
_xs_old.push_back(
&MathLib::GlobalVectorProvider<Vector>::provider.getVector(x));
} else {
- *_xs_old[_offset] = x;
+ BLAS::copy(x, *_xs_old[_offset]);
_offset = (_offset+1) % _num_steps; // treat _xs_old as a circular buffer
}
}
diff --git a/Tests/MathLib/TestGlobalMatrixInterface.cpp b/Tests/MathLib/TestGlobalMatrixInterface.cpp
index 2ab60c43aa03f3a99dc5aab4fc7ddebd4d69044b..71c17d97c1e995bed2d3f7794f537c7f4c959ccf 100644
--- a/Tests/MathLib/TestGlobalMatrixInterface.cpp
+++ b/Tests/MathLib/TestGlobalMatrixInterface.cpp
@@ -94,24 +94,12 @@ void checkGlobalMatrixInterfaceMPI(T_MATRIX &m, T_VECTOR &v)
MathLib::finalizeMatrixAssembly(m);
- // Multiply by a vector
- v = 1.;
- const bool deep_copy = false;
- T_VECTOR y(v, deep_copy);
- m.multiply(v, y);
-
- ASSERT_EQ(sqrt(3*(3*3 + 7*7)), y.getNorm());
-
// set a value
m.set(2 * mrank, 2 * mrank, 5.0);
MathLib::finalizeMatrixAssembly(m);
// add a value
m.add(2 * mrank+1, 2 * mrank+1, 5.0);
MathLib::finalizeMatrixAssembly(m);
- m.multiply(v, y);
-
- ASSERT_EQ(sqrt((3*7*7 + 3*12*12)), y.getNorm());
-
}
// Rectanglular matrix
@@ -154,12 +142,6 @@ void checkGlobalRectangularMatrixInterfaceMPI(T_MATRIX &m, T_VECTOR &v)
MathLib::finalizeMatrixAssembly(m);
- // Multiply by a vector
- v = 1.;
- T_VECTOR y(m.getNumberOfRows());
- m.multiply(v, y);
-
- ASSERT_NEAR(6.*sqrt(6.), y.getNorm(), 1.e-10);
}
#endif // end of: ifdef USE_PETSC // or MPI
diff --git a/Tests/MathLib/TestGlobalVectorInterface.cpp b/Tests/MathLib/TestGlobalVectorInterface.cpp
index aa7a0c7491bbe7453882ef0182226d1a3ab7b906..0a63772956d4057a6debb9490eadb966de4509a0 100644
--- a/Tests/MathLib/TestGlobalVectorInterface.cpp
+++ b/Tests/MathLib/TestGlobalVectorInterface.cpp
@@ -50,23 +50,15 @@ void checkGlobalVectorInterface()
x.add(0, 1.0);
ASSERT_EQ(2.0, x.get(0));
- T_VECTOR y(x);
- ASSERT_EQ(2.0, y.get(0));
- ASSERT_EQ(0.0, y.get(1));
- y += x;
-
- ASSERT_EQ(4.0, y.get(0));
- y -= x;
- ASSERT_EQ(2.0, y.get(0));
- y = 1.0;
- ASSERT_EQ(1.0, y.get(0));
- y = x;
- ASSERT_EQ(2.0, y.get(0));
-
std::vector<double> local_vec(2, 1.0);
std::vector<GlobalIndexType> vec_pos(2);
vec_pos[0] = 0;
vec_pos[1] = 3;
+
+ T_VECTOR y(x);
+ ASSERT_EQ(2.0, y.get(0));
+ ASSERT_EQ(0.0, y.get(1));
+
y.add(vec_pos, local_vec);
ASSERT_EQ(3.0, y.get(0));
ASSERT_EQ(0.0, y.get(1));
@@ -93,24 +85,11 @@ void checkGlobalVectorInterfacePETSc()
//x.get(0) is expensive, only get local value. Use it for test purpose
ASSERT_EQ(.0, x.get(r0));
- x = 10.;
-
// Value of x is not copied to y
const bool deep_copy = false;
T_VECTOR y(x, deep_copy);
ASSERT_EQ(0, y.get(r0));
- y = 10.0;
- ASSERT_EQ(10, y.get(r0));
-
- y += x;
- ASSERT_EQ(20, y.get(r0));
- ASSERT_EQ(80., y.getNorm());
-
- y -= x;
- ASSERT_EQ(10, y.get(r0));
- ASSERT_EQ(40., y.getNorm());
-
std::vector<double> local_vec(2, 10.0);
std::vector<GlobalIndexType> vec_pos(2);
@@ -119,10 +98,6 @@ void checkGlobalVectorInterfacePETSc()
y.add(vec_pos, local_vec);
- double normy = std::sqrt(6.0*400+10.0*100);
-
- ASSERT_NEAR(0.0, normy-y.getNorm(), 1.e-10);
-
double x0[16];
double z[] =
{
@@ -188,7 +163,6 @@ void checkGlobalVectorInterfacePETSc()
// Deep copy
MathLib::finalizeVectorAssembly(x_fixed_p);
T_VECTOR x_deep_copied(x_fixed_p);
- ASSERT_NEAR(sqrt(3.0*5), x_deep_copied.getNorm(), 1.e-10);
// -----------------------------------------------------------------
// Vector with ghost entries
diff --git a/Tests/MathLib/TestLinearSolver.cpp b/Tests/MathLib/TestLinearSolver.cpp
index af7cbfd8123bb3e878d0ed92b0dd882f543133fe..5707cf37638af5b76c5dd10950f0f7834b7e35bc 100644
--- a/Tests/MathLib/TestLinearSolver.cpp
+++ b/Tests/MathLib/TestLinearSolver.cpp
@@ -15,6 +15,7 @@
#include <gtest/gtest.h>
+#include "MathLib/LinAlg/BLAS.h"
#include "MathLib/LinAlg/Dense/DenseMatrix.h"
#include "MathLib/LinAlg/FinalizeMatrixAssembly.h"
#include "MathLib/LinAlg/ApplyKnownSolution.h"
@@ -169,7 +170,7 @@ void checkLinearSolverInterface(T_MATRIX& A, T_VECTOR& b,
double x1[6];
x.getGlobalVector(x0);
- A.multiply(x, b);
+ MathLib::BLAS::matMult(A, x, b);
// apply BC
std::vector<int> bc_id; // Type must be int to match Petsc_Int