sp
/
tempest
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

using maps for Eigen solver instead of copies of the vectors 

Former-commit-id: d53075ab36
tempestpy_adaptions
dehnert 9 years ago
parent
7944778b54
commit
ba43e23984
1 changed files with 33 additions and 96 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 129
      src/solver/EigenLinearEquationSolver.cpp

129
src/solver/EigenLinearEquationSolver.cpp
View File

@ -100,23 +100,13 @@ namespace storm {
template<typename ValueType>
void EigenLinearEquationSolver<ValueType>::solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
// Translate the vectors x and b into Eigen's format.
Eigen::Matrix<ValueType, Eigen::Dynamic, 1> eigenB(b.size());
for (uint64_t index = 0; index < b.size(); ++index) {
eigenB[index] = b[index];
}
Eigen::Matrix<ValueType, Eigen::Dynamic, 1> eigenX(x.size());
for (uint64_t index = 0; index < x.size(); ++index) {
eigenX[index] = x[index];
}
// Map the input vectors to Eigen's format.
auto eigenX = Eigen::Matrix<ValueType, Eigen::Dynamic, 1>::Map(x.data(), x.size());
auto eigenB = Eigen::Matrix<ValueType, Eigen::Dynamic, 1>::Map(b.data(), b.size());
if (this->getSettings().getSolutionMethod() == EigenLinearEquationSolverSettings<ValueType>::SolutionMethod::SparseLU) {
Eigen::SparseLU<Eigen::SparseMatrix<ValueType>, Eigen::COLAMDOrdering<int>> solver;
solver.compute(*this->eigenA);
if (solver.info() != Eigen::Success) {
std::cout << solver.lastErrorMessage() << std::endl;
}
solver._solve_impl(eigenB, eigenX);
} else {
if (this->getSettings().getPreconditioner() == EigenLinearEquationSolverSettings<ValueType>::Preconditioner::Ilu) {
@ -133,40 +123,27 @@ namespace storm {
solver.solveWithGuess(eigenB, eigenX);
}
}
// Translate the solution from Eigen's format into our representation.
for (uint64_t index = 0; index < eigenX.size(); ++index) {
x[index] = eigenX[index];
}
}
template<typename ValueType>
void EigenLinearEquationSolver<ValueType>::performMatrixVectorMultiplication(std::vector<ValueType>& x, std::vector<ValueType> const* b, uint_fast64_t n, std::vector<ValueType>* multiplyResult) const {
// Translate the vectors x and b into Eigen's format.
std::unique_ptr<Eigen::Matrix<ValueType, Eigen::Dynamic, 1>> eigenB;
// Typedef the map-type so we don't have to spell it out.
typedef decltype(Eigen::Matrix<ValueType, Eigen::Dynamic, 1>::Map(b->data(), b->size())) MapType;
// Map the input vectors x and b to Eigen's format.
std::unique_ptr<MapType> eigenB;
if (b != nullptr) {
eigenB = std::make_unique<Eigen::Matrix<ValueType, Eigen::Dynamic, 1>>(b->size());
for (uint64_t index = 0; index < b->size(); ++index) {
(*eigenB)[index] = (*b)[index];
}
}
Eigen::Matrix<ValueType, Eigen::Dynamic, 1> eigenX(x.size());
for (uint64_t index = 0; index < x.size(); ++index) {
eigenX[index] = x[index];
eigenB = std::make_unique<MapType>(Eigen::Matrix<ValueType, Eigen::Dynamic, 1>::Map(b->data(), b->size()));
}
auto eigenX = Eigen::Matrix<ValueType, Eigen::Dynamic, 1>::Map(x.data(), x.size());
// Perform n matrix-vector multiplications.
for (uint64_t iteration = 0; iteration < n; ++iteration) {
eigenX = *this->eigenA * eigenX;
eigenX *= *this->eigenA * eigenX;
if (eigenB != nullptr) {
eigenX += *eigenB;
}
}
// Translate the solution from Eigen's format into our representation.
for (uint64_t index = 0; index < eigenX.size(); ++index) {
x[index] = eigenX[index];
}
}
template<typename ValueType>
@ -183,41 +160,26 @@ namespace storm {
template<>
void EigenLinearEquationSolver<storm::RationalNumber>::solveEquationSystem(std::vector<storm::RationalNumber>& x, std::vector<storm::RationalNumber> const& b, std::vector<storm::RationalNumber>* multiplyResult) const {
// Translate the vectors x and b into Eigen's format.
Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1> eigenB(b.size());
for (uint64_t index = 0; index < b.size(); ++index) {
eigenB[index] = b[index];
}
Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1> eigenX(x.size());
for (uint64_t index = 0; index < x.size(); ++index) {
eigenX[index] = x[index];
}
// Map the input vectors to Eigen's format.
auto eigenX = Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1>::Map(x.data(), x.size());
auto eigenB = Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1>::Map(b.data(), b.size());
Eigen::SparseLU<Eigen::SparseMatrix<storm::RationalNumber>, Eigen::COLAMDOrdering<int>> solver;
solver.compute(*eigenA);
solver._solve_impl(eigenB, eigenX);
// Translate the solution from Eigen's format into our representation.
for (uint64_t index = 0; index < eigenX.size(); ++index) {
x[index] = eigenX[index];
}
}
template<>
void EigenLinearEquationSolver<storm::RationalNumber>::performMatrixVectorMultiplication(std::vector<storm::RationalNumber>& x, std::vector<storm::RationalNumber> const* b, uint_fast64_t n, std::vector<storm::RationalNumber>* multiplyResult) const {
// Translate the vectors x and b into Eigen's format.
std::unique_ptr<Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1>> eigenB;
// Typedef the map-type so we don't have to spell it out.
typedef decltype(Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1>::Map(b->data(), b->size())) MapType;
// Map the input vectors x and b to Eigen's format.
std::unique_ptr<MapType> eigenB;
if (b != nullptr) {
eigenB = std::make_unique<Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1>>(b->size());
for (uint64_t index = 0; index < b->size(); ++index) {
(*eigenB)[index] = (*b)[index];
}
}
Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1> eigenX(x.size());
for (uint64_t index = 0; index < x.size(); ++index) {
eigenX[index] = x[index];
eigenB = std::make_unique<MapType>(Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1>::Map(b->data(), b->size()));
}
Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1> eigenX = Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1>::Map(x.data(), x.size());
// Perform n matrix-vector multiplications.
for (uint64_t iteration = 0; iteration < n; ++iteration) {
@ -226,52 +188,32 @@ namespace storm {
eigenX += *eigenB;
}
}
// Translate the solution from Eigen's format into our representation.
for (uint64_t index = 0; index < eigenX.size(); ++index) {
x[index] = eigenX[index];
}
}
// Specialization form storm::RationalFunction
template<>
void EigenLinearEquationSolver<storm::RationalFunction>::solveEquationSystem(std::vector<storm::RationalFunction>& x, std::vector<storm::RationalFunction> const& b, std::vector<storm::RationalFunction>* multiplyResult) const {
// Translate the vectors x and b into Eigen's format.
Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1> eigenB(b.size());
for (uint64_t index = 0; index < b.size(); ++index) {
eigenB[index] = b[index];
}
Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1> eigenX(x.size());
for (uint64_t index = 0; index < x.size(); ++index) {
eigenX[index] = x[index];
}
// Map the input vectors to Eigen's format.
auto eigenX = Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1>::Map(x.data(), x.size());
auto eigenB = Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1>::Map(b.data(), b.size());
Eigen::SparseLU<Eigen::SparseMatrix<storm::RationalFunction>, Eigen::COLAMDOrdering<int>> solver;
solver.compute(*eigenA);
solver._solve_impl(eigenB, eigenX);
// Translate the solution from Eigen's format into our representation.
for (uint64_t index = 0; index < eigenX.size(); ++index) {
x[index] = eigenX[index];
}
}
template<>
void EigenLinearEquationSolver<storm::RationalFunction>::performMatrixVectorMultiplication(std::vector<storm::RationalFunction>& x, std::vector<storm::RationalFunction> const* b, uint_fast64_t n, std::vector<storm::RationalFunction>* multiplyResult) const {
// Translate the vectors x and b into Eigen's format.
std::unique_ptr<Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1>> eigenB;
// Typedef the map-type so we don't have to spell it out.
typedef decltype(Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1>::Map(b->data(), b->size())) MapType;
// Map the input vectors x and b to Eigen's format.
std::unique_ptr<MapType> eigenB;
if (b != nullptr) {
eigenB = std::make_unique<Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1>>(b->size());
for (uint64_t index = 0; index < b->size(); ++index) {
(*eigenB)[index] = (*b)[index];
}
}
Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1> eigenX(x.size());
for (uint64_t index = 0; index < x.size(); ++index) {
eigenX[index] = x[index];
eigenB = std::make_unique<MapType>(Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1>::Map(b->data(), b->size()));
}
Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1> eigenX = Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1>::Map(x.data(), x.size());
// Perform n matrix-vector multiplications.
for (uint64_t iteration = 0; iteration < n; ++iteration) {
@ -280,11 +222,6 @@ namespace storm {
eigenX += *eigenB;
}
}
// Translate the solution from Eigen's format into our representation.
for (uint64_t index = 0; index < eigenX.size(); ++index) {
x[index] = eigenX[index];
}
}
template<typename ValueType>

Write Preview
Loading…
Cancel
Save