system_solvers_split.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 "matrix_utilities.hpp"
#include "system_solvers_split.hpp"
 
#include "bitpit_common.hpp"
#include "bitpit_containers.hpp"
#include "bitpit_IO.hpp"
#include "bitpit_operators.hpp"
 
#include "petscksp.h"
#include "petscmat.h"
#include "petscvec.h"
 
#include <fstream>
#include <stdexcept>
#include <string>
#include <numeric>
#include <unordered_set>
 
#ifndef PETSC_NULLPTR
#define PETSC_NULLPTR PETSC_NULL
#endif
 
namespace bitpit {


SplitSystemMatrixAssembler::SplitSystemMatrixAssembler(SystemSplitStrategy splitStrategy, const std::vector<int> &splitSizes)
    : SystemMatrixAssembler(),
      m_splitStrategy(splitStrategy), m_splitSizes(splitSizes)
{
}

SystemSplitStrategy SplitSystemMatrixAssembler::getSplitStrategy() const
{
    return m_splitStrategy;
}

int SplitSystemMatrixAssembler::getSplitCount() const
{
    return m_splitSizes.size();
}

const std::vector<int> & SplitSystemMatrixAssembler::getSplitSizes() const
{
    return m_splitSizes;
}


SplitSystemSparseMatrixAssembler::SplitSystemSparseMatrixAssembler(const SparseMatrix *matrix, SystemSplitStrategy splitStrategy,
                                                                   const std::vector<int> &splitSizes)
    : SystemSparseMatrixAssembler(matrix),
      SplitSystemMatrixAssembler(splitStrategy, splitSizes)
{
}

#if BITPIT_ENABLE_MPI==1
#else
#endif
#if BITPIT_ENABLE_MPI==1
#else
#endif
#if BITPIT_ENABLE_MPI==1
#else
#endif

SplitSystemSolver::SplitSystemSolver(bool debug)
    : SplitSystemSolver("", false, false, debug)
{
}

SplitSystemSolver::SplitSystemSolver(bool transpose, bool debug)
    : SplitSystemSolver("", false, transpose, debug)
{
}

SplitSystemSolver::SplitSystemSolver(bool flatten, bool transpose, bool debug)
    : SplitSystemSolver("", flatten, transpose, debug)
{
}

SplitSystemSolver::SplitSystemSolver(const std::string &prefix, bool debug)
    : SplitSystemSolver(prefix, false, debug)
{
}

SplitSystemSolver::SplitSystemSolver(const std::string &prefix, bool transpose, bool debug)
    : SplitSystemSolver(prefix, false, transpose, debug)
{
}

SplitSystemSolver::SplitSystemSolver(const std::string &prefix, bool flatten, bool transpose, bool debug)
    : SystemSolver(prefix, flatten, transpose, debug),
      m_rhsSplitPermutation(PETSC_NULLPTR), m_solutionSplitPermutation(PETSC_NULLPTR)
{
}

int SplitSystemSolver::getBlockSize() const
{
    if (m_A == PETSC_NULLPTR) {
        return 0;
    }
 
    std::vector<int> splitBlockSizes = getSplitSizes();
    int blockSize = std::accumulate(splitBlockSizes.begin(), splitBlockSizes.end(), 0);
 
    return blockSize;
}

SystemSplitStrategy SplitSystemSolver::getSplitStrategy() const
{
    return m_splitStrategy;
}

int SplitSystemSolver::getSplitCount() const
{
    if (m_A == PETSC_NULLPTR) {
        return -1;
    }
 
    PetscInt nNestedRows;
    MatNestGetSize(m_A, &nNestedRows, PETSC_NULLPTR);
 
    return nNestedRows;
}

std::vector<int> SplitSystemSolver::getSplitSizes() const
{
    if (m_A == PETSC_NULLPTR) {
        return std::vector<int>();
    }
 
    int nSplits = getSplitCount();
    std::vector<int> splitSizes(nSplits, 0);
    for (int k = 0; k < nSplits; ++k) {
        int kk = getBlockSplitLinearIndex(k, k, nSplits);
 
        PetscInt splitBlockSize;
        MatGetBlockSize(m_splitAs[kk], &splitBlockSize);
        splitSizes[k] = splitBlockSize;
    }
 
    return splitSizes;
}

std::vector<int> SplitSystemSolver::getSplitOffsets() const
{
    std::vector<int> splitBlockSizes = getSplitSizes();
    std::vector<int> splitBlockOffsets(splitBlockSizes.size());
    splitBlockOffsets[0] = 0;
    std::partial_sum(splitBlockSizes.begin(), std::prev(splitBlockSizes.end()), std::next(splitBlockOffsets.begin()));
 
    return splitBlockOffsets;
}

void SplitSystemSolver::postKSPSolveActions()
{
    // Base actions
    SystemSolver::postKSPSolveActions();
 
    // Fill statuses of split KSP
    fillSplitKSPStatuses();
}

KSPOptions & SplitSystemSolver::getSplitKSPOptions(int split)
{
    if (!isAssembled()) {
        throw std::runtime_error("The options associated with the Krylov solver can only be accessed after assembling the system.");
    }
 
    return m_splitKSPOptions[split];
}

const KSPOptions & SplitSystemSolver::getSplitKSPOptions(int split) const
{
    if (!isAssembled()) {
        throw std::runtime_error("The options associated with the Krylov solver can only be accessed after assembling the system.");
    }
 
    return m_splitKSPOptions[split];
}

void SplitSystemSolver::initializeKSPOptions()
{
    SystemSolver::initializeKSPOptions();
    initializeSplitKSPOptions();
}

void SplitSystemSolver::initializeSplitKSPOptions()
{
    int nSplits = getSplitCount();
    m_splitKSPOptions.resize(nSplits);
    for (KSPOptions &splitOptions : m_splitKSPOptions) {
        resetKSPOptions(&splitOptions);
    }
}

void SplitSystemSolver::destroyKSPOptions()
{
    SystemSolver::destroyKSPOptions();
    destroySplitKSPOptions();
}

void SplitSystemSolver::destroySplitKSPOptions()
{
    m_splitKSPOptions.clear();
}

const KSPStatus & SplitSystemSolver::getSplitKSPStatus(int split) const
{
    if (!isAssembled()) {
        throw std::runtime_error("The status of the Krylov solver can only be accessed after assembling the system.");
    }
 
    return m_splitKSPStatuses[split];
}

void SplitSystemSolver::initializeKSPStatus()
{
    SystemSolver::initializeKSPStatus();
    initializeSplitKSPStatuses();
}

void SplitSystemSolver::initializeSplitKSPStatuses()
{
    int nSplits = getSplitCount();
    for (int k = static_cast<int>(m_splitKSPStatuses.size()); k < nSplits; ++k) {
        m_splitKSPStatuses.emplace_back();
        KSPStatus *splitStatus = &(m_splitKSPStatuses.back());
        resetKSPStatus(splitStatus);
    }
    m_splitKSPStatuses.resize(nSplits);
}

void SplitSystemSolver::fillKSPStatus()
{
    SystemSolver::fillKSPStatus();
    fillSplitKSPStatuses();
}

void SplitSystemSolver::fillSplitKSPStatuses()
{
    PC pc;
    KSPGetPC(m_KSP, &pc);
 
    PetscInt nSplits;
    KSP *splitKSPs;
    PCFieldSplitGetSubKSP(pc, &nSplits, &splitKSPs);
    for (PetscInt k = 0; k < nSplits; ++k) {
        KSPStatus *splitKSPStatus = &(m_splitKSPStatuses[k]);
        fillKSPStatus(splitKSPs[k], splitKSPStatus);
    }
}

void SplitSystemSolver::resetKSPStatus()
{
    SystemSolver::resetKSPStatus();
    resetSplitKSPStatuses();
}

void SplitSystemSolver::resetSplitKSPStatuses()
{
    for (KSPStatus &splitStatus : m_splitKSPStatuses) {
        resetKSPStatus(&splitStatus);
    }
}

void SplitSystemSolver::destroyKSPStatus()
{
    SystemSolver::destroyKSPStatus();
    destroySplitKSPStatuses();
}

void SplitSystemSolver::destroySplitKSPStatuses()
{
    m_splitKSPStatuses.clear();
}
 

void SplitSystemSolver::dumpInfo(std::ostream &systemStream) const
{
    SystemSolver::dumpInfo(systemStream);
 
    if (systemStream.good()) {
        utils::binary::write(systemStream, m_splitStrategy);
    }
}

void SplitSystemSolver::restoreInfo(std::istream &systemStream)
{
    SystemSolver::restoreInfo(systemStream);
 
    utils::binary::read(systemStream, m_splitStrategy);
}

void SplitSystemSolver::assembly(const SparseMatrix &matrix, SystemSplitStrategy splitStrategy,
                                 const std::vector<int> &splitSizes)
{
    assembly(matrix, splitStrategy, splitSizes, NaturalSystemMatrixOrdering());
}

void SplitSystemSolver::assembly(const SparseMatrix &matrix, SystemSplitStrategy splitStrategy,
                                 const std::vector<int> &splitSizes,
                                 const SystemMatrixOrdering &reordering)
{
    // Check if the matrix is assembled
    if (!matrix.isAssembled()) {
        throw std::runtime_error("Unable to assembly the system. The matrix is not yet assembled.");
    }
 
    // Update matrix
    SplitSystemSparseMatrixAssembler assembler(&matrix, splitStrategy, splitSizes);
    assembly<SplitSystemSolver>(static_cast<const Assembler &>(assembler), reordering);
}

void SplitSystemSolver::assembly(const Assembler &assembler)
{
    assembly(assembler, NaturalSystemMatrixOrdering());
}

void SplitSystemSolver::assembly(const Assembler &assembler, const SystemMatrixOrdering &reordering)
{
    assembly<SplitSystemSolver>(assembler, reordering);
}

void SplitSystemSolver::update(const SparseMatrix &elements)
{
    update(getRowCount(), nullptr, elements);
}

void SplitSystemSolver::update(long nRows, const long *rows, const SparseMatrix &elements)
{
    // Check if the element storage is assembled
    if (!elements.isAssembled()) {
        throw std::runtime_error("Unable to update the system. The element storage is not yet assembled.");
    }
 
    // Update matrix
    SplitSystemSparseMatrixAssembler assembler(&elements, getSplitStrategy(), getSplitSizes());
    update<SplitSystemSolver>(nRows, rows, static_cast<const Assembler &>(assembler));
}

void SplitSystemSolver::update(const Assembler &assembler)
{
    update(getRowCount(), nullptr, assembler);
}

void SplitSystemSolver::update(long nRows, const long *rows, const Assembler &assembler)
{
    update<SplitSystemSolver>(nRows, rows, assembler);
}

void SplitSystemSolver::matrixAssembly(const Assembler &assembler)
{
    const PetscInt *rowReordering = PETSC_NULLPTR;
    if (m_rowReordering) {
        ISGetIndices(m_rowReordering, &rowReordering);
    }
 
    // Matrix information
    long nRows = assembler.getRowCount();
 
    // Split information
    m_splitStrategy = assembler.getSplitStrategy();
 
    const std::vector<int> &splitBlockSizes = assembler.getSplitSizes();
    int nSplits = splitBlockSizes.size();
 
    // Create split matrices
    m_splitAs.assign(nSplits * nSplits, PETSC_NULLPTR);
    for (int splitRow = 0; splitRow < nSplits; ++splitRow) {
        for (int splitCol = 0; splitCol < nSplits; ++splitCol) {
            // Create matrix
            //
            // In LOWER mode, only the lower triangular portion of the splits need to
            // be created, whereas in DIAGONAL mode only the diagonal splits need to
            // be created.
            if (m_splitStrategy == SystemSplitStrategy::SYSTEM_SPLIT_STRATEGY_LOWER) {
                if (splitCol > splitRow) {
                    continue;
                }
            } else if (m_splitStrategy == SystemSplitStrategy::SYSTEM_SPLIT_STRATEGY_DIAGONAL) {
                if (splitCol != splitRow) {
                    continue;
                }
            }
 
            int splitIndex = getBlockSplitLinearIndex(splitRow, splitCol, nSplits);
            Mat *splitMatrix = m_splitAs.data() + splitIndex;
            createMatrix(splitBlockSizes[splitRow], splitBlockSizes[splitCol], splitMatrix);
 
            MatType splitMatrixType;
            MatGetType(*splitMatrix, &splitMatrixType);
 
            // Set matrix sizes
            long nSplitRowsElements = assembler.getRowCount() * splitBlockSizes[splitRow];
            long nSplitColsElements = assembler.getColCount() * splitBlockSizes[splitCol];
 
            long nGlobalSplitRowsElements;
            long nGlobalSplitColsElements;
#if BITPIT_ENABLE_MPI == 1
            nGlobalSplitRowsElements = assembler.getRowGlobalCount() * splitBlockSizes[splitRow];
            nGlobalSplitColsElements = assembler.getColGlobalCount() * splitBlockSizes[splitCol];
#else
            nGlobalSplitRowsElements = nSplitRowsElements;
            nGlobalSplitColsElements = nSplitColsElements;
#endif
 
            MatSetSizes(*splitMatrix, nSplitRowsElements, nSplitColsElements, nGlobalSplitRowsElements, nGlobalSplitColsElements);
 
            // Allocate matrix storage
            //
            // When the internal storage of the system matrix was created without taking into account
            // block information, preallocation information should be provided for each row of each
            // block.
            int rowAllocationExpansion;
            int colAllocationExpansion;
            if (strcmp(splitMatrixType, MATSEQAIJ) == 0) {
                rowAllocationExpansion = splitBlockSizes[splitRow];
                colAllocationExpansion = splitBlockSizes[splitCol];
#if BITPIT_ENABLE_MPI == 1
            } else if (strcmp(splitMatrixType, MATMPIAIJ) == 0) {
                rowAllocationExpansion = splitBlockSizes[splitRow];
                colAllocationExpansion = splitBlockSizes[splitCol];
#endif
            } else {
                rowAllocationExpansion = 1;
                colAllocationExpansion = 1;
            }
 
            long nAllocatedRowElements = rowAllocationExpansion * nRows;
 
            std::vector<int> d_nnz(nAllocatedRowElements, 0);
            for (long n = 0; n < nRows; ++n) {
                long matrixRow = n;
                if (rowReordering) {
                    matrixRow = rowReordering[matrixRow];
                }
 
                int nAssemblerRowNZ = assembler.getRowNZCount(n);
 
                long matrixRowOffset = matrixRow * rowAllocationExpansion;
                for (int n = 0; n < rowAllocationExpansion; ++n) {
                    d_nnz[matrixRowOffset + n] = colAllocationExpansion * nAssemblerRowNZ;
                }
            }
 
#if BITPIT_ENABLE_MPI == 1
            std::vector<int> o_nnz(nAllocatedRowElements, 0);
            if (isPartitioned()) {
                long nAssemblerCols = assembler.getColCount();
 
                long assemblerDiagonalBegin = assembler.getColGlobalOffset();
                long assemblerDiagonalEnd   = assemblerDiagonalBegin + nAssemblerCols;
 
                ConstProxyVector<long> assemblerRowPattern(static_cast<std::size_t>(0), assembler.getMaxRowNZCount());
                for (long n = 0; n < nRows; ++n) {
                    long matrixRow = n;
                    if (rowReordering) {
                        matrixRow = rowReordering[matrixRow];
                    }
 
                    assembler.getRowPattern(n, &assemblerRowPattern);
                    int nAssemblerRowNZ = assemblerRowPattern.size();
 
                    long matrixRowOffset = matrixRow * rowAllocationExpansion;
                    for (int k = 0; k < nAssemblerRowNZ; ++k) {
                        long id = assemblerRowPattern[k];
                        if (id < assemblerDiagonalBegin || id >= assemblerDiagonalEnd) {
                            for (int n = 0; n < rowAllocationExpansion; ++n) {
                                o_nnz[matrixRowOffset + n] += colAllocationExpansion;
                                d_nnz[matrixRowOffset + n] -= colAllocationExpansion;
                            }
                        }
                    }
                }
            }
#endif
 
            if (strcmp(splitMatrixType, MATSEQAIJ) == 0) {
                MatSeqAIJSetPreallocation(*splitMatrix, 0, d_nnz.data());
            } else if (strcmp(splitMatrixType, MATSEQBAIJ) == 0) {
                MatSeqBAIJSetPreallocation(*splitMatrix, splitBlockSizes[splitRow], 0, d_nnz.data());
#if BITPIT_ENABLE_MPI == 1
            } else if (strcmp(splitMatrixType, MATMPIAIJ) == 0) {
                MatMPIAIJSetPreallocation(*splitMatrix, 0, d_nnz.data(), 0, o_nnz.data());
            } else if (strcmp(splitMatrixType, MATMPIBAIJ) == 0) {
                MatMPIBAIJSetPreallocation(*splitMatrix, splitBlockSizes[splitRow], 0, d_nnz.data(), 0, o_nnz.data());
#endif
            } else {
                throw std::runtime_error("Matrix format not supported.");
            }
 
            // Each process will only set values for its own rows
            MatSetOption(*splitMatrix, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE);
 
#if PETSC_VERSION_GE(3, 12, 0)
            // The first assembly will set a superset of the off-process entries
            // required for all subsequent assemblies. This avoids a rendezvous
            // step in the MatAssembly functions.
            MatSetOption(*splitMatrix, MAT_SUBSET_OFF_PROC_ENTRIES, PETSC_TRUE);
#endif
        }
    }
 
    // Create nest matrix
    createMatrix(1, 1, nSplits, nSplits, m_splitAs.data(), &m_A);
 
    // Cleanup
    if (m_rowReordering) {
        ISRestoreIndices(m_rowReordering, &rowReordering);
    }
 
    // Fill matrix
    matrixUpdate(assembler.getRowCount(), nullptr, assembler);
 
    // No new allocations are now allowed
    //
    // When updating the matrix it will not be possible to alter the pattern,
    // it will be possible to change only the values.
    for (Mat splitMatrix : m_splitAs) {
        if (!splitMatrix) {
                continue;
        }
 
        MatSetOption(splitMatrix, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE);
    }
    MatSetOption(m_A, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE);
}

void SplitSystemSolver::matrixFill(const std::string &filePath)
{
    // Check if the matrix exists
    if (!m_A) {
        throw std::runtime_error("Matrix should be created before filling it.");
    }
 
    // Fill the matrix
    int nSplits = getSplitCount();
    for (PetscInt splitRow = 0; splitRow < nSplits; ++splitRow) {
        for (PetscInt splitCol = 0; splitCol < nSplits; ++splitCol) {
            int splitIndex = getBlockSplitLinearIndex(splitRow, splitCol, nSplits);
            std::string splitFilePath = generateSplitPath(filePath, splitRow, splitCol);
            fillMatrix(m_splitAs[splitIndex], filePath);
        }
    }
}

void SplitSystemSolver::matrixUpdate(long nRows, const long *rows, const Assembler &assembler)
{
    // Updating the matrix invalidates the KSP
    m_KSPDirty = true;
 
    // Get split information
    std::vector<int> splitBlockSizes = getSplitSizes();
    std::vector<int> splitBlockOffsets = getSplitOffsets();
    int nSplits = splitBlockSizes.size();
 
    // Get matrix information
    int blockSize = getBlockSize();
 
    std::vector<int> blockSplitRows(blockSize);
    for (int i = 0; i < nSplits; ++i) {
        for (int k = 0; k < splitBlockSizes[i]; ++k) {
            blockSplitRows[splitBlockOffsets[i] + k] = i;
        }
    }
 
    // Get assembler information
    int assemblerBlockSize = assembler.getBlockSize();
    if (assemblerBlockSize != blockSize) {
        std::string message = "Unable to update the matrix.";
        message += " The block size of the assembler is not equal to the block size of the system matrix.";
        throw std::runtime_error(message);
    }
 
    const long assemblerMaxRowNZ = assembler.getMaxRowNZCount();
 
    PetscInt colGlobalBegin;
    PetscInt colGlobalEnd;
    MatGetOwnershipRangeColumn(m_A, &colGlobalBegin, &colGlobalEnd);
    colGlobalBegin /= assemblerBlockSize;
    colGlobalEnd /= assemblerBlockSize;
 
    PetscInt rowGlobalOffset;
    MatGetOwnershipRange(m_A, &rowGlobalOffset, PETSC_NULLPTR);
    rowGlobalOffset /= assemblerBlockSize;
 
    // Initialize reordering
    const PetscInt *rowReordering = PETSC_NULLPTR;
    if (m_rowReordering) {
        ISGetIndices(m_rowReordering, &rowReordering);
    }
 
    const PetscInt *colReordering = PETSC_NULLPTR;
    if (m_colReordering) {
        ISGetIndices(m_colReordering, &colReordering);
    }
 
    // Generate the split indexes
    std::vector<std::vector<std::size_t>> scatterIndexes(nSplits);
    for (int split = 0; split < nSplits; ++split) {
        generateSplitIndexes(split, assemblerMaxRowNZ, scatterIndexes.data() + split);
    }
 
    // Get the options for assembling the matrix
    SystemMatrixAssembler::AssemblyOptions assemblyOptions = assembler.getOptions();
 
#if PETSC_VERSION_GE(3, 12, 0)
    // Check if it is possible to speedup insertion of values
    //
    // The option MAT_SORTED_FULL means that each process provides exactly its
    // local rows; all column indices for a given row are passed in a single call
    // to MatSetValues(), preallocation is perfect, row oriented, INSERT_VALUES
    // is used. If this options is set to PETSC_TRUE, the function MatSetValues
    // will be faster.
    //
    // This options needs at least PETSc 3.12.
    PetscBool matrixSortedFull = (assemblyOptions.full && assemblyOptions.sorted) ? PETSC_TRUE : PETSC_FALSE;
    MatSetOption(m_A, MAT_SORTED_FULL, matrixSortedFull);
#endif
 
    // Insert matrix values
    //
    // If the sizes of PETSc data types match the sizes of data types expected by bitpit
    // and no column reordering is needed, a direct update of the pattern can be performed,
    // otherwise an intermediate pattern storage is required.
    bool patternDirectUpdate = !colReordering && (sizeof(long) == sizeof(PetscInt));
 
    ConstProxyVector<long> rowPattern;
    std::vector<PetscInt> petscRowPatternStorage;
    const PetscInt *petscRowPattern;
    if (!patternDirectUpdate) {
        rowPattern.set(ConstProxyVector<long>::INTERNAL_STORAGE, 0, assemblerMaxRowNZ);
        petscRowPatternStorage.resize(assemblerMaxRowNZ);
        petscRowPattern = petscRowPatternStorage.data();
    }
 
    ConstProxyVector<double> rowValues;
    std::vector<std::vector<PetscScalar>> petscSplitRowValues(m_splitAs.size());
    for (PetscInt splitRow = 0; splitRow < nSplits; ++splitRow) {
        for (PetscInt splitCol = 0; splitCol < nSplits; ++splitCol) {
            int splitIndex = getBlockSplitLinearIndex(splitRow, splitCol, nSplits);
            if (!m_splitAs[splitIndex]) {
                continue;
            }
 
            int splitMatriBlockSize = splitBlockSizes[splitRow] * splitBlockSizes[splitCol];
            petscSplitRowValues[splitIndex].resize(splitMatriBlockSize * assemblerMaxRowNZ);
        }
    }
 
    for (long n = 0; n < nRows; ++n) {
        // Get row information
        long row;
        if (rows) {
            row = rows[n];
        } else {
            row = n;
        }
 
        if (rowReordering) {
            row = rowReordering[row];
        }
 
        const PetscInt globalRow = rowGlobalOffset + row;
 
        // Get row data
        assembler.getRowData(n, &rowPattern, &rowValues);
        if (rowValues.size() == 0) {
            continue;
        }
 
        // Get pattern in PETSc format
        const std::size_t rowPatternSize = rowPattern.size();
        if (patternDirectUpdate) {
            petscRowPattern = reinterpret_cast<const PetscInt *>(rowPattern.data());
        } else {
            for (std::size_t k = 0; k < rowPatternSize; ++k) {
                long globalCol = rowPattern[k];
                if (colReordering) {
                    if (globalCol >= colGlobalBegin && globalCol < colGlobalEnd) {
                        long col = globalCol - colGlobalBegin;
                        col = colReordering[col];
                        globalCol = colGlobalBegin + col;
                    }
                }
 
                petscRowPatternStorage[k] = globalCol;
            }
        }
 
        // Get values in PETSc format
        //
        // The assembler returns the block row values of the matrix in a row-major order.
        // We need to scatter these values to obtain the values of the block rows of the
        // split matrices.
        for (int blockRow = 0; blockRow < blockSize; ++blockRow) {
            int splitRow = blockSplitRows[blockRow];
            int splitBlockRow = blockRow - splitBlockOffsets[splitRow];
 
            std::size_t nRowElements = blockSize * rowPatternSize;
            std::size_t blockRowValuesOffset = nRowElements * blockRow;
            const double *blockRowValues = rowValues.data() + blockRowValuesOffset;
 
            for (int splitCol = 0; splitCol < nSplits; ++splitCol) {
                int splitIndex = getBlockSplitLinearIndex(splitRow, splitCol, nSplits);
                if (!m_splitAs[splitIndex]) {
                    continue;
                }
 
                const std::vector<std::size_t> &splitScatterIndexes = scatterIndexes[splitCol];
                std::size_t nSplitRowElements = splitBlockSizes[splitCol] * rowPatternSize;
                std::size_t splitBlockRowValuesOffset = nSplitRowElements * splitBlockRow;
                double *splitBlockRowValues = petscSplitRowValues[splitIndex].data() + splitBlockRowValuesOffset;
                for (std::size_t k = 0; k < nSplitRowElements; ++k) {
                    splitBlockRowValues[k] = blockRowValues[splitScatterIndexes[k]];
                }
            }
        }
 
        // Insert values
        for (std::size_t k = 0; k < m_splitAs.size(); ++k) {
            if (!m_splitAs[k]) {
                continue;
            }
 
            MatSetValuesBlocked(m_splitAs[k], 1, &globalRow, rowPatternSize, petscRowPattern,
                                petscSplitRowValues[k].data(), INSERT_VALUES);
        }
    }
 
    // Let PETSc assembly the matrix after the update
    for (Mat splitMatrix : m_splitAs) {
        if (!splitMatrix) {
            continue;
        }
 
        MatAssemblyBegin(splitMatrix, MAT_FINAL_ASSEMBLY);
    }
 
    for (Mat splitMatrix : m_splitAs) {
        if (!splitMatrix) {
            continue;
        }
 
        MatAssemblyEnd(splitMatrix, MAT_FINAL_ASSEMBLY);
    }
 
    MatAssemblyBegin(m_A, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(m_A, MAT_FINAL_ASSEMBLY);
 
    // Cleanup
    if (rowReordering) {
        ISRestoreIndices(m_rowReordering, &rowReordering);
    }
 
    if (colReordering) {
        ISRestoreIndices(m_colReordering, &colReordering);
    }
}

void SplitSystemSolver::matrixDump(std::ostream &systemStream, const std::string &directory,
                                   const std::string &prefix) const
{
    // Dump split matrices
    int nSplits = getSplitCount();
    if (systemStream.good()) {
        utils::binary::write(systemStream, nSplits);
    }
 
    for (PetscInt splitRow = 0; splitRow < nSplits; ++splitRow) {
        for (PetscInt splitCol = 0; splitCol < nSplits; ++splitCol) {
            int splitIndex = getBlockSplitLinearIndex(splitRow, splitCol, nSplits);
            std::string splitMatrixName = generateSplitPath(prefix + "A", splitRow, splitCol);
            dumpMatrix(m_splitAs[splitIndex], directory, splitMatrixName);
        }
    }
 
    // Dump information for creating the main matrix
    if (systemStream.good()) {
        PetscInt rowBlockSize;
        PetscInt colBlockSize;
        MatGetBlockSizes(m_A, &rowBlockSize, &colBlockSize);
        utils::binary::write(systemStream, static_cast<int>(rowBlockSize));
        utils::binary::write(systemStream, static_cast<int>(colBlockSize));
    }
}

#if BITPIT_ENABLE_MPI==1
void SplitSystemSolver::matrixRestore(std::istream &systemStream, const std::string &directory,
                                      const std::string &prefix, bool redistribute)
#else
void SplitSystemSolver::matrixRestore(std::istream &systemStream, const std::string &directory,
                                      const std::string &prefix)
#endif
{
    // Restore split matrices
    int nSplits;
    utils::binary::read(systemStream, nSplits);
 
    m_splitAs.assign(nSplits * nSplits, PETSC_NULLPTR);
    for (PetscInt splitRow = 0; splitRow < nSplits; ++splitRow) {
        for (PetscInt splitCol = 0; splitCol < nSplits; ++splitCol) {
            int splitIndex = getBlockSplitLinearIndex(splitRow, splitCol, nSplits);
            std::string splitMatrixName = generateSplitPath(prefix + "A", splitRow, splitCol);
#if BITPIT_ENABLE_MPI==1
            restoreMatrix(directory, splitMatrixName, redistribute, m_splitAs.data() + splitIndex);
#else
            restoreMatrix(directory, splitMatrixName, m_splitAs.data() + splitIndex);
#endif
        }
    }
 
    // Restore main matrix
    int rowBlockSize;
    utils::binary::read(systemStream, rowBlockSize);
    int colBlockSize;
    utils::binary::read(systemStream, colBlockSize);
 
    createMatrix(rowBlockSize, colBlockSize, nSplits, nSplits, m_splitAs.data(), &m_A);
}

void SplitSystemSolver::matrixDestroy()
{
    for (Mat &splitMatrix : m_splitAs) {
        destroyMatrix(&splitMatrix);
    }
 
    destroyMatrix(&m_A);
}

void SplitSystemSolver::vectorsCreate()
{
    // Create vectors
    SystemSolver::vectorsCreate();
 
    // Create the split permutations
    vectorsCreateSplitPermutations();
}

void SplitSystemSolver::vectorsReorder(bool invert)
{
    SystemSolver::vectorsReorder(invert);
 
    reorderVector(m_rhs, m_solutionSplitPermutation, !invert);
    reorderVector(m_solution, m_rhsSplitPermutation, !invert);
}

void SplitSystemSolver::vectorsRestore(std::istream &systemStream, const std::string &directory,
                                       const std::string &prefix)
{
    // Create vectors
    SystemSolver::vectorsRestore(systemStream, directory, prefix);
 
    // Create the split permutations
    vectorsCreateSplitPermutations();
}

void SplitSystemSolver::vectorsCreateSplitPermutations()
{
    int blockSize = getBlockSize();
 
#if BITPIT_ENABLE_MPI == 1
    PetscInt rhsGlobalBegin;
    PetscInt rhsGlobalEnd;
    VecGetOwnershipRange(m_rhs, &rhsGlobalBegin, &rhsGlobalEnd);
    rhsGlobalBegin /= blockSize;
    rhsGlobalEnd /= blockSize;
    PetscInt rhsSize = rhsGlobalEnd - rhsGlobalBegin;
 
    PetscInt solutionGlobalBegin;
    PetscInt solutionGlobalEnd;
    VecGetOwnershipRange(m_solution, &solutionGlobalBegin, &solutionGlobalEnd);
    solutionGlobalBegin /= blockSize;
    solutionGlobalEnd /= blockSize;
    PetscInt solutionSize = solutionGlobalEnd - solutionGlobalBegin;
 
    generateSplitPermutation(rhsSize, static_cast<long>(rhsGlobalBegin), &m_rhsSplitPermutation);
    generateSplitPermutation(solutionSize, static_cast<long>(solutionGlobalBegin), &m_solutionSplitPermutation);
#else
    PetscInt rhsSize;
    VecGetSize(m_rhs, &rhsSize);
    rhsSize /= blockSize;
 
    PetscInt solutionSize;
    VecGetSize(m_solution, &solutionSize);
    solutionSize /= blockSize;
 
    generateSplitPermutation(rhsSize, &m_rhsSplitPermutation);
    generateSplitPermutation(solutionSize, &m_solutionSplitPermutation);
#endif
}

void SplitSystemSolver::vectorsDestroy()
{
    // Destroy split permutations
    if (m_rhsSplitPermutation) {
        ISDestroy(&m_rhsSplitPermutation);
    }
 
    if (m_solutionSplitPermutation) {
        ISDestroy(&m_solutionSplitPermutation);
    }
 
    // Destroy RHS and solution vectors
    SystemSolver::vectorsDestroy();
}

void SplitSystemSolver::exportMatrix(const std::string &filePath, FileFormat fileFormat) const
{
    // Export main matrix
    exportMatrix(m_A, filePath, fileFormat);
 
    // Export split matrices
    int nSplits = getSplitCount();
    for (int splitRow = 0; splitRow < nSplits; ++splitRow) {
        for (int splitCol = 0; splitCol < nSplits; ++splitCol) {
            int splitIndex = getBlockSplitLinearIndex(splitRow, splitCol, nSplits);
            std::string splitFilePath = generateSplitPath(filePath, splitRow, splitCol);
            exportMatrix(m_splitAs[splitIndex], splitFilePath, fileFormat);
        }
    }
}

void SplitSystemSolver::setupPreconditioner()
{
    PC pc;
    KSPGetPC(m_KSP, &pc);
 
    // Setup main preconditioner
    PCSetType(pc, PCFIELDSPLIT);
    if (m_splitStrategy == SystemSplitStrategy::SYSTEM_SPLIT_STRATEGY_LOWER) {
        PCFieldSplitSetType(pc, PC_COMPOSITE_MULTIPLICATIVE);
    } else if (m_splitStrategy == SystemSplitStrategy::SYSTEM_SPLIT_STRATEGY_DIAGONAL) {
        PCFieldSplitSetType(pc, PC_COMPOSITE_ADDITIVE);
    } else {
        PCFieldSplitSetType(pc, PC_COMPOSITE_ADDITIVE);
    }
 
    int nSplits = getSplitCount();
    std::vector<IS> splitIndexSets(nSplits);
    MatNestGetISs(m_A, splitIndexSets.data(), PETSC_NULLPTR);
    for (int k = 0; k < nSplits; ++k) {
        PCFieldSplitSetIS(pc, PETSC_NULLPTR, splitIndexSets[k]);
    }
 
    PCSetUp(pc);
 
    // Set preconditioners of the split
    setupSplitPreconditioners();
}

void SplitSystemSolver::setupSplitPreconditioners()
{
    PC pc;
    KSPGetPC(m_KSP, &pc);
 
    PetscInt nSplits;
    KSP *splitKSPs;
    PCFieldSplitGetSubKSP(pc, &nSplits, &splitKSPs);
 
    for (PetscInt k = 0; k < nSplits; ++k) {
        PC splitPc;
        KSPGetPC(splitKSPs[k], &splitPc);
        const KSPOptions &splitOptions = getSplitKSPOptions(k);
 
        setupPreconditioner(splitPc, splitOptions);
    }
}

void SplitSystemSolver::setupKrylov()
{
    // Setup main Krylov subspace method
    //
    // In LOWER and DIAGONAL modes, the solution of the system is performed by the Krylov subspace
    // methods of the single splits. Thus, in these modes, the main Krylov subspace method only
    // needs to apply the preconditioner.
    if (m_splitStrategy == SystemSplitStrategy::SYSTEM_SPLIT_STRATEGY_LOWER) {
        KSPSetType(m_KSP, KSPPREONLY);
    } else if (m_splitStrategy == SystemSplitStrategy::SYSTEM_SPLIT_STRATEGY_DIAGONAL) {
        KSPSetType(m_KSP, KSPPREONLY);
    } else {
        setupKrylov(m_KSP, getKSPOptions());
    }
 
    // Setup Krylov subspace methods of the splits
    setupSplitKrylovs();
}

void SplitSystemSolver::setupSplitKrylovs()
{
    PC pc;
    KSPGetPC(m_KSP, &pc);
 
    PetscInt nSplits;
    KSP *splitKSPs;
    PCFieldSplitGetSubKSP(pc, &nSplits, &splitKSPs);
 
    for (PetscInt k = 0; k < nSplits; ++k) {
        // In FULL mode, the solution of the system is performed by the Krylov subspace
        // methods of the single splits. Thus, in these modes, the Krylov subspaces of
        // the single blocks only need to apply the preconditioner.
        const KSPOptions &splitOptions = getSplitKSPOptions(k);
        if (m_splitStrategy == SystemSplitStrategy::SYSTEM_SPLIT_STRATEGY_LOWER) {
            setupKrylov(splitKSPs[k], splitOptions);
        } else if (m_splitStrategy == SystemSplitStrategy::SYSTEM_SPLIT_STRATEGY_DIAGONAL) {
            setupKrylov(splitKSPs[k], splitOptions);
        } else {
            KSPSetType(splitKSPs[k], KSPPREONLY);
        }
    }
}

void SplitSystemSolver::preKrylovSetupActions()
{
    // Execute base actions
    SystemSolver::preKrylovSetupActions();
 
    // Enable convergence monitor
    if (m_convergenceMonitorEnabled) {
        int nSplits = getSplitCount();
        for (int split = 0; split < nSplits; ++split) {
            std::string prefix = "fieldsplit_" + std::to_string(split) + "_";
 
#if (PETSC_VERSION_MAJOR >= 3 && PETSC_VERSION_MINOR >= 7)
            PetscOptionsSetValue(PETSC_NULLPTR, ("-" + prefix + "ksp_monitor_true_residual").c_str(), "");
            PetscOptionsSetValue(PETSC_NULLPTR, ("-" + prefix + "ksp_converged_reason").c_str(), "");
            PetscOptionsSetValue(PETSC_NULLPTR, ("-" + prefix + "ksp_monitor_singular_value").c_str(), "");
#else
            PetscOptionsSetValue(("-" + prefix + "ksp_monitor_true_residual").c_str(), "");
            PetscOptionsSetValue(("-" + prefix + "ksp_converged_reason").c_str(), "");
            PetscOptionsSetValue(("-" + prefix + "ksp_monitor_singular_value").c_str(), "");
#endif
        }
    }
}

#if BITPIT_ENABLE_MPI == 1
#endif
#if BITPIT_ENABLE_MPI == 1
void SplitSystemSolver::generateSplitPermutation(long nItems, long itemsGlobalOffset, IS *permutation) const
#else
void SplitSystemSolver::generateSplitPermutation(long nItems, IS *permutation) const
#endif
{
   std::vector<std::size_t> workspace;
 
    // Get split information
    std::vector<int> splitBlockSizes = getSplitSizes();
    std::vector<int> splitBlockOffsets = getSplitOffsets();
    int nSplits = splitBlockSizes.size();
 
    // Matrix information
    int blockSize = getBlockSize();
    std::size_t nElements = blockSize * nItems;
 
    // Create permutations
    PetscInt elementsGlobalOffset;
#if BITPIT_ENABLE_MPI == 1
    elementsGlobalOffset = static_cast<PetscInt>(blockSize * itemsGlobalOffset);
#else
    elementsGlobalOffset = 0;
#endif
 
    PetscInt *permutationStorage;
    PetscMalloc(nElements * sizeof(PetscInt), &permutationStorage);
 
    PetscInt *splitPermutationStorage = permutationStorage;
    for (int split = 0; split < nSplits; ++split) {
        int splitBlockSize = splitBlockSizes[split];
        std::size_t nSplitElements = splitBlockSize * nItems;
 
        workspace.resize(nSplitElements);
        generateSplitIndexes(split, nItems, &workspace);
        for (std::size_t k = 0; k < nSplitElements; ++k) {
            splitPermutationStorage[k] = elementsGlobalOffset + static_cast<PetscInt>(workspace[k]);
 
        }
        splitPermutationStorage += nSplitElements;
    }
 
#if BITPIT_ENABLE_MPI == 1
    ISCreateGeneral(getCommunicator(), nItems, permutationStorage, PETSC_OWN_POINTER, permutation);
#else
    ISCreateGeneral(PETSC_COMM_SELF, nItems, permutationStorage, PETSC_OWN_POINTER, permutation);
#endif
    ISSetPermutation(*permutation);
}

void SplitSystemSolver::generateSplitIndexes(int split, long nItems, std::vector<std::size_t> *indexes) const
{
    // Get split information
    int splitBlockSize   = getSplitSizes()[split];
    int splitBlockOffset = getSplitOffsets()[split];
 
    // Get matrix information
    int blockSize = getBlockSize();
 
    // Initialize index storage
    indexes->resize(nItems * splitBlockSize);
 
    // Evaluate the indexes
    for (long k = 0; k < nItems; ++k) {
        for (int n = 0; n < splitBlockSize; ++n) {
            (*indexes)[splitBlockSize * k + n] = blockSize * k + splitBlockOffset + n;
        }
    }
}

std::string SplitSystemSolver::generateSplitPath(const std::string &path, int i) const
{
    return generateSplitPath(path, std::to_string(i));
}

std::string SplitSystemSolver::generateSplitPath(const std::string &path, int i, int j) const
{
    return generateSplitPath(path, std::to_string(i) + std::to_string(j));
}

std::string SplitSystemSolver::generateSplitPath(const std::string &path, const std::string &index) const
{
    std::string splitLabel = "_split_" + index;
 
    std::size_t dotPosition = path.find_last_of(".");
 
    std::string splitPath;
    std::string extension;
    std::string basePath;
    if (dotPosition != std::string::npos) {
        extension = path.substr(dotPosition + 1);
        basePath  = path.substr(0, dotPosition);
        splitPath = basePath + splitLabel + "." + extension;
    } else {
        splitPath = path + splitLabel;
    }
 
    return splitPath;
}

int SplitSystemSolver::getBlockSplitLinearIndex(int i, int j) const
{
    int nSplits = getSplitCount();
 
    return getBlockSplitLinearIndex(i, j, nSplits);
}

int SplitSystemSolver::getBlockSplitLinearIndex(int i, int j, int nSplits) const
{
    return linearalgebra::linearIndexRowMajor(i, j, nSplits, nSplits);
}
 
}