From b2c0cfc57cb3e03ab599694f4f2d32d68c726b78 Mon Sep 17 00:00:00 2001 From: PBerger Date: Tue, 26 Feb 2013 04:16:44 +0100 Subject: [PATCH] Added a conversion routine GmmXX -> Storm Sparse Matrix Added Jacobi to possible LE Solvers in the GMM Model Checker --- src/adapters/GmmxxAdapter.h | 42 ++++++++++- src/modelchecker/GmmxxDtmcPrctlModelChecker.h | 73 ++++++++++++++----- src/storage/JacobiDecomposition.h | 5 ++ src/storage/SparseMatrix.h | 13 +++- 4 files changed, 110 insertions(+), 23 deletions(-) diff --git a/src/adapters/GmmxxAdapter.h b/src/adapters/GmmxxAdapter.h index 5c5faa670..74f0e6254 100644 --- a/src/adapters/GmmxxAdapter.h +++ b/src/adapters/GmmxxAdapter.h @@ -22,7 +22,7 @@ namespace adapters { class GmmxxAdapter { public: /*! - * Converts a sparse matrix into the sparse matrix in the gmm++ format. + * Converts a sparse matrix into a sparse matrix in the gmm++ format. * @return A pointer to a row-major sparse matrix in gmm++ format. */ template @@ -46,6 +46,46 @@ public: return result; } + + /*! + * Converts a sparse matrix in the gmm++ format to Storm Sparse Matrix format. + * @return A pointer to a row-major sparse matrix in our format. + */ + template + static storm::storage::SparseMatrix* fromGmmxxSparseMatrix(gmm::csr_matrix const& matrix) { + uint_fast64_t realNonZeros = gmm::nnz(matrix); + LOG4CPLUS_DEBUG(logger, "Converting matrix with " << realNonZeros << " non-zeros from gmm++ format into Storm."); + + // Prepare the resulting matrix. + storm::storage::SparseMatrix* result = new storm::storage::SparseMatrix(matrix.nrows(), matrix.ncols()); + + // Set internal NonZero Counter + result->nonZeroEntryCount = realNonZeros; + result->setState(result->Initialized); + + if (!result->prepareInternalStorage(false)) { + LOG4CPLUS_ERROR(logger, "Unable to allocate internal storage while converting GMM++ Matrix to Storm."); + delete result; + return nullptr; + } else { + + // Copy Row Indications + std::copy(matrix.jc.begin(), matrix.jc.end(), std::back_inserter(result->rowIndications)); + // Copy Columns Indications + std::copy(matrix.ir.begin(), matrix.ir.end(), std::back_inserter(result->columnIndications)); + // And do the same thing with the actual values. + std::copy(matrix.pr.begin(), matrix.pr.end(), std::back_inserter(result->valueStorage)); + + result->currentSize = realNonZeros; + result->lastRow = matrix.nrows() - 1; + } + + result->finalize(); + + LOG4CPLUS_DEBUG(logger, "Done converting matrix to storm format."); + + return result; + } }; } //namespace adapters diff --git a/src/modelchecker/GmmxxDtmcPrctlModelChecker.h b/src/modelchecker/GmmxxDtmcPrctlModelChecker.h index d4aeef764..a5a3b567d 100644 --- a/src/modelchecker/GmmxxDtmcPrctlModelChecker.h +++ b/src/modelchecker/GmmxxDtmcPrctlModelChecker.h @@ -353,7 +353,7 @@ public: * Registers all options associated with the gmm++ matrix library. */ static void putOptions(boost::program_options::options_description* desc) { - desc->add_options()("lemethod", boost::program_options::value()->default_value("bicgstab")->notifier(&validateLeMethod), "Sets the method used for linear equation solving. Must be in {bicgstab, qmr}."); + desc->add_options()("lemethod", boost::program_options::value()->default_value("bicgstab")->notifier(&validateLeMethod), "Sets the method used for linear equation solving. Must be in {bicgstab, qmr, jacobi}."); desc->add_options()("maxiter", boost::program_options::value()->default_value(10000), "Sets the maximal number of iterations for iterative equation solving."); desc->add_options()("precision", boost::program_options::value()->default_value(1e-6), "Sets the precision for iterative equation solving."); desc->add_options()("precond", boost::program_options::value()->default_value("ilu")->notifier(&validatePreconditioner), "Sets the preconditioning technique for linear equation solving. Must be in {ilu, diagonal, ildlt, none}."); @@ -365,7 +365,7 @@ public: * Throws an exception of type InvalidSettings in case the selected method is illegal. */ static void validateLeMethod(const std::string& lemethod) { - if ((lemethod != "bicgstab") && (lemethod != "qmr")) { + if ((lemethod != "bicgstab") && (lemethod != "qmr") && (lemethod != "jacobi")) { throw exceptions::InvalidSettingsException() << "Argument " << lemethod << " for option 'lemethod' is invalid."; } } @@ -423,6 +423,13 @@ private: } else if (precond == "none") { gmm::bicgstab(A, x, b, gmm::identity_matrix(), iter); } + + // Check if the solver converged and issue a warning otherwise. + if (iter.converged()) { + LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iter.get_iteration() << " iterations."); + } else { + LOG4CPLUS_WARN(logger, "Iterative solver did not converge."); + } // FIXME: gmres has been disabled, because it triggers gmm++ compilation errors /* } else if (s->getString("lemethod").compare("gmres") == 0) { LOG4CPLUS_INFO(logger, "Using GMRES method."); @@ -446,13 +453,16 @@ private: } else if (precond == "none") { gmm::qmr(A, x, b, gmm::identity_matrix(), iter); } - } - // Check if the solver converged and issue a warning otherwise. - if (iter.converged()) { - LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iter.get_iteration() << " iterations."); - } else { - LOG4CPLUS_WARN(logger, "Iterative solver did not converge."); + // Check if the solver converged and issue a warning otherwise. + if (iter.converged()) { + LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iter.get_iteration() << " iterations."); + } else { + LOG4CPLUS_WARN(logger, "Iterative solver did not converge."); + } + } else if (s->getString("lemethod") == "jacobi") { + LOG4CPLUS_INFO(logger, "Using Jacobi method."); + solveLinearEquationSystemWithJacobi(A, x, b); } } @@ -467,7 +477,7 @@ private: * @return The solution of the system of linear equations in form of the elements of the vector * x. */ - void solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b) const { + void solveLinearEquationSystemWithJacobi(gmm::csr_matrix const& A, std::vector& x, std::vector const& b) const { // Get the settings object to customize linear solving. storm::settings::Settings* s = storm::settings::instance(); @@ -475,41 +485,68 @@ private: if (precision <= 0) { LOG4CPLUS_ERROR(logger, "Selected precision for linear equation solving must be strictly greater than zero for Jacobi method."); } - + + // Convert the Source Matrix to Storm format for Decomposition + storm::storage::SparseMatrix* stormFormatA = storm::adapters::GmmxxAdapter::fromGmmxxSparseMatrix(A); + // Get a Jacobi Decomposition of the Input Matrix A - storm::storage::JacobiDecomposition* jacobiDecomposition = A.getJacobiDecomposition(); + storm::storage::JacobiDecomposition* jacobiDecomposition = stormFormatA->getJacobiDecomposition(); + + // The Storm Version is not needed after the decomposition step + delete stormFormatA; // Convert the Diagonal matrix to GMM format - gmm::csr_matrix* gmmxxDiagonalInverted = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(jacobiDecomposition->getJacobiDInv()); + gmm::csr_matrix* gmmxxDiagonalInverted = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(jacobiDecomposition->getJacobiDInvReference()); // Convert the LU Matrix to GMM format - gmm::csr_matrix* gmmxxLU = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(jacobiDecomposition->getJacobiLU()); + gmm::csr_matrix* gmmxxLU = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(jacobiDecomposition->getJacobiLUReference()); - delete jacobiDecomposition; LOG4CPLUS_INFO(logger, "Starting iterative Jacobi Solver."); // x_(k + 1) = D^-1 * (b - R * x_k) + // In x we keep a copy of the result for swapping in the loop (e.g. less copy-back) std::vector* xNext = new std::vector(x.size()); const std::vector* xCopy = xNext; std::vector* xCurrent = &x; + // Target vector for precision calculation + std::vector* residuum = new std::vector(x.size()); + uint_fast64_t iterationCount = 0; do { - // R * x_k -> xCurrent + // R * x_k (xCurrent is x_k) -> xNext gmm::mult(*gmmxxLU, *xCurrent, *xNext); - // b - R * x_k + // b - R * x_k (xNext contains R * x_k) -> xNext gmm::add(b, gmm::scaled(*xNext, -1.0), *xNext); - // D^-1 * (b - R * x_k) + // D^-1 * (b - R * x_k) -> xNext gmm::mult(*gmmxxDiagonalInverted, *xNext, *xNext); + // swap xNext with xCurrent so that the next iteration can use xCurrent again without having to copy the vector std::vector* swap = xNext; xNext = xCurrent; xCurrent = swap; ++iterationCount; - } while (gmm::vect_norminf(*xCurrent) > precision); + // Precision calculation via ||A * x_k - b|| < precision + gmm::mult(A, *xCurrent, *residuum); + gmm::add(gmm::scaled(*residuum, -1.0), b, *residuum); + } while (gmm::vect_norminf(*residuum) > precision); + + // If the last iteration did not write to the original x + // we have to swith them + if (xCurrent == xCopy) { + x.swap(*xCurrent); + } + // xCopy always points to the Swap-Copy of x we created delete xCopy; + // Delete the residuum vector + delete residuum; + // Delete the decompositions + delete jacobiDecomposition; + // and the GMM Matrices + delete gmmxxDiagonalInverted; + delete gmmxxLU; LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iterationCount << " iterations."); } diff --git a/src/storage/JacobiDecomposition.h b/src/storage/JacobiDecomposition.h index 9aec6f061..ef4f8bb99 100644 --- a/src/storage/JacobiDecomposition.h +++ b/src/storage/JacobiDecomposition.h @@ -28,6 +28,11 @@ template 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 * const jacobiLuMatrix, storm::storage::SparseMatrix * const jacobiDInvMatrix) : jacobiLuMatrix(jacobiLuMatrix), jacobiDInvMatrix(jacobiDInvMatrix) { } diff --git a/src/storage/SparseMatrix.h b/src/storage/SparseMatrix.h index 970cbcd07..503d7ccc8 100644 --- a/src/storage/SparseMatrix.h +++ b/src/storage/SparseMatrix.h @@ -932,14 +932,19 @@ 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 */ storm::storage::JacobiDecomposition* getJacobiDecomposition() const { uint_fast64_t rowCount = this->getRowCount(); - SparseMatrix *resultLU = new SparseMatrix(this); - SparseMatrix *resultDinv = new SparseMatrix(rowCount); - // no entries apart from those on the diagonal - resultDinv->initialize(0); + uint_fast64_t colCount = this->getColumnCount(); + if (rowCount != colCount) { + throw storm::exceptions::InvalidArgumentException() << "SparseMatrix::getJacobiDecomposition requires the Matrix to be square!"; + } + storm::storage::SparseMatrix *resultLU = new storm::storage::SparseMatrix(*this); + storm::storage::SparseMatrix *resultDinv = new storm::storage::SparseMatrix(rowCount, colCount); + // no entries apart from those on the diagonal (rowCount) + resultDinv->initialize(rowCount); // constant 1 for diagonal inversion T constOne = storm::utility::constGetOne();