diff --git a/src/adapters/GmmxxAdapter.h b/src/adapters/GmmxxAdapter.h
index 43a3f8302..3fd990a46 100644
--- a/src/adapters/GmmxxAdapter.h
+++ b/src/adapters/GmmxxAdapter.h
@@ -5,10 +5,10 @@
  *      Author: Christian Dehnert
  */
 
-#ifndef GMMXXADAPTER_H_
-#define GMMXXADAPTER_H_
+#ifndef STORM_ADAPTERS_GMMXXADAPTER_H_
+#define STORM_ADAPTERS_GMMXXADAPTER_H_
 
-#include "src/storage/SquareSparseMatrix.h"
+#include "src/storage/SparseMatrix.h"
 
 #include "log4cplus/logger.h"
 #include "log4cplus/loggingmacros.h"
@@ -26,74 +26,20 @@ public:
 	 * @return A pointer to a column-major sparse matrix in gmm++ format.
 	 */
 	template<class T>
-	static gmm::csr_matrix<T>* toGmmxxSparseMatrix(storm::storage::SquareSparseMatrix<T> const& matrix) {
-		uint_fast64_t realNonZeros = matrix.getNonZeroEntryCount() + matrix.getDiagonalNonZeroEntryCount();
+	static gmm::csr_matrix<T>* toGmmxxSparseMatrix(storm::storage::SparseMatrix<T> const& matrix) {
+		uint_fast64_t realNonZeros = matrix.getNonZeroEntryCount();
 		LOG4CPLUS_DEBUG(logger, "Converting matrix with " << realNonZeros << " non-zeros to gmm++ format.");
 
 		// Prepare the resulting matrix.
 		gmm::csr_matrix<T>* result = new gmm::csr_matrix<T>(matrix.rowCount, matrix.rowCount);
 
-		// Reserve enough elements for the row indications.
-		result->jc.reserve(matrix.rowCount + 1);
-
-		// For the column indications and the actual values, we have to gather
-		// the values in a temporary array first, as we have to integrate
-		// the values from the diagonal. For the row indications, we can just count the number of
-		// inserted diagonal elements and add it to the previous value.
-		uint_fast64_t* tmpColumnIndicationsArray = new uint_fast64_t[realNonZeros];
-		T* tmpValueArray = new T[realNonZeros];
-		T zero(0);
-		uint_fast64_t currentPosition = 0;
-		uint_fast64_t insertedDiagonalElements = 0;
-		for (uint_fast64_t i = 0; i < matrix.rowCount; ++i) {
-			// Compute correct start index of row.
-			result->jc[i] = matrix.rowIndications[i] + insertedDiagonalElements;
-
-			// If the current row has no non-zero which is not on the diagonal, we have to check the
-			// diagonal element explicitly.
-			if (matrix.rowIndications[i + 1] - matrix.rowIndications[i] == 0) {
-				if (matrix.diagonalStorage[i] != zero) {
-					tmpColumnIndicationsArray[currentPosition] = i;
-					tmpValueArray[currentPosition] = matrix.diagonalStorage[i];
-					++currentPosition; ++insertedDiagonalElements;
-				}
-			} else {
-				// Otherwise, we can just enumerate the non-zeros which are not on the diagonal
-				// and fit in the diagonal element where appropriate.
-				bool includedDiagonal = false;
-				for (uint_fast64_t j = matrix.rowIndications[i]; j < matrix.rowIndications[i + 1]; ++j) {
-					if (matrix.diagonalStorage[i] != zero && !includedDiagonal && matrix.columnIndications[j] > i) {
-						includedDiagonal = true;
-						tmpColumnIndicationsArray[currentPosition] = i;
-						tmpValueArray[currentPosition] = matrix.diagonalStorage[i];
-						++currentPosition; ++insertedDiagonalElements;
-					}
-					tmpColumnIndicationsArray[currentPosition] = matrix.columnIndications[j];
-					tmpValueArray[currentPosition] = matrix.valueStorage[j];
-					++currentPosition;
-				}
-
-				// If the diagonal element is non-zero and was not inserted until now (i.e. all
-				// off-diagonal elements in the row are before the diagonal element.
-				if (!includedDiagonal && matrix.diagonalStorage[i] != zero) {
-					tmpColumnIndicationsArray[currentPosition] = i;
-					tmpValueArray[currentPosition] = matrix.diagonalStorage[i];
-					++currentPosition; ++insertedDiagonalElements;
-				}
-			}
-		}
-		// Fill in sentinel element at the end.
-		result->jc[matrix.rowCount] = static_cast<unsigned int>(realNonZeros);
-
-		// Now, we can copy the temporary array to the GMMXX format.
+		// Copy Row Indications
+		std::copy(matrix.rowIndications.begin(), matrix.rowIndications.end(), std::back_inserter(result->jc));
+		// Copy Columns Indications
 		result->ir.resize(realNonZeros);
-		std::copy(tmpColumnIndicationsArray, tmpColumnIndicationsArray + realNonZeros, result->ir.begin());
-		delete[] tmpColumnIndicationsArray;
-
+		std::copy(matrix.columnIndications.begin(), matrix.columnIndications.end(), std::back_inserter(result->ir));
 		// And do the same thing with the actual values.
-		result->pr.resize(realNonZeros);
-		std::copy(tmpValueArray, tmpValueArray + realNonZeros, result->pr.begin());
-		delete[] tmpValueArray;
+		std::copy(matrix.valueStorage.begin(), matrix.valueStorage.end(), std::back_inserter(result->pr));
 
 		LOG4CPLUS_DEBUG(logger, "Done converting matrix to gmm++ format.");
 
@@ -105,4 +51,4 @@ public:
 
 } //namespace storm
 
-#endif /* GMMXXADAPTER_H_ */
+#endif /* STORM_ADAPTERS_GMMXXADAPTER_H_ */
diff --git a/src/formula/BoundedEventually.h b/src/formula/BoundedEventually.h
index 58973456d..0a66f6b83 100644
--- a/src/formula/BoundedEventually.h
+++ b/src/formula/BoundedEventually.h
@@ -120,7 +120,7 @@ public:
 		BoundedEventually<T>* result = new BoundedEventually<T>();
 		result->setBound(bound);
 		if (child != nullptr) {
-			result->setRight(child->clone());
+			result->setChild(child->clone());
 		}
 		return result;
 	}
diff --git a/src/modelChecker/EigenDtmcPrctlModelChecker.h b/src/modelChecker/EigenDtmcPrctlModelChecker.h
index 3c387bec5..e6ea58725 100644
--- a/src/modelChecker/EigenDtmcPrctlModelChecker.h
+++ b/src/modelChecker/EigenDtmcPrctlModelChecker.h
@@ -48,7 +48,7 @@ public:
 		storm::storage::BitVector* rightStates = this->checkStateFormula(formula.getRight());
 
 		// Copy the matrix before we make any changes.
-		storm::storage::SquareSparseMatrix<Type> tmpMatrix(*this->getModel().getTransitionProbabilityMatrix());
+		storm::storage::SparseMatrix<Type> tmpMatrix(*this->getModel().getTransitionProbabilityMatrix());
 
 		// Make all rows absorbing that violate both sub-formulas or satisfy the second sub-formula.
 		tmpMatrix.makeRowsAbsorbing((~*leftStates | *rightStates) | *rightStates);
@@ -148,7 +148,7 @@ public:
 			typedef Eigen::Map<VectorType> MapType;
 
 			// Now we can eliminate the rows and columns from the original transition probability matrix.
-			storm::storage::SquareSparseMatrix<double>* submatrix = this->getModel().getTransitionProbabilityMatrix()->getSubmatrix(maybeStates);
+			storm::storage::SparseMatrix<double>* submatrix = this->getModel().getTransitionProbabilityMatrix()->getSubmatrix(maybeStates);
 			// Converting the matrix to the form needed for the equation system. That is, we go from
 			// x = A*x + b to (I-A)x = b.
 			submatrix->convertToEquationSystem();
diff --git a/src/modelChecker/GmmxxDtmcPrctlModelChecker.h b/src/modelChecker/GmmxxDtmcPrctlModelChecker.h
index 2d4d65d49..c90efc3b0 100644
--- a/src/modelChecker/GmmxxDtmcPrctlModelChecker.h
+++ b/src/modelChecker/GmmxxDtmcPrctlModelChecker.h
@@ -48,7 +48,7 @@ public:
 		storm::storage::BitVector* rightStates = this->checkStateFormula(formula.getRight());
 
 		// Copy the matrix before we make any changes.
-		storm::storage::SquareSparseMatrix<Type> tmpMatrix(*this->getModel().getTransitionProbabilityMatrix());
+		storm::storage::SparseMatrix<Type> tmpMatrix(*this->getModel().getTransitionProbabilityMatrix());
 
 		// Make all rows absorbing that violate both sub-formulas or satisfy the second sub-formula.
 		tmpMatrix.makeRowsAbsorbing(~(*leftStates | *rightStates) | *rightStates);
