system_matrix.cpp

/*---------------------------------------------------------------------------*\
 *
 *  bitpit
 *
 *  Copyright (C) 2015-2021 OPTIMAD engineering Srl
 *
 *  -------------------------------------------------------------------------
 *  License
 *  This file is part of bitpit.
 *
 *  bitpit is free software: you can redistribute it and/or modify it
 *  under the terms of the GNU Lesser General Public License v3 (LGPL)
 *  as published by the Free Software Foundation.
 *
 *  bitpit 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 bitpit. If not, see <http://www.gnu.org/licenses/>.
 *
\*---------------------------------------------------------------------------*/
 
#include "system_matrix.hpp"
 
namespace bitpit {

#if BITPIT_ENABLE_MPI==1
#endif
 
#if BITPIT_ENABLE_MPI==1
SparseMatrix::SparseMatrix()
    : SparseMatrix(MPI_COMM_SELF)
{
}

SparseMatrix::SparseMatrix(MPI_Comm communicator)
    : m_blockSize(0), m_nRows(0), m_nCols(0), m_nNZ(0), m_maxRowNZ(0), m_lastRow(-1),
      m_assembled(false),
      m_partitioned(false),
      m_global_nRows(0), m_global_nCols(0), m_global_nNZ(0),
      m_global_maxRowNZ(0), m_global_rowOffset(0), m_global_colOffset(0)
{
    // Set the communicator
    setCommunicator(communicator);
}
#else
SparseMatrix::SparseMatrix()
    : m_blockSize(0), m_nRows(0), m_nCols(0), m_nNZ(0), m_maxRowNZ(0), m_lastRow(-1),
      m_assembled(false)
{
}
#endif
 
#if BITPIT_ENABLE_MPI==1
SparseMatrix::SparseMatrix(long nRows, long nCols, long nNZ)
    : SparseMatrix(MPI_COMM_SELF, false, 1, nRows, nCols, nNZ)
{
}

SparseMatrix::SparseMatrix(int blockSize, long nRows, long nCols, long nNZ)
    : SparseMatrix(MPI_COMM_SELF, false, blockSize, nRows, nCols, nNZ)
{
}

SparseMatrix::SparseMatrix(MPI_Comm communicator, bool partitioned, long nRows, long nCols, long nNZ)
    : SparseMatrix(communicator)
{
    _initialize(partitioned, 1, nRows, nCols, nNZ);
}

SparseMatrix::SparseMatrix(MPI_Comm communicator, bool partitioned, int blockSize, long nRows, long nCols, long nNZ)
    : SparseMatrix(communicator)
{
    _initialize(partitioned, blockSize, nRows, nCols, nNZ);
}
#else
SparseMatrix::SparseMatrix(long nRows, long nCols, long nNZ)
    : SparseMatrix(1, nRows, nCols, nNZ)
{
}

SparseMatrix::SparseMatrix(int blockSize, long nRows, long nCols, long nNZ)
    : SparseMatrix()
{
    _initialize(blockSize, nRows, nCols, nNZ);
}
#endif
 
#if BITPIT_ENABLE_MPI==1
SparseMatrix::SparseMatrix(const SparseMatrix &other)
    : m_blockSize(other.m_blockSize),
      m_nRows(other.m_nRows),
      m_nCols(other.m_nCols),
      m_nNZ(other.m_nNZ),
      m_maxRowNZ(other.m_maxRowNZ),
      m_lastRow(other.m_lastRow),
      m_assembled(other.m_assembled),
      m_partitioned(other.m_partitioned),
      m_global_nRows(other.m_global_nRows),
      m_global_nCols(other.m_global_nCols),
      m_global_nNZ(other.m_global_nNZ),
      m_global_maxRowNZ(other.m_global_maxRowNZ),
      m_global_rowOffset(other.m_global_rowOffset),
      m_global_colOffset(other.m_global_colOffset),
      m_pattern(other.m_pattern),
      m_values (other.m_values)
{
    // Set the communicator
    setCommunicator(other.m_communicator);
}
#endif

SparseMatrix::~SparseMatrix()
{
#if BITPIT_ENABLE_MPI==1
    // Free the MPI communicator
    freeCommunicator();
#endif
}
 
#if BITPIT_ENABLE_MPI==1
void SparseMatrix::initialize(bool partitioned, long nRows, long nCols, long nNZ)
{
    _initialize(partitioned, 1, nRows, nCols, nNZ);
}

void SparseMatrix::initialize(bool partitioned, int blockSize, long nRows, long nCols, long nNZ)
{
    _initialize(partitioned, blockSize, nRows, nCols, nNZ);
}
#endif

void SparseMatrix::initialize(long nRows, long nCols, long nNZ)
{
    initialize(1, nRows, nCols, nNZ);
}

void SparseMatrix::initialize(int blockSize, long nRows, long nCols, long nNZ)
{
#if BITPIT_ENABLE_MPI==1
    _initialize(false, blockSize, nRows, nCols, nNZ);
#else
    _initialize(blockSize, nRows, nCols, nNZ);
#endif
}
 
#if BITPIT_ENABLE_MPI==1
void SparseMatrix::_initialize(bool partitioned, int blockSize, long nRows, long nCols, long nNZ)
{
    // Serial initialization
    _initialize(blockSize, nRows, nCols, nNZ);
 
    // Parallel initialization
    m_partitioned = partitioned;
 
    m_global_nRows = m_nRows;
    m_global_nCols = m_nCols;
    if (m_partitioned) {
        MPI_Allreduce(MPI_IN_PLACE, &m_global_nRows, 1, MPI_LONG, MPI_SUM, m_communicator);
        MPI_Allreduce(MPI_IN_PLACE, &m_global_nCols, 1, MPI_LONG, MPI_SUM, m_communicator);
    }
 
    m_global_rowOffset = 0;
    m_global_colOffset = 0;
    if (m_partitioned) {
        int nProcessors;
        MPI_Comm_size(m_communicator, &nProcessors);
 
        std::vector<long> nGlobalRows(nProcessors);
        MPI_Allgather(&nRows, 1, MPI_LONG, nGlobalRows.data(), 1, MPI_LONG, m_communicator);
 
        std::vector<long> nGlobalCols(nProcessors);
        MPI_Allgather(&nCols, 1, MPI_LONG, nGlobalCols.data(), 1, MPI_LONG, m_communicator);
 
        int rank;
        MPI_Comm_rank(m_communicator, &rank);
        for (int i = 0; i < rank; ++i) {
            m_global_rowOffset += nGlobalRows[i];
            m_global_colOffset += nGlobalCols[i];
        }
    }
}
#endif

void SparseMatrix::_initialize(int blockSize, long nRows, long nCols, long nNZ)
{
    assert(blockSize >= 0);
    assert(nRows >= 0);
    assert(nCols >= 0);
    assert(nNZ >= 0);
 
    clear();
 
    m_blockSize = blockSize;
    m_nRows     = nRows;
    m_nCols     = nCols;
 
    initializePatternStorage(nNZ);
    initializeValueStorage(nNZ);
}

void SparseMatrix::clear(bool release)
{
    clearPatternStorage(release);
    clearValueStorage(release);
 
    m_blockSize =  0;
    m_nRows     =  0;
    m_nCols     =  0;
    m_nNZ       =  0;
    m_maxRowNZ  =  0;
    m_lastRow   = -1;
 
#if BITPIT_ENABLE_MPI==1
    m_partitioned = false;
 
    m_global_nRows     = 0;
    m_global_nCols     = 0;
    m_global_nNZ       = 0;
    m_global_maxRowNZ  = 0;
    m_global_rowOffset = 0;
    m_global_colOffset = 0;
#endif
 
    m_assembled = false;
}

void SparseMatrix::squeeze()
{
    squeezePatternStorage();
    squeezeValueStorage();
}

void SparseMatrix::display(std::ostream &stream, double negligiblity, int indent) const
{
    // Initialize padding
    std::string padding(indent, ' ');
 
    // Display base information
    int blockSize = getBlockSize();
    int nBlockElements = blockSize * blockSize;
 
    stream << padding << "General information" << std::endl;
    stream << padding << "  Block size ............................................ " << blockSize << std::endl;
 
    // Gather local information
    long nRows = getRowCount();
    long nCols = getColCount();
 
    long nRowElements = getRowElementCount();
    long nColElements = getColElementCount();
 
    long nNZBlocks   = getNZCount();
    long nNZElements = getNZElementCount();
 
    long maxRowNZBlocks   = getMaxRowNZCount();
    long maxRowNZElements = getMaxRowNZElementCount();
 
    // Count local non-negligible blocks and elements
    int nZeroStreak       = 0;
    long nNZBlockValues   = 0;
    long nNZElementValues = 0;
    for (long k = 0; k < nNZElements; ++k) {
        double value = m_values[k];
        if (!bitpit::utils::DoubleFloatingEqual()(value, 0., negligiblity, negligiblity)) {
            ++nNZElementValues;
            nZeroStreak = 0;
        } else {
            ++nZeroStreak;
        }
 
        if ((k % nBlockElements) == 0) {
            if (nZeroStreak < nBlockElements) {
                ++nNZBlockValues;
            }
            nZeroStreak = 0;
        }
    }
 
    // Evaluate local sparisty
    std::size_t nBlocks   = nRows * nCols;
    std::size_t nElements = nRowElements * nColElements;
 
    double blockSparsity   = static_cast<double>(nBlocks - nNZBlocks) / nBlocks;
    double elementSparsity = static_cast<double>(nElements - nNZElements) / nElements;
 
    double blockValueSparsity   = static_cast<double>(nBlocks - nNZBlockValues) / nBlocks;
    double elementValueSparsity = static_cast<double>(nElements - nNZElementValues) / nElements;
 
    // Evaluate local memory usage
    double valueStorageMemory = static_cast<double>(nNZElements * sizeof(decltype(m_values)::value_type));
 
    // Display local information
    stream << padding << "Local information: " << std::endl;
    stream << padding << "  Maximum number of non-zero blocks per row ............. " << maxRowNZBlocks << std::endl;
    stream << padding << "  Maximum number of non-zero elements per row ........... " << maxRowNZElements << std::endl;
    stream << padding << "  Number of block rows .................................. " << nRows << std::endl;
    stream << padding << "  Number of block columns ............................... " << nCols << std::endl;
    stream << padding << "  Number of non-zero blocks (pattern) ................... " << nNZBlocks << std::endl;
    stream << padding << "  Number of non-zero elements (pattern) ................. " << nNZElements << std::endl;
    stream << padding << "  Number of non-zero blocks (non-neglibile values) ...... " << nNZBlockValues << std::endl;
    stream << padding << "  Number of non-zero elements (non-neglibile values) .... " << nNZElementValues << std::endl;
    stream << padding << "  Sparsity of the blocks (pattern) ...................... " << blockSparsity << std::endl;
    stream << padding << "  Sparsity of the elements (pattern) .................... " << elementSparsity << std::endl;
    stream << padding << "  Sparsity of the blocks (non-neglibile values) ......... " << blockValueSparsity << std::endl;
    stream << padding << "  Sparsity of the elements (non-neglibile values) ....... " << elementValueSparsity << std::endl;
 
    stream << padding << "  Memory used by the value storage ...................... ";
    if (valueStorageMemory > 1024 * 1024 * 1024) {
        stream  << valueStorageMemory / (1024 * 1024 * 1024) << " GB";
    } else if (valueStorageMemory > 1024 * 1024) {
        stream  << valueStorageMemory / (1024 * 1024) << " MB";
    } else if (valueStorageMemory > 1024) {
        stream  << valueStorageMemory / (1024) << " KB";
    }
    stream << std::endl;
 
#if BITPIT_ENABLE_MPI==1
    // Display global information
    if (isPartitioned()) {
        // Get MPI data types
        MPI_Datatype valueMPIDataType;
        if (std::is_same<decltype(m_values)::value_type, double>::value) {
            valueMPIDataType = MPI_DOUBLE;
        } else if (std::is_same<decltype(m_values)::value_type, float>::value) {
            valueMPIDataType = MPI_FLOAT;
        } else {
            throw std::runtime_error("Unable to identify the MPI data type of the matrix values.");
        }
 
        // Gather global information
        long nGlobalRows = getRowGlobalCount();
        long nGlobalCols = getColGlobalCount();
 
        long nGlobalRowElements = getRowGlobalElementCount();
        long nGlobalColElements = getColGlobalElementCount();
 
        long nGlobalNZBlocks   = getNZGlobalCount();
        long nGlobalNZElements = getNZGlobalElementCount();
 
        long maxGlobalRowNZBlocks;
        long maxGlobalRowNZElements;
        MPI_Allreduce(&maxRowNZBlocks,   &maxGlobalRowNZBlocks,   1, MPI_LONG, MPI_MAX, getCommunicator());
        MPI_Allreduce(&maxRowNZElements, &maxGlobalRowNZElements, 1, MPI_LONG, MPI_MAX, getCommunicator());
 
        // Count global non-negligible blocks and elements
        long nGlobalNZBlockValues;
        long nGlobalNZElementValues;
        MPI_Allreduce(&nNZBlockValues,   &nGlobalNZBlockValues,   1, MPI_LONG, MPI_SUM, getCommunicator());
        MPI_Allreduce(&nNZElementValues, &nGlobalNZElementValues, 1, MPI_LONG, MPI_SUM, getCommunicator());
 
        // Evaluate global sparisty
        std::size_t nGlobalBlocks   = nGlobalRows * nGlobalCols;
        std::size_t nGlobalElements = nGlobalRowElements * nGlobalColElements;
 
        double globalBlockSparsity   = static_cast<double>(nGlobalBlocks - nGlobalNZBlocks) / nGlobalBlocks;
        double globalElementSparsity = static_cast<double>(nGlobalElements - nGlobalNZElements) / nGlobalElements;
 
        double globalBlockValueSparsity   = static_cast<double>(nGlobalBlocks - nGlobalNZBlockValues) / nGlobalBlocks;
        double globalElementValueSparsity = static_cast<double>(nGlobalElements - nGlobalNZElementValues) / nGlobalElements;
 
        // Evaluate global memory usage
        double valueStorageGlobalMemory;
        MPI_Allreduce(&valueStorageMemory, &valueStorageGlobalMemory, 1, valueMPIDataType, MPI_SUM, getCommunicator());
 
        // Display information
        stream << padding << "Global information: " << std::endl;
        stream << padding << "  Maximum number of non-zero blocks per row ............. " << maxGlobalRowNZBlocks << std::endl;
        stream << padding << "  Maximum number of non-zero elements per row ........... " << maxGlobalRowNZElements << std::endl;
        stream << padding << "  Number of block columns ............................... " << nGlobalCols << std::endl;
        stream << padding << "  Number of block rows .................................. " << nGlobalRows << std::endl;
        stream << padding << "  Number of non-zero blocks (pattern) ................... " << nGlobalNZBlocks << std::endl;
        stream << padding << "  Number of non-zero elements (pattern) ................. " << nGlobalNZElements << std::endl;
        stream << padding << "  Number of non-zero blocks (non-neglibile values) ...... " << nGlobalNZBlockValues << std::endl;
        stream << padding << "  Number of non-zero elements (non-neglibile values) .... " << nGlobalNZElementValues << std::endl;
        stream << padding << "  Sparsity of the blocks (pattern) ...................... " << globalBlockSparsity << std::endl;
        stream << padding << "  Sparsity of the elements (pattern) .................... " << globalElementSparsity << std::endl;
        stream << padding << "  Sparsity of the blocks (non-neglibile values) ......... " << globalBlockValueSparsity << std::endl;
        stream << padding << "  Sparsity of the elements (non-neglibile values) ....... " << globalElementValueSparsity << std::endl;
 
        stream << padding << "  Memory used by the value storage ...................... ";
        if (valueStorageGlobalMemory > 1024 * 1024 * 1024) {
            stream  << valueStorageGlobalMemory / (1024 * 1024 * 1024) << " GB";
        } else if (valueStorageGlobalMemory > 1024 * 1024) {
            stream  << valueStorageGlobalMemory / (1024 * 1024) << " MB";
        } else if (valueStorageGlobalMemory > 1024) {
            stream  << valueStorageGlobalMemory / (1024) << " KB";
        }
        stream << std::endl;
    }
#endif
 
    // Display information
    stream << padding << "Display information: " << std::endl;
    stream << padding << "  Negligibility threshold ............................... " << negligiblity << std::endl;
}

void SparseMatrix::initializePatternStorage()
{
    long nNZ = getNZCount();
 
    m_pattern.reserve(nNZ);
}

void SparseMatrix::initializePatternStorage(long nNZ)
{
    long nRows = getRowCount();
 
    m_pattern.reserve(nRows, nNZ);
}

void SparseMatrix::squeezePatternStorage()
{
    m_pattern.shrinkToFit();
}

void SparseMatrix::clearPatternStorage(bool release)
{
    m_pattern.clear(release);
}

void SparseMatrix::initializeValueStorage()
{
    long nNZ = getNZCount();
 
    initializeValueStorage(nNZ);
}

void SparseMatrix::initializeValueStorage(long nNZ)
{
    long nNZElements = getNZElementCount(nNZ);
 
    m_values.reserve(nNZElements);
}

void SparseMatrix::squeezeValueStorage()
{
    m_values.shrink_to_fit();
}

void SparseMatrix::clearValueStorage(bool release)
{
    m_values.clear();
 
    if (release) {
        squeezeValueStorage();
    }
}

#if BITPIT_ENABLE_MPI==1
#endif
void SparseMatrix::assembly()
{
    // Early return if the matrix is already assembled.
    if (isAssembled()) {
        return;
    }
 
    // Assembly can be called only after adding all the rows
    if (countMissingRows() != 0) {
        throw std::runtime_error("Assembly can be called only after adding all the rows.");
    }
 
#if BITPIT_ENABLE_MPI==1
    // Update global information of the non-zero elements
    if (m_partitioned) {
        MPI_Allreduce(&m_maxRowNZ, &m_global_maxRowNZ, 1, MPI_LONG, MPI_MAX, m_communicator);
        MPI_Allreduce(&m_nNZ, &m_global_nNZ, 1, MPI_LONG, MPI_SUM, m_communicator);
    } else {
        m_global_maxRowNZ = m_maxRowNZ;
        m_global_nNZ      = m_nNZ;
    }
#endif
 
    // The function is now assembled
    m_assembled = true;
}

bool SparseMatrix::isAssembled() const
{
    return m_assembled;
}

long SparseMatrix::countMissingRows() const
{
    long nRows      = getRowCount();
    long nAddedRows = countAddedRows();
 
    return (nRows - nAddedRows);
}

long SparseMatrix::countAddedRows() const
{
    return (m_lastRow + 1);
}

int SparseMatrix::getBlockSize() const
{
    return m_blockSize;
}

long SparseMatrix::getRowCount() const
{
    return m_nRows;
}

long SparseMatrix::getColCount() const
{
    return m_nCols;
}

long SparseMatrix::getRowElementCount() const
{
    long nElements = getBlockSize() * getRowCount();
 
    return nElements;
}

long SparseMatrix::getColElementCount() const
{
    long nElements = getBlockSize() * getColCount();
 
    return nElements;
}

long SparseMatrix::getNZCount() const
{
    return m_nNZ;
}

long SparseMatrix::getRowNZCount(long row) const
{
    return m_pattern.getItemCount(row);
}

long SparseMatrix::getMaxRowNZCount() const
{
    return m_maxRowNZ;
}

long SparseMatrix::getNZElementCount() const
{
    long nNZ = getNZCount();
 
    return getNZElementCount(nNZ);
}

long SparseMatrix::getNZElementCount(long nNZ) const
{
    int blockSize  = getBlockSize();
    long nElements = blockSize * blockSize * nNZ;
 
    return nElements;
}

long SparseMatrix::getRowNZElementCount(long row) const
{
    long nElements = getBlockSize() * getRowNZCount(row);
 
    return nElements;
}

long SparseMatrix::getMaxRowNZElementCount() const
{
    long nElements = getBlockSize() * getMaxRowNZCount();
 
    return nElements;
}
 
#if BITPIT_ENABLE_MPI==1
bool SparseMatrix::isPartitioned() const
{
    return m_partitioned;
}

const MPI_Comm & SparseMatrix::getCommunicator() const
{
    return m_communicator;
}

void SparseMatrix::setCommunicator(MPI_Comm communicator)
{
    if ((communicator != MPI_COMM_NULL) && (communicator != MPI_COMM_SELF)) {
        MPI_Comm_dup(communicator, &m_communicator);
    } else {
        m_communicator = MPI_COMM_SELF;
    }
}

void SparseMatrix::freeCommunicator()
{
    if (m_communicator != MPI_COMM_SELF) {
        int finalizedCalled;
        MPI_Finalized(&finalizedCalled);
        if (!finalizedCalled) {
            MPI_Comm_free(&m_communicator);
        }
    }
}

long SparseMatrix::getRowGlobalCount() const
{
    return m_global_nRows;
}

long SparseMatrix::getRowGlobalOffset() const
{
    return m_global_rowOffset;
}

long SparseMatrix::getRowGlobalElementCount() const
{
    long nElements = getBlockSize() * getRowGlobalCount();
 
    return nElements;
}

long SparseMatrix::getRowGlobalElementOffset() const
{
    long offset = getBlockSize() * getRowGlobalOffset();
 
    return offset;
}

long SparseMatrix::getColGlobalCount() const
{
    return m_global_nCols;
}

long SparseMatrix::getColGlobalOffset() const
{
    return m_global_colOffset;
}

long SparseMatrix::getColGlobalElementCount() const
{
    long nElements = getBlockSize() * getColGlobalCount();
 
    return nElements;
}

long SparseMatrix::getColGlobalElementOffset() const
{
    long offset = getBlockSize() * getColGlobalOffset();
 
    return offset;
}

long SparseMatrix::getNZGlobalCount() const
{
    return m_global_nNZ;
}

long SparseMatrix::getMaxRowNZGlobalCount() const
{
    return m_global_maxRowNZ;
}

long SparseMatrix::getNZGlobalElementCount() const
{
    long nGlobalNZ = getNZGlobalCount();
 
    return getNZElementCount(nGlobalNZ);
}

long SparseMatrix::getMaxRowNZGlobalElementCount() const
{
    long nElement = getBlockSize() * getMaxRowNZGlobalCount();
 
    return nElement;
}

std::vector<long> SparseMatrix::extractLocalGlobalRows() const
{
    const std::size_t *rowExtents = m_pattern.indices();
 
    std::vector<long> localGlobalRows;
    localGlobalRows.reserve(m_nRows);
    for (long i = 0; i < m_nRows; ++i) {
        std::size_t nRowNZ = rowExtents[i + 1] - rowExtents[i];
        if (nRowNZ > 0) {
            localGlobalRows.push_back(m_global_rowOffset + i);
        }
    }
 
    return localGlobalRows;
}

std::vector<long> SparseMatrix::extractGhostGlobalRows() const
{
    return std::vector<long>();
}

std::vector<long> SparseMatrix::extractLocalGlobalCols() const
{
    long firstGlobalCol = m_global_colOffset;
    long lastGlobalCol  = m_global_colOffset + m_nCols - 1;
 
    const long *globalCols = m_pattern.data();
 
    std::vector<long> localGlobalCols;
    localGlobalCols.reserve(m_nNZ);
    for (long k = 0; k < m_nNZ; ++k) {
        long globalCol = globalCols[k];
        if (globalCol < firstGlobalCol || globalCol >= lastGlobalCol) {
            utils::addToOrderedVector<long>(globalCol, localGlobalCols);
        }
    }
 
    return localGlobalCols;
}

std::vector<long> SparseMatrix::extractGhostGlobalCols() const
{
    long firstGlobalCol = m_global_colOffset;
    long lastGlobalCol  = m_global_colOffset + m_nCols - 1;
 
    const long *globalCols = m_pattern.data();
 
    std::vector<long> ghostGlobalCols;
    if (m_nNZ > 0) {
        ghostGlobalCols.reserve(m_nNZ);
        for (long k = 0; k < m_nNZ; ++k) {
            long globalCol = globalCols[k];
            if (globalCol >= firstGlobalCol || globalCol < lastGlobalCol) {
                utils::addToOrderedVector<long>(globalCol, ghostGlobalCols);
            }
        }
    }
 
    return ghostGlobalCols;
}
#endif

void SparseMatrix::addRow(const std::vector<long> &rowPattern, const std::vector<double> &rowValues)
{
    addRow(rowPattern.size(), rowPattern.data(), rowValues.data());
}

void SparseMatrix::addRow(long nRowNZ, const long *rowPattern, const double *rowValues)
{
    if (countMissingRows() == 0) {
        throw std::runtime_error("Unable to add another row: all rows have already been defined.");
    }
 
    // Add the row pattern
    m_pattern.pushBack(nRowNZ, rowPattern);
 
    // Add row values
    const int blockSize = getBlockSize();
    const int nBlockElements = blockSize * blockSize;
 
    m_values.insert(m_values.end(), rowValues, rowValues + nBlockElements * nRowNZ);
 
    // Update the non-zero element counters
    m_maxRowNZ = std::max(nRowNZ, m_maxRowNZ);
    m_nNZ += nRowNZ;
 
    // Update the index of the last row
    m_lastRow++;
}

ConstProxyVector<long> SparseMatrix::getRowPattern(long row) const
{
    ConstProxyVector<long> pattern;
    getRowPattern(row, &pattern);
 
    return pattern;
}

void SparseMatrix::getRowPattern(long row, ConstProxyVector<long> *pattern) const
{
    const std::size_t rowPatternSize  = getRowNZCount(row);
    const long *rowPattern = getRowPatternData(row);
    pattern->set(rowPattern, rowPatternSize);
}

ConstProxyVector<double> SparseMatrix::getRowValues(long row) const
{
    ConstProxyVector<double> values;
    getRowValues(row, &values);
 
    return values;
}

void SparseMatrix::getRowValues(long row, ConstProxyVector<double> *values) const
{
    const int blockSize = getBlockSize();
    const int nBlockElements = blockSize * blockSize;
 
    const std::size_t nRowValues = nBlockElements * getRowNZCount(row);
    const double *rowValues = getRowValuesData(row);
 
    values->set(rowValues, nRowValues);
}

long * SparseMatrix::getRowPatternData(long row)
{
    return const_cast<long *>(const_cast<const SparseMatrix *>(this)->getRowPatternData(row));
}

const long * SparseMatrix::getRowPatternData(long row) const
{
    const std::size_t *rowExtent      = m_pattern.indices(row);
    const std::size_t rowPatternBegin = rowExtent[0];
 
    const long *rowPattern = m_pattern.data() + rowPatternBegin;
 
    return rowPattern;
}

double * SparseMatrix::getRowValuesData(long row)
{
    return const_cast<double *>(const_cast<const SparseMatrix *>(this)->getRowValuesData(row));
}

const double * SparseMatrix::getRowValuesData(long row) const
{
    const int blockSize = getBlockSize();
    const int nBlockElements = blockSize * blockSize;
 
    const std::size_t *rowExtent  = m_pattern.indices(row);
    const std::size_t valuesBegin = nBlockElements * rowExtent[0];
 
    const double *rowValues = m_values.data() + valuesBegin;
 
    return rowValues;
}

std::unique_ptr<SparseMatrix> SparseMatrix::computeTranspose() const
{
    // Early return if the matrix is not assembled.
    if (!isAssembled()) {
        return nullptr;
    }
 
    // Create the transpose matrix
    std::unique_ptr<SparseMatrix> transpose;
#if BITPIT_ENABLE_MPI==1
    transpose = std::unique_ptr<SparseMatrix>(new SparseMatrix(m_communicator, m_partitioned, m_nCols, m_nRows, m_nNZ));
#else
    transpose = std::unique_ptr<SparseMatrix>(new SparseMatrix(m_nCols, m_nRows, m_nNZ));
#endif
 
    // Create an empty transpose pattern
    std::vector<std::size_t> transposeRowSizes(transpose->m_nRows, 0);
    for (int i = 0; i < m_nNZ; ++i) {
        long column = *(m_pattern.data() + i);
        ++transposeRowSizes[column];
    }
 
    transpose->m_pattern.initialize(transpose->m_nRows, transposeRowSizes.data(), static_cast<long>(0));
 
    // Create the empty storage for the values
    transpose->m_values.resize(m_blockSize * m_nNZ);
 
    // Set non-zero information
    transpose->m_nNZ = m_nNZ;
 
    transpose->m_maxRowNZ = 0;
    for (int i = 0; i < transpose->m_nRows; ++i) {
        transpose->m_maxRowNZ = std::max(static_cast<long>(transposeRowSizes[i]), transpose->m_maxRowNZ);
    }
 
    // Fill patter and values of the transpose matrix
    const std::size_t *rowExtents = m_pattern.indices();
    std::fill(transposeRowSizes.begin(), transposeRowSizes.end(), 0);
    for (long row = 0; row < m_nRows; ++row) {
        const std::size_t rowPatternSize = rowExtents[row + 1] - rowExtents[row];
        const long *rowPattern = m_pattern.data() + rowExtents[row];
        const double *rowValues = m_values.data() + rowExtents[row];
        for (std::size_t k = 0; k < rowPatternSize; ++k) {
            long column = rowPattern[k];
            const std::size_t transposeRowLastIndex = *(transpose->m_pattern.indices(column)) + transposeRowSizes[column];
 
            long *transposeRowLastPattern = transpose->m_pattern.data() + transposeRowLastIndex;
            *transposeRowLastPattern = row;
 
            double *transposeRowLastValue = transpose->m_values.data() + transposeRowLastIndex;
            *transposeRowLastValue = rowValues[k];
 
            ++transposeRowSizes[column];
        }
    }
    transpose->m_lastRow = transpose->m_nRows - 1;
 
    // Assembly the transpose matrix
    transpose->assembly();
 
    return transpose;
}
 
}