Refactored the Jacobi Decomposition

Former-commit-id: 55d5d38475

resources/3rdparty/gmm-4.2/include/gmm/gmm_matrix.h

@@ -491,7 +491,7 @@ namespace gmm
 
   template <typename T, int shift = 0>
   struct csc_matrix {
-    typedef unsigned int IND_TYPE;
+    typedef unsigned long long IND_TYPE;
 
     std::vector<T> pr;
     std::vector<IND_TYPE> ir;
@@ -639,7 +639,7 @@ namespace gmm
   template <typename T, int shift = 0>
   struct csr_matrix {
 
-    typedef unsigned int IND_TYPE;
+	  typedef unsigned long long IND_TYPE;
 
     std::vector<T> pr;        // values.
     std::vector<IND_TYPE> ir; // col indices.

src/modelchecker/prctl/GmmxxDtmcPrctlModelChecker.h

@@ -10,7 +10,7 @@
 
 #include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h"
 #include "src/adapters/GmmxxAdapter.h"
-#include "src/storage/JacobiDecomposition.h"
+#include "src/storage/SparseMatrix.h"
 #include "src/utility/ConstTemplates.h"
 #include "src/utility/Settings.h"
 
@@ -241,12 +241,12 @@ private:
 		bool relative = s->get<bool>("relative");
 
 		// Get a Jacobi decomposition of the matrix A.
-		storm::storage::JacobiDecomposition<Type>* jacobiDecomposition = A.getJacobiDecomposition();
+		storm::storage::SparseMatrix<Type>::SparseJacobiDecomposition_t jacobiDecomposition = A.getJacobiDecomposition();
 
 		// Convert the (inverted) diagonal matrix to gmm++'s format.
-		gmm::csr_matrix<Type>* gmmxxDiagonalInverted = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(jacobiDecomposition->getJacobiDInvReference());
+		gmm::csr_matrix<Type>* gmmxxDiagonalInverted = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(std::move(jacobiDecomposition.second));
 		// Convert the LU matrix to gmm++'s format.
-		gmm::csr_matrix<Type>* gmmxxLU = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(jacobiDecomposition->getJacobiLUReference());
+		gmm::csr_matrix<Type>* gmmxxLU = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(std::move(jacobiDecomposition.first));
 
 		LOG4CPLUS_INFO(logger, "Starting iterative Jacobi Solver.");
 
@@ -295,7 +295,6 @@ private:
 
 		// Also delete the other dynamically allocated variables.
 		delete residuum;
-		delete jacobiDecomposition;
 		delete gmmxxDiagonalInverted;
 		delete gmmxxLU;
 

src/storage/JacobiDecomposition.h

@@ -1,105 +0,0 @@
-#ifndef STORM_STORAGE_JACOBIDECOMPOSITION_H_
-#define STORM_STORAGE_JACOBIDECOMPOSITION_H_
-
-#include "boost/integer/integer_mask.hpp"
-
-#include "log4cplus/logger.h"
-#include "log4cplus/loggingmacros.h"
-
-#include "src/exceptions/InvalidAccessException.h"
-
-extern log4cplus::Logger logger;
-
-namespace storm {
-
-namespace storage {
-
-/*!
- * Forward declaration against Cycle
- */
-template <class T>
-class SparseMatrix;
-
-
-/*!
- * A simple container for a sparse Jacobi decomposition
- */
-template <class T>
-class JacobiDecomposition {
-
-public:
-	/*!
-	 *	Standard constructor
-	 * Initializes this object with the two given sparse matrices
-	 * Ownership of both matrices stay with THIS object.
-	 */
-	JacobiDecomposition(storm::storage::SparseMatrix<T> * const jacobiLuMatrix, storm::storage::SparseMatrix<T> * const jacobiDInvMatrix) : jacobiLuMatrix(jacobiLuMatrix), jacobiDInvMatrix(jacobiDInvMatrix) {
-	}
-
-	~JacobiDecomposition() {
-		delete this->jacobiDInvMatrix;
-		delete this->jacobiLuMatrix;
-	}
-
-	/*!
-	 * Accessor for the internal LU Matrix.
-	 * Ownership stays with this class.
-	 * @return A reference to the Jacobi LU Matrix
-	 */
-	storm::storage::SparseMatrix<T>& getJacobiLUReference() {
-		return *(this->jacobiLuMatrix);
-	}
-
-	/*!
-	 * Accessor for the internal D^{-1} Matrix.
-	 * Ownership stays with this class.
-	 * @return A reference to the Jacobi D^{-1} Matrix
-	 */
-	storm::storage::SparseMatrix<T>& getJacobiDInvReference() {
-		return *(this->jacobiDInvMatrix);
-	}
-
-	/*!
-	 * Accessor for the internal LU Matrix.
-	 * Ownership stays with this class.
-	 * @return A pointer to the Jacobi LU Matrix
-	 */
-	storm::storage::SparseMatrix<T>* getJacobiLU() {
-		return this->jacobiLuMatrix;
-	}
-
-	/*!
-	 * Accessor for the internal D^{-1} Matrix.
-	 * Ownership stays with this class.
-	 * @return A pointer to the Jacobi D^{-1} Matrix
-	 */
-	storm::storage::SparseMatrix<T>* getJacobiDInv() {
-		return this->jacobiDInvMatrix;
-	}
-
-private:
-
-	/*!
-	 * The copy constructor is disabled for this class.
-	 */
-	//JacobiDecomposition(const JacobiDecomposition<T>& that) = delete; // not possible in VS2012
-	JacobiDecomposition(const JacobiDecomposition<T>& that) {
-		throw storm::exceptions::InvalidAccessException() << "The copy constructor of JacobiDecomposition is explicitly disabled.";
-	}
-
-	/*!
-	 * Pointer to the LU Matrix
-	 */
-	storm::storage::SparseMatrix<T> *jacobiLuMatrix;
-	
-	/*!
-	 * Pointer to the D^{-1} Matrix
-	 */
-	storm::storage::SparseMatrix<T> *jacobiDInvMatrix;
-};
-
-} // namespace storage
-
-} // namespace storm
-
-#endif // STORM_STORAGE_JACOBIDECOMPOSITION_H_

src/storage/SparseMatrix.h

@@ -14,7 +14,6 @@
 #include "src/exceptions/OutOfRangeException.h"
 #include "src/exceptions/FileIoException.h"
 #include "src/storage/BitVector.h"
-#include "src/storage/JacobiDecomposition.h"
 
 #include "src/utility/ConstTemplates.h"
 #include "src/utility/Hash.h"
@@ -46,6 +45,11 @@ namespace storage {
 template<class T>
 class SparseMatrix {
 public:
+	/*!
+	 * Return Type of the Jacobi Decompostion
+	 */
+	typedef std::pair<storm::storage::SparseMatrix<T>, storm::storage::SparseMatrix<T>> SparseJacobiDecomposition_t;
+
 	/*!
 	 * Declare adapter classes as friends to use internal data.
 	 */
@@ -938,29 +942,30 @@ public:
 	/*!
 	 * Calculates the Jacobi-Decomposition of this sparse matrix.
 	 * The source Sparse Matrix must be square.
-	 * @return A pointer to a class containing the matrix L+U and the inverted diagonal matrix D^-1
+	 * @return A std::pair containing the matrix L+U and the inverted diagonal matrix D^-1
 	 */
-	storm::storage::JacobiDecomposition<T>* getJacobiDecomposition() const {
+	SparseJacobiDecomposition_t getJacobiDecomposition() const {
 		uint_fast64_t rowCount = this->getRowCount();
 		uint_fast64_t colCount = this->getColumnCount();
 		if (rowCount != colCount) {
 			throw storm::exceptions::InvalidArgumentException() << "SparseMatrix::getJacobiDecomposition requires the Matrix to be square.";
 		}
-		storm::storage::SparseMatrix<T> *resultLU = new storm::storage::SparseMatrix<T>(*this);
-		storm::storage::SparseMatrix<T> *resultDinv = new storm::storage::SparseMatrix<T>(rowCount, colCount);
+		storm::storage::SparseMatrix<T> resultLU(*this);
+		storm::storage::SparseMatrix<T> resultDinv(rowCount, colCount);
 		// no entries apart from those on the diagonal (rowCount)
-		resultDinv->initialize(rowCount);
+		resultDinv.initialize(rowCount);
 
 		// constant 1 for diagonal inversion
 		T constOne = storm::utility::constGetOne<T>();
 
 		// copy diagonal entries to other matrix
-		for (unsigned int i = 0; i < rowCount; ++i) {
-			resultDinv->addNextValue(i, i, constOne / resultLU->getValue(i, i));
-			resultLU->getValue(i, i) = storm::utility::constGetZero<T>();
+		for (uint_fast64_t i = 0; i < rowCount; ++i) {
+			resultDinv.addNextValue(i, i, constOne / resultLU.getValue(i, i));
+			resultLU.getValue(i, i) = storm::utility::constGetZero<T>();
 		}
+		resultDinv.finalize();
 
-		return new storm::storage::JacobiDecomposition<T>(resultLU, resultDinv);
+		return std::make_pair(std::move(resultLU), std::move(resultDinv));
 	}
 
 	/*!