diff --git a/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp b/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp
index 412c22625..4f24cde22 100644
--- a/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp
+++ b/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp
@@ -300,14 +300,14 @@ namespace storm {
                 storm::utility::vector::scaleVectorInPlace(result, std::get<3>(foxGlynnResult)[0]);
                 std::function<ValueType(ValueType const&, ValueType const&)> addAndScale = [&foxGlynnResult] (ValueType const& a, ValueType const& b) { return a + std::get<3>(foxGlynnResult)[0] * b; };
                 if (addVector != nullptr) {
-                    storm::utility::vector::applyPointwiseInPlace(result, *addVector, addAndScale);
+                    storm::utility::vector::applyPointwise(result, *addVector, result, addAndScale);
                 }
                 ++startingIteration;
             } else {
                 if (computeCumulativeReward) {
                     result = std::vector<ValueType>(values.size());
                     std::function<ValueType (ValueType const&)> scaleWithUniformizationRate = [&uniformizationRate] (ValueType const& a) -> ValueType { return a / uniformizationRate; };
-                    storm::utility::vector::applyPointwiseInPlace(result, scaleWithUniformizationRate);
+                    storm::utility::vector::applyPointwise(result, result, scaleWithUniformizationRate);
                 } else {
                     result = std::vector<ValueType>(values.size());
                 }
@@ -325,7 +325,7 @@ namespace storm {
                 // For the iterations below the left truncation point, we need to add and scale the result with the uniformization rate.
                 for (uint_fast64_t index = 1; index < startingIteration; ++index) {
                     solver->performMatrixVectorMultiplication(values, nullptr, 1, &multiplicationResult);
-                    storm::utility::vector::applyPointwiseInPlace(result, values, addAndScale);
+                    storm::utility::vector::applyPointwise(result, values, result, addAndScale);
                 }
             }
             
@@ -337,7 +337,7 @@ namespace storm {
                 solver->performMatrixVectorMultiplication(values, addVector, 1, &multiplicationResult);
                 
                 weight = std::get<3>(foxGlynnResult)[index - std::get<0>(foxGlynnResult)];
-                storm::utility::vector::applyPointwiseInPlace(result, values, addAndScale);
+                storm::utility::vector::applyPointwise(result, values, result, addAndScale);
             }
             
             return result;
@@ -409,7 +409,7 @@ namespace storm {
             if (this->getModel().hasTransitionRewards()) {
                 totalRewardVector = this->getModel().getTransitionMatrix().getPointwiseProductRowSumVector(this->getModel().getTransitionRewardMatrix());
                 if (this->getModel().hasStateRewards()) {
-                    storm::utility::vector::addVectorsInPlace(totalRewardVector, this->getModel().getStateRewardVector());
+                    storm::utility::vector::addVectors(totalRewardVector, this->getModel().getStateRewardVector(), totalRewardVector);
                 }
             } else {
                 totalRewardVector = std::vector<ValueType>(this->getModel().getStateRewardVector());
diff --git a/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp b/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp
index 04b36d823..abec6a929 100644
--- a/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp
+++ b/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp
@@ -103,22 +103,22 @@ namespace storm {
             for (uint_fast64_t currentStep = 0; currentStep < numberOfSteps; ++currentStep) {
                 // Start by (re-)computing bProbabilistic = bProbabilisticFixed + aProbabilisticToMarkovian * vMarkovian.
                 aProbabilisticToMarkovian.multiplyWithVector(markovianNonGoalValues, bProbabilistic);
-                storm::utility::vector::addVectorsInPlace(bProbabilistic, bProbabilisticFixed);
+                storm::utility::vector::addVectors(bProbabilistic, bProbabilisticFixed, bProbabilistic);
                 
                 // Now perform the inner value iteration for probabilistic states.
                 solver->solveEquationSystem(min, probabilisticNonGoalValues, bProbabilistic, &multiplicationResultScratchMemory, &aProbabilisticScratchMemory);
                 
                 // (Re-)compute bMarkovian = bMarkovianFixed + aMarkovianToProbabilistic * vProbabilistic.
                 aMarkovianToProbabilistic.multiplyWithVector(probabilisticNonGoalValues, bMarkovian);
-                storm::utility::vector::addVectorsInPlace(bMarkovian, bMarkovianFixed);
+                storm::utility::vector::addVectors(bMarkovian, bMarkovianFixed, bMarkovian);
                 aMarkovian.multiplyWithVector(markovianNonGoalValues, markovianNonGoalValuesSwap);
                 std::swap(markovianNonGoalValues, markovianNonGoalValuesSwap);
-                storm::utility::vector::addVectorsInPlace(markovianNonGoalValues, bMarkovian);
+                storm::utility::vector::addVectors(markovianNonGoalValues, bMarkovian, markovianNonGoalValues);
             }
             
             // After the loop, perform one more step of the value iteration for PS states.
             aProbabilisticToMarkovian.multiplyWithVector(markovianNonGoalValues, bProbabilistic);
-            storm::utility::vector::addVectorsInPlace(bProbabilistic, bProbabilisticFixed);
+            storm::utility::vector::addVectors(bProbabilistic, bProbabilisticFixed, bProbabilistic);
             solver->solveEquationSystem(min, probabilisticNonGoalValues, bProbabilistic, &multiplicationResultScratchMemory, &aProbabilisticScratchMemory);
         }
         
@@ -217,7 +217,7 @@ namespace storm {
             if (model.hasTransitionRewards()) {
                 totalRewardVector = model.getTransitionMatrix().getPointwiseProductRowSumVector(model.getTransitionRewardMatrix());
                 if (model.hasStateRewards()) {
-                    storm::utility::vector::addVectorsInPlace(totalRewardVector, model.getStateRewardVector());
+                    storm::utility::vector::addVectors(totalRewardVector, model.getStateRewardVector(), totalRewardVector);
                 }
             } else {
                 totalRewardVector = std::vector<ValueType>(model.getStateRewardVector());
diff --git a/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.cpp
index d71d206a4..b927031e0 100644
--- a/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.cpp
+++ b/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.cpp
@@ -69,18 +69,26 @@ namespace storm {
                     submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs());
                     submatrix = (model.getRowColumnIdentity() * maybeStatesAdd) - submatrix;
                     
+                    submatrix.exportToDot("submatrix.dot");
+                    
                     // Create the solution vector.
                     std::vector<ValueType> x(maybeStates.getNonZeroCount(), ValueType(0.5));
                     
                     // Translate the symbolic matrix/vector to their explicit representations and solve the equation system.
+                    STORM_LOG_DEBUG("Converting the symbolic matrix to a sparse matrix.");
                     storm::storage::SparseMatrix<ValueType> explicitSubmatrix = submatrix.toMatrix(odd, odd);
+                    
+                    STORM_LOG_DEBUG("Converting the symbolic vector to an explicit vector.");
                     std::vector<ValueType> b = subvector.template toVector<ValueType>(odd);
                     
+                    STORM_LOG_DEBUG("Solving explicit linear equation system.");
                     std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(explicitSubmatrix);
                     solver->solveEquationSystem(x, b);
                     
                     // Now that we have the explicit solution of the system, we need to transform it to a symbolic representation.
-                    storm::dd::Add<DdType> numericResult; // = storm::dd::Add<DdType>(x, odd);
+                    STORM_LOG_DEBUG("Converting the explicit result to a symbolic form.");
+                    storm::dd::Add<DdType> numericResult(model.getManager().asSharedPointer(), x, odd, model.getRowVariables());
+                    
                     return statesWithProbability01.second.toAdd() + numericResult;
                 } else {
                     return statesWithProbability01.second.toAdd();
diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp
index 1b543b5fc..46ba68183 100644
--- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp
+++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp
@@ -171,7 +171,7 @@ namespace storm {
             if (this->getModel().hasTransitionRewards()) {
                 totalRewardVector = this->getModel().getTransitionMatrix().getPointwiseProductRowSumVector(this->getModel().getTransitionRewardMatrix());
                 if (this->getModel().hasStateRewards()) {
-                    storm::utility::vector::addVectorsInPlace(totalRewardVector, this->getModel().getStateRewardVector());
+                    storm::utility::vector::addVectors(totalRewardVector, this->getModel().getStateRewardVector(), totalRewardVector);
                 }
             } else {
                 totalRewardVector = std::vector<ValueType>(this->getModel().getStateRewardVector());
@@ -272,7 +272,7 @@ namespace storm {
                         // first.
                         std::vector<ValueType> subStateRewards(b.size());
                         storm::utility::vector::selectVectorValues(subStateRewards, maybeStates, stateRewardVector.get());
-                        storm::utility::vector::addVectorsInPlace(b, subStateRewards);
+                        storm::utility::vector::addVectors(b, subStateRewards, b);
                     }
                 } else {
                     // If only a state-based reward model is  available, we take this vector as the
diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp
index 516090698..664a15fdd 100644
--- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp
+++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp
@@ -183,7 +183,7 @@ namespace storm {
             if (this->getModel().hasTransitionRewards()) {
                 totalRewardVector = this->getModel().getTransitionMatrix().getPointwiseProductRowSumVector(this->getModel().getTransitionRewardMatrix());
                 if (this->getModel().hasStateRewards()) {
-                    storm::utility::vector::addVectorsInPlace(totalRewardVector, this->getModel().getStateRewardVector());
+                    storm::utility::vector::addVectors(totalRewardVector, this->getModel().getStateRewardVector(), totalRewardVector);
                 }
             } else {
                 totalRewardVector = std::vector<ValueType>(this->getModel().getStateRewardVector());
@@ -285,7 +285,7 @@ namespace storm {
                         // first.
                         std::vector<ValueType> subStateRewards(b.size());
                         storm::utility::vector::selectVectorValuesRepeatedly(subStateRewards, maybeStates, transitionMatrix.getRowGroupIndices(), stateRewardVector.get());
-                        storm::utility::vector::addVectorsInPlace(b, subStateRewards);
+                        storm::utility::vector::addVectors(b, subStateRewards, b);
                     }
                 } else {
                     // If only a state-based reward model is  available, we take this vector as the
diff --git a/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp b/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp
index 72dfe070e..e4087ad52 100644
--- a/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp
+++ b/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp
@@ -176,7 +176,7 @@ namespace storm {
                     // first.
                     std::vector<ValueType> subStateRewards(stateRewards.size());
                     storm::utility::vector::selectVectorValues(subStateRewards, maybeStates, model.getStateRewardVector());
-                    storm::utility::vector::addVectorsInPlace(stateRewards, subStateRewards);
+                    storm::utility::vector::addVectors(stateRewards, subStateRewards, stateRewards);
                 }
             } else {
                 // If only a state-based reward model is  available, we take this vector as the
diff --git a/src/modelchecker/results/SymbolicQuantitativeCheckResult.cpp b/src/modelchecker/results/SymbolicQuantitativeCheckResult.cpp
index c265efd5c..fc18df5c6 100644
--- a/src/modelchecker/results/SymbolicQuantitativeCheckResult.cpp
+++ b/src/modelchecker/results/SymbolicQuantitativeCheckResult.cpp
@@ -49,14 +49,18 @@ namespace storm {
         template<storm::dd::DdType Type>
         std::ostream& SymbolicQuantitativeCheckResult<Type>::writeToStream(std::ostream& out) const {
             out << "[";
-            bool first = true;
-            for (auto valuationValuePair : this->values) {
-                if (!first) {
-                    out << ", ";
-                } else {
-                    first = false;
+            if (this->values.isZero()) {
+                out << "0";
+            } else {
+                bool first = true;
+                for (auto valuationValuePair : this->values) {
+                    if (!first) {
+                        out << ", ";
+                    } else {
+                        first = false;
+                    }
+                    out << valuationValuePair.second;
                 }
-                out << valuationValuePair.second;
             }
             out << "]";
             return out;
diff --git a/src/models/sparse/Model.cpp b/src/models/sparse/Model.cpp
index 6541d517f..e1230d051 100644
--- a/src/models/sparse/Model.cpp
+++ b/src/models/sparse/Model.cpp
@@ -142,7 +142,7 @@ namespace storm {
             void Model<ValueType>::convertTransitionRewardsToStateRewards() {
                 STORM_LOG_THROW(this->hasTransitionRewards(), storm::exceptions::InvalidOperationException, "Cannot reduce non-existant transition rewards to state rewards.");
                 if (this->hasStateRewards()) {
-                    storm::utility::vector::addVectorsInPlace(stateRewardVector.get(), transitionMatrix.getPointwiseProductRowSumVector(transitionRewardMatrix.get()));
+                    storm::utility::vector::addVectors(stateRewardVector.get(), transitionMatrix.getPointwiseProductRowSumVector(transitionRewardMatrix.get()), stateRewardVector.get());
                 } else {
                     this->stateRewardVector = transitionMatrix.getPointwiseProductRowSumVector(transitionRewardMatrix.get());
                 }
diff --git a/src/solver/GmmxxLinearEquationSolver.cpp b/src/solver/GmmxxLinearEquationSolver.cpp
index 0bc5a3c46..789ee6e55 100644
--- a/src/solver/GmmxxLinearEquationSolver.cpp
+++ b/src/solver/GmmxxLinearEquationSolver.cpp
@@ -148,10 +148,8 @@ namespace storm {
         template<typename ValueType>
         uint_fast64_t GmmxxLinearEquationSolver<ValueType>::solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
             // Get a Jacobi decomposition of the matrix A.
-            std::pair<storm::storage::SparseMatrix<ValueType>, storm::storage::SparseMatrix<ValueType>> jacobiDecomposition = A.getJacobiDecomposition();
+            std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> jacobiDecomposition = A.getJacobiDecomposition();
             
-            // Convert the (inverted) diagonal matrix to gmm++'s format.
-            std::unique_ptr<gmm::csr_matrix<ValueType>> gmmDinv = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(std::move(jacobiDecomposition.second));
             // Convert the LU matrix to gmm++'s format.
             std::unique_ptr<gmm::csr_matrix<ValueType>> gmmLU = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(std::move(jacobiDecomposition.first));
         
@@ -176,7 +174,7 @@ namespace storm {
                 // Compute D^-1 * (b - LU * x) and store result in nextX.
                 gmm::mult(*gmmLU, *currentX, tmpX);
                 gmm::add(b, gmm::scaled(tmpX, -storm::utility::one<ValueType>()), tmpX);
-                gmm::mult(*gmmDinv, tmpX, *nextX);
+                storm::utility::vector::multiplyVectorsPointwise(jacobiDecomposition.second, tmpX, *nextX);
                 
                 // Now check if the process already converged within our precision.
                 converged = storm::utility::vector::equalModuloPrecision(*currentX, *nextX, precision, relative);
diff --git a/src/solver/NativeLinearEquationSolver.cpp b/src/solver/NativeLinearEquationSolver.cpp
index 47c72dae8..9206ae301 100644
--- a/src/solver/NativeLinearEquationSolver.cpp
+++ b/src/solver/NativeLinearEquationSolver.cpp
@@ -34,7 +34,9 @@ namespace storm {
         template<typename ValueType>
         void NativeLinearEquationSolver<ValueType>::solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
             // Get a Jacobi decomposition of the matrix A.
-            std::pair<storm::storage::SparseMatrix<ValueType>, storm::storage::SparseMatrix<ValueType>> jacobiDecomposition = A.getJacobiDecomposition();
+            std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> jacobiDecomposition = A.getJacobiDecomposition();
+            
+            std::cout << "LU has " << jacobiDecomposition.first.getNonzeroEntryCount() << " nonzeros." << std::endl;
             
             // To avoid copying the contents of the vector in the loop, we create a temporary x to swap with.
             bool multiplyResultProvided = true;
@@ -56,9 +58,8 @@ namespace storm {
             while (!converged && iterationCount < maximalNumberOfIterations) {
                 // Compute D^-1 * (b - LU * x) and store result in nextX.
                 jacobiDecomposition.first.multiplyWithVector(*currentX, tmpX);
-                storm::utility::vector::scaleVectorInPlace(tmpX, -storm::utility::one<ValueType>());
-                storm::utility::vector::addVectorsInPlace(tmpX, b);
-                jacobiDecomposition.second.multiplyWithVector(tmpX, *nextX);
+                storm::utility::vector::subtractVectors(b, tmpX, tmpX);
+                storm::utility::vector::multiplyVectorsPointwise(jacobiDecomposition.second, tmpX, *nextX);
                 
                 // Swap the two pointers as a preparation for the next iteration.
                 std::swap(nextX, currentX);
@@ -103,7 +104,7 @@ namespace storm {
                 
                 // If requested, add an offset to the current result vector.
                 if (b != nullptr) {
-                    storm::utility::vector::addVectorsInPlace(*currentX, *b);
+                    storm::utility::vector::addVectors(*currentX, *b, *currentX);
                 }
             }
             
diff --git a/src/solver/NativeMinMaxLinearEquationSolver.cpp b/src/solver/NativeMinMaxLinearEquationSolver.cpp
index 48ef291c7..51da3e692 100644
--- a/src/solver/NativeMinMaxLinearEquationSolver.cpp
+++ b/src/solver/NativeMinMaxLinearEquationSolver.cpp
@@ -50,7 +50,7 @@ namespace storm {
             while (!converged && iterations < maximalNumberOfIterations) {
                 // Compute x' = A*x + b.
                 A.multiplyWithVector(*currentX, *multiplyResult);
-                storm::utility::vector::addVectorsInPlace(*multiplyResult, b);
+                storm::utility::vector::addVectors(*multiplyResult, b, *multiplyResult);
                 
                 // Reduce the vector x' by applying min/max for all non-deterministic choices as given by the topmost
                 // element of the min/max operator stack.
@@ -106,7 +106,7 @@ namespace storm {
                 
                 // Add b if it is non-null.
                 if (b != nullptr) {
-                    storm::utility::vector::addVectorsInPlace(*multiplyResult, *b);
+                    storm::utility::vector::addVectors(*multiplyResult, *b, *multiplyResult);
                 }
                 
                 // Reduce the vector x' by applying min/max for all non-deterministic choices as given by the topmost
diff --git a/src/solver/TopologicalMinMaxLinearEquationSolver.cpp b/src/solver/TopologicalMinMaxLinearEquationSolver.cpp
index 66de10ed9..8ea3478b0 100644
--- a/src/solver/TopologicalMinMaxLinearEquationSolver.cpp
+++ b/src/solver/TopologicalMinMaxLinearEquationSolver.cpp
@@ -254,7 +254,7 @@ namespace storm {
 						while (!converged && localIterations < this->maximalNumberOfIterations) {
 							// Compute x' = A*x + b.
 							sccSubmatrix.multiplyWithVector(*currentX, sccMultiplyResult);
-							storm::utility::vector::addVectorsInPlace<ValueType>(sccMultiplyResult, sccSubB);
+							storm::utility::vector::addVectors<ValueType>(sccMultiplyResult, sccSubB, sccMultiplyResult);
 
 							//A.multiplyWithVector(scc, nondeterministicChoiceIndices, *currentX, multiplyResult);
 							//storm::utility::addVectors(scc, nondeterministicChoiceIndices, multiplyResult, b);
diff --git a/src/storage/SparseMatrix.cpp b/src/storage/SparseMatrix.cpp
index 335e79910..9a8affbd1 100644
--- a/src/storage/SparseMatrix.cpp
+++ b/src/storage/SparseMatrix.cpp
@@ -11,6 +11,7 @@
 #include "src/adapters/CarlAdapter.h"
 
 #include "src/exceptions/InvalidStateException.h"
+#include "src/exceptions/NotImplementedException.h"
 #include "src/utility/macros.h"
 
 #include "log4cplus/logger.h"
@@ -701,40 +702,31 @@ namespace storm {
         }
         
         template<typename ValueType>
-        typename std::pair<storm::storage::SparseMatrix<ValueType>, storm::storage::SparseMatrix<ValueType>> SparseMatrix<ValueType>::getJacobiDecomposition() const {
-            if (rowCount != columnCount) {
-                throw storm::exceptions::InvalidArgumentException() << "Illegal call to SparseMatrix::invertDiagonal: matrix is non-square.";
-            }
-            storm::storage::SparseMatrix<ValueType> resultLU(*this);
-            resultLU.deleteDiagonalEntries();
+        typename std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> SparseMatrix<ValueType>::getJacobiDecomposition() const {
+            STORM_LOG_THROW(this->getRowCount() == this->getColumnCount(), storm::exceptions::InvalidArgumentException, "Canno compute Jacobi decomposition of non-square matrix.");
             
-            SparseMatrixBuilder<ValueType> dInvBuilder(rowCount, columnCount, rowCount);
+            // Prepare the resulting data structures.
+            SparseMatrixBuilder<ValueType> luBuilder(this->getRowCount(), this->getColumnCount());
+            std::vector<ValueType> invertedDiagonal(rowCount);
             
             // Copy entries to the appropriate matrices.
             for (index_type rowNumber = 0; rowNumber < rowCount; ++rowNumber) {
-                
-                // Because the matrix may have several entries on the diagonal, we need to sum them before we are able
-                // to invert the entry.
-                ValueType diagonalValue = storm::utility::zero<ValueType>();
                 for (const_iterator it = this->begin(rowNumber), ite = this->end(rowNumber); it != ite; ++it) {
                     if (it->getColumn() == rowNumber) {
-                        diagonalValue += it->getValue();
-                    } else if (it->getColumn() > rowNumber) {
-                        break;
+                        invertedDiagonal[rowNumber] = storm::utility::one<ValueType>() / it->getValue();
+                    } else {
+                        luBuilder.addNextValue(rowNumber, it->getColumn(), it->getValue());
                     }
                 }
-                dInvBuilder.addNextValue(rowNumber, rowNumber, storm::utility::one<ValueType>() / diagonalValue);
             }
             
-            return std::make_pair(std::move(resultLU), dInvBuilder.build());
+            return std::make_pair(luBuilder.build(), std::move(invertedDiagonal));
         }
         
 #ifdef STORM_HAVE_CARL
         template<>
-        typename std::pair<storm::storage::SparseMatrix<RationalFunction>, storm::storage::SparseMatrix<RationalFunction>> SparseMatrix<RationalFunction>::getJacobiDecomposition() const {
-            // NOT SUPPORTED
-            // TODO do whatever storm does in such cases.
-            assert(false);
+        typename std::pair<storm::storage::SparseMatrix<RationalFunction>, std::vector<RationalFunction>> SparseMatrix<RationalFunction>::getJacobiDecomposition() const {
+            STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "This operation is not supported.");
         }
 #endif
         
diff --git a/src/storage/SparseMatrix.h b/src/storage/SparseMatrix.h
index 1bac30a6c..8eded966b 100644
--- a/src/storage/SparseMatrix.h
+++ b/src/storage/SparseMatrix.h
@@ -605,9 +605,9 @@ namespace storm {
             /*!
              * Calculates the Jacobi decomposition of this sparse matrix. For this operation, the matrix must be square.
              *
-             * @return A pair (L+U, D^-1) containing the matrix L+U and the inverted diagonal matrix D^-1.
+             * @return A pair (L+U, D^-1) containing the matrix L+U and the inverted diagonal D^-1 (as a vector).
              */
-            std::pair<storm::storage::SparseMatrix<value_type>, storm::storage::SparseMatrix<value_type>> getJacobiDecomposition() const;
+            std::pair<storm::storage::SparseMatrix<value_type>, std::vector<value_type>> getJacobiDecomposition() const;
             
             /*!
              * Performs a pointwise matrix multiplication of the matrix with the given matrix and returns a vector
diff --git a/src/storage/dd/CuddAdd.cpp b/src/storage/dd/CuddAdd.cpp
index 72787075b..71ec05758 100644
--- a/src/storage/dd/CuddAdd.cpp
+++ b/src/storage/dd/CuddAdd.cpp
@@ -18,6 +18,10 @@ namespace storm {
             // Intentionally left empty.
         }
         
+        Add<DdType::CUDD>::Add(std::shared_ptr<DdManager<DdType::CUDD> const> ddManager, std::vector<double> const& values, Odd<DdType::CUDD> const& odd, std::set<storm::expressions::Variable> const& metaVariables) : Dd<DdType::CUDD>(ddManager, metaVariables) {
+            cuddAdd = fromVector(ddManager, values, odd, metaVariables);
+        }
+        
         Bdd<DdType::CUDD> Add<DdType::CUDD>::toBdd() const {
             return Bdd<DdType::CUDD>(this->getDdManager(), this->getCuddAdd().BddPattern(), this->getContainedMetaVariables());
         }
@@ -155,19 +159,19 @@ namespace storm {
             std::set_union(this->getContainedMetaVariables().begin(), this->getContainedMetaVariables().end(), other.getContainedMetaVariables().begin(), other.getContainedMetaVariables().end(), std::inserter(metaVariables, metaVariables.begin()));
             return Add<DdType::CUDD>(this->getDdManager(), this->getCuddAdd().Pow(other.getCuddAdd()), metaVariables);
         }
-
+        
         Add<DdType::CUDD> Add<DdType::CUDD>::mod(Add<DdType::CUDD> const& other) const {
             std::set<storm::expressions::Variable> metaVariables;
             std::set_union(this->getContainedMetaVariables().begin(), this->getContainedMetaVariables().end(), other.getContainedMetaVariables().begin(), other.getContainedMetaVariables().end(), std::inserter(metaVariables, metaVariables.begin()));
             return Add<DdType::CUDD>(this->getDdManager(), this->getCuddAdd().Mod(other.getCuddAdd()), metaVariables);
         }
-
+        
         Add<DdType::CUDD> Add<DdType::CUDD>::logxy(Add<DdType::CUDD> const& other) const {
             std::set<storm::expressions::Variable> metaVariables;
             std::set_union(this->getContainedMetaVariables().begin(), this->getContainedMetaVariables().end(), other.getContainedMetaVariables().begin(), other.getContainedMetaVariables().end(), std::inserter(metaVariables, metaVariables.begin()));
             return Add<DdType::CUDD>(this->getDdManager(), this->getCuddAdd().LogXY(other.getCuddAdd()), metaVariables);
         }
-
+        
         Add<DdType::CUDD> Add<DdType::CUDD>::floor() const {
             return Add<DdType::CUDD>(this->getDdManager(), this->getCuddAdd().Floor(), this->getContainedMetaVariables());
         }
@@ -464,17 +468,20 @@ namespace storm {
                     ddRowVariableIndices.push_back(ddVariable.getIndex());
                 }
             }
+            std::sort(ddRowVariableIndices.begin(), ddRowVariableIndices.end());
+            
             for (auto const& variable : columnMetaVariables) {
                 DdMetaVariable<DdType::CUDD> const& metaVariable = this->getDdManager()->getMetaVariable(variable);
                 for (auto const& ddVariable : metaVariable.getDdVariables()) {
                     ddColumnVariableIndices.push_back(ddVariable.getIndex());
                 }
             }
+            std::sort(ddColumnVariableIndices.begin(), ddColumnVariableIndices.end());
             
             // Prepare the vectors that represent the matrix.
             std::vector<uint_fast64_t> rowIndications(rowOdd.getTotalOffset() + 1);
             std::vector<storm::storage::MatrixEntry<uint_fast64_t, double>> columnsAndValues(this->getNonZeroCount());
-
+            
             // Create a trivial row grouping.
             std::vector<uint_fast64_t> trivialRowGroupIndices(rowIndications.size());
             uint_fast64_t i = 0;
@@ -526,6 +533,7 @@ namespace storm {
                 }
                 rowAndColumnMetaVariables.insert(variable);
             }
+            std::sort(ddRowVariableIndices.begin(), ddRowVariableIndices.end());
             for (auto const& variable : columnMetaVariables) {
                 DdMetaVariable<DdType::CUDD> const& metaVariable = this->getDdManager()->getMetaVariable(variable);
                 for (auto const& ddVariable : metaVariable.getDdVariables()) {
@@ -533,12 +541,14 @@ namespace storm {
                 }
                 rowAndColumnMetaVariables.insert(variable);
             }
+            std::sort(ddColumnVariableIndices.begin(), ddColumnVariableIndices.end());
             for (auto const& variable : groupMetaVariables) {
                 DdMetaVariable<DdType::CUDD> const& metaVariable = this->getDdManager()->getMetaVariable(variable);
                 for (auto const& ddVariable : metaVariable.getDdVariables()) {
                     ddGroupVariableIndices.push_back(ddVariable.getIndex());
                 }
             }
+            std::sort(ddGroupVariableIndices.begin(), ddGroupVariableIndices.end());
             
             // TODO: assert that the group variables are at the very top of the variable ordering?
             
@@ -705,7 +715,67 @@ namespace storm {
                 addToVectorRec(Cudd_T(dd), currentLevel + 1, maxLevel, currentOffset + odd.getElseOffset(), odd.getThenSuccessor(), ddVariableIndices, targetVector);
             }
         }
-
+        
+        template<typename ValueType>
+        ADD Add<DdType::CUDD>::fromVector(std::shared_ptr<DdManager<DdType::CUDD> const> ddManager, std::vector<ValueType> const& values, Odd<DdType::CUDD> const& odd, std::set<storm::expressions::Variable> const& metaVariables) {
+            std::vector<uint_fast64_t> ddVariableIndices = getSortedVariableIndices(*ddManager, metaVariables);
+            uint_fast64_t offset = 0;
+            return ADD(ddManager->getCuddManager(), fromVectorRec(ddManager->getCuddManager().getManager(), offset, 0, ddVariableIndices.size(), values, odd, ddVariableIndices));
+        }
+        
+        template<typename ValueType>
+        DdNode* Add<DdType::CUDD>::fromVectorRec(::DdManager* manager, uint_fast64_t& currentOffset, uint_fast64_t currentLevel, uint_fast64_t maxLevel, std::vector<ValueType> const& values, Odd<DdType::CUDD> const& odd, std::vector<uint_fast64_t> const& ddVariableIndices) {
+            if (currentLevel == maxLevel) {
+                // If we are in a terminal node of the ODD, we need to check whether the then-offset of the ODD is one
+                // (meaning the encoding is a valid one) or zero (meaning the encoding is not valid). Consequently, we
+                // need to copy the next value of the vector iff the then-offset is greater than zero.
+                if (odd.getThenOffset() > 0) {
+                    return Cudd_addConst(manager, values[currentOffset++]);
+                } else {
+                    return Cudd_ReadZero(manager);
+                }
+            } else {
+                // If the total offset is zero, we can just return the constant zero DD.
+                if (odd.getThenOffset() + odd.getElseOffset() == 0) {
+                    return Cudd_ReadZero(manager);
+                }
+                
+                // Determine the new else-successor.
+                DdNode* elseSuccessor = nullptr;
+                if (odd.getElseOffset() > 0) {
+                    elseSuccessor = fromVectorRec(manager, currentOffset, currentLevel + 1, maxLevel, values, odd.getElseSuccessor(), ddVariableIndices);
+                } else {
+                    elseSuccessor = Cudd_ReadZero(manager);
+                }
+                Cudd_Ref(elseSuccessor);
+                
+                // Determine the new then-successor.
+                DdNode* thenSuccessor = nullptr;
+                if (odd.getThenOffset() > 0) {
+                    thenSuccessor = fromVectorRec(manager, currentOffset, currentLevel + 1, maxLevel, values, odd.getThenSuccessor(), ddVariableIndices);
+                } else {
+                    thenSuccessor = Cudd_ReadZero(manager);
+                }
+                Cudd_Ref(thenSuccessor);
+                
+                // Create a node representing ITE(currentVar, thenSuccessor, elseSuccessor);
+                DdNode* result = Cudd_addIthVar(manager, static_cast<int>(ddVariableIndices[currentLevel]));
+                Cudd_Ref(result);
+                DdNode* newResult = Cudd_addIte(manager, result, thenSuccessor, elseSuccessor);
+                Cudd_Ref(newResult);
+                
+                // Dispose of the intermediate results
+                Cudd_RecursiveDeref(manager, result);
+                Cudd_RecursiveDeref(manager, thenSuccessor);
+                Cudd_RecursiveDeref(manager, elseSuccessor);
+                
+                // Before returning, we remove the protection imposed by the previous call to Cudd_Ref.
+                Cudd_Deref(newResult);
+                
+                return newResult;
+            }
+        }
+        
         void Add<DdType::CUDD>::exportToDot(std::string const& filename) const {
             if (filename.empty()) {
                 std::vector<ADD> cuddAddVector = { this->getCuddAdd() };
@@ -751,7 +821,7 @@ namespace storm {
         DdForwardIterator<DdType::CUDD> Add<DdType::CUDD>::end(bool enumerateDontCareMetaVariables) const {
             return DdForwardIterator<DdType::CUDD>(this->getDdManager(), nullptr, nullptr, 0, true, nullptr, enumerateDontCareMetaVariables);
         }
-                
+        
         std::ostream& operator<<(std::ostream& out, const Add<DdType::CUDD>& add) {
             add.exportToDot();
             return out;
diff --git a/src/storage/dd/CuddAdd.h b/src/storage/dd/CuddAdd.h
index 07cfcc20e..65de2dd10 100644
--- a/src/storage/dd/CuddAdd.h
+++ b/src/storage/dd/CuddAdd.h
@@ -27,6 +27,16 @@ namespace storm {
             friend class Bdd<DdType::CUDD>;
             friend class Odd<DdType::CUDD>;
             
+            /*!
+             * Creates an ADD from the given explicit vector.
+             *
+             * @param ddManager The manager responsible for this DD.
+             * @param values The vector that is to be represented by the ADD.
+             * @param odd An ODD that is used to do the translation.
+             * @param metaVariables The meta variables to use to encode the vector.
+             */
+            Add(std::shared_ptr<DdManager<DdType::CUDD> const> ddManager, std::vector<double> const& values, Odd<DdType::CUDD> const& odd, std::set<storm::expressions::Variable> const& metaVariables);
+            
             // Instantiate all copy/move constructors/assignments with the default implementation.
             Add() = default;
             Add(Add<DdType::CUDD> const& other) = default;
@@ -675,6 +685,33 @@ namespace storm {
             template<typename ValueType>
             void addToVectorRec(DdNode const* dd, uint_fast64_t currentLevel, uint_fast64_t maxLevel, uint_fast64_t currentOffset, Odd<DdType::CUDD> const& odd, std::vector<uint_fast64_t> const& ddVariableIndices, std::vector<ValueType>& targetVector) const;
             
+            /*!
+             * Builds an ADD representing the given vector.
+             *
+             * @param ddManager The manager responsible for the ADD.
+             * @param values The vector that is to be represented by the ADD.
+             * @param odd The ODD used for the translation.
+             * @param metaVariables The meta variables used for the translation.
+             * @return The resulting (CUDD) ADD.
+             */
+            template<typename ValueType>
+            static ADD fromVector(std::shared_ptr<DdManager<DdType::CUDD> const> ddManager, std::vector<ValueType> const& values, Odd<DdType::CUDD> const& odd, std::set<storm::expressions::Variable> const& metaVariables);
+            
+            /*!
+             * Builds an ADD representing the given vector.
+             *
+             * @param manager The manager responsible for the ADD.
+             * @param currentOffset The current offset in the vector.
+             * @param currentLevel The current level in the DD.
+             * @param maxLevel The maximal level in the DD.
+             * @param values The vector that is to be represented by the ADD.
+             * @param odd The ODD used for the translation.
+             * @param ddVariableIndices The (sorted) list of DD variable indices to use.
+             * @return The resulting (CUDD) ADD node.
+             */
+            template<typename ValueType>
+            static DdNode* fromVectorRec(::DdManager* manager, uint_fast64_t& currentOffset, uint_fast64_t currentLevel, uint_fast64_t maxLevel, std::vector<ValueType> const& values, Odd<DdType::CUDD> const& odd, std::vector<uint_fast64_t> const& ddVariableIndices);
+            
             // The ADD created by CUDD.
             ADD cuddAdd;
         };
diff --git a/src/storage/dd/CuddDd.cpp b/src/storage/dd/CuddDd.cpp
index 43e2afd3f..d56589824 100644
--- a/src/storage/dd/CuddDd.cpp
+++ b/src/storage/dd/CuddDd.cpp
@@ -50,9 +50,13 @@ namespace storm {
         }
         
         std::vector<uint_fast64_t> Dd<DdType::CUDD>::getSortedVariableIndices() const {
+            return getSortedVariableIndices(*this->getDdManager(), this->getContainedMetaVariables());
+        }
+        
+        std::vector<uint_fast64_t> Dd<DdType::CUDD>::getSortedVariableIndices(DdManager<DdType::CUDD> const& manager, std::set<storm::expressions::Variable> const& metaVariables) {
             std::vector<uint_fast64_t> ddVariableIndices;
-            for (auto const& metaVariableName : this->getContainedMetaVariables()) {
-                auto const& metaVariable = this->getDdManager()->getMetaVariable(metaVariableName);
+            for (auto const& metaVariableName : metaVariables) {
+                auto const& metaVariable = manager.getMetaVariable(metaVariableName);
                 for (auto const& ddVariable : metaVariable.getDdVariables()) {
                     ddVariableIndices.push_back(ddVariable.getIndex());
                 }
diff --git a/src/storage/dd/CuddDd.h b/src/storage/dd/CuddDd.h
index fdeb085ef..b80bfc671 100644
--- a/src/storage/dd/CuddDd.h
+++ b/src/storage/dd/CuddDd.h
@@ -105,6 +105,8 @@ namespace storm {
              */
             std::shared_ptr<DdManager<DdType::CUDD> const> getDdManager() const;
             
+        protected:
+            
             /*!
              * Retrieves the (sorted) list of the variable indices of DD variables contained in this DD.
              *
@@ -112,7 +114,14 @@ namespace storm {
              */
             std::vector<uint_fast64_t> getSortedVariableIndices() const;
             
-        protected:
+            /*!
+             * Retrieves the (sorted) list of the variable indices of the DD variables given by the meta variable set.
+             *
+             * @param manager The manager responsible for the DD.
+             * @param metaVariable The set of meta variables for which to retrieve the index list.
+             * @return The sorted list of variable indices.
+             */
+            static std::vector<uint_fast64_t> getSortedVariableIndices(DdManager<DdType::CUDD> const& manager, std::set<storm::expressions::Variable> const& metaVariables);
             
             /*!
              * Retrieves the set of all meta variables contained in the DD.
diff --git a/src/storage/dd/CuddDdManager.cpp b/src/storage/dd/CuddDdManager.cpp
index 63a84a58b..d7adf1a6c 100644
--- a/src/storage/dd/CuddDdManager.cpp
+++ b/src/storage/dd/CuddDdManager.cpp
@@ -278,5 +278,9 @@ namespace storm {
         void DdManager<DdType::CUDD>::triggerReordering() {
             this->getCuddManager().ReduceHeap(this->reorderingTechnique, 0);
         }
+        
+        std::shared_ptr<DdManager<DdType::CUDD> const> DdManager<DdType::CUDD>::asSharedPointer() const {
+            return this->shared_from_this();
+        }
     }
 }
\ No newline at end of file
diff --git a/src/storage/dd/CuddDdManager.h b/src/storage/dd/CuddDdManager.h
index 07c11eb19..522be1fba 100644
--- a/src/storage/dd/CuddDdManager.h
+++ b/src/storage/dd/CuddDdManager.h
@@ -165,6 +165,13 @@ namespace storm {
              */
             DdMetaVariable<DdType::CUDD> const& getMetaVariable(storm::expressions::Variable const& variable) const;
             
+            /*!
+             * Retrieves the manager as a shared pointer.
+             *
+             * @return A shared pointer to the manager.
+             */
+            std::shared_ptr<DdManager<DdType::CUDD> const> asSharedPointer() const;
+            
         private:
             /*!
              * Retrieves a list of names of the DD variables in the order of their index.
diff --git a/src/storage/dd/CuddOdd.cpp b/src/storage/dd/CuddOdd.cpp
index 40805d9e9..31bc42dd6 100644
--- a/src/storage/dd/CuddOdd.cpp
+++ b/src/storage/dd/CuddOdd.cpp
@@ -91,6 +91,10 @@ namespace storm {
             }
         }
         
+        bool Odd<DdType::CUDD>::isTerminalNode() const {
+            return this->elseNode == nullptr && this->thenNode == nullptr;
+        }
+        
         std::shared_ptr<Odd<DdType::CUDD>> Odd<DdType::CUDD>::buildOddFromAddRec(DdNode* dd, Cudd const& manager, uint_fast64_t currentLevel, uint_fast64_t maxLevel, std::vector<uint_fast64_t> const& ddVariableIndices, std::vector<std::unordered_map<DdNode*, std::shared_ptr<Odd<DdType::CUDD>>>>& uniqueTableForLevels) {
             // Check whether the ODD for this node has already been computed (for this level) and if so, return this instead.
             auto const& iterator = uniqueTableForLevels[currentLevel].find(dd);
diff --git a/src/storage/dd/CuddOdd.h b/src/storage/dd/CuddOdd.h
index 2d64275a3..9269e4ab1 100644
--- a/src/storage/dd/CuddOdd.h
+++ b/src/storage/dd/CuddOdd.h
@@ -96,6 +96,13 @@ namespace storm {
              */
             uint_fast64_t getNodeCount() const;
             
+            /*!
+             * Checks whether the given ODD node is a terminal node, i.e. has no successors.
+             *
+             * @return True iff the node is terminal.
+             */
+            bool isTerminalNode() const;
+            
         private:
             // Declare a hash functor that is used for the unique tables in the construction process.
             class HashFunctor {
diff --git a/src/utility/cli.h b/src/utility/cli.h
index 202786ba2..780201cf3 100644
--- a/src/utility/cli.h
+++ b/src/utility/cli.h
@@ -501,7 +501,7 @@ namespace storm {
                     std::shared_ptr<storm::models::symbolic::Dtmc<DdType>> dtmc = model->template as<storm::models::symbolic::Dtmc<DdType>>();
                     storm::modelchecker::HybridDtmcPrctlModelChecker<DdType, double> modelchecker(*dtmc);
                     if (modelchecker.canHandle(*formula.get())) {
-                        modelchecker.check(*formula.get());
+                        result = modelchecker.check(*formula.get());
                     }
                 } else {
                     STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "This functionality is not yet implemented.");
diff --git a/src/utility/vector.h b/src/utility/vector.h
index bfc60b3a6..b5567ae5b 100644
--- a/src/utility/vector.h
+++ b/src/utility/vector.h
@@ -130,67 +130,78 @@ namespace storm {
             }
             
             /*!
-             * Applies the given operation pointwise on the two given vectors and writes the result into the first
-             * vector.
+             * Applies the given operation pointwise on the two given vectors and writes the result to the third vector.
+             * To botain an in-place operation, the third vector may be equal to any of the other two vectors.
              *
-             * @param target The first operand and target vector.
+             * @param firstOperand The first operand.
              * @param secondOperand The second operand.
+             * @param target The target vector.
              */
             template<class T>
-            void applyPointwiseInPlace(std::vector<T>& target, std::vector<T> const& secondOperand, std::function<T (T const&, T const&)> function) {
-#ifdef DEBUG
-                if (target.size() != summand.size()) {
-                    throw storm::exceptions::InvalidArgumentException() << "Invalid call to storm::utility::vector::applyPointwiseInPlace: operand lengths mismatch.";
-                }
-#endif
+            void applyPointwise(std::vector<T> const& firstOperand, std::vector<T> const& secondOperand, std::vector<T>& target, std::function<T (T const&, T const&)> function) {
 #ifdef STORM_HAVE_INTELTBB
                 tbb::parallel_for(tbb::blocked_range<uint_fast64_t>(0, target.size()),
                                   [&](tbb::blocked_range<uint_fast64_t> const& range) {
-                                      std::transform(target.begin() + range.begin(), target.begin() + range.end(), secondOperand.begin() + range.begin(), target.begin() + range.begin(), function);
+                                      std::transform(firstOperand.begin() + range.begin(), firstOperand.begin() + range.end(), secondOperand.begin() + range.begin(), target.begin() + range.begin(), function);
                                   });
 #else
-                std::transform(target.begin(), target.end(), secondOperand.begin(), target.begin(), function);
+                std::transform(firstOperand.begin(), firstOperand.end(), secondOperand.begin(), target.begin(), function);
 #endif
             }
             
             /*!
              * Applies the given function pointwise on the given vector.
              *
-             * @param target The vector to which to apply the function.
+             * @param operand The vector to which to apply the function.
+             * @param target The target vector.
              * @param function The function to apply.
              */
             template<class T>
-            void applyPointwiseInPlace(std::vector<T>& target, std::function<T (T const&)> const& function) {
+            void applyPointwise(std::vector<T> const& operand, std::vector<T>& target, std::function<T (T const&)> const& function) {
 #ifdef STORM_HAVE_INTELTBB
                 tbb::parallel_for(tbb::blocked_range<uint_fast64_t>(0, target.size()),
                                   [&](tbb::blocked_range<uint_fast64_t> const& range) {
-                                      std::transform(target.begin() + range.begin(), target.begin() + range.end(), target.begin() + range.begin(), function);
+                                      std::transform(operand.begin() + range.begin(), operand.begin() + range.end(), target.begin() + range.begin(), function);
                                   });
 #else
-                std::transform(target.begin(), target.end(), target.begin(), function);
+                std::transform(operand.begin(), operand.end(), target.begin(), function);
 #endif
             }
             
             /*!
-             * Adds the two given vectors and writes the result into the first operand.
+             * Adds the two given vectors and writes the result to the target vector.
              *
-             * @param target The first summand and target vector.
-             * @param summand The second summand.
+             * @param firstOperand The first operand.
+             * @param secondOperand The second operand
+             * @param target The target vector.
              */
             template<class T>
-            void addVectorsInPlace(std::vector<T>& target, std::vector<T> const& summand) {
-                applyPointwiseInPlace<T>(target, summand, std::plus<T>());
+            void addVectors(std::vector<T> const& firstOperand, std::vector<T> const& secondOperand, std::vector<T>& target) {
+                applyPointwise<T>(firstOperand, secondOperand, target, std::plus<T>());
             }
             
             /*!
-             * Subtracts the two given vectors and writes the result into the first operand.
+             * Subtracts the two given vectors and writes the result to the target vector.
              *
-             * @param target The first summand and target vector.
-             * @param summand The second summand.
+             * @param firstOperand The first operand.
+             * @param secondOperand The second operand
+             * @param target The target vector.
+             */
+            template<class T>
+            void subtractVectors(std::vector<T> const& firstOperand, std::vector<T> const& secondOperand, std::vector<T>& target) {
+                applyPointwise<T>(firstOperand, secondOperand, target, std::minus<T>());
+            }
+            
+            /*!
+             * Multiplies the two given vectors (pointwise) and writes the result to the target vector.
+             *
+             * @param firstOperand The first operand.
+             * @param secondOperand The second operand
+             * @param target The target vector.
              */
             template<class T>
-            void subtractVectorsInPlace(std::vector<T>& target, std::vector<T> const& summand) {
-                applyPointwiseInPlace<T>(target, summand, std::minus<T>());
+            void multiplyVectorsPointwise(std::vector<T> const& firstOperand, std::vector<T> const& secondOperand, std::vector<T>& target) {
+                applyPointwise<T>(firstOperand, secondOperand, target, std::multiplies<T>());
             }
             
             /*!
@@ -201,7 +212,7 @@ namespace storm {
              */
             template<class T>
             void scaleVectorInPlace(std::vector<T>& target, T const& factor) {
-                applyPointwiseInPlace<T>(target, [&] (T const& argument) { return argument * factor; });
+                applyPointwise<T>(target, target, [&] (T const& argument) { return argument * factor; });
             }
             
             /*!
diff --git a/test/functional/storage/CuddDdTest.cpp b/test/functional/storage/CuddDdTest.cpp
index 71efcadaf..6255e912c 100644
--- a/test/functional/storage/CuddDdTest.cpp
+++ b/test/functional/storage/CuddDdTest.cpp
@@ -376,6 +376,9 @@ TEST(CuddDd, BddOddTest) {
         EXPECT_TRUE(i+1 == ddAsVector[i]);
     }
     
+    storm::dd::Add<storm::dd::DdType::CUDD> vectorAdd(manager, ddAsVector, odd, {x.first});
+    EXPECT_EQ(dd, vectorAdd);
+    
     // Create a non-trivial matrix.
     dd = manager->getIdentity(x.first).equals(manager->getIdentity(x.second)) * manager->getRange(x.first).toAdd();
     dd += manager->getEncoding(x.first, 1).toAdd() * manager->getRange(x.second).toAdd() + manager->getEncoding(x.second, 1).toAdd() * manager->getRange(x.first).toAdd();
diff --git a/test/functional/storage/SparseMatrixTest.cpp b/test/functional/storage/SparseMatrixTest.cpp
index b677f6f3c..20716916a 100644
--- a/test/functional/storage/SparseMatrixTest.cpp
+++ b/test/functional/storage/SparseMatrixTest.cpp
@@ -447,7 +447,7 @@ TEST(SparseMatrix, JacobiDecomposition) {
     ASSERT_NO_THROW(matrix = matrixBuilder.build());
     
     ASSERT_NO_THROW(matrix.getJacobiDecomposition());
-    std::pair<storm::storage::SparseMatrix<double>, storm::storage::SparseMatrix<double>> jacobiDecomposition = matrix.getJacobiDecomposition();
+    std::pair<storm::storage::SparseMatrix<double>, std::vector<double>> jacobiDecomposition = matrix.getJacobiDecomposition();
     
     storm::storage::SparseMatrixBuilder<double> luBuilder(4, 4, 3);
     ASSERT_NO_THROW(luBuilder.addNextValue(0, 1, 1.2));
@@ -456,14 +456,7 @@ TEST(SparseMatrix, JacobiDecomposition) {
     storm::storage::SparseMatrix<double> lu;
     ASSERT_NO_THROW(lu = luBuilder.build());
     
-    storm::storage::SparseMatrixBuilder<double> dinvBuilder(4, 4, 4);
-    ASSERT_NO_THROW(dinvBuilder.addNextValue(0, 0, 1 / 1.1));
-    ASSERT_NO_THROW(dinvBuilder.addNextValue(1, 1, 1 / 0.5));
-    ASSERT_NO_THROW(dinvBuilder.addNextValue(2, 2, 1 / 0.99));
-    ASSERT_NO_THROW(dinvBuilder.addNextValue(3, 3, 1 / 0.11));
-    storm::storage::SparseMatrix<double> dinv;
-    ASSERT_NO_THROW(dinv = dinvBuilder.build());
-    
+    std::vector<double> dinv = {1/1.1, 1/0.5, 1/0.99, 1/0.11};    
     ASSERT_TRUE(lu == jacobiDecomposition.first);
     ASSERT_TRUE(dinv == jacobiDecomposition.second);
 }