diff --git a/MathLib/LinAlg/Lis/LisMatrix.cpp b/MathLib/LinAlg/Lis/LisMatrix.cpp
index 8aae74c5ccbb5dbcaaf72def19b17f5bfacfb18e..1bd85c440d84ead22b1fecb4819f20f98fe49ff0 100644
--- a/MathLib/LinAlg/Lis/LisMatrix.cpp
+++ b/MathLib/LinAlg/Lis/LisMatrix.cpp
@@ -23,105 +23,104 @@ namespace MathLib
{
LisMatrix::LisMatrix(std::size_t n_rows, LisOption::MatrixType mat_type)
- : _n_rows(n_rows), _mat_type(mat_type), _is_assembled(false)
+ : _n_rows(n_rows), _mat_type(mat_type), _is_assembled(false)
{
int ierr = lis_matrix_create(0, &_AA);
checkLisError(ierr);
ierr = lis_matrix_set_size(_AA, 0, n_rows);
checkLisError(ierr);
lis_matrix_get_range(_AA, &_is, &_ie);
+ ierr = lis_vector_duplicate(_AA, &_diag);
+ checkLisError(ierr);
}
LisMatrix::~LisMatrix()
{
int ierr = lis_matrix_destroy(_AA);
checkLisError(ierr);
+ ierr = lis_vector_destroy(_diag);
+ checkLisError(ierr);
}
void LisMatrix::setZero()
{
- // A matrix has to be destroyed and created again because Lis doesn't provide a
- // function to set matrix entries to zero
- int ierr = lis_matrix_destroy(_AA);
- checkLisError(ierr);
- ierr = lis_matrix_create(0, &_AA);
- checkLisError(ierr);
- ierr = lis_matrix_set_size(_AA, 0, _n_rows);
- checkLisError(ierr);
-
- _is_assembled = false;
+ // A matrix has to be destroyed and created again because Lis doesn't provide a
+ // function to set matrix entries to zero
+ int ierr = lis_matrix_destroy(_AA);
+ checkLisError(ierr);
+ ierr = lis_matrix_create(0, &_AA);
+ checkLisError(ierr);
+ ierr = lis_matrix_set_size(_AA, 0, _n_rows);
+ checkLisError(ierr);
+ ierr = lis_vector_set_all(0.0, _diag);
+ checkLisError(ierr);
+
+ _is_assembled = false;
}
int LisMatrix::setValue(std::size_t rowId, std::size_t colId, double v)
{
- lis_matrix_set_value(LIS_INS_VALUE, rowId, colId, v, _AA);
- _is_assembled = false;
- return 0;
+ lis_matrix_set_value(LIS_INS_VALUE, rowId, colId, v, _AA);
+ if (rowId==colId)
+ lis_vector_set_value(LIS_INS_VALUE, rowId, v, _diag);
+ _is_assembled = false;
+ return 0;
}
int LisMatrix::addValue(std::size_t rowId, std::size_t colId, double v)
{
- lis_matrix_set_value(LIS_ADD_VALUE, rowId, colId, v, _AA);
- _is_assembled = false;
- return 0;
+ lis_matrix_set_value(LIS_ADD_VALUE, rowId, colId, v, _AA);
+ if (rowId==colId)
+ lis_vector_set_value(LIS_ADD_VALUE, rowId, v, _diag);
+ _is_assembled = false;
+ return 0;
}
void LisMatrix::write(const std::string &filename) const
{
- if (!_is_assembled)
- throw std::logic_error("LisMatrix::write(): matrix not assembled.");
+ if (!_is_assembled)
+ throw std::logic_error("LisMatrix::write(): matrix not assembled.");
lis_output_matrix(_AA, LIS_FMT_MM, const_cast<char*>(filename.c_str()));
}
double LisMatrix::getMaxDiagCoeff()
{
- if (!_is_assembled)
- throw std::logic_error("LisMatrix::getMaxDiagCoeff(): matrix not assembled.");
- LIS_VECTOR diag;
- int ierr = lis_vector_create(0, &diag);
- checkLisError(ierr);
- ierr = lis_vector_set_size(diag, 0, _n_rows);
- checkLisError(ierr);
-
- ierr = lis_matrix_get_diagonal(_AA, diag);
- checkLisError(ierr);
-
- double abs_max_entry;
- ierr = lis_vector_get_value(diag, 0, &abs_max_entry);
- checkLisError(ierr);
- abs_max_entry = std::abs(abs_max_entry);
- for (std::size_t k(1); k<_n_rows; ++k) {
- double tmp;
- ierr = lis_vector_get_value(diag, k, &tmp);
- checkLisError(ierr);
- if (abs_max_entry < std::abs(tmp)) {
- abs_max_entry = std::abs(tmp);
- }
- }
-
- return abs_max_entry;
+ double abs_max_entry;
+ int ierr = lis_vector_get_value(_diag, 0, &abs_max_entry);
+ checkLisError(ierr);
+ abs_max_entry = std::abs(abs_max_entry);
+ for (std::size_t k(1); k<_n_rows; ++k) {
+ double tmp;
+ ierr = lis_vector_get_value(_diag, k, &tmp);
+ checkLisError(ierr);
+ if (abs_max_entry < std::abs(tmp)) {
+ abs_max_entry = std::abs(tmp);
+ }
+ }
+
+ return abs_max_entry;
}
void LisMatrix::matvec (const LisVector &x, LisVector &y) const
{
- if (!_is_assembled)
- throw std::logic_error("LisMatrix::matvec(): matrix not assembled.");
+ if (!_is_assembled)
+ throw std::logic_error("LisMatrix::matvec(): 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();
-
- if (!mat.isAssembled()) {
- int ierr = lis_matrix_set_type(A, static_cast<int>(mat.getMatrixType()));
- checkLisError(ierr);
- ierr = lis_matrix_assemble(A);
- checkLisError(ierr);
- mat._is_assembled = true;
- }
- return true;
+ LIS_MATRIX &A = mat.getRawMatrix();
+
+ if (!mat.isAssembled()) {
+ int ierr = lis_matrix_set_type(A, static_cast<int>(mat.getMatrixType()));
+ checkLisError(ierr);
+ ierr = lis_matrix_assemble(A);
+ checkLisError(ierr);
+ mat._is_assembled = true;
+ }
+ return true;
}
diff --git a/MathLib/LinAlg/Lis/LisMatrix.h b/MathLib/LinAlg/Lis/LisMatrix.h
index b6a61662ed02f53967f0d3201bf1083947e5d6b4..110d43543ee82419624c56d86051fb785001f06a 100644
--- a/MathLib/LinAlg/Lis/LisMatrix.h
+++ b/MathLib/LinAlg/Lis/LisMatrix.h
@@ -98,6 +98,7 @@ private:
std::size_t const _n_rows;
LisOption::MatrixType const _mat_type;
LIS_MATRIX _AA;
+ LIS_VECTOR _diag;
bool _is_assembled;
int _is; ///< location where the partial matrix _AA starts in global matrix.
int _ie; ///< location where the partial matrix _AA ends in global matrix.
diff --git a/MathLib/LinAlg/Lis/LisVector.cpp b/MathLib/LinAlg/Lis/LisVector.cpp
index 7965f654679b99054e52fa2048d423da016c5acf..4b176e74d65d8e48e197bb4471d1cf91f3d1961c 100644
--- a/MathLib/LinAlg/Lis/LisVector.cpp
+++ b/MathLib/LinAlg/Lis/LisVector.cpp
@@ -13,19 +13,18 @@
*/
#include "LisVector.h"
+#include "LisCheck.h"
namespace MathLib
{
LisVector::LisVector(std::size_t length)
-: _length(length)
{
lis_vector_create(0, &_vec);
lis_vector_set_size(_vec, 0, length);
}
LisVector::LisVector(LisVector const &src)
-: _length(src.size())
{
lis_vector_duplicate(src._vec, &_vec);
lis_vector_copy(src._vec, _vec);
@@ -58,6 +57,15 @@ LisVector& LisVector::operator= (double v)
return *this;
}
+std::size_t LisVector::size() const
+{
+ LIS_INT dummy;
+ LIS_INT size;
+ int const ierr = lis_vector_get_size(_vec, &dummy, &size);
+ checkLisError(ierr);
+ return size;
+}
+
void LisVector::write (const std::string &filename) const
{
lis_output_vector(_vec, LIS_FMT_PLAIN, const_cast<char*>(filename.c_str()));
diff --git a/MathLib/LinAlg/Lis/LisVector.h b/MathLib/LinAlg/Lis/LisVector.h
index b8a258cebb696c0e2883f7fcdbf755b485cba4b2..99c8c25b351de1a0c03ab238eda4c2d046acceb9 100644
--- a/MathLib/LinAlg/Lis/LisVector.h
+++ b/MathLib/LinAlg/Lis/LisVector.h
@@ -45,13 +45,13 @@ public:
virtual ~LisVector();
/// return a vector length
- std::size_t size() const {return _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 _length; }
+ std::size_t getRangeEnd() const { return this->size(); }
/// set all values in this vector
LisVector& operator= (double v);
@@ -103,7 +103,6 @@ public:
}
}
private:
- std::size_t const _length;
LIS_VECTOR _vec;
};