From d24fb0cf9ae59e112128ba1d87a75bbcabd305de Mon Sep 17 00:00:00 2001
From: dehnert <dehnert@cs.rwth-aachen.de>
Date: Fri, 24 Jun 2016 19:53:32 +0200
Subject: [PATCH] avoid temporary in Eigen solver by providing .noalias().
 slightly rewrote matrix-vector expression to benefit more from Eigen's
 optimization capabilities

Former-commit-id: 838eac1449c9c7ec3c9a40e506f138f93072ac3d
---
 src/solver/EigenLinearEquationSolver.cpp | 23 ++++++++++++-----------
 1 file changed, 12 insertions(+), 11 deletions(-)

diff --git a/src/solver/EigenLinearEquationSolver.cpp b/src/solver/EigenLinearEquationSolver.cpp
index 7a9d8867c..34eda29cc 100644
--- a/src/solver/EigenLinearEquationSolver.cpp
+++ b/src/solver/EigenLinearEquationSolver.cpp
@@ -147,9 +147,10 @@ namespace storm {
             auto currentX = &eigenX;
             auto nextX = &eigenMultiplyResult;
             for (uint64_t iteration = 0; iteration < n; ++iteration) {
-                *nextX = *eigenA * *currentX;
-                if (eigenB != nullptr) {
-                    *nextX += *eigenB;
+                if (eigenB) {
+                    nextX->noalias() = *eigenA * *currentX + *eigenB;
+                } else {
+                    nextX->noalias() = *eigenA * *currentX;
                 }
                 std::swap(nextX, currentX);
             }
@@ -209,11 +210,11 @@ namespace storm {
             auto currentX = &eigenX;
             auto nextX = &eigenMultiplyResult;
             for (uint64_t iteration = 0; iteration < n; ++iteration) {
-                *nextX = *eigenA * *currentX;
-                if (eigenB != nullptr) {
-                    *nextX += *eigenB;
+                if (eigenB) {
+                    nextX->noalias() = *eigenA * *currentX + *eigenB;
+                } else {
+                    nextX->noalias() = *eigenA * *currentX;
                 }
-                std::swap(currentX, nextX);
             }
             
             // If the last result we obtained is not the one in the input vector x, we swap the result there.
@@ -261,11 +262,11 @@ namespace storm {
             auto currentX = &eigenX;
             auto nextX = &eigenMultiplyResult;
             for (uint64_t iteration = 0; iteration < n; ++iteration) {
-                *nextX = *eigenA * *currentX;
-                if (eigenB != nullptr) {
-                    *nextX += *eigenB;
+                if (eigenB) {
+                    nextX->noalias() = *eigenA * *currentX + *eigenB;
+                } else {
+                    nextX->noalias() = *eigenA * *currentX;
                 }
-                std::swap(nextX, currentX);
             }
             
             // If the last result we obtained is not the one in the input vector x, we swap the result there.