PARP Research Group

Universidad de Murcia

---

src/qvmath/qvsparseblockmatrix.h Go to the documentation of this file. /* * Copyright (C) 2010, 2011, 2012. PARP Research Group. * <http://perception.inf.um.es> * University of Murcia, Spain. * * This file is part of the QVision library. * * QVision is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as * published by the Free Software Foundation, version 3 of the License. * * QVision is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU Lesser General Public License for more details. * * You should have received a copy of the GNU Lesser General Public * License along with QVision. If not, see <http://www.gnu.org/licenses/>. */ #ifndef QVSPARSEBLOCKMATRIX_H #define QVSPARSEBLOCKMATRIX_H #include <QMap> #include <qvmath/qvmatrix.h> #ifdef MKL_AVAILABLE #include "mkl_dss.h" #endif class QVSparseBlockMatrix: public QMap<int, QMap<int, QVMatrix> > { protected: int majorRows, majorCols, minorRows, minorCols; public: int blockCount; QVSparseBlockMatrix(): QMap<int, QMap<int, QVMatrix> >(), majorRows(0), majorCols(0), minorRows(0), minorCols(0), blockCount(0) { } QVSparseBlockMatrix (const QVSparseBlockMatrix &other): QMap<int, QMap<int, QVMatrix> >(other), majorRows(other.majorRows), majorCols(other.majorCols), minorRows(other.minorRows), minorCols(other.minorCols), blockCount(other.blockCount) { } QVSparseBlockMatrix(const int majorRows, const int majorCols, const int minorRows, const int minorCols): QMap<int, QMap<int, QVMatrix> >(), majorRows(majorRows), majorCols(majorCols), minorRows(minorRows), minorCols(minorCols), blockCount(0) { } QVSparseBlockMatrix(const int majorRows, const int majorCols, const QVMatrix &other): QMap<int, QMap<int, QVMatrix> >(), majorRows(majorRows), majorCols(majorCols), minorRows(other.getRows() / majorRows), minorCols(other.getCols() / majorCols) { if( (minorRows*majorRows != other.getRows()) || (minorCols*majorCols != other.getCols()) ) { std::cout << "[QVSparseBlockMatrix] Error: tried to construct a sparse block matrix of incompatible sizes with original dense matrix." << std::endl << "\tSparse matrix number of blocks:\t" << majorRows << "x" << majorCols << std::endl << "\tDense matrix number of blocks:\t" << other.getRows() << "x" << other.getCols() << std::endl; exit(1); } for(int i = 0; i < majorRows; i++) for(int j = 0; j < majorCols; j++) { QVMatrix M = other.getSubmatrix(i*minorRows, (i+1)*minorRows-1, j*minorCols, (j+1)*minorCols-1); if (M != QVMatrix(minorRows, minorCols, 0.0)) operator[](i)[j] = M; } } inline int getMajorRows() const {return majorRows;}; inline int getMajorCols() const {return majorCols;}; inline int getMinorRows() const {return minorRows;}; inline int getMinorCols() const {return minorCols;}; operator QVMatrix () const { QVMatrix result(majorRows * minorRows, majorCols * minorCols, 0.0); /*for(int i = 0; i < majorRows; i++) for(int j = 0; j < majorCols; j++) if (operator[](i).contains(j)) result.setSubmatrix(i*minorRows, j*minorCols, operator[](i)[j]);*/ foreach(int i, keys()) foreach(int j, operator[](i).keys()) result.setSubmatrix(i*minorRows, j*minorCols, operator[](i)[j]); return result; } QVSparseBlockMatrix transpose() const; void setBlock(const int majorRow, const int majorCol, const QVMatrix &M) { if (majorRow < 0 || majorRow >= majorRows || majorCol < 0 || majorCol >= majorCols) { std::cout << "[QVSparseBlockMatrix::setBlock]: accessing a block outside the sparse matrix." << std::endl << "\tProvided block index:\t" << majorRow << "," << majorCol << std::endl << "\tValid block ranges for the matrix:\t" << majorRows << "x" << majorCols << std::endl; exit(1); } if (minorRows != M.getRows() || minorCols != M.getCols()) { std::cout << "[QVSparseBlockMatrix::setBlock] Error: tried to assing a matrix of incorrect dimensions to a block." << std::endl << "\tMatrix size:\t" << M.getRows() << "x" << M.getCols() << std::endl << "\tSparse matrix block size:\t" << minorRows << "x" << minorCols << std::endl; exit(1); } QMap<int, QVMatrix> & row = operator[](majorRow); if (!row.contains(majorCol)) blockCount ++; operator[](majorRow)[majorCol] = M; } QVMatrix & getBlock(const int majorRow, const int majorCol) { if (majorRow < 0 || majorRow >= majorRows || majorCol < 0 || majorCol >= majorCols) { std::cout << "[QVSparseBlockMatrix::getBlock] Error: accessing a block outside the sparse matrix." << std::endl << "\tProvided block index:\t" << majorRow << "," << majorCol << std::endl << "\tValid block ranges for the matrix:\t" << majorRows << "x" << majorCols << std::endl; exit(1); } QMap<int, QVMatrix> & row = operator[](majorRow); // If no submatrix is at the indexed block, create one filled with NULL values. if (!row.contains(majorCol)) row[majorCol] = QVMatrix(minorRows, minorCols, 0.0); return row[majorCol]; } bool isNullBlock(const int majorRow, const int majorCol) const { if (majorRow < 0 || majorRow >= majorRows || majorCol < 0 || majorCol >= majorCols) { std::cout << "[QVSparseBlockMatrix::isNullBlock] Error: accessing a block outside the sparse matrix." << std::endl << "\tProvided block index:\t" << majorRow << "," << majorCol << std::endl << "\tValid block ranges for the matrix:\t" << majorRows << "x" << majorCols << std::endl; exit(1); } return !operator[](majorRow).contains(majorCol); } /*const QVMatrix & getBlock(const int majorRow, const int majorCol) const { if (majorRow < 0 || majorRow >= majorRows || majorCol < 0 || majorCol >= majorCols) { std::cout << "ERROR: accessing a block outside the sparse matrix." << std::endl << "\tProvided block index:\t" << majorRow << "," << majorCol << std::endl << "\tValid block ranges for the matrix:\t" << majorRows << "x" << majorCols << std::endl; exit(1); } const QMap<int, QVMatrix> & row = operator[](majorRow); // @todo This operation should create a block of adequate dimensions in the sparse matrix if no block is found. if (!row.contains(majorCol)) { std::cout << "ERROR: accessing a block which is NULL in the sparse matrix." << std::endl; exit(1); } return row[majorCol]; }*/ void printBlocks() const { std::cout << "Blocks for sparse matrix:" << std::endl; for(int i = 0; i < majorRows; i++) { std::cout << "\t"; for(int j = 0; j < majorCols; j++) if (operator[](i).contains(j)) std::cout << "[]"; else std::cout << "<>"; std::cout << std::endl; } } QList<int> getBlockRowIndexes(const int majorRow) const { return operator[](majorRow).keys(); } QVSparseBlockMatrix operator*(const QVSparseBlockMatrix &other) const { return dotProduct(other, false, false); }; QVVector operator*(const QVVector &vector) const { return dotProduct(vector, false); } QVSparseBlockMatrix dotProduct(const QVSparseBlockMatrix &other, const bool transposeFirstOperand = false, const bool transposeSecondOperand = false) const; QVVector dotProduct(const QVVector &vector, const bool transposeMatrix = false) const; QVSparseBlockMatrix & operator=(const QVSparseBlockMatrix &other); void clear() { QMap<int, QMap<int, QVMatrix> >::clear(); } // @brief Gets the trace of the matrix. // @note this function works only on sparse matrices with square blocks. #ifndef DOXYGEN_IGNORE_THIS double trace() { // Evaluate the trace. double trace = 0.0; int count = 0; for(int k = 0; k < MIN(getMajorRows(),getMajorCols()); k++) { if (isNullBlock(k,k)) continue; const QVMatrix M = getBlock(k,k); if (M.getCols() < getMinorCols() || M.getRows() < getMinorRows()) { std::cout << "ERROR: block in the diagonal is NULL." << std::endl; exit(1); } for (int j = 0; j < MIN(getMinorRows(), getMinorCols()); j++) { trace += M(j,j); count++; } } return trace; } #endif static QVSparseBlockMatrix randomSquare(const int NB,const int N, const double NZProb,const bool symmetric=true,const bool positive=true); }; #ifndef DOXYGEN_IGNORE_THIS // Internal class: Pardiso sparse matrix format (for solving sparse systems using MKL DSS): class MKLPardisoSparseFormat { public: MKLPardisoSparseFormat(const QVSparseBlockMatrix qvspmatrix, bool isSymmetric); MKLPardisoSparseFormat(): majorRows(0), majorCols(0), minorRows(0), minorCols(0), nRows(0), nCols(0), nNonZeros(0), rowIndex(NULL), columns(NULL), values(NULL) {}; ~MKLPardisoSparseFormat(); int majorRows, majorCols, minorRows, minorCols; int nRows; // Total number of rows. int nCols; // Total number of cols. int nNonZeros; // Total number of non zero elements. #ifdef MKL_AVAILABLE _INTEGER_t *rowIndex; // rowIndex array (see MKL Pardiso sparse format). _INTEGER_t *columns; // columns arrays (see MKL Pardiso sparse format). _DOUBLE_PRECISION_t *values;// values arrays (see MKL Pardiso sparse format). #else int *rowIndex; // rowIndex array (see MKL Pardiso sparse format). int *columns; // columns arrays (see MKL Pardiso sparse format). double *values;// values arrays (see MKL Parddiso sparse format). #endif inline int getMajorRows() const {return majorRows;}; inline int getMajorCols() const {return majorCols;}; inline int getMinorRows() const {return minorRows;}; inline int getMinorCols() const {return minorCols;}; }; #ifdef MKL_AVAILABLE void squareSymmetricSparseMatrixToPardisoFormat(const QVSparseBlockMatrix &qvspmatrix, MKLPardisoSparseFormat &pardiso); #endif #endif #endif

---

QVision framework. PARP research group. Copyright © 2007, 2008, 2009, 2010, 2011.