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

Merge pull request #172 from norihiro-w/fix-LisMatrix-getMaxDiagCoeff 

fix getMaxDiagCoeff() in LisMatrix. prepare a diagonal vector by ourself...
parents 6a41cb0f db8f308d
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
MathLib/LinAlg/Lis/LisMatrix.cpp
+ 55
− 56
View file @ ca07a0fa
...@@ -23,105 +23,104 @@ namespace MathLib ...@@ -23,105 +23,104 @@ namespace MathLib
{ {
LisMatrix::LisMatrix(std::size_t n_rows, LisOption::MatrixType mat_type) 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); int ierr = lis_matrix_create(0, &_AA);
checkLisError(ierr); checkLisError(ierr);
ierr = lis_matrix_set_size(_AA, 0, n_rows); ierr = lis_matrix_set_size(_AA, 0, n_rows);
checkLisError(ierr); checkLisError(ierr);
lis_matrix_get_range(_AA, &_is, &_ie); lis_matrix_get_range(_AA, &_is, &_ie);
ierr = lis_vector_duplicate(_AA, &_diag);
checkLisError(ierr);
} }
LisMatrix::~LisMatrix() LisMatrix::~LisMatrix()
{ {
int ierr = lis_matrix_destroy(_AA); int ierr = lis_matrix_destroy(_AA);
checkLisError(ierr); checkLisError(ierr);
ierr = lis_vector_destroy(_diag);
checkLisError(ierr);
} }
void LisMatrix::setZero() void LisMatrix::setZero()
{ {
// A matrix has to be destroyed and created again because Lis doesn't provide a // A matrix has to be destroyed and created again because Lis doesn't provide a
// function to set matrix entries to zero // function to set matrix entries to zero
int ierr = lis_matrix_destroy(_AA); int ierr = lis_matrix_destroy(_AA);
checkLisError(ierr); checkLisError(ierr);
ierr = lis_matrix_create(0, &_AA); ierr = lis_matrix_create(0, &_AA);
checkLisError(ierr); checkLisError(ierr);
ierr = lis_matrix_set_size(_AA, 0, _n_rows); ierr = lis_matrix_set_size(_AA, 0, _n_rows);
checkLisError(ierr); checkLisError(ierr);
ierr = lis_vector_set_all(0.0, _diag);
_is_assembled = false; checkLisError(ierr);
_is_assembled = false;
} }
int LisMatrix::setValue(std::size_t rowId, std::size_t colId, double v) int LisMatrix::setValue(std::size_t rowId, std::size_t colId, double v)
{ {
lis_matrix_set_value(LIS_INS_VALUE, rowId, colId, v, _AA); lis_matrix_set_value(LIS_INS_VALUE, rowId, colId, v, _AA);
_is_assembled = false; if (rowId==colId)
return 0; 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) int LisMatrix::addValue(std::size_t rowId, std::size_t colId, double v)
{ {
lis_matrix_set_value(LIS_ADD_VALUE, rowId, colId, v, _AA); lis_matrix_set_value(LIS_ADD_VALUE, rowId, colId, v, _AA);
_is_assembled = false; if (rowId==colId)
return 0; lis_vector_set_value(LIS_ADD_VALUE, rowId, v, _diag);
_is_assembled = false;
return 0;
} }
void LisMatrix::write(const std::string &filename) const void LisMatrix::write(const std::string &filename) const
{ {
if (!_is_assembled) if (!_is_assembled)
throw std::logic_error("LisMatrix::write(): matrix not assembled."); throw std::logic_error("LisMatrix::write(): matrix not assembled.");
lis_output_matrix(_AA, LIS_FMT_MM, const_cast<char*>(filename.c_str())); lis_output_matrix(_AA, LIS_FMT_MM, const_cast<char*>(filename.c_str()));
} }
double LisMatrix::getMaxDiagCoeff() double LisMatrix::getMaxDiagCoeff()
{ {
if (!_is_assembled) double abs_max_entry;
throw std::logic_error("LisMatrix::getMaxDiagCoeff(): matrix not assembled."); int ierr = lis_vector_get_value(_diag, 0, &abs_max_entry);
LIS_VECTOR diag; checkLisError(ierr);
int ierr = lis_vector_create(0, &diag); abs_max_entry = std::abs(abs_max_entry);
checkLisError(ierr); for (std::size_t k(1); k<_n_rows; ++k) {
ierr = lis_vector_set_size(diag, 0, _n_rows); double tmp;
checkLisError(ierr); ierr = lis_vector_get_value(_diag, k, &tmp);
checkLisError(ierr);
ierr = lis_matrix_get_diagonal(_AA, diag); if (abs_max_entry < std::abs(tmp)) {
checkLisError(ierr); abs_max_entry = std::abs(tmp);
}
double abs_max_entry; }
ierr = lis_vector_get_value(diag, 0, &abs_max_entry);
checkLisError(ierr); return abs_max_entry;
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 void LisMatrix::matvec (const LisVector &x, LisVector &y) const
{ {
if (!_is_assembled) if (!_is_assembled)
throw std::logic_error("LisMatrix::matvec(): matrix not assembled."); throw std::logic_error("LisMatrix::matvec(): matrix not assembled.");
int ierr = lis_matvec(_AA, const_cast<LisVector*>(&x)->getRawVector(), y.getRawVector()); int ierr = lis_matvec(_AA, const_cast<LisVector*>(&x)->getRawVector(), y.getRawVector());
checkLisError(ierr); checkLisError(ierr);
} }
bool finalizeMatrixAssembly(LisMatrix &mat) bool finalizeMatrixAssembly(LisMatrix &mat)
{ {
LIS_MATRIX &A = mat.getRawMatrix(); LIS_MATRIX &A = mat.getRawMatrix();
if (!mat.isAssembled()) { if (!mat.isAssembled()) {
int ierr = lis_matrix_set_type(A, static_cast<int>(mat.getMatrixType())); int ierr = lis_matrix_set_type(A, static_cast<int>(mat.getMatrixType()));
checkLisError(ierr); checkLisError(ierr);
ierr = lis_matrix_assemble(A); ierr = lis_matrix_assemble(A);
checkLisError(ierr); checkLisError(ierr);
mat._is_assembled = true; mat._is_assembled = true;
} }
return true; return true;
} }
......
MathLib/LinAlg/Lis/LisMatrix.h
+ 1
− 0
View file @ ca07a0fa
...@@ -98,6 +98,7 @@ private: ...@@ -98,6 +98,7 @@ private:
std::size_t const _n_rows; std::size_t const _n_rows;
LisOption::MatrixType const _mat_type; LisOption::MatrixType const _mat_type;
LIS_MATRIX _AA; LIS_MATRIX _AA;
LIS_VECTOR _diag;
bool _is_assembled; bool _is_assembled;
int _is; ///< location where the partial matrix _AA starts in global matrix. int _is; ///< location where the partial matrix _AA starts in global matrix.
int _ie; ///< location where the partial matrix _AA ends in global matrix. int _ie; ///< location where the partial matrix _AA ends in global matrix.
......
MathLib/LinAlg/Lis/LisVector.cpp
+ 10
− 2
View file @ ca07a0fa
...@@ -13,19 +13,18 @@ ...@@ -13,19 +13,18 @@
*/ */
#include "LisVector.h" #include "LisVector.h"
#include "LisCheck.h"
namespace MathLib namespace MathLib
{ {
LisVector::LisVector(std::size_t length) LisVector::LisVector(std::size_t length)
: _length(length)
{ {
lis_vector_create(0, &_vec); lis_vector_create(0, &_vec);
lis_vector_set_size(_vec, 0, length); lis_vector_set_size(_vec, 0, length);
} }
LisVector::LisVector(LisVector const &src) LisVector::LisVector(LisVector const &src)
: _length(src.size())
{ {
lis_vector_duplicate(src._vec, &_vec); lis_vector_duplicate(src._vec, &_vec);
lis_vector_copy(src._vec, _vec); lis_vector_copy(src._vec, _vec);
...@@ -58,6 +57,15 @@ LisVector& LisVector::operator= (double v) ...@@ -58,6 +57,15 @@ LisVector& LisVector::operator= (double v)
return *this; 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 void LisVector::write (const std::string &filename) const
{ {
lis_output_vector(_vec, LIS_FMT_PLAIN, const_cast<char*>(filename.c_str())); lis_output_vector(_vec, LIS_FMT_PLAIN, const_cast<char*>(filename.c_str()));
......
MathLib/LinAlg/Lis/LisVector.h
+ 2
− 3
View file @ ca07a0fa
...@@ -45,13 +45,13 @@ public: ...@@ -45,13 +45,13 @@ public:
virtual ~LisVector(); virtual ~LisVector();
/// return a vector length /// 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 /// return a start index of the active data range
std::size_t getRangeBegin() const { return 0;} std::size_t getRangeBegin() const { return 0;}
/// return an end index of the active data range /// 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 /// set all values in this vector
LisVector& operator= (double v); LisVector& operator= (double v);
...@@ -103,7 +103,6 @@ public: ...@@ -103,7 +103,6 @@ public:
} }
} }
private: private:
std::size_t const _length;
LIS_VECTOR _vec; LIS_VECTOR _vec;
}; };
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment