From 95b000436b1f9035b2b8086b82fb8fab4b93bfa0 Mon Sep 17 00:00:00 2001 From: PBerger Date: Fri, 14 Dec 2012 01:32:21 +0100 Subject: [PATCH] Added a JacobiDecomposition container and conversion function. Added const where possible. --- src/storage/JacobiDecomposition.h | 96 +++++++++++++++++++++++++++++++ src/storage/SquareSparseMatrix.h | 91 +++++++++++++++++++++++++---- src/utility/const_templates.h | 20 ++++++- 3 files changed, 194 insertions(+), 13 deletions(-) create mode 100644 src/storage/JacobiDecomposition.h diff --git a/src/storage/JacobiDecomposition.h b/src/storage/JacobiDecomposition.h new file mode 100644 index 000000000..5ff33b0aa --- /dev/null +++ b/src/storage/JacobiDecomposition.h @@ -0,0 +1,96 @@ +#ifndef MRMC_STORAGE_JACOBIDECOMPOSITION_H_ +#define MRMC_STORAGE_JACOBIDECOMPOSITION_H_ + +#include "boost/integer/integer_mask.hpp" + +#include "log4cplus/logger.h" +#include "log4cplus/loggingmacros.h" + +extern log4cplus::Logger logger; + +namespace mrmc { + +namespace storage { + +/*! + * Forward declaration against Cycle + */ +template +class SquareSparseMatrix; + + +/*! + * A simple container for a sparse Jacobi decomposition + */ +template +class JacobiDecomposition { + +public: + JacobiDecomposition(mrmc::storage::SquareSparseMatrix * const jacobiLuMatrix, mrmc::storage::SquareSparseMatrix * const jacobiDInvMatrix) : this->jacobiLuMatrix(jacobiLuMatrix), this->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 + */ + mrmc::storage::SquareSparseMatrix& 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 + */ + mrmc::storage::SquareSparseMatrix& getJacobiDInvReference() { + return *(this->jacobiDInvMatrix); + } + + /*! + * Accessor for the internal LU Matrix. + * Ownership stays with this class. + * @return A pointer to the Jacobi LU Matrix + */ + mrmc::storage::SquareSparseMatrix* 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 + */ + mrmc::storage::SquareSparseMatrix* getJacobiDInv() { + return this->jacobiDInvMatrix; + } + +private: + + /*! + * The copy constructor is disabled for this class. + */ + //JacobiDecomposition(const JacobiDecomposition& that) = delete; // not possible in VS2012 + JacobiDecomposition(const JacobiDecomposition& that) {} + + /*! + * Pointer to the LU Matrix + */ + mrmc::storage::SquareSparseMatrix *jacobiLuMatrix; + + /*! + * Pointer to the D^{-1} Matrix + */ + mrmc::storage::SquareSparseMatrix *jacobiDInvMatrix; +}; + +} // namespace storage + +} // namespace mrmc + +#endif // MRMC_STORAGE_JACOBIDECOMPOSITION_H_ diff --git a/src/storage/SquareSparseMatrix.h b/src/storage/SquareSparseMatrix.h index d57d9ddb6..4fe87daed 100644 --- a/src/storage/SquareSparseMatrix.h +++ b/src/storage/SquareSparseMatrix.h @@ -1,5 +1,5 @@ -#ifndef MRMC_SPARSE_STATIC_SPARSE_MATRIX_H_ -#define MRMC_SPARSE_STATIC_SPARSE_MATRIX_H_ +#ifndef MRMC_STORAGE_SQUARESPARSEMATRIX_H_ +#define MRMC_STORAGE_SQUARESPARSEMATRIX_H_ #include #include @@ -12,6 +12,7 @@ #include "src/exceptions/out_of_range.h" #include "src/exceptions/file_IO_exception.h" #include "src/storage/BitVector.h" +#include "src/storage/JacobiDecomposition.h" #include "src/utility/const_templates.h" #include "Eigen/Sparse" @@ -340,7 +341,7 @@ public: * @return True iff the value is set in the matrix, false otherwise. * On false, 0 will be written to *target. */ - inline bool getValue(uint_fast64_t row, uint_fast64_t col, T* const target) { + 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)) { LOG4CPLUS_ERROR(logger, "Trying to read a value from illegal position (" << row << ", " << col << ")."); @@ -379,6 +380,52 @@ public: return false; } + /*! + * Gets the matrix element at the given row and column to the given value. + * NOTE: This function does not check the internal status for errors for performance reasons. + * WARNING: It is possible to modify through this function. Usage only valid + * for elements EXISTING in the sparse matrix! If the requested value does not exist, + * an exception will be thrown. + * @param row The row in which the element is to be read. + * @param col The column in which the element is to be read. + * + * @return A reference to the value + */ + inline T& getValue(uint_fast64_t row, uint_fast64_t col) { + // Check for illegal access indices. + if ((row > rowCount) || (col > rowCount)) { + LOG4CPLUS_ERROR(logger, "Trying to read a value from illegal position (" << row << ", " << col << ")."); + throw mrmc::exceptions::out_of_range("Trying to read a value from illegal position."); + return false; + } + + // 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 + // over the accessed row to find the element. + uint_fast64_t rowStart = rowIndications[row]; + uint_fast64_t rowEnd = rowIndications[row + 1]; + while (rowStart < rowEnd) { + // If the lement is found, write the content to the specified + // position and return true. + if (columnIndications[rowStart] == col) { + return valueStorage[rowStart]; + } + // If the column of the current element is already larger than the + // requested column, the requested element cannot be contained + // in the matrix and we may therefore stop searching. + if (columnIndications[rowStart] > col) { + break; + } + ++rowStart; + } + + throw mrmc::exceptions::invalid_argument("Trying to get a reference to a non-existant value."); + } + /*! * Returns the number of rows of the matrix. */ @@ -713,7 +760,7 @@ public: * @return The sum of the elements in the given row whose column bits * are set to one on the given constraint. */ - T getConstrainedRowSum(const uint_fast64_t row, const mrmc::storage::BitVector& constraint) { + T getConstrainedRowSum(const uint_fast64_t row, const mrmc::storage::BitVector& constraint) const { T result(0); for (uint_fast64_t i = rowIndications[row]; i < rowIndications[row + 1]; ++i) { if (constraint.get(columnIndications[i])) { @@ -731,7 +778,7 @@ public: * @param resultVector A pointer to the resulting vector that has at least * as many elements as there are bits set to true in the constraint. */ - void getConstrainedRowCountVector(const mrmc::storage::BitVector& rowConstraint, const mrmc::storage::BitVector& columnConstraint, std::vector* resultVector) { + void getConstrainedRowCountVector(const mrmc::storage::BitVector& rowConstraint, const mrmc::storage::BitVector& columnConstraint, std::vector* resultVector) const { uint_fast64_t currentRowCount = 0; for (auto row : rowConstraint) { (*resultVector)[currentRowCount++] = getConstrainedRowSum(row, columnConstraint); @@ -744,7 +791,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(mrmc::storage::BitVector& constraint) { + SquareSparseMatrix* getSubmatrix(mrmc::storage::BitVector& constraint) const { LOG4CPLUS_DEBUG(logger, "Creating a sub-matrix with " << constraint.getNumberOfSetBits() << " rows."); // Check for valid constraint. @@ -828,11 +875,31 @@ public: } } + /*! + * Calculates the Jacobi-Decomposition of this sparse matrix. + * @return A pointer to a class containing the matrix L+U and the inverted diagonal matrix D^-1 + */ + mrmc::storage::JacobiDecomposition* getJacobiDecomposition() const { + uint_fast64_t rowCount = this->getRowCount(); + SquareSparseMatrix *resultLU = new SquareSparseMatrix(this); + SquareSparseMatrix *resultDinv = new SquareSparseMatrix(rowCount); + // no entries apart from those on the diagonal + resultDinv->initialize(0); + // copy diagonal entries to other matrix + T dummy; + for (int i = 0; i < rowCount; ++i) { + resultDinv->addNextValue(i, i, resultLU->getValue(i, i)); + resultLU->getValue(i, i) = mrmc::utility::constGetZero(&dummy); + } + + return new mrmc::storage::JacobiDecomposition(resultLU, resultDinv); + } + /*! * Returns the size of the matrix in memory measured in bytes. * @return The size of the matrix in memory measured in bytes. */ - uint_fast64_t getSizeInMemory() { + uint_fast64_t getSizeInMemory() const { uint_fast64_t size = sizeof(*this); // Add value_storage size. size += sizeof(T) * nonZeroEntryCount; @@ -865,7 +932,7 @@ public: return this->columnIndications + this->rowIndications[row + 1]; } - void print() { + 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; @@ -999,7 +1066,7 @@ private: * the given Eigen matrix. */ template - _Index getEigenSparseMatrixCorrectNonZeroEntryCount(const Eigen::SparseMatrix<_Scalar, _Options, _Index>& eigen_sparse_matrix) { + _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(); @@ -1040,8 +1107,10 @@ private: }; -} // namespace sparse + + +} // namespace storage } // namespace mrmc -#endif // MRMC_SPARSE_STATIC_SPARSE_MATRIX_H_ +#endif // MRMC_STORAGE_SQUARESPARSEMATRIX_H_ diff --git a/src/utility/const_templates.h b/src/utility/const_templates.h index 0261cdecc..ea2231cba 100644 --- a/src/utility/const_templates.h +++ b/src/utility/const_templates.h @@ -35,8 +35,15 @@ inline int_fast32_t constGetZero(int_fast32_t&) { return 0; } -/*! @internal - * Specialization of constGetZero for int_fast32_t +/*! @cond TEMPLATE_SPECIALIZATION + * (exclude the specializations from the documentation) */ +template <> +inline uint_fast64_t constGetZero(uint_fast64_t&) { + return 0; +} + +/*! @cond TEMPLATE_SPECIALIZATION + * Specialization of constGetZero for double */ template <> inline double constGetZero(double&) { @@ -66,6 +73,15 @@ inline int_fast32_t constGetOne(int_fast32_t&) { return 1; } +/*! @cond TEMPLATE_SPECIALIZATION + * (exclude the specializations from the documentation) */ +template<> +inline uint_fast64_t constGetOne(uint_fast64_t&) { + return 1; +} + +/*! @cond TEMPLATE_SPECIALIZATION + * (exclude the specializations from the documentation) */ template<> inline double constGetOne(double&) { return 1.0;