@@ -130,7 +130,7 @@ public:
 		// Only try to solve system if there are states for which the probability is unknown.
 		if (maybeStates.getNumberOfSetBits() > 0) {
 			// Now we can eliminate the rows and columns from the original transition probability matrix.
-			storm::storage::SquareSparseMatrix<Type>* submatrix = this->getModel().getTransitionProbabilityMatrix()->getSubmatrix(maybeStates);
+			storm::storage::SparseMatrix<Type>* submatrix = this->getModel().getTransitionProbabilityMatrix()->getSubmatrix(maybeStates);
 			// Converting the matrix from the fixpoint notation to the form needed for the equation
 			// system. That is, we go from x = A*x + b to (I-A)x = b.
 			submatrix->convertToEquationSystem();
@@ -261,7 +261,7 @@ public:
 		storm::storage::BitVector maybeStates = ~(*targetStates) & ~infinityStates;
 		if (maybeStates.getNumberOfSetBits() > 0) {
 			// Now we can eliminate the rows and columns from the original transition probability matrix.
-			storm::storage::SquareSparseMatrix<Type>* submatrix = this->getModel().getTransitionProbabilityMatrix()->getSubmatrix(maybeStates);
+			storm::storage::SparseMatrix<Type>* submatrix = this->getModel().getTransitionProbabilityMatrix()->getSubmatrix(maybeStates);
 			// Converting the matrix from the fixpoint notation to the form needed for the equation
 			// system. That is, we go from x = A*x + b to (I-A)x = b.
 			submatrix->convertToEquationSystem();
diff --git a/src/models/Ctmc.h b/src/models/Ctmc.h
index 7f21528f0..413350124 100644
--- a/src/models/Ctmc.h
+++ b/src/models/Ctmc.h
@@ -15,7 +15,7 @@
 
 #include "AtomicPropositionsLabeling.h"
 #include "GraphTransitions.h"
-#include "src/storage/SquareSparseMatrix.h"
+#include "src/storage/SparseMatrix.h"
 #include "src/exceptions/InvalidArgumentException.h"
 
 namespace storm {
@@ -39,10 +39,10 @@ public:
 	 * @param stateLabeling The labeling that assigns a set of atomic
 	 * propositions to each state.
 	 */
-	Ctmc(std::shared_ptr<storm::storage::SquareSparseMatrix<T>> rateMatrix,
+	Ctmc(std::shared_ptr<storm::storage::SparseMatrix<T>> rateMatrix,
 			std::shared_ptr<storm::models::AtomicPropositionsLabeling> stateLabeling,
 			std::shared_ptr<std::vector<T>> stateRewards = nullptr,
-			std::shared_ptr<storm::storage::SquareSparseMatrix<T>> transitionRewardMatrix = nullptr)
+			std::shared_ptr<storm::storage::SparseMatrix<T>> transitionRewardMatrix = nullptr)
 			: rateMatrix(rateMatrix), stateLabeling(stateLabeling),
 			  stateRewards(stateRewards), transitionRewardMatrix(transitionRewardMatrix),
 			  backwardTransitions(nullptr) {
@@ -104,7 +104,7 @@ public:
 	 * @return A pointer to the matrix representing the transition probability
 	 * function.
 	 */
-	std::shared_ptr<storm::storage::SquareSparseMatrix<T>> getTransitionRateMatrix() const {
+	std::shared_ptr<storm::storage::SparseMatrix<T>> getTransitionRateMatrix() const {
 		return this->rateMatrix;
 	}
 
@@ -112,7 +112,7 @@ public:
 	 * Returns a pointer to the matrix representing the transition rewards.
 	 * @return A pointer to the matrix representing the transition rewards.
 	 */
-	std::shared_ptr<storm::storage::SquareSparseMatrix<T>> getTransitionRewardMatrix() const {
+	std::shared_ptr<storm::storage::SparseMatrix<T>> getTransitionRewardMatrix() const {
 		return this->transitionRewardMatrix;
 	}
 
@@ -164,7 +164,7 @@ public:
 private:
 
 	/*! A matrix representing the transition rate function of the CTMC. */
-	std::shared_ptr<storm::storage::SquareSparseMatrix<T>> rateMatrix;
+	std::shared_ptr<storm::storage::SparseMatrix<T>> rateMatrix;
 
 	/*! The labeling of the states of the CTMC. */
 	std::shared_ptr<storm::models::AtomicPropositionsLabeling> stateLabeling;
@@ -173,7 +173,7 @@ private:
 	std::shared_ptr<std::vector<T>> stateRewards;
 
 	/*! The transition-based rewards of the CTMC. */
-	std::shared_ptr<storm::storage::SquareSparseMatrix<T>> transitionRewardMatrix;
+	std::shared_ptr<storm::storage::SparseMatrix<T>> transitionRewardMatrix;
 
 	/*!
 	 * A data structure that stores the predecessors for all states. This is
diff --git a/src/models/Dtmc.h b/src/models/Dtmc.h
index cb75afc28..55abacffe 100644
--- a/src/models/Dtmc.h
+++ b/src/models/Dtmc.h
@@ -15,7 +15,7 @@
 
 #include "AtomicPropositionsLabeling.h"
 #include "GraphTransitions.h"
-#include "src/storage/SquareSparseMatrix.h"
+#include "src/storage/SparseMatrix.h"
 #include "src/exceptions/InvalidArgumentException.h"
 #include "src/utility/CommandLine.h"
 
@@ -40,10 +40,10 @@ public:
 	 * @param stateLabeling The labeling that assigns a set of atomic
 	 * propositions to each state.
 	 */
-	Dtmc(std::shared_ptr<storm::storage::SquareSparseMatrix<T>> probabilityMatrix,
+	Dtmc(std::shared_ptr<storm::storage::SparseMatrix<T>> probabilityMatrix,
 			std::shared_ptr<storm::models::AtomicPropositionsLabeling> stateLabeling,
 			std::shared_ptr<std::vector<T>> stateRewards = nullptr,
-			std::shared_ptr<storm::storage::SquareSparseMatrix<T>> transitionRewardMatrix = nullptr)
+			std::shared_ptr<storm::storage::SparseMatrix<T>> transitionRewardMatrix = nullptr)
 			: probabilityMatrix(probabilityMatrix), stateLabeling(stateLabeling),
 			  stateRewards(stateRewards), transitionRewardMatrix(transitionRewardMatrix),
 			  backwardTransitions(nullptr) {
@@ -111,7 +111,7 @@ public:
 	 * @return A pointer to the matrix representing the transition probability
 	 * function.
 	 */
-	std::shared_ptr<storm::storage::SquareSparseMatrix<T>> getTransitionProbabilityMatrix() const {
+	std::shared_ptr<storm::storage::SparseMatrix<T>> getTransitionProbabilityMatrix() const {
 		return this->probabilityMatrix;
 	}
 
@@ -119,7 +119,7 @@ public:
 	 * Returns a pointer to the matrix representing the transition rewards.
 	 * @return A pointer to the matrix representing the transition rewards.
 	 */
-	std::shared_ptr<storm::storage::SquareSparseMatrix<T>> getTransitionRewardMatrix() const {
+	std::shared_ptr<storm::storage::SparseMatrix<T>> getTransitionRewardMatrix() const {
 		return this->transitionRewardMatrix;
 	}
 
@@ -210,7 +210,7 @@ private:
 	}
 
 	/*! A matrix representing the transition probability function of the DTMC. */
-	std::shared_ptr<storm::storage::SquareSparseMatrix<T>> probabilityMatrix;
+	std::shared_ptr<storm::storage::SparseMatrix<T>> probabilityMatrix;
 
 	/*! The labeling of the states of the DTMC. */
 	std::shared_ptr<storm::models::AtomicPropositionsLabeling> stateLabeling;
@@ -219,7 +219,7 @@ private:
 	std::shared_ptr<std::vector<T>> stateRewards;
 
 	/*! The transition-based rewards of the DTMC. */
-	std::shared_ptr<storm::storage::SquareSparseMatrix<T>> transitionRewardMatrix;
+	std::shared_ptr<storm::storage::SparseMatrix<T>> transitionRewardMatrix;
 
 	/*!
 	 * A data structure that stores the predecessors for all states. This is
diff --git a/src/models/GraphTransitions.h b/src/models/GraphTransitions.h
index 3ccc6bed5..576d99f8c 100644
--- a/src/models/GraphTransitions.h
+++ b/src/models/GraphTransitions.h
@@ -8,7 +8,7 @@
 #ifndef STORM_MODELS_GRAPHTRANSITIONS_H_
 #define STORM_MODELS_GRAPHTRANSITIONS_H_
 
-#include "src/storage/SquareSparseMatrix.h"
+#include "src/storage/SparseMatrix.h"
 
 #include <algorithm>
 #include <memory>
@@ -39,7 +39,7 @@ public:
 	 * @param forward If set to true, this objects will store the graph structure
 	 * of the backwards transition relation.
 	 */
-	GraphTransitions(std::shared_ptr<storm::storage::SquareSparseMatrix<T>> transitionMatrix, bool forward)
+	GraphTransitions(std::shared_ptr<storm::storage::SparseMatrix<T>> transitionMatrix, bool forward)
 			: successorList(nullptr), stateIndications(nullptr), numberOfStates(transitionMatrix->getRowCount()), numberOfNonZeroTransitions(transitionMatrix->getNonZeroEntryCount()) {
 		if (forward) {
 			this->initializeForward(transitionMatrix);
@@ -87,18 +87,18 @@ private:
 	 * Initializes this graph transitions object using the forward transition
 	 * relation given by means of a sparse matrix.
 	 */
-	void initializeForward(std::shared_ptr<storm::storage::SquareSparseMatrix<T>> transitionMatrix) {
+	void initializeForward(std::shared_ptr<storm::storage::SparseMatrix<T>> transitionMatrix) {
 		this->successorList = new uint_fast64_t[numberOfNonZeroTransitions];
 		this->stateIndications = new uint_fast64_t[numberOfStates + 1];
 
 		// First, we copy the index list from the sparse matrix as this will
 		// stay the same.
-		std::copy(transitionMatrix->getRowIndicationsPointer(), transitionMatrix->getRowIndicationsPointer() + numberOfStates + 1, this->stateIndications);
+		std::copy(transitionMatrix->getRowIndicationsPointer().begin(), transitionMatrix->getRowIndicationsPointer().end(), this->stateIndications);
 
 		// Now we can iterate over all rows of the transition matrix and record
 		// the target state.
 		for (uint_fast64_t i = 0, currentNonZeroElement = 0; i < numberOfStates; i++) {
-			for (auto rowIt = transitionMatrix->beginConstColumnNoDiagIterator(i); rowIt != transitionMatrix->endConstColumnNoDiagIterator(i); ++rowIt) {
+			for (auto rowIt = transitionMatrix->beginConstColumnIterator(i); rowIt != transitionMatrix->endConstColumnIterator(i); ++rowIt) {
 				this->stateIndications[currentNonZeroElement++] = *rowIt;
 			}
 		}
@@ -109,7 +109,7 @@ private:
 	 * relation, whose forward transition relation is given by means of a sparse
 	 * matrix.
 	 */
-	void initializeBackward(std::shared_ptr<storm::storage::SquareSparseMatrix<T>> transitionMatrix) {
+	void initializeBackward(std::shared_ptr<storm::storage::SparseMatrix<T>> transitionMatrix) {
 		this->successorList = new uint_fast64_t[numberOfNonZeroTransitions]();
 		this->stateIndications = new uint_fast64_t[numberOfStates + 1]();
 
@@ -117,7 +117,7 @@ private:
 		// NOTE: We disregard the diagonal here, as we only consider "true"
 		// predecessors.
 		for (uint_fast64_t i = 0; i < numberOfStates; i++) {
-			for (auto rowIt = transitionMatrix->beginConstColumnNoDiagIterator(i); rowIt != transitionMatrix->endConstColumnNoDiagIterator(i); ++rowIt) {
+			for (auto rowIt = transitionMatrix->beginConstColumnIterator(i); rowIt != transitionMatrix->endConstColumnIterator(i); ++rowIt) {
 				this->stateIndications[*rowIt + 1]++;
 			}
 		}
@@ -140,7 +140,7 @@ private:
 		// Now we are ready to actually fill in the list of predecessors for
 		// every state. Again, we start by considering all but the last row.
 		for (uint_fast64_t i = 0; i < numberOfStates; i++) {
-			for (auto rowIt = transitionMatrix->beginConstColumnNoDiagIterator(i); rowIt != transitionMatrix->endConstColumnNoDiagIterator(i); ++rowIt) {
+			for (auto rowIt = transitionMatrix->beginConstColumnIterator(i); rowIt != transitionMatrix->endConstColumnIterator(i); ++rowIt) {
 				this->successorList[nextIndicesList[*rowIt]++] = i;
 			}
 		}
diff --git a/src/parser/DeterministicSparseTransitionParser.cpp b/src/parser/DeterministicSparseTransitionParser.cpp
index d82e4e3d9..662adc1b7 100644
--- a/src/parser/DeterministicSparseTransitionParser.cpp
+++ b/src/parser/DeterministicSparseTransitionParser.cpp
@@ -32,13 +32,10 @@ namespace parser{
  *	non-zero cells and maximum node id.
  *
  *	This method does the first pass through the .tra file and computes
- *	the number of non-zero elements that are not diagonal elements,
- *	which correspondents to the number of transitions that are not
- *	self-loops.
- *	(Diagonal elements are treated in a special way).
+ *	the number of non-zero elements.
  *	It also calculates the maximum node id and stores it in maxnode.
  *
- *	@return The number of non-zero elements that are not on the diagonal
+ *	@return The number of non-zero elements
  *	@param buf Data to scan. Is expected to be some char array.
  *	@param maxnode Is set to highest id of all nodes.
  */
@@ -88,7 +85,8 @@ uint_fast64_t DeterministicSparseTransitionParser::firstPass(char* buf, uint_fas
 			LOG4CPLUS_ERROR(logger, "Expected a positive probability but got \"" << val << "\".");
 			return 0;
 		}
-		if (row == col) non_zero--;
+		// not needed anymore
+		//if (row == col) non_zero--;
 		buf = trimWhitespaces(buf);
 	}
 
@@ -154,7 +152,7 @@ DeterministicSparseTransitionParser::DeterministicSparseTransitionParser(std::st
 	 *	non-zero elements has to be specified (which is non_zero, computed by make_first_pass)
 	 */
 	LOG4CPLUS_INFO(logger, "Attempting to create matrix of size " << (maxnode+1) << " x " << (maxnode+1) << ".");
-	this->matrix = std::shared_ptr<storm::storage::SquareSparseMatrix<double>>(new storm::storage::SquareSparseMatrix<double>(maxnode + 1));
+	this->matrix = std::shared_ptr<storm::storage::SparseMatrix<double>>(new storm::storage::SparseMatrix<double>(maxnode + 1));
 	if (this->matrix == NULL)
 	{
 		LOG4CPLUS_ERROR(logger, "Could not create matrix of size " << (maxnode+1) << " x " << (maxnode+1) << ".");
diff --git a/src/parser/DeterministicSparseTransitionParser.h b/src/parser/DeterministicSparseTransitionParser.h
index a5b8560de..1d699d0c8 100644
--- a/src/parser/DeterministicSparseTransitionParser.h
+++ b/src/parser/DeterministicSparseTransitionParser.h
@@ -1,7 +1,7 @@
 #ifndef STORM_PARSER_TRAPARSER_H_
 #define STORM_PARSER_TRAPARSER_H_
 
-#include "src/storage/SquareSparseMatrix.h"
+#include "src/storage/SparseMatrix.h"
 
 #include "src/parser/Parser.h"
 #include "src/utility/OsDetection.h"
@@ -19,12 +19,12 @@ class DeterministicSparseTransitionParser : public Parser {
 	public:
 		DeterministicSparseTransitionParser(std::string const &filename);
 		
-		std::shared_ptr<storm::storage::SquareSparseMatrix<double>> getMatrix() {
+		std::shared_ptr<storm::storage::SparseMatrix<double>> getMatrix() {
 			return this->matrix;
 		}
 	
 	private:
-		std::shared_ptr<storm::storage::SquareSparseMatrix<double>> matrix;
+		std::shared_ptr<storm::storage::SparseMatrix<double>> matrix;
 		
 		uint_fast64_t firstPass(char* buf, uint_fast64_t &maxnode);
 	
diff --git a/src/parser/DtmcParser.cpp b/src/parser/DtmcParser.cpp
index 4f4ed2645..1002c021f 100644
--- a/src/parser/DtmcParser.cpp
+++ b/src/parser/DtmcParser.cpp
@@ -29,7 +29,7 @@ DtmcParser::DtmcParser(std::string const & transitionSystemFile, std::string con
 	uint_fast64_t stateCount = tp.getMatrix()->getRowCount();
 
 	std::shared_ptr<std::vector<double>> stateRewards = nullptr;
-	std::shared_ptr<storm::storage::SquareSparseMatrix<double>> transitionRewards = nullptr;
+	std::shared_ptr<storm::storage::SparseMatrix<double>> transitionRewards = nullptr;
 
 	storm::parser::AtomicPropositionLabelingParser lp(stateCount, labelingFile);
 	if (stateRewardFile != "") {
diff --git a/src/parser/NonDeterministicSparseTransitionParser.cpp b/src/parser/NonDeterministicSparseTransitionParser.cpp
index eaf72bc7d..ad2654f98 100644
--- a/src/parser/NonDeterministicSparseTransitionParser.cpp
+++ b/src/parser/NonDeterministicSparseTransitionParser.cpp
@@ -167,7 +167,7 @@ NonDeterministicSparseTransitionParser::NonDeterministicSparseTransitionParser(s
 	 *	non-zero elements has to be specified (which is non_zero, computed by make_first_pass)
 	 */
 	LOG4CPLUS_INFO(logger, "Attempting to create matrix of size " << (maxnode+1) << " x " << (maxnode+1) << ".");
-	this->matrix = std::shared_ptr<storm::storage::SquareSparseMatrix<double>>(new storm::storage::SquareSparseMatrix<double>(maxnode + 1));
+	this->matrix = std::shared_ptr<storm::storage::SparseMatrix<double>>(new storm::storage::SparseMatrix<double>(maxnode + 1));
 	if (this->matrix == NULL)
 	{
 		LOG4CPLUS_ERROR(logger, "Could not create matrix of size " << (maxnode+1) << " x " << (maxnode+1) << ".");
diff --git a/src/parser/NonDeterministicSparseTransitionParser.h b/src/parser/NonDeterministicSparseTransitionParser.h
index bc144b160..bd46b2593 100644
--- a/src/parser/NonDeterministicSparseTransitionParser.h
+++ b/src/parser/NonDeterministicSparseTransitionParser.h
@@ -1,7 +1,7 @@
 #ifndef STORM_PARSER_NONDETTRAPARSER_H_
 #define STORM_PARSER_NONDETTRAPARSER_H_
 
-#include "src/storage/SquareSparseMatrix.h"
+#include "src/storage/SparseMatrix.h"
 
 #include "src/parser/Parser.h"
 #include "src/utility/OsDetection.h"
@@ -20,12 +20,12 @@ class NonDeterministicSparseTransitionParser : public Parser {
 	public:
 		NonDeterministicSparseTransitionParser(std::string const &filename);
 		
-		std::shared_ptr<storm::storage::SquareSparseMatrix<double>> getMatrix() {
+		std::shared_ptr<storm::storage::SparseMatrix<double>> getMatrix() {
 			return this->matrix;
 		}
 	
 	private:
-		std::shared_ptr<storm::storage::SquareSparseMatrix<double>> matrix;
+		std::shared_ptr<storm::storage::SparseMatrix<double>> matrix;
 		
 		std::unique_ptr<std::vector<uint_fast64_t>> firstPass(char* buf, uint_fast64_t &maxnode, uint_fast64_t &maxchoice);
 	
diff --git a/src/storage/JacobiDecomposition.h b/src/storage/JacobiDecomposition.h
index d62f36712..f6de68e98 100644
--- a/src/storage/JacobiDecomposition.h
+++ b/src/storage/JacobiDecomposition.h
@@ -16,7 +16,7 @@ namespace storage {
  * Forward declaration against Cycle
  */
 template <class T>
-class SquareSparseMatrix;
+class SparseMatrix;
 
 
 /*!
@@ -26,7 +26,7 @@ template <class T>
 class JacobiDecomposition {
 
 public:
-	JacobiDecomposition(storm::storage::SquareSparseMatrix<T> * const jacobiLuMatrix, storm::storage::SquareSparseMatrix<T> * const jacobiDInvMatrix) : jacobiLuMatrix(jacobiLuMatrix), jacobiDInvMatrix(jacobiDInvMatrix) {
+	JacobiDecomposition(storm::storage::SparseMatrix<T> * const jacobiLuMatrix, storm::storage::SparseMatrix<T> * const jacobiDInvMatrix) : jacobiLuMatrix(jacobiLuMatrix), jacobiDInvMatrix(jacobiDInvMatrix) {
 	}
 
 	~JacobiDecomposition() {
@@ -39,7 +39,7 @@ public:
 	 * Ownership stays with this class.
 	 * @return A reference to the Jacobi LU Matrix
 	 */
-	storm::storage::SquareSparseMatrix<T>& getJacobiLUReference() {
+	storm::storage::SparseMatrix<T>& getJacobiLUReference() {
 		return *(this->jacobiLuMatrix);
 	}
 
@@ -48,7 +48,7 @@ public:
 	 * Ownership stays with this class.
 	 * @return A reference to the Jacobi D^{-1} Matrix
 	 */
-	storm::storage::SquareSparseMatrix<T>& getJacobiDInvReference() {
+	storm::storage::SparseMatrix<T>& getJacobiDInvReference() {
 		return *(this->jacobiDInvMatrix);
 	}
 
@@ -57,7 +57,7 @@ public:
 	 * Ownership stays with this class.
 	 * @return A pointer to the Jacobi LU Matrix
 	 */
-	storm::storage::SquareSparseMatrix<T>* getJacobiLU() {
+	storm::storage::SparseMatrix<T>* getJacobiLU() {
 		return this->jacobiLuMatrix;
 	}
 
@@ -66,7 +66,7 @@ public:
 	 * Ownership stays with this class.
 	 * @return A pointer to the Jacobi D^{-1} Matrix
 	 */
-	storm::storage::SquareSparseMatrix<T>* getJacobiDInv() {
+	storm::storage::SparseMatrix<T>* getJacobiDInv() {
 		return this->jacobiDInvMatrix;
 	}
 
@@ -81,12 +81,12 @@ private:
 	/*!
 	 * Pointer to the LU Matrix
 	 */
-	storm::storage::SquareSparseMatrix<T> *jacobiLuMatrix;
+	storm::storage::SparseMatrix<T> *jacobiLuMatrix;
 	
 	/*!
 	 * Pointer to the D^{-1} Matrix
 	 */
-	storm::storage::SquareSparseMatrix<T> *jacobiDInvMatrix;
+	storm::storage::SparseMatrix<T> *jacobiDInvMatrix;
 };
 
 } // namespace storage
diff --git a/src/storage/SquareSparseMatrix.h b/src/storage/SparseMatrix.h
similarity index 74%
rename from src/storage/SquareSparseMatrix.h
rename to src/storage/SparseMatrix.h
index ae48cf4ed..6e4ca5066 100644
--- a/src/storage/SquareSparseMatrix.h
+++ b/src/storage/SparseMatrix.h
@@ -1,10 +1,11 @@
-#ifndef STORM_STORAGE_SQUARESPARSEMATRIX_H_
-#define STORM_STORAGE_SQUARESPARSEMATRIX_H_
+#ifndef STORM_STORAGE_SPARSEMATRIX_H_
+#define STORM_STORAGE_SPARSEMATRIX_H_
 
 #include <exception>
 #include <new>
 #include <algorithm>
 #include <iostream>
+#include <iterator>
 #include "boost/integer/integer_mask.hpp"
 
 #include "src/exceptions/InvalidStateException.h"
@@ -31,13 +32,12 @@ namespace storm {
 namespace storage {
 
 /*!
- * A sparse matrix class with a constant number of non-zero entries on the non-diagonal fields
- * and a separate dense storage for the diagonal elements.
+ * A sparse matrix class with a constant number of non-zero entries.
  * NOTE: Addressing *is* zero-based, so the valid range for getValue and addNextValue is 0..(rows - 1)
  * where rows is the first argument to the constructor.
  */
 template<class T>
-class SquareSparseMatrix {
+class SparseMatrix {
 public:
 	/*!
 	 * Declare adapter classes as friends to use internal data.
@@ -73,9 +73,25 @@ public:
 	 * Constructs a sparse matrix object with the given number of rows.
 	 * @param rows The number of rows of the matrix
 	 */
-	SquareSparseMatrix(uint_fast64_t rows)
-			: rowCount(rows), nonZeroEntryCount(0), valueStorage(nullptr),
-			  diagonalStorage(nullptr),columnIndications(nullptr), rowIndications(nullptr),
+	SparseMatrix(uint_fast64_t rows, uint_fast64_t cols)
+			: rowCount(rows), colCount(cols), nonZeroEntryCount(0),
+			  internalStatus(MatrixStatus::UnInitialized), currentSize(0), lastRow(0) { }
+
+	/* Sadly, Delegate Constructors are not yet available with MSVC2012
+	//! Constructor
+	/*!
+	 * Constructs a square sparse matrix object with the given number rows
+	 * @param size The number of rows and cols in the matrix
+	 */ /*
+	SparseMatrix(uint_fast64_t size) : SparseMatrix(size, size) { }
+	*/
+
+	//! Constructor
+	/*!
+	 * Constructs a square sparse matrix object with the given number rows
+	 * @param size The number of rows and cols in the matrix
+	 */
+	SparseMatrix(uint_fast64_t size) : rowCount(size), colCount(size), nonZeroEntryCount(0),
 			  internalStatus(MatrixStatus::UnInitialized), currentSize(0), lastRow(0) { }
 
 	//! Copy Constructor
@@ -83,8 +99,8 @@ public:
 	 * Copy Constructor. Performs a deep copy of the given sparse matrix.
 	 * @param ssm A reference to the matrix to be copied.
 	 */
-	SquareSparseMatrix(const SquareSparseMatrix<T> &ssm)
-			: rowCount(ssm.rowCount), nonZeroEntryCount(ssm.nonZeroEntryCount),
+	SparseMatrix(const SparseMatrix<T> &ssm)
+			: rowCount(ssm.rowCount), colCount(ssm.colCount), nonZeroEntryCount(ssm.nonZeroEntryCount),
 			  internalStatus(ssm.internalStatus), currentSize(ssm.currentSize), lastRow(ssm.lastRow) {
 		LOG4CPLUS_WARN(logger, "Invoking copy constructor.");
 		// Check whether copying the matrix is safe.
@@ -98,24 +114,12 @@ public:
 				LOG4CPLUS_ERROR(logger, "Unable to allocate internal storage.");
 				throw std::bad_alloc();
 			} else {
-				// Now that all storages have been prepared, copy over all
-				// elements. Start by copying the elements of type value and
-				// copy them seperately in order to invoke copy their copy
-				// constructor. This may not be necessary, but it is safer to
-				// do so in any case.
-				for (uint_fast64_t i = 0; i < nonZeroEntryCount; ++i) {
-					// use T() to force use of the copy constructor for complex T types
-					valueStorage[i] = T(ssm.valueStorage[i]);
-				}
-				for (uint_fast64_t i = 0; i < rowCount; ++i) {
-					// use T() to force use of the copy constructor for complex T types
-					diagonalStorage[i] = T(ssm.diagonalStorage[i]);
-				}
+				std::copy(ssm.valueStorage.begin(), ssm.valueStorage.end(), std::back_inserter(valueStorage));
 
 				// The elements that are not of the value type but rather the
 				// index type may be copied directly.
-				std::copy(ssm.columnIndications, ssm.columnIndications + nonZeroEntryCount, columnIndications);
-				std::copy(ssm.rowIndications, ssm.rowIndications + rowCount + 1, rowIndications);
+				std::copy(ssm.columnIndications.begin(), ssm.columnIndications.end(), std::back_inserter(columnIndications));
+				std::copy(ssm.rowIndications.begin(), ssm.rowIndications.end(), std::back_inserter(rowIndications));
 			}
 		}
 	}
@@ -124,20 +128,11 @@ public:
 	/*!
 	 * Destructor. Performs deletion of the reserved storage arrays.
 	 */
-	~SquareSparseMatrix() {
+	~SparseMatrix() {
 		setState(MatrixStatus::UnInitialized);
-		if (valueStorage != nullptr) {
-			delete[] valueStorage;
-		}
-		if (columnIndications != nullptr) {
-			delete[] columnIndications;
-		}
-		if (rowIndications != nullptr) {
-			delete[] rowIndications;
-		}
-		if (diagonalStorage != nullptr) {
-			delete[] diagonalStorage;
-		}
+		valueStorage.resize(0);
+		columnIndications.resize(0);
+		rowIndications.resize(0);
 	}
 
 	/*!
@@ -155,11 +150,11 @@ public:
 			triggerErrorState();
 			LOG4CPLUS_ERROR(logger, "Trying to initialize matrix that is not uninitialized.");
 			throw storm::exceptions::InvalidStateException("Trying to initialize matrix that is not uninitialized.");
-		} else if (rowCount == 0) {
+		} else if ((rowCount == 0) || (colCount == 0)) {
 			triggerErrorState();
-			LOG4CPLUS_ERROR(logger, "Trying to create initialize a matrix with 0 rows.");
-			throw storm::exceptions::InvalidArgumentException("Trying to create initialize a matrix with 0 rows.");
-		} else if (((rowCount * rowCount) - rowCount) < nonZeroEntries) {
+			LOG4CPLUS_ERROR(logger, "Trying to create initialize a matrix with 0 rows or 0 columns.");
+			throw storm::exceptions::InvalidArgumentException("Trying to create initialize a matrix with 0 rows or 0 columns.");
+		} else if ((rowCount * colCount) < nonZeroEntries) {
 			triggerErrorState();
 			LOG4CPLUS_ERROR(logger, "Trying to initialize a matrix with more non-zero entries than there can be.");
 			throw storm::exceptions::InvalidArgumentException("Trying to initialize a matrix with more non-zero entries than there can be.");
@@ -196,37 +191,72 @@ public:
 			throw storm::exceptions::InvalidArgumentException("Trying to initialize from an Eigen matrix that is not in compressed form.");
 		}
 
-		// Compute the actual (i.e. non-diagonal) number of non-zero entries.
-		nonZeroEntryCount = getEigenSparseMatrixCorrectNonZeroEntryCount(eigenSparseMatrix);
+		if (eigenSparseMatrix.rows() > this->rowCount) {
+			triggerErrorState();
+			LOG4CPLUS_ERROR(logger, "Trying to initialize from an Eigen matrix that has more rows than the target matrix.");
+			throw storm::exceptions::InvalidArgumentException("Trying to initialize from an Eigen matrix that has more rows than the target matrix.");
+		}
+		if (eigenSparseMatrix.cols() > this->colCount) {
+			triggerErrorState();
+			LOG4CPLUS_ERROR(logger, "Trying to initialize from an Eigen matrix that has more columns than the target matrix.");
+			throw storm::exceptions::InvalidArgumentException("Trying to initialize from an Eigen matrix that has more columns than the target matrix.");
+		}
+
+		const _Index entryCount = eigenSparseMatrix.nonZeros();
+		nonZeroEntryCount = entryCount;
 		lastRow = 0;
 
 		// Try to prepare the internal storage and throw an error in case of
 		// failure.
-		if (!prepareInternalStorage()) {
-			triggerErrorState();
-			LOG4CPLUS_ERROR(logger, "Unable to allocate internal storage.");
-			throw std::bad_alloc();
-		} else {
-			// Get necessary pointers to the contents of the Eigen matrix.
-			const T* valuePtr = eigenSparseMatrix.valuePtr();
-			const _Index* indexPtr = eigenSparseMatrix.innerIndexPtr();
-			const _Index* outerPtr = eigenSparseMatrix.outerIndexPtr();
-
-			// If the given matrix is in RowMajor format, copying can simply
-			// be done by adding all values in order.
-			// Direct copying is, however, prevented because we have to
-			// separate the diagonal entries from others.
-			if (isEigenRowMajor(eigenSparseMatrix)) {
-				// Because of the RowMajor format outerSize evaluates to the
-				// number of rows.
-				const _Index rowCount = eigenSparseMatrix.outerSize();
-				for (_Index row = 0; row < rowCount; ++row) {
-					for (_Index col = outerPtr[row]; col < outerPtr[row + 1]; ++col) {
-						addNextValue(row, indexPtr[col], valuePtr[col]);
-					}
+		
+		// Get necessary pointers to the contents of the Eigen matrix.
+		const T* valuePtr = eigenSparseMatrix.valuePtr();
+		const _Index* indexPtr = eigenSparseMatrix.innerIndexPtr();
+		const _Index* outerPtr = eigenSparseMatrix.outerIndexPtr();
+
+		// If the given matrix is in RowMajor format, copying can simply
+		// be done by adding all values in order.
+		// Direct copying is, however, prevented because we have to
+		// separate the diagonal entries from others.
+		if (isEigenRowMajor(eigenSparseMatrix)) {
+			// Because of the RowMajor format outerSize evaluates to the
+			// number of rows.
+			if (!prepareInternalStorage(false)) {
+				triggerErrorState();
+				LOG4CPLUS_ERROR(logger, "Unable to allocate internal storage.");
+				throw std::bad_alloc();
+			} else {
+				if ((eigenSparseMatrix.innerSize() > nonZeroEntryCount) || (entryCount > nonZeroEntryCount)) {
+					triggerErrorState();
+					LOG4CPLUS_ERROR(logger, "Invalid internal composition of Eigen Sparse Matrix");
+					throw storm::exceptions::InvalidArgumentException("Invalid internal composition of Eigen Sparse Matrix");
+				}
+				std::vector<uint_fast64_t> eigenColumnTemp;
+				std::vector<uint_fast64_t> eigenRowTemp;
+				std::vector<T> eigenValueTemp;
+				uint_fast64_t outerSize = eigenSparseMatrix.outerSize() + 1;
+
+				for (uint_fast64_t i = 0; i < entryCount; ++i) {
+					eigenColumnTemp.push_back(indexPtr[i]);
+					eigenValueTemp.push_back(valuePtr[i]);
 				}
+				for (uint_fast64_t i = 0; i < outerSize; ++i) {
+					eigenRowTemp.push_back(outerPtr[i]);
+				}
+
+				std::copy(eigenRowTemp.begin(), eigenRowTemp.end(), std::back_inserter(this->rowIndications));
+				std::copy(eigenColumnTemp.begin(), eigenColumnTemp.end(), std::back_inserter(this->columnIndications));
+				std::copy(eigenValueTemp.begin(), eigenValueTemp.end(), std::back_inserter(this->valueStorage));
+
+				currentSize = entryCount;
+				lastRow = rowCount;
+			}
+		} else {
+			if (!prepareInternalStorage()) {
+				triggerErrorState();
+				LOG4CPLUS_ERROR(logger, "Unable to allocate internal storage.");
+				throw std::bad_alloc();
 			} else {
-				const _Index entryCount = eigenSparseMatrix.nonZeros();
 				// Because of the ColMajor format outerSize evaluates to the
 				// number of columns.
 				const _Index colCount = eigenSparseMatrix.outerSize();
@@ -250,8 +280,7 @@ public:
 					// add it in case it is also in the current row.
 					if ((positions[currentColumn] < outerPtr[currentColumn + 1])
 							&& (indexPtr[positions[currentColumn]] == currentRow)) {
-						addNextValue(currentRow, currentColumn,
-								valuePtr[positions[currentColumn]]);
+						addNextValue(currentRow, currentColumn,	valuePtr[positions[currentColumn]]);
 						// Remember that we found one more non-zero element.
 						++i;
 						// Mark this position as "used".
@@ -268,8 +297,8 @@ public:
 				}
 				delete[] positions;
 			}
-			setState(MatrixStatus::Initialized);
 		}
+		setState(MatrixStatus::Initialized);
 	}
 
 	/*!
@@ -283,30 +312,27 @@ public:
 	void addNextValue(const uint_fast64_t row, const uint_fast64_t col,	const T& value) {
 		// Check whether the given row and column positions are valid and throw
 		// error otherwise.
-		if ((row > rowCount) || (col > rowCount)) {
+		if ((row > rowCount) || (col > colCount)) {
 			triggerErrorState();
 			LOG4CPLUS_ERROR(logger, "Trying to add a value at illegal position (" << row << ", " << col << ").");
 			throw storm::exceptions::OutOfRangeException("Trying to add a value at illegal position.");
 		}
 
-		if (row == col) { // Set a diagonal element.
-			diagonalStorage[row] = value;
-		} else { // Set a non-diagonal element.
-			// If we switched to another row, we have to adjust the missing
-			// entries in the row_indications array.
-			if (row != lastRow) {
-				for (uint_fast64_t i = lastRow + 1; i <= row; ++i) {
-					rowIndications[i] = currentSize;
-				}
-				lastRow = row;
+		
+		// If we switched to another row, we have to adjust the missing
+		// entries in the row_indications array.
+		if (row != lastRow) {
+			for (uint_fast64_t i = lastRow + 1; i <= row; ++i) {
+				rowIndications[i] = currentSize;
 			}
+			lastRow = row;
+		}
 
-			// Finally, set the element and increase the current size.
-			valueStorage[currentSize] = value;
-			columnIndications[currentSize] = col;
+		// Finally, set the element and increase the current size.
+		valueStorage[currentSize] = value;
+		columnIndications[currentSize] = col;
 
-			++currentSize;
-		}
+		++currentSize;
 	}
 
 	/*
@@ -355,18 +381,12 @@ public:
 	 */
 	inline bool getValue(uint_fast64_t row, uint_fast64_t col, T* const target) const {
 		// Check for illegal access indices.
-		if ((row > rowCount) || (col > rowCount)) {
+		if ((row > rowCount) || (col > colCount)) {
 			LOG4CPLUS_ERROR(logger, "Trying to read a value from illegal position (" << row << ", " << col << ").");
 			throw storm::exceptions::OutOfRangeException("Trying to read a value from illegal position.");
 			return false;
 		}
 
-		// Read elements on the diagonal directly.
-		if (row == col) {
-			*target = diagonalStorage[row];
-			return true;
-		}
-
 		// In case the element is not on the diagonal, we have to iterate
 		// over the accessed row to find the element.
 		uint_fast64_t rowStart = rowIndications[row];
@@ -405,17 +425,12 @@ public:
 	 */
 	inline T& getValue(uint_fast64_t row, uint_fast64_t col) {
 		// Check for illegal access indices.
-		if ((row > rowCount) || (col > rowCount)) {
+		if ((row > rowCount) || (col > colCount)) {
 			LOG4CPLUS_ERROR(logger, "Trying to read a value from illegal position (" << row << ", " << col << ").");
 			throw storm::exceptions::OutOfRangeException("Trying to read a value from illegal position.");
 		}
 
-		// Read elements on the diagonal directly.
-		if (row == col) {
-			return diagonalStorage[row];
-		}
-
-		// In case the element is not on the diagonal, we have to iterate
+		// we have to iterate
 		// over the accessed row to find the element.
 		uint_fast64_t rowStart = rowIndications[row];
 		uint_fast64_t rowEnd = rowIndications[row + 1];
@@ -445,20 +460,18 @@ public:
 	}
 
 	/*!
-	 * Returns a pointer to the value storage of the matrix. This storage does
-	 * *not* include elements on the diagonal.
-	 * @return A pointer to the value storage of the matrix.
+	 * Returns the number of columns of the matrix.
 	 */
-	T* getStoragePointer() const {
-		return valueStorage;
+	uint_fast64_t getColumnCount() const {
+		return colCount;
 	}
 
 	/*!
-	 * Returns a pointer to the storage of elements on the diagonal.
-	 * @return A pointer to the storage of elements on the diagonal.
+	 * Returns a pointer to the value storage of the matrix.
+	 * @return A pointer to the value storage of the matrix.
 	 */
-	T* getDiagonalStoragePointer() const {
-		return diagonalStorage;
+	std::vector<T> const & getStoragePointer() const {
+		return valueStorage;
 	}
 
 	/*!
@@ -467,17 +480,17 @@ public:
 	 * @return A pointer to the array that stores the start indices of non-zero
 	 * entries in the value storage for each row.
 	 */
-	uint_fast64_t* getRowIndicationsPointer() const {
+	std::vector<uint_fast64_t> const & getRowIndicationsPointer() const {
 		return rowIndications;
 	}
 
 	/*!
 	 * Returns a pointer to an array that stores the column of each non-zero
-	 * element that is not on the diagonal.
+	 * element.
 	 * @return A pointer to an array that stores the column of each non-zero
-	 * element that is not on the diagonal.
+	 * element.
 	 */
-	uint_fast64_t* getColumnIndicationsPointer() const {
+	std::vector<uint_fast64_t> const & getColumnIndicationsPointer() const {
 		return columnIndications;
 	}
 
@@ -548,10 +561,6 @@ public:
 #define STORM_USE_TRIPLETCONVERT
 #			ifdef STORM_USE_TRIPLETCONVERT
 
-			// FIXME: Wouldn't it be more efficient to add the elements in
-			// order including the diagonal elements? Otherwise, Eigen has to
-			// perform some sorting.
-
 			// Prepare the triplet storage.
 			typedef Eigen::Triplet<T> IntTriplet;
 			std::vector<IntTriplet> tripletList;
@@ -572,12 +581,6 @@ public:
 				}
 			}
 
-			// Then add the elements on the diagonal.
-			for (uint_fast64_t i = 0; i < rowCount; ++i) {
-				if (diagonalStorage[i] == 0) zeroCount++;
-				tripletList.push_back(IntTriplet(static_cast<int_fast32_t>(i), static_cast<int_fast32_t>(i), diagonalStorage[i]));
-			}
-
 			// Let Eigen create a matrix from the given list of triplets.
 			mat->setFromTriplets(tripletList.begin(), tripletList.end());
 
@@ -596,10 +599,6 @@ public:
 				rowStart = rowIndications[row];
 				rowEnd = rowIndications[row + 1];
 
-				// Insert the element on the diagonal.
-				mat->insert(row, row) = diagonalStorage[row];
-				count++;
-
 				// Insert the elements that are not on the diagonal
 				while (rowStart < rowEnd) {
 					mat->insert(row, columnIndications[rowStart]) = valueStorage[rowStart];
@@ -628,19 +627,6 @@ public:
 		return nonZeroEntryCount;
 	}
 
-	/*!
-	 * Returns the number of non-zero entries on the diagonal.
-	 * @return The number of non-zero entries on the diagonal.
-	 */
-	uint_fast64_t getDiagonalNonZeroEntryCount() const {
-		uint_fast64_t result = 0;
-		T zero(0);
-		for (uint_fast64_t i = 0; i < rowCount; ++i) {
-			if (diagonalStorage[i] != zero) ++result;
-		}
-		return result;
-	}
-
 	/*!
 	 * This function makes the rows given by the bit vector absorbing.
 	 * @param rows A bit vector indicating which rows to make absorbing.
@@ -658,7 +644,7 @@ public:
 	/*!
 	 * This function makes the given row absorbing. This means that all
 	 * entries in will be set to 0 and the value 1 will be written
-	 * to the element on the diagonal.
+	 * to the element on the (pseudo-) diagonal.
 	 * @param row The row to be made absorbing.
 	 * @returns True iff the operation was successful.
 	 */
@@ -675,13 +661,31 @@ public:
 		uint_fast64_t rowStart = rowIndications[row];
 		uint_fast64_t rowEnd = rowIndications[row + 1];
 
+		if (rowStart >= rowEnd) {
+			LOG4CPLUS_ERROR(logger, "The row " << row << " can not be made absorbing, no state in row, would have to recreate matrix!");
+			throw storm::exceptions::InvalidStateException("A row can not be made absorbing, no state in row, would have to recreate matrix!");
+		}
+		uint_fast64_t pseudoDiagonal = row % colCount;
+
+		bool foundDiagonal = false;
 		while (rowStart < rowEnd) {
-			valueStorage[rowStart] = storm::utility::constGetZero<T>();
+			if (!foundDiagonal && columnIndications[rowStart] >= pseudoDiagonal) {
+				foundDiagonal = true;
+				// insert/replace the diagonal entry
+				columnIndications[rowStart] = pseudoDiagonal;
+				valueStorage[rowStart] = storm::utility::constGetOne<T>();
+			} else {
+				valueStorage[rowStart] = storm::utility::constGetZero<T>();
+			}
 			++rowStart;
 		}
 
-		// Set the element on the diagonal to one.
-		diagonalStorage[row] = storm::utility::constGetOne<T>();
+		if (!foundDiagonal) {
+			--rowStart;
+			columnIndications[rowStart] = pseudoDiagonal;
+			valueStorage[rowStart] = storm::utility::constGetOne<T>();
+		}
+
 		return true;
 	}
 
@@ -724,7 +728,7 @@ public:
 	 * @param constraint A bit vector indicating which rows and columns to drop.
 	 * @return A pointer to a sparse matrix that is a sub-matrix of the current one.
 	 */
-	SquareSparseMatrix* getSubmatrix(storm::storage::BitVector& constraint) const {
+	SparseMatrix* getSubmatrix(storm::storage::BitVector& constraint) const {
 		LOG4CPLUS_DEBUG(logger, "Creating a sub-matrix with " << constraint.getNumberOfSetBits() << " rows.");
 
 		// Check for valid constraint.
@@ -745,7 +749,7 @@ public:
 		}
 
 		// Create and initialize resulting matrix.
-		SquareSparseMatrix* result = new SquareSparseMatrix(constraint.getNumberOfSetBits());
+		SparseMatrix* result = new SparseMatrix(constraint.getNumberOfSetBits());
 		result->initialize(subNonZeroEntries);
 
 		// Create a temporary array that stores for each index whose bit is set
@@ -763,8 +767,6 @@ public:
 		// Copy over selected entries.
 		uint_fast64_t rowCount = 0;
 		for (auto rowIndex : constraint) {
-			result->addNextValue(rowCount, rowCount, diagonalStorage[rowIndex]);
-
 			for (uint_fast64_t i = rowIndications[rowIndex]; i < rowIndications[rowIndex + 1]; ++i) {
 				if (constraint.get(columnIndications[i])) {
 					result->addNextValue(rowCount, bitsSetBeforeIndex[columnIndications[i]], valueStorage[i]);
@@ -794,8 +796,17 @@ public:
 	 */
 	void invertDiagonal() {
 		T one(1);
-		for (uint_fast64_t i = 0; i < rowCount; ++i) {
-			diagonalStorage[i] = one - diagonalStorage[i];
+		for (uint_fast64_t row = 0; row < rowCount; ++row) {
+			uint_fast64_t rowStart = rowIndications[row];
+			uint_fast64_t rowEnd = rowIndications[row + 1];
+			uint_fast64_t pseudoDiagonal = row % colCount;
+			while (rowStart < rowEnd) {
+				if (columnIndications[rowStart] == pseudoDiagonal) {
+					valueStorage[rowStart] = one - valueStorage[rowStart];
+					break;
+				}
+				++rowStart;
+			}
 		}
 	}
 
@@ -803,8 +814,16 @@ public:
 	 * Negates all non-zero elements that are not on the diagonal.
 	 */
 	void negateAllNonDiagonalElements() {
-		for (uint_fast64_t i = 0; i < nonZeroEntryCount; ++i) {
-			valueStorage[i] = - valueStorage[i];
+		for (uint_fast64_t row = 0; row < rowCount; ++row) {
+			uint_fast64_t rowStart = rowIndications[row];
+			uint_fast64_t rowEnd = rowIndications[row + 1];
+			uint_fast64_t pseudoDiagonal = row % colCount;
+			while (rowStart < rowEnd) {
+				if (columnIndications[rowStart] != pseudoDiagonal) {
+					valueStorage[rowStart] = - valueStorage[rowStart];
+				}
+				++rowStart;
+			}
 		}
 	}
 
@@ -814,8 +833,8 @@ public:
 	 */
 	storm::storage::JacobiDecomposition<T>* getJacobiDecomposition() const {
 		uint_fast64_t rowCount = this->getRowCount();
-		SquareSparseMatrix<T> *resultLU = new SquareSparseMatrix<T>(this);
-		SquareSparseMatrix<T> *resultDinv = new SquareSparseMatrix<T>(rowCount);
+		SparseMatrix<T> *resultLU = new SparseMatrix<T>(this);
+		SparseMatrix<T> *resultDinv = new SparseMatrix<T>(rowCount);
 		// no entries apart from those on the diagonal
 		resultDinv->initialize(0);
 		// copy diagonal entries to other matrix
@@ -836,7 +855,7 @@ public:
 	 * @return A vector containing the sum of the elements in each row of the matrix resulting from
 	 * pointwise multiplication of the current matrix with the given matrix.
 	 */
-	std::vector<T>* getPointwiseProductRowSumVector(storm::storage::SquareSparseMatrix<T> const& otherMatrix) {
+	std::vector<T>* getPointwiseProductRowSumVector(storm::storage::SparseMatrix<T> const& otherMatrix) {
 		// Prepare result.
 		std::vector<T>* result = new std::vector<T>(rowCount);
 
@@ -844,7 +863,6 @@ public:
 		// in case the given matrix does not have a non-zero element at this column position, or
 		// multiply the two entries and add the result to the corresponding position in the vector.
 		for (uint_fast64_t row = 0; row < rowCount && row < otherMatrix.rowCount; ++row) {
-			(*result)[row] += diagonalStorage[row] * otherMatrix.diagonalStorage[row];
 			for (uint_fast64_t element = rowIndications[row], nextOtherElement = otherMatrix.rowIndications[row]; element < rowIndications[row + 1] && nextOtherElement < otherMatrix.rowIndications[row + 1]; ++element) {
 				if (columnIndications[element] < otherMatrix.columnIndications[nextOtherElement]) {
 					continue;
@@ -868,25 +886,23 @@ public:
 	uint_fast64_t getSizeInMemory() const {
 		uint_fast64_t size = sizeof(*this);
 		// Add value_storage size.
-		size += sizeof(T) * nonZeroEntryCount;
-		// Add diagonal_storage size.
-		size += sizeof(T) * (rowCount + 1);
+		size += sizeof(T) * valueStorage.capacity();
 		// Add column_indications size.
-		size += sizeof(uint_fast64_t) * nonZeroEntryCount;
+		size += sizeof(uint_fast64_t) * columnIndications.capacity();
 		// Add row_indications size.
-		size += sizeof(uint_fast64_t) * (rowCount + 1);
+		size += sizeof(uint_fast64_t) * rowIndications.capacity();
 		return size;
 	}
 
 	/*!
 	 * Returns an iterator to the columns of the non-zero entries of the given
-	 * row that are not on the diagonal.
+	 * row.
 	 * @param row The row whose columns the iterator will return.
 	 * @return An iterator to the columns of the non-zero entries of the given
-	 * row that are not on the diagonal.
+	 * row.
 	 */
-	constIndexIterator beginConstColumnNoDiagIterator(uint_fast64_t row) const {
-		return this->columnIndications + this->rowIndications[row];
+	constIndexIterator beginConstColumnIterator(uint_fast64_t row) const {
+		return &(this->columnIndications[0]) + this->rowIndications[row];
 	}
 
 	/*!
@@ -894,18 +910,18 @@ public:
 	 * @param row The row for which the iterator should point to the past-the-end
 	 * element.
 	 */
-	constIndexIterator endConstColumnNoDiagIterator(uint_fast64_t row) const {
-		return this->columnIndications + this->rowIndications[row + 1];
+	constIndexIterator endConstColumnIterator(uint_fast64_t row) const {
+		return &(this->columnIndications[0]) + this->rowIndications[row + 1];
 	}
 	
 	/*!
 	 *	Returns an iterator over the elements of the given row. The iterator
-	 *	will include neither the diagonal element nor zero entries.
+	 *	will include no zero entries.
 	 *	@param row The row whose elements the iterator will return.
 	 *	@return An iterator over the elements of the given row.
 	 */
-	constIterator beginConstNoDiagIterator(uint_fast64_t row) const {
-		return this->valueStorage + this->rowIndications[row];
+	constIterator beginConstIterator(uint_fast64_t row) const {
+		return &(this->valueStorage[0]) + this->rowIndications[row];
 	}
 	/*!
 	 *	Returns an iterator pointing to the first element after the given
@@ -914,32 +930,28 @@ public:
 	 *	past-the-end element.
 	 *	@return An iterator to the element after the given row.
 	 */
-	constIterator endConstNoDiagIterator(uint_fast64_t row) const {
-		return this->valueStorage + this->rowIndications[row + 1];
+	constIterator endConstIterator(uint_fast64_t row) const {
+		return &(this->valueStorage[0]) + this->rowIndications[row + 1];
 	}
 	
 	/*!
 	 *	@brief Calculate sum of all entries in given row.
 	 *
-	 *	Adds up all values in the given row (including the diagonal value)
+	 *	Adds up all values in the given row
 	 *	and returns the sum.
 	 *	@param row The row that should be added up.
 	 *	@return Sum of the row.
 	 */
 	T getRowSum(uint_fast64_t row) const {
-		T sum = this->diagonalStorage[row];
-		for (auto it = this->beginConstNoDiagIterator(row); it != this->endConstNoDiagIterator(row); it++) {
+		T sum = storm::utility::constGetZero<T>();
+		for (auto it = this->beginConstIterator(row); it != this->endConstIterator(row); it++) {
 			sum += *it;
 		}
 		return sum;
 	}
 
 	void print() const {
-		std::cout << "diag: --------------------------------" << std::endl;
-		for (uint_fast64_t i = 0; i < rowCount; ++i) {
-			std::cout << "(" << i << "," << i << ") = " << diagonalStorage[i] << std::endl;
-		}
-		std::cout << "non diag: ----------------------------" << std::endl;
+		std::cout << "entries: ----------------------------" << std::endl;
 		for (uint_fast64_t i = 0; i < rowCount; ++i) {
 			for (uint_fast64_t j = rowIndications[i]; j < rowIndications[i + 1]; ++j) {
 				std::cout << "(" << i << "," << columnIndications[j] << ") = " << valueStorage[j] << std::endl;
@@ -955,31 +967,31 @@ private:
 	uint_fast64_t rowCount;
 
 	/*!
-	 * The number of non-zero elements that are not on the diagonal.
+	 * The number of columns of the matrix.
 	 */
-	uint_fast64_t nonZeroEntryCount;
+	uint_fast64_t colCount;
 
 	/*!
-	 * Stores all non-zero values that are not on the diagonal.
+	 * The number of non-zero elements.
 	 */
-	T* valueStorage;
+	uint_fast64_t nonZeroEntryCount;
 
 	/*!
-	 * Stores all elements on the diagonal, even the ones that are zero.
+	 * Stores all non-zero values.
 	 */
-	T* diagonalStorage;
+	std::vector<T> valueStorage;
 
 	/*!
-	 * Stores the column for each non-zero element that is not on the diagonal.
+	 * Stores the column for each non-zero element.
 	 */
-	uint_fast64_t* columnIndications;
+	std::vector<uint_fast64_t> columnIndications;
 
 	/*!
-	 * Array containing the boundaries (indices) in the value_storage array
+	 * Vector containing the boundaries (indices) in the value_storage array
 	 * for each row. All elements of value_storage with indices between the
 	 * i-th and the (i+1)-st element of this array belong to row i.
 	 */
-	uint_fast64_t* rowIndications;
+	std::vector<uint_fast64_t> rowIndications;
 
 	/*!
 	 * The internal status of the matrix.
@@ -1017,24 +1029,37 @@ private:
 	/*!
 	 * Prepares the internal CSR storage. For this, it requires
 	 * non_zero_entry_count and row_count to be set correctly.
+	 * @param alsoPerformAllocation If set to true, all entries are pre-allocated. This is the default.
 	 * @return True on success, false otherwise (allocation failed).
 	 */
-	bool prepareInternalStorage() {
-		// Set up the arrays for the elements that are not on the diagonal.
-		valueStorage = new (std::nothrow) T[nonZeroEntryCount]();
-		columnIndications = new (std::nothrow) uint_fast64_t[nonZeroEntryCount]();
-
-		// Set up the row_indications array and reserve one element more than
-		// there are rows in order to put a sentinel element at the end,
-		// which eases iteration process.
-		rowIndications = new (std::nothrow) uint_fast64_t[rowCount + 1]();
-
-		// Set up the array for the elements on the diagonal.
-		diagonalStorage = new (std::nothrow) T[rowCount]();
+	bool prepareInternalStorage(const bool alsoPerformAllocation) {
+		if (alsoPerformAllocation) {
+			// Set up the arrays for the elements that are not on the diagonal.
+			valueStorage.resize(nonZeroEntryCount, storm::utility::constGetZero<T>());
+			columnIndications.resize(nonZeroEntryCount, 0);
+
+			// Set up the row_indications vector and reserve one element more than
+			// there are rows in order to put a sentinel element at the end,
+			// which eases iteration process.
+			rowIndications.resize(rowCount + 1, 0);
+
+			// Return whether all the allocations could be made without error.
+			return ((valueStorage.capacity() >= nonZeroEntryCount) && (columnIndications.capacity() >= nonZeroEntryCount)
+				&& (rowIndications.capacity() >= (rowCount + 1)));
+		} else {
+			valueStorage.reserve(nonZeroEntryCount);
+			columnIndications.reserve(nonZeroEntryCount);
+			rowIndications.reserve(rowCount + 1);
+			return true;
+		}
+	}
 
-		// Return whether all the allocations could be made without error.
-		return ((valueStorage != NULL) && (columnIndications != NULL)
-				&& (rowIndications != NULL) && (diagonalStorage != NULL));
+	/*!
+	 * Shorthand for prepareInternalStorage(true)
+	 * @see prepareInternalStorage(const bool)
+	 */
+	bool prepareInternalStorage() {
+		return this->prepareInternalStorage(true);
 	}
 
 	/*!
@@ -1060,57 +1085,10 @@ private:
 		return false;
 	}
 
-	/*!
-	 * Helper function to determine the number of non-zero elements that are
-	 * not on the diagonal of the given Eigen matrix.
-	 * @param eigen_sparse_matrix The Eigen matrix to analyze.
-	 * @return The number of non-zero elements that are not on the diagonal of
-	 * the given Eigen matrix.
-	 */
-	template<typename _Scalar, int _Options, typename _Index>
-	_Index getEigenSparseMatrixCorrectNonZeroEntryCount(const Eigen::SparseMatrix<_Scalar, _Options, _Index>& eigen_sparse_matrix) const {
-		const _Index* indexPtr = eigen_sparse_matrix.innerIndexPtr();
-		const _Index* outerPtr = eigen_sparse_matrix.outerIndexPtr();
-
-		const _Index entryCount = eigen_sparse_matrix.nonZeros();
-		const _Index outerCount = eigen_sparse_matrix.outerSize();
-
-		uint_fast64_t diagNonZeros = 0;
-
-		// For RowMajor, row is the current row and col the column and for
-		// ColMajor, row is the current column and col the row, but this is
-		// not important as we are only looking for elements on the diagonal.
-		_Index innerStart = 0;
-		_Index innerEnd = 0;
-		_Index innerMid = 0;
-		for (_Index row = 0; row < outerCount; ++row) {
-			innerStart = outerPtr[row];
-			innerEnd = outerPtr[row + 1] - 1;
-
-			// Now use binary search (but defer equality detection).
-			while (innerStart < innerEnd) {
-				innerMid = innerStart + ((innerEnd - innerStart) / 2);
-
-				if (indexPtr[innerMid] < row) {
-					innerStart = innerMid + 1;
-				} else {
-					innerEnd = innerMid;
-				}
-			}
-
-			// Check whether we have found an element on the diagonal.
-			if ((innerStart == innerEnd) && (indexPtr[innerStart] == row)) {
-				++diagNonZeros;
-			}
-		}
-
-		return static_cast<_Index>(entryCount - diagNonZeros);
-	}
-
 };
 
 } // namespace storage
 
 } // namespace storm
 
-#endif // STORM_STORAGE_SQUARESPARSEMATRIX_H_
+#endif // STORM_STORAGE_SPARSEMATRIX_H_
diff --git a/src/storm.cpp b/src/storm.cpp
index 0c5951961..cd7fda7e0 100644
--- a/src/storm.cpp
+++ b/src/storm.cpp
@@ -20,7 +20,7 @@
 
 #include "storm-config.h"
 #include "src/models/Dtmc.h"
-#include "src/storage/SquareSparseMatrix.h"
+#include "src/storage/SparseMatrix.h"
 #include "src/models/AtomicPropositionsLabeling.h"
 #include "src/modelChecker/EigenDtmcPrctlModelChecker.h"
 #include "src/modelChecker/GmmxxDtmcPrctlModelChecker.h"
diff --git a/src/utility/IoUtility.cpp b/src/utility/IoUtility.cpp
index e3fc7dbaf..d2f78a25a 100644
--- a/src/utility/IoUtility.cpp
+++ b/src/utility/IoUtility.cpp
@@ -1,9 +1,9 @@
 /*
- * IoUtility.cpp
- *
- *  Created on: 17.10.2012
- *      Author: Thomas Heinemann
- */
+* IoUtility.cpp
+*
+*  Created on: 17.10.2012
+*      Author: Thomas Heinemann
+*/
 
 #include "src/utility/IoUtility.h"
 #include "src/parser/DeterministicSparseTransitionParser.h"
@@ -13,63 +13,58 @@
 
 namespace storm {
 
-namespace utility {
-
-void dtmcToDot(storm::models::Dtmc<double> const &dtmc, std::string filename) {
-   std::shared_ptr<storm::storage::SquareSparseMatrix<double>> matrix(dtmc.getTransitionProbabilityMatrix());
-   double* diagonal_storage = matrix->getDiagonalStoragePointer();
-
-   std::ofstream file;
-   file.open(filename);
-
-   file << "digraph dtmc {\n";
-
-   //Specify the nodes and their labels
-   for (uint_fast64_t i = 1; i < dtmc.getNumberOfStates(); i++) {
-      file << "\t" << i << "[label=\"" << i << "\\n{";
-      char komma=' ';
-      std::set<std::string> propositions = dtmc.getPropositionsForState(i);
-      for(auto it = propositions.begin();
-               it != propositions.end();
-               it++) {
-      	file << komma << *it;
-         komma=',';
-      }
-
-      file << " }\"];\n";
-
-   }
-
-   for (uint_fast64_t row = 0; row < dtmc.getNumberOfStates(); row++ ) {
-   	//write diagonal entry/self loop first
-   	if (diagonal_storage[row] != 0) {
-   	            file << "\t" << row << " -> " << row << " [label=" << diagonal_storage[row] <<"]\n";
-   	         }
-   	//Then, iterate through the row and write each non-diagonal value into the file
-   	for (	auto it = matrix->beginConstColumnNoDiagIterator(row);
-   			it !=  matrix->endConstColumnNoDiagIterator(row);
-   			it++) {
-   		double value = 0;
-   		matrix->getValue(row,*it,&value);
-   		file << "\t" << row << " -> " << *it << " [label=" << value << "]\n";
-   	}
-   }
-
-   file << "}\n";
-   file.close();
-}
+	namespace utility {
 
-//TODO: Should this stay here or be integrated in the new parser structure?
-/*storm::models::Dtmc<double>* parseDTMC(std::string const &tra_file, std::string const &lab_file) {
-	storm::parser::DeterministicSparseTransitionParser tp(tra_file);
-	uint_fast64_t node_count = tp.getMatrix()->getRowCount();
+		void dtmcToDot(storm::models::Dtmc<double> const &dtmc, std::string filename) {
+			std::shared_ptr<storm::storage::SparseMatrix<double>> matrix(dtmc.getTransitionProbabilityMatrix());
+			std::ofstream file;
+			file.open(filename);
 
-	storm::parser::AtomicPropositionLabelingParser lp(node_count, lab_file);
+			file << "digraph dtmc {\n";
 
-	storm::models::Dtmc<double>* result = new storm::models::Dtmc<double>(tp.getMatrix(), lp.getLabeling());
-	return result;
-}*/
+			//Specify the nodes and their labels
+			for (uint_fast64_t i = 1; i < dtmc.getNumberOfStates(); i++) {
+				file << "\t" << i << "[label=\"" << i << "\\n{";
+				char komma=' ';
+				std::set<std::string> propositions = dtmc.getPropositionsForState(i);
+				for(auto it = propositions.begin();
+					it != propositions.end();
+					it++) {
+						file << komma << *it;
+						komma=',';
+				}
 
-}
+				file << " }\"];\n";
+
+			}
+
+			for (uint_fast64_t row = 0; row < dtmc.getNumberOfStates(); row++ ) {
+
+				//Then, iterate through the row and write each non-diagonal value into the file
+				for (	auto it = matrix->beginConstColumnIterator(row);
+					it !=  matrix->endConstColumnIterator(row);
+					it++) {
+						double value = 0;
+						matrix->getValue(row,*it,&value);
+						file << "\t" << row << " -> " << *it << " [label=" << value << "]\n";
+				}
+			}
+
+			file << "}\n";
+			file.close();
+		}
+
+		//TODO: Should this stay here or be integrated in the new parser structure?
+		/*storm::models::Dtmc<double>* parseDTMC(std::string const &tra_file, std::string const &lab_file) {
+		storm::parser::DeterministicSparseTransitionParser tp(tra_file);
+		uint_fast64_t node_count = tp.getMatrix()->getRowCount();
+
+		storm::parser::AtomicPropositionLabelingParser lp(node_count, lab_file);
+
+		storm::models::Dtmc<double>* result = new storm::models::Dtmc<double>(tp.getMatrix(), lp.getLabeling());
+		return result;
+		}*/
+
+	}
 
 }
diff --git a/test/parser/ReadTraFileTest.cpp b/test/parser/ReadTraFileTest.cpp
index 77377d857..2bda7561e 100644
--- a/test/parser/ReadTraFileTest.cpp
+++ b/test/parser/ReadTraFileTest.cpp
@@ -7,7 +7,7 @@
 
 #include "gtest/gtest.h"
 #include "storm-config.h"
-#include "src/storage/SquareSparseMatrix.h"
+#include "src/storage/SparseMatrix.h"
 #include "src/parser/DeterministicSparseTransitionParser.h"
 #include "src/exceptions/FileIoException.h"
 #include "src/exceptions/WrongFileFormatException.h"
@@ -24,7 +24,7 @@ TEST(ReadTraFileTest, NonExistingFileTest) {
 TEST(ReadTraFileTest, ParseFileTest1) {
 	storm::parser::DeterministicSparseTransitionParser* parser;
 	ASSERT_NO_THROW(parser = new storm::parser::DeterministicSparseTransitionParser(STORM_CPP_TESTS_BASE_PATH "/parser/tra_files/csl_general_input_01.tra"));
-	std::shared_ptr<storm::storage::SquareSparseMatrix<double>> result = parser->getMatrix();
+	std::shared_ptr<storm::storage::SparseMatrix<double>> result = parser->getMatrix();
 
 	if (result != NULL) {
 		double val = 0;
@@ -53,13 +53,13 @@ TEST(ReadTraFileTest, ParseFileTest1) {
 		ASSERT_TRUE(result->getValue(3,2,&val));
 		ASSERT_EQ(val,0.0806451612903225806451612903225812);
 
-		ASSERT_TRUE(result->getValue(3,3,&val));
+		ASSERT_FALSE(result->getValue(3,3,&val));
 		ASSERT_EQ(val,0);
 
 		ASSERT_TRUE(result->getValue(3,4,&val));
 		ASSERT_EQ(val,0.080645161290322580645161290322581);
 
-		ASSERT_TRUE(result->getValue(4,4,&val));
+		ASSERT_FALSE(result->getValue(4,4,&val));
 		ASSERT_EQ(val,0);
 
 		delete parser;
diff --git a/test/storage/SquareSparseMatrixTest.cpp b/test/storage/SparseMatrixTest.cpp
similarity index 58%
rename from test/storage/SquareSparseMatrixTest.cpp
rename to test/storage/SparseMatrixTest.cpp
index 61d64c933..2a29fb736 100644
--- a/test/storage/SquareSparseMatrixTest.cpp
+++ b/test/storage/SparseMatrixTest.cpp
@@ -1,73 +1,73 @@
 #include "gtest/gtest.h"
-#include "src/storage/SquareSparseMatrix.h"
+#include "src/storage/SparseMatrix.h"
 #include "src/exceptions/InvalidArgumentException.h"
 #include "src/exceptions/OutOfRangeException.h"
 
-TEST(SquareSparseMatrixTest, ZeroRowsTest) {
-	storm::storage::SquareSparseMatrix<int> *ssm = new storm::storage::SquareSparseMatrix<int>(0);
-	ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::UnInitialized);
+TEST(SparseMatrixTest, ZeroRowsTest) {
+	storm::storage::SparseMatrix<int> *ssm = new storm::storage::SparseMatrix<int>(0);
+	ASSERT_EQ(ssm->getState(), storm::storage::SparseMatrix<int>::MatrixStatus::UnInitialized);
 
 	ASSERT_THROW(ssm->initialize(50), storm::exceptions::InvalidArgumentException);
-	ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::Error);
+	ASSERT_EQ(ssm->getState(), storm::storage::SparseMatrix<int>::MatrixStatus::Error);
 
 	delete ssm;
 }
 
-TEST(SquareSparseMatrixTest, TooManyEntriesTest) {
-	storm::storage::SquareSparseMatrix<int> *ssm = new storm::storage::SquareSparseMatrix<int>(2);
-	ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::UnInitialized);
+TEST(SparseMatrixTest, TooManyEntriesTest) {
+	storm::storage::SparseMatrix<int> *ssm = new storm::storage::SparseMatrix<int>(2);
+	ASSERT_EQ(ssm->getState(), storm::storage::SparseMatrix<int>::MatrixStatus::UnInitialized);
 
 	ASSERT_THROW(ssm->initialize(10), storm::exceptions::InvalidArgumentException);
-	ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::Error);
+	ASSERT_EQ(ssm->getState(), storm::storage::SparseMatrix<int>::MatrixStatus::Error);
 
 	delete ssm;
 }
 
-TEST(SquareSparseMatrixTest, addNextValueTest) {
-	storm::storage::SquareSparseMatrix<int> *ssm = new storm::storage::SquareSparseMatrix<int>(5);
-	ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::UnInitialized);
+TEST(SparseMatrixTest, addNextValueTest) {
+	storm::storage::SparseMatrix<int> *ssm = new storm::storage::SparseMatrix<int>(5);
+	ASSERT_EQ(ssm->getState(), storm::storage::SparseMatrix<int>::MatrixStatus::UnInitialized);
 
 	ASSERT_NO_THROW(ssm->initialize(1));
-	ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::Initialized);
+	ASSERT_EQ(ssm->getState(), storm::storage::SparseMatrix<int>::MatrixStatus::Initialized);
 
 	ASSERT_THROW(ssm->addNextValue(-1, 1, 1), storm::exceptions::OutOfRangeException);
-	ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::Error);
+	ASSERT_EQ(ssm->getState(), storm::storage::SparseMatrix<int>::MatrixStatus::Error);
 
 	ASSERT_THROW(ssm->addNextValue(1, -1, 1), storm::exceptions::OutOfRangeException);
-	ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::Error);
+	ASSERT_EQ(ssm->getState(), storm::storage::SparseMatrix<int>::MatrixStatus::Error);
 
 	ASSERT_THROW(ssm->addNextValue(6, 1, 1), storm::exceptions::OutOfRangeException);
-	ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::Error);
+	ASSERT_EQ(ssm->getState(), storm::storage::SparseMatrix<int>::MatrixStatus::Error);
 
 	ASSERT_THROW(ssm->addNextValue(1, 6, 1), storm::exceptions::OutOfRangeException);
-	ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::Error);
+	ASSERT_EQ(ssm->getState(), storm::storage::SparseMatrix<int>::MatrixStatus::Error);
 
 	delete ssm;
 }
 
-TEST(SquareSparseMatrixTest, finalizeTest) {
-	storm::storage::SquareSparseMatrix<int> *ssm = new storm::storage::SquareSparseMatrix<int>(5);
-	ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::UnInitialized);
+TEST(SparseMatrixTest, finalizeTest) {
+	storm::storage::SparseMatrix<int> *ssm = new storm::storage::SparseMatrix<int>(5);
+	ASSERT_EQ(ssm->getState(), storm::storage::SparseMatrix<int>::MatrixStatus::UnInitialized);
 
 	ASSERT_NO_THROW(ssm->initialize(5));
-	ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::Initialized);
+	ASSERT_EQ(ssm->getState(), storm::storage::SparseMatrix<int>::MatrixStatus::Initialized);
 
 	ASSERT_NO_THROW(ssm->addNextValue(1, 2, 1));
 	ASSERT_NO_THROW(ssm->addNextValue(1, 3, 1));
 	ASSERT_NO_THROW(ssm->addNextValue(1, 4, 1));
 	ASSERT_NO_THROW(ssm->addNextValue(1, 5, 1));
-	ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::Initialized);
+	ASSERT_EQ(ssm->getState(), storm::storage::SparseMatrix<int>::MatrixStatus::Initialized);
 
 	ASSERT_THROW(ssm->finalize(), storm::exceptions::InvalidStateException);
-	ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::Error);
+	ASSERT_EQ(ssm->getState(), storm::storage::SparseMatrix<int>::MatrixStatus::Error);
 
 	delete ssm;
 }
 
-TEST(SquareSparseMatrixTest, Test) {
+TEST(SparseMatrixTest, Test) {
 	// 25 rows, 50 non zero entries
-	storm::storage::SquareSparseMatrix<int> *ssm = new storm::storage::SquareSparseMatrix<int>(25);
-	ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::UnInitialized);
+	storm::storage::SparseMatrix<int> *ssm = new storm::storage::SparseMatrix<int>(25);
+	ASSERT_EQ(ssm->getState(), storm::storage::SparseMatrix<int>::MatrixStatus::UnInitialized);
 
 	int values[50] = {
 		0, 1, 2, 3, 4, 5, 6, 7, 8, 9,
@@ -96,15 +96,15 @@ TEST(SquareSparseMatrixTest, Test) {
 	};
 
 	ASSERT_NO_THROW(ssm->initialize(50));
-	ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::Initialized);
+	ASSERT_EQ(ssm->getState(), storm::storage::SparseMatrix<int>::MatrixStatus::Initialized);
 
 	for (int i = 0; i < 50; ++i) {
 		ASSERT_NO_THROW(ssm->addNextValue(position_row[i], position_col[i], values[i]));
 	}
-	ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::Initialized);
+	ASSERT_EQ(ssm->getState(), storm::storage::SparseMatrix<int>::MatrixStatus::Initialized);
 
 	ASSERT_NO_THROW(ssm->finalize());
-	ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::ReadReady);
+	ASSERT_EQ(ssm->getState(), storm::storage::SparseMatrix<int>::MatrixStatus::ReadReady);
 
 	int target;
 	for (int i = 0; i < 50; ++i) {
@@ -116,24 +116,20 @@ TEST(SquareSparseMatrixTest, Test) {
 	for (int row = 15; row < 24; ++row) {
 		for (int col = 1; col <= 25; ++col) {
 		   target = 1;
-		   if (row != col) {
-		      ASSERT_FALSE(ssm->getValue(row, col, &target));
-		   } else {
-		      ASSERT_TRUE(ssm->getValue(row, col, &target));
-		   }
+		   ASSERT_FALSE(ssm->getValue(row, col, &target));
 
 		   ASSERT_EQ(0, target);
 		}
 	}
-	ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::ReadReady);
+	ASSERT_EQ(ssm->getState(), storm::storage::SparseMatrix<int>::MatrixStatus::ReadReady);
 
 	delete ssm;
 }
 
-TEST(SquareSparseMatrixTest, ConversionFromDenseEigen_ColMajor_SparseMatrixTest) {
+TEST(SparseMatrixTest, ConversionFromDenseEigen_ColMajor_SparseMatrixTest) {
 	// 10 rows, 100 non zero entries
-	storm::storage::SquareSparseMatrix<int> *ssm = new storm::storage::SquareSparseMatrix<int>(10);
-	ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::UnInitialized);
+	storm::storage::SparseMatrix<int> *ssm = new storm::storage::SparseMatrix<int>(10);
+	ASSERT_EQ(ssm->getState(), storm::storage::SparseMatrix<int>::MatrixStatus::UnInitialized);
 
 	Eigen::SparseMatrix<int> esm(10, 10);
 	for (int row = 0; row < 10; ++row) {
@@ -149,7 +145,7 @@ TEST(SquareSparseMatrixTest, ConversionFromDenseEigen_ColMajor_SparseMatrixTest)
 
 	ASSERT_NO_THROW(ssm->finalize());
 
-	ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::ReadReady);
+	ASSERT_EQ(ssm->getState(), storm::storage::SparseMatrix<int>::MatrixStatus::ReadReady);
 
 	int target = -1;
 	for (int row = 0; row < 10; ++row) {
@@ -162,10 +158,10 @@ TEST(SquareSparseMatrixTest, ConversionFromDenseEigen_ColMajor_SparseMatrixTest)
 	delete ssm;
 }
 
-TEST(SquareSparseMatrixTest, ConversionFromDenseEigen_RowMajor_SparseMatrixTest) {
+TEST(SparseMatrixTest, ConversionFromDenseEigen_RowMajor_SparseMatrixTest) {
 	// 10 rows, 100 non zero entries
-	storm::storage::SquareSparseMatrix<int> *ssm = new storm::storage::SquareSparseMatrix<int>(10);
-	ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::UnInitialized);
+	storm::storage::SparseMatrix<int> *ssm = new storm::storage::SparseMatrix<int>(10);
+	ASSERT_EQ(ssm->getState(), storm::storage::SparseMatrix<int>::MatrixStatus::UnInitialized);
 
 	Eigen::SparseMatrix<int, Eigen::RowMajor> esm(10, 10);
 	for (int row = 0; row < 10; ++row) {
@@ -181,7 +177,7 @@ TEST(SquareSparseMatrixTest, ConversionFromDenseEigen_RowMajor_SparseMatrixTest)
 
 	ASSERT_NO_THROW(ssm->finalize());
 
-	ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::ReadReady);
+	ASSERT_EQ(ssm->getState(), storm::storage::SparseMatrix<int>::MatrixStatus::ReadReady);
 
 	int target = -1;
 	for (int row = 0; row < 10; ++row) {
@@ -194,10 +190,10 @@ TEST(SquareSparseMatrixTest, ConversionFromDenseEigen_RowMajor_SparseMatrixTest)
 	delete ssm;
 }
 
-TEST(SquareSparseMatrixTest, ConversionFromSparseEigen_ColMajor_SparseMatrixTest) {
+TEST(SparseMatrixTest, ConversionFromSparseEigen_ColMajor_SparseMatrixTest) {
 	// 10 rows, 15 non zero entries
-	storm::storage::SquareSparseMatrix<int> *ssm = new storm::storage::SquareSparseMatrix<int>(10);
-	ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::UnInitialized);
+	storm::storage::SparseMatrix<int> *ssm = new storm::storage::SparseMatrix<int>(10);
+	ASSERT_EQ(ssm->getState(), storm::storage::SparseMatrix<int>::MatrixStatus::UnInitialized);
 
 	Eigen::SparseMatrix<int> esm(10, 10);
 	
@@ -231,7 +227,7 @@ TEST(SquareSparseMatrixTest, ConversionFromSparseEigen_ColMajor_SparseMatrixTest
 
 	ASSERT_NO_THROW(ssm->finalize());
 
-	ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::ReadReady);
+	ASSERT_EQ(ssm->getState(), storm::storage::SparseMatrix<int>::MatrixStatus::ReadReady);
 
 	int target = -1;
 
@@ -243,17 +239,17 @@ TEST(SquareSparseMatrixTest, ConversionFromSparseEigen_ColMajor_SparseMatrixTest
 	delete ssm;
 }
 
-TEST(SquareSparseMatrixTest, ConversionFromSparseEigen_RowMajor_SparseMatrixTest) {
+TEST(SparseMatrixTest, ConversionFromSparseEigen_RowMajor_SparseMatrixTest) {
 	// 10 rows, 15 non zero entries
-	storm::storage::SquareSparseMatrix<int> *ssm = new storm::storage::SquareSparseMatrix<int>(10);
-	ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::UnInitialized);
+	storm::storage::SparseMatrix<int> *ssm = new storm::storage::SparseMatrix<int>(10, 10);
+	ASSERT_EQ(ssm->getState(), storm::storage::SparseMatrix<int>::MatrixStatus::UnInitialized);
 
 	Eigen::SparseMatrix<int, Eigen::RowMajor> esm(10, 10);
 	
 	typedef Eigen::Triplet<int> IntTriplet;
 	std::vector<IntTriplet> tripletList;
 	tripletList.reserve(15);
-	tripletList.push_back(IntTriplet(1, 0, 0));
+	tripletList.push_back(IntTriplet(1, 0, 15));
 	tripletList.push_back(IntTriplet(1, 1, 1));
 	tripletList.push_back(IntTriplet(1, 2, 2));
 	tripletList.push_back(IntTriplet(1, 3, 3));
@@ -280,38 +276,42 @@ TEST(SquareSparseMatrixTest, ConversionFromSparseEigen_RowMajor_SparseMatrixTest
 
 	ASSERT_NO_THROW(ssm->finalize());
 
-	ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::ReadReady);
+	ASSERT_EQ(ssm->getState(), storm::storage::SparseMatrix<int>::MatrixStatus::ReadReady);
+	
+	const std::vector<uint_fast64_t> rowP = ssm->getRowIndicationsPointer();
+	const std::vector<uint_fast64_t> colP = ssm->getColumnIndicationsPointer();
+	const std::vector<int> valP = ssm->getStoragePointer();
 
 	int target = -1;
-
 	for (auto &coeff: tripletList) {
-		ASSERT_TRUE(ssm->getValue(coeff.row(), coeff.col(), &target));
+		bool retVal = ssm->getValue(coeff.row(), coeff.col(), &target);
+		ASSERT_TRUE(retVal);
 		ASSERT_EQ(target, coeff.value());
 	}
 	
 	delete ssm;
 }
 
-TEST(SquareSparseMatrixTest, ConversionToSparseEigen_RowMajor_SparseMatrixTest) {
+TEST(SparseMatrixTest, ConversionToSparseEigen_RowMajor_SparseMatrixTest) {
 	int values[100];
-	storm::storage::SquareSparseMatrix<int> *ssm = new storm::storage::SquareSparseMatrix<int>(10);
+	storm::storage::SparseMatrix<int> *ssm = new storm::storage::SparseMatrix<int>(10);
 
 	for (uint_fast32_t i = 0; i < 100; ++i) {
 		values[i] = i;
 	}
 
-	ASSERT_NO_THROW(ssm->initialize(100 - 10));
-	ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::Initialized);
+	ASSERT_NO_THROW(ssm->initialize(100));
+	ASSERT_EQ(ssm->getState(), storm::storage::SparseMatrix<int>::MatrixStatus::Initialized);
 
 	for (uint_fast32_t row = 0; row < 10; ++row) {
 		for (uint_fast32_t col = 0; col < 10; ++col) {
 			ASSERT_NO_THROW(ssm->addNextValue(row, col, values[row * 10 + col]));
 		}
 	}
-	ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::Initialized);
+	ASSERT_EQ(ssm->getState(), storm::storage::SparseMatrix<int>::MatrixStatus::Initialized);
 
 	ASSERT_NO_THROW(ssm->finalize());
-	ASSERT_EQ(ssm->getState(), storm::storage::SquareSparseMatrix<int>::MatrixStatus::ReadReady);
+	ASSERT_EQ(ssm->getState(), storm::storage::SparseMatrix<int>::MatrixStatus::ReadReady);
 
 	Eigen::SparseMatrix<int, Eigen::RowMajor, int_fast32_t>* esm = ssm->toEigenSparseMatrix();