fixed some more problems with Eigen solver

Former-commit-id: c6ed18c4ab

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<ValueType, Eigen::Dynamic, 1>::Map(b->data(), b->size())) MapType;
             
+            bool multiplyResultProvided = multiplyResult != nullptr;
+            if (!multiplyResult) {
+                multiplyResult = new std::vector<ValueType>(eigenA->cols());
+            }
+            auto eigenMultiplyResult = Eigen::Matrix<ValueType, Eigen::Dynamic, 1>::Map(multiplyResult->data(), multiplyResult->size());
+            
             // Map the input vectors x and b to Eigen's format.
             std::unique_ptr<MapType> eigenB;
             if (b != nullptr) {
@@ -138,11 +144,23 @@ namespace storm {
             auto eigenX = Eigen::Matrix<ValueType, Eigen::Dynamic, 1>::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<storm::RationalNumber, Eigen::Dynamic, 1>::Map(b->data(), b->size())) MapType;
 
+            bool multiplyResultProvided = multiplyResult != nullptr;
+            if (!multiplyResult) {
+                multiplyResult = new std::vector<storm::RationalNumber>(eigenA->cols());
+            }
+            auto eigenMultiplyResult = Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1>::Map(multiplyResult->data(), multiplyResult->size());
+            
             // Map the input vectors x and b to Eigen's format.
             std::unique_ptr<MapType> eigenB;
             if (b != nullptr) {
                 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());
+            auto eigenX = Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1>::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<storm::RationalFunction, Eigen::Dynamic, 1>::Map(b->data(), b->size())) MapType;
 
+            bool multiplyResultProvided = multiplyResult != nullptr;
+            if (!multiplyResult) {
+                multiplyResult = new std::vector<storm::RationalFunction>(eigenA->cols());
+            }
+            auto eigenMultiplyResult = Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1>::Map(multiplyResult->data(), multiplyResult->size());
+            
             // Map the input vectors x and b to Eigen's format.
             std::unique_ptr<MapType> eigenB;
             if (b != nullptr) {
                 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());
+            auto eigenX = Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1>::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;
             }
         }
         

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;

test/functional/solver/EigenLinearEquationSolverTest.cpp

@@ -207,7 +207,7 @@ TEST(DISABLED_EigenLinearEquationSolver, BiCGSTAB_Diagonal) {
     ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
 }
 
-TEST(EigenLinearEquationSolver, MatrixVectorMultplication) {
+TEST(EigenLinearEquationSolver, MatrixVectorMultiplication) {
     ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder<double> builder);
     storm::storage::SparseMatrixBuilder<double> builder;
     ASSERT_NO_THROW(builder.addNextValue(0, 1, 0.5));

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<double> builder);
     storm::storage::SparseMatrixBuilder<double> builder;
     ASSERT_NO_THROW(builder.addNextValue(0, 1, 0.5));

test/functional/solver/GmmxxLinearEquationSolverTest.cpp

@@ -227,7 +227,7 @@ TEST(GmmxxLinearEquationSolver, gmresdiag) {
     ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
 }
 
-TEST(GmmxxLinearEquationSolver, MatrixVectorMultplication) {
+TEST(GmmxxLinearEquationSolver, MatrixVectorMultiplication) {
     ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder<double> builder);
     storm::storage::SparseMatrixBuilder<double> builder;
     ASSERT_NO_THROW(builder.addNextValue(0, 1, 0.5));

test/functional/solver/NativeLinearEquationSolverTest.cpp

@@ -34,7 +34,7 @@ TEST(NativeLinearEquationSolver, SolveWithStandardOptions) {
     ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
 }
 
-TEST(NativeLinearEquationSolver, MatrixVectorMultplication) {
+TEST(NativeLinearEquationSolver, MatrixVectorMultiplication) {
     ASSERT_NO_THROW(storm::storage::SparseMatrixBuilder<double> builder);
     storm::storage::SparseMatrixBuilder<double> builder;
     ASSERT_NO_THROW(builder.addNextValue(0, 1, 0.5));