From 69da4ff147f824439080d413a109969d58f5f77d Mon Sep 17 00:00:00 2001 From: dehnert Date: Fri, 24 Jun 2016 19:44:04 +0200 Subject: [PATCH] fixed some more problems with Eigen solver Former-commit-id: c6ed18c4ab97eeb096a12fe5bd45fe63007f980e --- src/solver/EigenLinearEquationSolver.cpp | 70 ++++++++++++++++--- test/functional/builder/parametric_die.pm | 34 +++++++++ .../solver/EigenLinearEquationSolverTest.cpp | 2 +- .../EliminationLinearEquationSolverTest.cpp | 2 +- .../solver/GmmxxLinearEquationSolverTest.cpp | 2 +- .../solver/NativeLinearEquationSolverTest.cpp | 2 +- 6 files changed, 100 insertions(+), 12 deletions(-) create mode 100644 test/functional/builder/parametric_die.pm diff --git a/src/solver/EigenLinearEquationSolver.cpp b/src/solver/EigenLinearEquationSolver.cpp index 409c4427d..7a9d8867c 100644 --- a/src/solver/EigenLinearEquationSolver.cpp +++ b/src/solver/EigenLinearEquationSolver.cpp @@ -130,6 +130,12 @@ namespace storm { // Typedef the map-type so we don't have to spell it out. typedef decltype(Eigen::Matrix::Map(b->data(), b->size())) MapType; + bool multiplyResultProvided = multiplyResult != nullptr; + if (!multiplyResult) { + multiplyResult = new std::vector(eigenA->cols()); + } + auto eigenMultiplyResult = Eigen::Matrix::Map(multiplyResult->data(), multiplyResult->size()); + // Map the input vectors x and b to Eigen's format. std::unique_ptr eigenB; if (b != nullptr) { @@ -138,11 +144,23 @@ namespace storm { auto eigenX = Eigen::Matrix::Map(x.data(), x.size()); // Perform n matrix-vector multiplications. + auto currentX = &eigenX; + auto nextX = &eigenMultiplyResult; for (uint64_t iteration = 0; iteration < n; ++iteration) { - eigenX = *this->eigenA * eigenX; + *nextX = *eigenA * *currentX; if (eigenB != nullptr) { - eigenX += *eigenB; + *nextX += *eigenB; } + std::swap(nextX, currentX); + } + + // If the last result we obtained is not the one in the input vector x, we swap the result there. + if (currentX != &eigenX) { + std::swap(*nextX, *currentX); + } + + if (!multiplyResultProvided) { + delete multiplyResult; } } @@ -174,19 +192,37 @@ namespace storm { // Typedef the map-type so we don't have to spell it out. typedef decltype(Eigen::Matrix::Map(b->data(), b->size())) MapType; + bool multiplyResultProvided = multiplyResult != nullptr; + if (!multiplyResult) { + multiplyResult = new std::vector(eigenA->cols()); + } + auto eigenMultiplyResult = Eigen::Matrix::Map(multiplyResult->data(), multiplyResult->size()); + // Map the input vectors x and b to Eigen's format. std::unique_ptr eigenB; if (b != nullptr) { eigenB = std::make_unique(Eigen::Matrix::Map(b->data(), b->size())); } - Eigen::Matrix eigenX = Eigen::Matrix::Map(x.data(), x.size()); + auto eigenX = Eigen::Matrix::Map(x.data(), x.size()); // Perform n matrix-vector multiplications. + auto currentX = &eigenX; + auto nextX = &eigenMultiplyResult; for (uint64_t iteration = 0; iteration < n; ++iteration) { - eigenX = *eigenA * eigenX; + *nextX = *eigenA * *currentX; if (eigenB != nullptr) { - eigenX += *eigenB; + *nextX += *eigenB; } + std::swap(currentX, nextX); + } + + // If the last result we obtained is not the one in the input vector x, we swap the result there. + if (currentX != &eigenX) { + std::swap(*nextX, *currentX); + } + + if (!multiplyResultProvided) { + delete multiplyResult; } } @@ -208,19 +244,37 @@ namespace storm { // Typedef the map-type so we don't have to spell it out. typedef decltype(Eigen::Matrix::Map(b->data(), b->size())) MapType; + bool multiplyResultProvided = multiplyResult != nullptr; + if (!multiplyResult) { + multiplyResult = new std::vector(eigenA->cols()); + } + auto eigenMultiplyResult = Eigen::Matrix::Map(multiplyResult->data(), multiplyResult->size()); + // Map the input vectors x and b to Eigen's format. std::unique_ptr eigenB; if (b != nullptr) { eigenB = std::make_unique(Eigen::Matrix::Map(b->data(), b->size())); } - Eigen::Matrix eigenX = Eigen::Matrix::Map(x.data(), x.size()); + auto eigenX = Eigen::Matrix::Map(x.data(), x.size()); // Perform n matrix-vector multiplications. + auto currentX = &eigenX; + auto nextX = &eigenMultiplyResult; for (uint64_t iteration = 0; iteration < n; ++iteration) { - eigenX = *eigenA * eigenX; + *nextX = *eigenA * *currentX; if (eigenB != nullptr) { - eigenX += *eigenB; + *nextX += *eigenB; } + std::swap(nextX, currentX); + } + + // If the last result we obtained is not the one in the input vector x, we swap the result there. + if (currentX != &eigenX) { + std::swap(*nextX, *currentX); + } + + if (!multiplyResultProvided) { + delete multiplyResult; } } diff --git a/test/functional/builder/parametric_die.pm b/test/functional/builder/parametric_die.pm new file mode 100644 index 000000000..bd4e8b103 --- /dev/null +++ b/test/functional/builder/parametric_die.pm @@ -0,0 +1,34 @@ +// Knuth's model of a fair die using only fair coins +dtmc + +const double p; + +module die + + // local state + s : [0..7] init 0; + // value of the dice + d : [0..6] init 0; + + [] s=0 -> p : (s'=1) + (1-p) : (s'=2); + [] s=1 -> p : (s'=3) + (1-p) : (s'=4); + [] s=2 -> p : (s'=5) + (1-p) : (s'=6); + [] s=3 -> p : (s'=1) + (1-p) : (s'=7) & (d'=1); + [] s=4 -> p : (s'=7) & (d'=2) + (1-p) : (s'=7) & (d'=3); + [] s=5 -> p : (s'=7) & (d'=4) + (1-p) : (s'=7) & (d'=5); + [] s=6 -> p : (s'=2) + (1-p) : (s'=7) & (d'=6); + [] s=7 -> 1: (s'=7); + +endmodule + +rewards "coin_flips" + [] s<7 : 1; +endrewards + +label "one" = s=7&d=1; +label "two" = s=7&d=2; +label "three" = s=7&d=3; +label "four" = s=7&d=4; +label "five" = s=7&d=5; +label "six" = s=7&d=6; +label "done" = s=7; diff --git a/test/functional/solver/EigenLinearEquationSolverTest.cpp b/test/functional/solver/EigenLinearEquationSolverTest.cpp index 2bf575c2a..5bc1b4a42 100644 --- a/test/functional/solver/EigenLinearEquationSolverTest.cpp +++ b/test/functional/solver/EigenLinearEquationSolverTest.cpp @@ -207,7 +207,7 @@ TEST(DISABLED_EigenLinearEquationSolver, BiCGSTAB_Diagonal) { ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule().getPrecision()); } -TEST(EigenLinearEquationSolver, MatrixVectorMultplication) { +TEST(EigenLinearEquationSolver, MatrixVectorMultiplication) { ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder builder); storm::storage::SparseMatrixBuilder builder; ASSERT_NO_THROW(builder.addNextValue(0, 1, 0.5)); diff --git a/test/functional/solver/EliminationLinearEquationSolverTest.cpp b/test/functional/solver/EliminationLinearEquationSolverTest.cpp index 39007a5d7..7a964a1d6 100644 --- a/test/functional/solver/EliminationLinearEquationSolverTest.cpp +++ b/test/functional/solver/EliminationLinearEquationSolverTest.cpp @@ -34,7 +34,7 @@ TEST(EliminationLinearEquationSolver, Solve) { ASSERT_LT(std::abs(x[2] - (-1)), 1e-15); } -TEST(EliminationLinearEquationSolver, MatrixVectorMultplication) { +TEST(EliminationLinearEquationSolver, MatrixVectorMultiplication) { ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder builder); storm::storage::SparseMatrixBuilder builder; ASSERT_NO_THROW(builder.addNextValue(0, 1, 0.5)); diff --git a/test/functional/solver/GmmxxLinearEquationSolverTest.cpp b/test/functional/solver/GmmxxLinearEquationSolverTest.cpp index 181d3429a..44282f522 100644 --- a/test/functional/solver/GmmxxLinearEquationSolverTest.cpp +++ b/test/functional/solver/GmmxxLinearEquationSolverTest.cpp @@ -227,7 +227,7 @@ TEST(GmmxxLinearEquationSolver, gmresdiag) { ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule().getPrecision()); } -TEST(GmmxxLinearEquationSolver, MatrixVectorMultplication) { +TEST(GmmxxLinearEquationSolver, MatrixVectorMultiplication) { ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder builder); storm::storage::SparseMatrixBuilder builder; ASSERT_NO_THROW(builder.addNextValue(0, 1, 0.5)); diff --git a/test/functional/solver/NativeLinearEquationSolverTest.cpp b/test/functional/solver/NativeLinearEquationSolverTest.cpp index a1313a797..cfe61741c 100644 --- a/test/functional/solver/NativeLinearEquationSolverTest.cpp +++ b/test/functional/solver/NativeLinearEquationSolverTest.cpp @@ -34,7 +34,7 @@ TEST(NativeLinearEquationSolver, SolveWithStandardOptions) { ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule().getPrecision()); } -TEST(NativeLinearEquationSolver, MatrixVectorMultplication) { +TEST(NativeLinearEquationSolver, MatrixVectorMultiplication) { ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder builder); storm::storage::SparseMatrixBuilder builder; ASSERT_NO_THROW(builder.addNextValue(0, 1, 0.5));