system_solvers_large.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_solvers_large.hpp"
 
#include "bitpit_common.hpp"
#include "bitpit_containers.hpp"
#include "bitpit_IO.hpp"
#include "bitpit_operators.hpp"
 
#include "petscksp.h"
#include "petscmat.h"
#include "petscsystypes.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 {



long NaturalSystemMatrixOrdering::getRowPermutationRank(long row) const
{
    return row;
}

long NaturalSystemMatrixOrdering::getColPermutationRank(long col) const
{
    return col;
}



SystemSparseMatrixAssembler::SystemSparseMatrixAssembler(const SparseMatrix *matrix)
    : SystemMatrixAssembler(), m_matrix(matrix)
{
}
 
#if BITPIT_ENABLE_MPI==1
bool SystemSparseMatrixAssembler::isPartitioned() const
{
    return m_matrix->isPartitioned();
}

const MPI_Comm & SystemSparseMatrixAssembler::getCommunicator() const
{
    return m_matrix->getCommunicator();
}
#endif

SystemMatrixAssembler::AssemblyOptions SystemSparseMatrixAssembler::getOptions() const
{
    AssemblyOptions options;
    options.full   = true;
    options.sorted = false;
 
    return options;
}

bool SystemSolver::getTranspose() const
{
    return m_transpose;
}

void SystemSolver::setTranspose(bool transpose)
{
    // Early return if the transpose flag is already correctly set.
    if (m_transpose == transpose) {
        return;
    }
 
    // Clear the workspace
    clearWorkspace();
 
    // Re-create the vectors
    if (isAssembled()) {
        VecDestroy(&m_rhs);
        VecDestroy(&m_solution);
 
        vectorsCreate();
    }
 
    // Set the flag
    m_transpose = transpose;
}

int SystemSparseMatrixAssembler::getBlockSize() const
{
    return m_matrix->getBlockSize();
}

long SystemSparseMatrixAssembler::getRowCount() const
{
    return m_matrix->getRowCount();
}

long SystemSparseMatrixAssembler::getColCount() const
{
    return m_matrix->getColCount();
}

long SystemSparseMatrixAssembler::getRowElementCount() const
{
    long nRowElements = getBlockSize() * getRowCount();
 
    return nRowElements;
}

long SystemSparseMatrixAssembler::getColElementCount() const
{
    long nColElements = getBlockSize() * getColCount();
 
    return nColElements;
}
 
#if BITPIT_ENABLE_MPI==1
/**
 * Get number of global (block) rows handled by the assembler.
 *
 * If the matrix is a block matrix (i.e., the block size is greater than one),
 * this function will return the global number of block rows, where a block row
 * is defined as a group of block-size matrix rows.
 *
 * \result The number of global rows handled by the assembler.
 */
long SystemSparseMatrixAssembler::getRowGlobalCount() const
{
    return m_matrix->getRowGlobalCount();
}

/**
 * Get number of global (block) columns handled by the assembler.
 *
 * If the matrix is a block matrix (i.e., the block size is greater than one),
 * this function will return the global number of block columns, where a block
 * column is defined as a group of block-size matrix columns.
 *
 * \result The number of global (block) columns handled by the assembler.
 */
long SystemSparseMatrixAssembler::getColGlobalCount() const
{
    return m_matrix->getColGlobalCount();
}

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

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

long SystemSparseMatrixAssembler::getRowGlobalOffset() const
{
    return m_matrix->getRowGlobalOffset();
}

long SystemSparseMatrixAssembler::getColGlobalOffset() const
{
    return m_matrix->getColGlobalOffset();
}

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

long SystemSparseMatrixAssembler::getColGlobalElementOffset() const
{
    long offset = getBlockSize() * getColGlobalOffset();
 
    return offset;
}
#endif

long SystemSparseMatrixAssembler::getRowNZCount(long rowIndex) const
{
    return m_matrix->getRowNZCount(rowIndex);
}

long SystemSparseMatrixAssembler::getMaxRowNZCount() const
{
    return m_matrix->getMaxRowNZCount();
}

 */
void SystemSparseMatrixAssembler::getRowPattern(long rowIndex, ConstProxyVector<long> *pattern) const
{
    m_matrix->getRowPattern(rowIndex, pattern);
}

 *
 * \param rowIndex is the index of the row in the assembler
 * \param values on output will contain the values of the specified (block) row.
 * If the block size is greater than one, values will be stored in a logically
 * two-dimensional array that uses a col-major order
 */
void SystemSparseMatrixAssembler::getRowValues(long rowIndex, ConstProxyVector<double> *values) const
{
    m_matrix->getRowValues(rowIndex, values);
}

 *
 * \param rowIndex is the index of the row in the assembler
 * \param pattern on output will contain the pattern of the specified (block) row
 * \param values on output will contain the values of the specified (block) row.
 * If the block size is greater than one, values will be stored in a logically
 * two-dimensional array that uses a col-major order
 */
void SystemSparseMatrixAssembler::getRowData(long rowIndex, ConstProxyVector<long> *pattern, ConstProxyVector<double> *values) const
{
    m_matrix->getRowPattern(rowIndex, pattern);
    m_matrix->getRowValues(rowIndex, values);
}


PetscErrorCode PetscManager::displayLogView()
{
    return PetscLogView(PETSC_VIEWER_STDOUT_WORLD);
}

void PetscManager::addInitOption(const std::string &option)
{
    if (!areOptionsEditable()) {
        throw std::runtime_error("Initialization options can be set only before initializing the solver.");
    }
 
    m_options.push_back(option);
}

void PetscManager::addInitOptions(int argc, char **argv)
{
    if (!areOptionsEditable()) {
        throw std::runtime_error("Initialization options can be set only before initializing the solver.");
    }
 
    for (int i = 1; i < argc; ++i) {
        m_options.push_back(argv[i]);
    }
}

void PetscManager::addInitOptions(const std::vector<std::string> &options)
{
    if (!areOptionsEditable()) {
        throw std::runtime_error("Initialization options can be set only before initializing the solver.");
    }
 
    for (const std::string &option : options) {
        m_options.push_back(option);
    }
}

void PetscManager::clearInitOptions()
{
    m_options.clear();
}

void PetscManager::enableLogView()
{
    if (m_logViewEnabled) {
        return;
    }
 
#if (PETSC_VERSION_MAJOR >= 3 && PETSC_VERSION_MINOR >= 7)
    PetscLogDefaultBegin();
#else
    PetscLogBegin();
#endif
    PetscRegisterFinalize(&(PetscManager::displayLogView));
    m_logViewEnabled = true;
}

PetscManager::PetscManager()
    : m_externalMPIInitialization(true), m_externalPETScInitialization(true),
      m_options(std::vector<std::string>(1, "bitpit")),  m_logViewEnabled(false)
{
}

PetscManager::~PetscManager()
{
    finalize();
}

bool PetscManager::areOptionsEditable() const
{
    PetscBool isPetscInitialized;
    PetscInitialized(&isPetscInitialized);
 
    return (!isPetscInitialized);
}

#if BITPIT_ENABLE_MPI==1
#endif
bool PetscManager::initialize(bool debug)
{
    // Early return if PETSc is already initialized
    PetscBool isPetscInitialized;
    PetscInitialized(&isPetscInitialized);
    if (isPetscInitialized) {
        return false;
    }
 
#if BITPIT_ENABLE_MPI==1
    // Early return if MPI is already finalized
    int isMPIFinalized;
    MPI_Finalized(&isMPIFinalized);
    if (isMPIFinalized) {
        throw std::runtime_error("PETSc finalization cannot be called after MPI finaliation.");
    }
#endif
 
    // Generate command line arguments
    //
    // The first argument is the executable name and it is set to a
    // dummy value.
    std::string help        = "None";
    std::string programName = "bitpit_petsc_manager";
 
    int argc = 1 + m_options.size();
    char **argv = new char*[argc + 1];
    argv[0] = strdup(programName.data());
    for (std::size_t i = 0; i < m_options.size(); i++) {
        argv[1 + i] = strdup(m_options[i].data());
    }
    argv[argc] = nullptr;
 
#if BITPIT_ENABLE_MPI==1
    // Initialize MPI
    int isMPIInitialized;
    MPI_Initialized(&isMPIInitialized);
    if (!isMPIInitialized) {
        MPI_Init(&argc, &argv);
        m_externalMPIInitialization = false;
    }
#endif
 
    // Initialize PETSc
    PetscInitialize(&argc, &argv, 0, help.data());
    m_externalPETScInitialization = false;
 
    // Enable log view
    if (debug) {
        enableLogView();
    }
 
    // Clean-up command line arguments
    for (int i = 0; i < argc; ++i) {
        free(argv[i]);
    }
 
    delete[] argv;
 
    // Initialization completed
    return true;
}

#if BITPIT_ENABLE_MPI==1
#else
#endif
bool PetscManager::finalize()
{
#if BITPIT_ENABLE_MPI==1
    // Early return if MPI is already finalized
    int isMPIFinalized;
    MPI_Finalized(&isMPIFinalized);
    if (isMPIFinalized) {
        return false;
    }
#endif
 
    // Early return if PETSc was initialized externally
    if (m_externalPETScInitialization) {
        return false;
    }
 
    // Finalize PETSc
    PetscBool isPetscFinalized;
    PetscFinalized(&isPetscFinalized);
    if (!isPetscFinalized) {
        PetscFinalize();
    }
 
#if BITPIT_ENABLE_MPI==1
    // Finalize MPI
    if (!m_externalMPIInitialization) {
        MPI_Finalize();
    }
#endif
 
    return (!isPetscFinalized);
}

#if BITPIT_ENABLE_MPI==1
#else
#endif
 
PetscManager SystemSolver::m_petscManager = PetscManager();
 
int SystemSolver::m_nInstances = 0;

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

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

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

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

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

SystemSolver::SystemSolver(const std::string &prefix, bool flatten, bool transpose, bool debug)
    : m_flatten(flatten), m_transpose(transpose),
      m_A(PETSC_NULLPTR), m_rhs(PETSC_NULLPTR), m_solution(PETSC_NULLPTR),
      m_rowReordering(PETSC_NULLPTR), m_colReordering(PETSC_NULLPTR),
      m_convergenceMonitorEnabled(debug),
      m_KSP(PETSC_NULLPTR), m_KSPDirty(true),
      m_prefix(prefix), m_assembled(false),
#if BITPIT_ENABLE_MPI==1
      m_communicator(MPI_COMM_SELF), m_partitioned(false),
#endif
      m_forceConsistency(false)
{
    // Initialize PETSc
    if (m_nInstances == 0) {
        m_petscManager.initialize(debug);
    }
 
    // Increase the number of instances
    ++m_nInstances;
}

SystemSolver::~SystemSolver()
{
    // Clear the solver
    clear();
 
    // Decrease the number of instances
    --m_nInstances;
}

void SystemSolver::clear()
{
    clearWorkspace();
 
    destroyKSPOptions();
    destroyKSPStatus();
 
    vectorsDestroy();
    matrixDestroy();
 
    clearReordering();
 
#if BITPIT_ENABLE_MPI==1
    freeCommunicator();
#endif
 
    m_assembled = false;
}

void SystemSolver::clearWorkspace()
{
    // Early return if the KSP has not been created
    if (!m_KSP) {
        return;
    }
 
    // Destroy the KSP
    destroyKSP();
}

void SystemSolver::assembly(const SparseMatrix &matrix)
{
    assembly(matrix, NaturalSystemMatrixOrdering());
}

void SystemSolver::assembly(const SparseMatrix &matrix, 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.");
    }
 
    // Assembly the system matrix
    SystemSparseMatrixAssembler assembler(&matrix);
    assembly<SystemSolver>(static_cast<const Assembler &>(assembler), reordering);
}

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

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

void SystemSolver::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
    SystemSparseMatrixAssembler assembler(&elements);
    update<SystemSolver>(nRows, rows, static_cast<const Assembler &>(assembler));
}

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

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

int SystemSolver::getBlockSize() const
{
    if (m_A == PETSC_NULLPTR) {
        return 0;
    }
 
    PetscInt blockSize;
    MatGetBlockSize(m_A, &blockSize);
 
    return blockSize;
}

long SystemSolver::getRowCount() const
{
    if (m_A == PETSC_NULLPTR) {
        return 0;
    }
 
    PetscInt nRows;
    MatGetLocalSize(m_A, &nRows, NULL);
    nRows /= getBlockSize();
 
    return nRows;
}

long SystemSolver::getColCount() const
{
    if (m_A == PETSC_NULLPTR) {
        return 0;
    }
 
    PetscInt nCols;
    MatGetLocalSize(m_A, NULL, &nCols);
    nCols /= getBlockSize();
 
    return nCols;
}

long SystemSolver::getRowElementCount() const
{
    if (m_A == PETSC_NULLPTR) {
        return 0;
    }
 
    PetscInt nRows;
    MatGetLocalSize(m_A, &nRows, NULL);
 
    return nRows;
}

long SystemSolver::getColElementCount() const
{
    if (m_A == PETSC_NULLPTR) {
        return 0;
    }
 
    PetscInt nCols;
    MatGetLocalSize(m_A, NULL, &nCols);
 
    return nCols;
}
 
#if BITPIT_ENABLE_MPI==1
long SystemSolver::getRowGlobalCount() const
{
    if (m_A == PETSC_NULLPTR) {
        return 0;
    }
 
    PetscInt nRows;
    MatGetSize(m_A, &nRows, NULL);
    nRows /= getBlockSize();
 
    return nRows;
}

long SystemSolver::getColGlobalCount() const
{
    if (m_A == PETSC_NULLPTR) {
        return 0;
    }
 
    PetscInt nCols;
    MatGetSize(m_A, NULL, &nCols);
    nCols /= getBlockSize();
 
    return nCols;
}

long SystemSolver::getRowGlobalElementCount() const
{
    if (m_A == PETSC_NULLPTR) {
        return 0;
    }
 
    PetscInt nRows;
    MatGetSize(m_A, &nRows, NULL);
 
    return nRows;
}

long SystemSolver::getColGlobalElementCount() const
{
    if (m_A == PETSC_NULLPTR) {
        return 0;
    }
 
    PetscInt nCols;
    MatGetSize(m_A, NULL, &nCols);
 
    return nCols;
}

bool SystemSolver::isPartitioned() const
{
    return m_partitioned;
}
#endif

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

void SystemSolver::solve()
{
    // Check if the system is assembled
    if (!isAssembled()) {
        throw std::runtime_error("Unable to solve the system. The system is not yet assembled.");
    }
 
    // Perform actions before KSP solution
    preKSPSolveActions();
 
    // Solve KSP
    solveKSP();
 
    // Perform actions after KSP solution
    postKSPSolveActions();
}

void SystemSolver::solve(const std::vector<double> &rhs, std::vector<double> *solution)
{
    // Fills the vectors
    vectorsFill(rhs, *solution);
 
    // Solve the system
    solve();
 
    // Export the solution
    exportVector(m_solution, solution);
}

void SystemSolver::solveKSP()
{
    PetscErrorCode solverError;
    if (!m_transpose) {
        solverError = KSPSolve(m_KSP, m_rhs, m_solution);
    } else {
        solverError = KSPSolveTranspose(m_KSP, m_rhs, m_solution);
    }
 
    if (solverError) {
        const char *petscMessage = nullptr;
        PetscErrorMessage(solverError, &petscMessage, PETSC_NULLPTR);
        std::string message = "Unable to solver the system. " + std::string(petscMessage);
        throw std::runtime_error(message);
    }
}

void SystemSolver::preKSPSolveActions()
{
    // Prepare KSP
    prepareKSP();
 
    // Reorder vectors
    vectorsReorder(true);
 
    // Force consistency
    if (m_forceConsistency) {
        removeNullSpaceFromRHS();
    }
}

void SystemSolver::postKSPSolveActions()
{
    // Fill status of KSP
    fillKSPStatus();
 
    // Reorder vectors
    vectorsReorder(false);
 
    // Finalize KSP
    finalizeKSP();
}

int SystemSolver::getDumpVersion() const
{
    const int DUMP_VERSION = 1;
 
    return DUMP_VERSION;
}

void SystemSolver::matrixAssembly(const Assembler &assembler)
{
    const PetscInt *rowReordering = PETSC_NULLPTR;
    if (m_rowReordering) {
        ISGetIndices(m_rowReordering, &rowReordering);
    }
 
    // Create the matrix
    int blockSize = assembler.getBlockSize();
    createMatrix(blockSize, blockSize, &m_A);
 
    MatType matrixType;
    MatGetType(m_A, &matrixType);
 
    // Get sizes
    long nAssemblerRows = assembler.getRowCount();
 
    long nRowsElements = assembler.getRowElementCount();
    long nColsElements = assembler.getColElementCount();
 
    long nGlobalRowsElements;
    long nGlobalColsElements;
#if BITPIT_ENABLE_MPI == 1
    nGlobalRowsElements = assembler.getRowGlobalElementCount();
    nGlobalColsElements = assembler.getColGlobalElementCount();
#else
    nGlobalRowsElements = nRowsElements;
    nGlobalColsElements = nColsElements;
#endif
 
    MatSetSizes(m_A, nRowsElements, nColsElements, nGlobalRowsElements, nGlobalColsElements);
 
    // Allocate 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 allocationExpansion;
    if (strcmp(matrixType, MATSEQAIJ) == 0) {
        allocationExpansion = blockSize;
#if BITPIT_ENABLE_MPI == 1
    } else if (strcmp(matrixType, MATMPIAIJ) == 0) {
        allocationExpansion = blockSize;
#endif
    } else {
        allocationExpansion = 1;
    }
 
    long nAllocatedRows = allocationExpansion * nAssemblerRows;
 
    std::vector<int> d_nnz(nAllocatedRows, 0);
    for (long n = 0; n < nAssemblerRows; ++n) {
        long matrixRow = n;
        if (rowReordering) {
            matrixRow = rowReordering[matrixRow];
        }
 
        int nAssemblerRowNZ = assembler.getRowNZCount(n);
 
        long matrixRowOffset = matrixRow * allocationExpansion;
        for (int i = 0; i < allocationExpansion; ++i) {
            d_nnz[matrixRowOffset + i] = allocationExpansion * nAssemblerRowNZ;
        }
    }
 
#if BITPIT_ENABLE_MPI == 1
    std::vector<int> o_nnz(nAllocatedRows, 0);
    if (m_partitioned) {
        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 < nAssemblerRows; ++n) {
            long matrixRow = n;
            if (rowReordering) {
                matrixRow = rowReordering[matrixRow];
            }
 
            assembler.getRowPattern(n, &assemblerRowPattern);
            int nAssemblerRowNZ = assemblerRowPattern.size();
 
            long matrixRowOffset = matrixRow * allocationExpansion;
            for (int k = 0; k < nAssemblerRowNZ; ++k) {
                long id = assemblerRowPattern[k];
                if (id < assemblerDiagonalBegin || id >= assemblerDiagonalEnd) {
                    for (int i = 0; i < allocationExpansion; ++i) {
                        o_nnz[matrixRowOffset + i] += allocationExpansion;
                        d_nnz[matrixRowOffset + i] -= allocationExpansion;
                    }
                }
            }
        }
    }
#endif
 
    if (strcmp(matrixType, MATSEQAIJ) == 0) {
        MatSeqAIJSetPreallocation(m_A, 0, d_nnz.data());
    } else if (strcmp(matrixType, MATSEQBAIJ) == 0) {
        MatSeqBAIJSetPreallocation(m_A, blockSize, 0, d_nnz.data());
#if BITPIT_ENABLE_MPI == 1
    } else if (strcmp(matrixType, MATMPIAIJ) == 0) {
        MatMPIAIJSetPreallocation(m_A, 0, d_nnz.data(), 0, o_nnz.data());
    } else if (strcmp(matrixType, MATMPIBAIJ) == 0) {
        MatMPIBAIJSetPreallocation(m_A, blockSize, 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(m_A, 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(m_A, MAT_SUBSET_OFF_PROC_ENTRIES, PETSC_TRUE);
#endif
 
    // 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.
    MatSetOption(m_A, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE);
}

void SystemSolver::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
    fillMatrix(m_A, filePath);
}

void SystemSolver::matrixUpdate(long nRows, const long *rows, const Assembler &assembler)
{
    // Updating the matrix invalidates the KSP
    m_KSPDirty = true;
 
    // Get block size
    const int blockSize = getBlockSize();
    if (assembler.getBlockSize() != 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);
    }
 
    // 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);
    }
 
    // Global information
    PetscInt colGlobalBegin;
    PetscInt colGlobalEnd;
    MatGetOwnershipRangeColumn(m_A, &colGlobalBegin, &colGlobalEnd);
    colGlobalBegin /= blockSize;
    colGlobalEnd /= blockSize;
 
    PetscInt rowGlobalOffset;
    MatGetOwnershipRange(m_A, &rowGlobalOffset, PETSC_NULLPTR);
    rowGlobalOffset /= blockSize;
 
    // 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
 
    // Check if it possible to perform a fast update
    //
    // A fast update allows to set all the values of a row at once (without
    // the need to get the row pattern), it can be performed if:
    //  - the system matrix has already been assembled;
    //  - the system matrix has a unitary block size;
    //  - the assembler is providing all the values of the row;
    //  - values provided by the assembler are sorted by ascending column.
    //
    // If fast update is used, row values will be set using a special PETSc
    // function (MatSetValuesRow) that allows to set all the values of a
    // row at once, without requiring the pattern of the row.
    //
    // Fast update is not related to the option MAT_SORTED_FULL, that option
    // is used to speedup the standard function MatSetValues (which still
    // requires the pattern of the row).
    bool fastUpdate = isAssembled() && (blockSize == 1) && assemblyOptions.full && assemblyOptions.sorted;
 
    // Update element values
    //
    // If the sizes of PETSc data types match the sizes of data types expected by
    // bitpit a direct update can be performed, otherwise the matrix is updated
    // using intermediate data storages.
    const long assemblerMaxRowNZ = std::max(assembler.getMaxRowNZCount(), 0L);
 
    bool patternDirectUpdate = !colReordering && (sizeof(long) == sizeof(PetscInt));
    bool valuesDirectUpdate  = (sizeof(double) == sizeof(PetscScalar));
 
    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<PetscScalar> petscRowValuesStorage;
    const PetscScalar *petscRowValues;
    if (!valuesDirectUpdate) {
        long assemblerMaxRowNZElements = blockSize * blockSize * assemblerMaxRowNZ;
        rowValues.set(ConstProxyVector<double>::INTERNAL_STORAGE, 0, assemblerMaxRowNZElements);
        petscRowValuesStorage.resize(assemblerMaxRowNZElements);
        petscRowValues = petscRowValuesStorage.data();
    }
 
    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
        if (fastUpdate) {
            assembler.getRowValues(n, &rowValues);
        } else {
            assembler.getRowData(n, &rowPattern, &rowValues);
        }
 
        if (rowValues.size() == 0) {
            continue;
        }
 
        // Get values in PETSc format
        const std::size_t rowPatternSize = rowPattern.size();
 
        if (valuesDirectUpdate) {
            petscRowValues = reinterpret_cast<const PetscScalar *>(rowValues.data());
        } else {
            std::copy(rowValues.cbegin(), rowValues.cend(), petscRowValuesStorage.begin());
        }
 
        if (fastUpdate) {
            MatSetValuesRow(m_A, globalRow, petscRowValues);
        } else {
            // Get pattern in PETSc format
            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;
                }
            }
 
            // Set data
            if (blockSize > 1) {
                MatSetValuesBlocked(m_A, 1, &globalRow, rowPatternSize, petscRowPattern, petscRowValues, INSERT_VALUES);
            } else {
                MatSetValues(m_A, 1, &globalRow, rowPatternSize, petscRowPattern, petscRowValues, INSERT_VALUES);
            }
        }
    }
 
    // Let petsc assembly the matrix after the update
    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 SystemSolver::matrixDump(std::ostream &systemStream, const std::string &directory,
                              const std::string &prefix) const
{
    BITPIT_UNUSED(systemStream);
 
    dumpMatrix(m_A, directory, prefix + "A");
}

#if BITPIT_ENABLE_MPI==1
void SystemSolver::matrixRestore(std::istream &systemStream, const std::string &directory,
                                 const std::string &prefix, bool redistribute)
#else
void SystemSolver::matrixRestore(std::istream &systemStream, const std::string &directory,
                                 const std::string &prefix)
#endif
{
    BITPIT_UNUSED(systemStream);
 
#if BITPIT_ENABLE_MPI==1
    restoreMatrix(directory, prefix + "A", redistribute, &m_A);
#else
    restoreMatrix(directory, prefix + "A", &m_A);
#endif
}

void SystemSolver::matrixDestroy()
{
    destroyMatrix(&m_A);
}

void SystemSolver::vectorsCreate()
{
    if (!m_transpose) {
        MatCreateVecs(m_A, &m_solution, &m_rhs);
    } else {
        MatCreateVecs(m_A, &m_rhs, &m_solution);
    }
}

void SystemSolver::vectorsFill(const std::vector<double> &rhs, const std::vector<double> &solution)
{
    if (!rhs.empty()) {
        fillVector(m_rhs, rhs);
    }
 
    if (!solution.empty()) {
        fillVector(m_solution, solution);
    }
}

void SystemSolver::vectorsFill(const std::string &rhsFilePath, const std::string &solutionFilePath)
{
    if (!rhsFilePath.empty()) {
        fillVector(m_rhs, rhsFilePath);
    }
 
    if (!solutionFilePath.empty()) {
        fillVector(m_solution, solutionFilePath);
    }
}

void SystemSolver::vectorsReorder(bool invert)
{
    if (!m_transpose) {
        reorderVector(m_rhs, m_colReordering, invert);
        reorderVector(m_solution, m_rowReordering, invert);
    } else {
        reorderVector(m_rhs, m_rowReordering, invert);
        reorderVector(m_solution, m_colReordering, invert);
    }
}

void SystemSolver::vectorsDump(std::ostream &systemStream, const std::string &directory,
                               const std::string &prefix) const
{
    BITPIT_UNUSED(systemStream);
 
    dumpVector(m_rhs, directory, prefix + "rhs");
    dumpVector(m_solution, directory, prefix + "solution");
}

void SystemSolver::vectorsRestore(std::istream &systemStream, const std::string &directory,
                                  const std::string &prefix)
{
    BITPIT_UNUSED(systemStream);
 
    restoreVector(directory, prefix + "rhs", m_A, VECTOR_SIDE_RIGHT, &m_rhs);
    restoreVector(directory, prefix + "solution", m_A, VECTOR_SIDE_LEFT, &m_solution);
}

void SystemSolver::vectorsDestroy()
{
    destroyVector(&m_rhs);
    destroyVector(&m_solution);
}

double * SystemSolver::getRHSRawPtr()
{
    PetscScalar *raw_rhs;
    VecGetArray(m_rhs, &raw_rhs);
 
    return raw_rhs;
}

const double * SystemSolver::getRHSRawPtr() const
{
    return getRHSRawReadPtr();
}

const double * SystemSolver::getRHSRawReadPtr() const
{
    const PetscScalar *raw_rhs;
    VecGetArrayRead(m_rhs, &raw_rhs);
 
    return raw_rhs;
}

void SystemSolver::restoreRHSRawPtr(double *raw_rhs)
{
    VecRestoreArray(m_rhs, &raw_rhs);
}

void SystemSolver::restoreRHSRawReadPtr(const double *raw_rhs) const
{
    VecRestoreArrayRead(m_rhs, &raw_rhs);
}

double * SystemSolver::getSolutionRawPtr()
{
    PetscScalar *raw_solution;
    VecGetArray(m_solution, &raw_solution);
 
    return raw_solution;
}

const double * SystemSolver::getSolutionRawPtr() const
{
    return getSolutionRawReadPtr();
}

const double * SystemSolver::getSolutionRawReadPtr() const
{
    const PetscScalar *raw_solution;
    VecGetArrayRead(m_solution, &raw_solution);
 
    return raw_solution;
}

void SystemSolver::restoreSolutionRawPtr(double *raw_solution)
{
    VecRestoreArray(m_solution, &raw_solution);
}

void SystemSolver::restoreSolutionRawReadPtr(const double *raw_solution) const
{
    VecRestoreArrayRead(m_solution, &raw_solution);
}

void SystemSolver::exportMatrix(const std::string &filePath, FileFormat fileFormat) const
{
    exportMatrix(m_A, filePath, fileFormat);
}

void SystemSolver::importMatrix(const std::string &filePath)
{
    // Check if the matrix exists
    if (!m_A) {
        throw std::runtime_error("Matrix should be created before filling it.");
    }
 
    // Clear workspace
    clearWorkspace();
 
    // Clear reordering
    clearReordering();
 
    // Fill the matrix
    matrixFill(filePath);
 
    // Re-create vectors
    vectorsDestroy();
    vectorsCreate();
}

void SystemSolver::exportRHS(const std::string &filePath, FileFormat fileFormat) const
{
    exportVector(m_rhs, filePath, fileFormat);
}

void SystemSolver::importRHS(const std::string &filePath)
{
    // Check if the system is assembled
    if (!m_assembled) {
        throw std::runtime_error("The RHS vector can be loaded only after assembling the system.");
    }
 
    // Fill the RHS vector
    vectorsFill(filePath, "");
 
    // Check if the imported RHS is compatible with the matrix
    PetscInt size;
    VecGetLocalSize(m_rhs, &size);
 
    PetscInt expectedSize;
    if (!m_transpose) {
        expectedSize = getColCount();
    } else {
        expectedSize = getRowCount();
    }
    expectedSize *= getBlockSize();
 
    if (size != expectedSize) {
        log::cout() << "The imported RHS vector is not compatible with the matrix" << std::endl;
        log::cout() << "The size of the imported RHS vector is " << size << std::endl;
        log::cout() << "The expected size of RHS vector is " << expectedSize << std::endl;
        throw std::runtime_error("The imported RHS vector is not compatible with the matrix");
    }
}

void SystemSolver::exportSolution(const std::string &filePath, FileFormat fileFormat) const
{
    exportVector(m_solution, filePath, fileFormat);
}

void SystemSolver::importSolution(const std::string &filePath)
{
    // Check if the system is assembled
    if (!m_assembled) {
        throw std::runtime_error("The solution vector can be loaded only after assembling the system.");
    }
 
    // Fill the solution vector
    vectorsFill("", filePath);
 
    // Check if the imported solution is compatible with the matrix
    PetscInt size;
    VecGetLocalSize(m_solution, &size);
 
    PetscInt expectedSize;
    if (!m_transpose) {
        expectedSize = getRowCount();
    } else {
        expectedSize = getColCount();
    }
    expectedSize *= getBlockSize();
 
    if (size != expectedSize) {
        log::cout() << "The imported solution vector is not compatible with the matrix" << std::endl;
        log::cout() << "The size of the imported solution vector is " << size << std::endl;
        log::cout() << "The expected size of solution vector is " << expectedSize << std::endl;
        throw std::runtime_error("The imported solution vector is not compatible with the matrix");
    }
}

void SystemSolver::dumpSystem(const std::string &header, const std::string &directory,
                              const std::string &prefix) const
{
    int partitioningBlock = getBinaryArchiveBlock();
 
    // Open stream that will contain system information
    OBinaryArchive systemArchive;
    if (partitioningBlock <= 0) {
        openBinaryArchive(header, directory, prefix + "info", -1, &systemArchive);
    }
    std::ostream &systemStream = systemArchive.getStream();
 
    // Dump system information
    dumpInfo(systemStream);
 
    // Dump matrix
    matrixDump(systemStream, directory, prefix);
 
    // Dump vectors
    vectorsDump(systemStream, directory, prefix);
 
    // Open stream with system information
    closeBinaryArchive(&systemArchive);
}

#if BITPIT_ENABLE_MPI==1
#endif
#if BITPIT_ENABLE_MPI==1
void SystemSolver::restoreSystem(MPI_Comm communicator, bool redistribute,
                                 const std::string &directory, const std::string &prefix)
#else
void SystemSolver::restoreSystem(const std::string &directory, const std::string &prefix)
#endif
{
    // Clear the system
    clear();
 
#if BITPIT_ENABLE_MPI == 1
    // Set the communicator
    setCommunicator(communicator);
#endif
 
    // Open stream with system information
    IBinaryArchive systemArchive;
    openBinaryArchive(directory, prefix + "info", -1, &systemArchive);
    std::istream &systemStream = systemArchive.getStream();
 
    // Dump system information
    restoreInfo(systemStream);
 
    // Restore the matrix
#if BITPIT_ENABLE_MPI==1
    matrixRestore(systemStream, directory, prefix, redistribute);
#else
    matrixRestore(systemStream, directory, prefix);
#endif
 
    // Restore RHS and solution vectors
    vectorsRestore(systemStream, directory, prefix);
 
    // Close stream with system information
    closeBinaryArchive(&systemArchive);
 
    // The system is now assembled
    m_assembled = true;
 
    // Initialize KSP options
    initializeKSPOptions();
 
    // Initialize KSP statuses
    initializeKSPStatus();
}

void SystemSolver::dumpInfo(std::ostream &systemStream) const
{
    if (systemStream.good()) {
#if BITPIT_ENABLE_MPI==1
        utils::binary::write(systemStream, m_partitioned);
#endif
        utils::binary::write(systemStream, m_transpose);
    }
}

void SystemSolver::restoreInfo(std::istream &systemStream)
{
#if BITPIT_ENABLE_MPI == 1
    // Detect if the system is partitioned
    utils::binary::read(systemStream, m_partitioned);
#endif
 
    // Set transpose flag
    bool transpose;
    utils::binary::read(systemStream, transpose);
    setTranspose(transpose);
 
}

void SystemSolver::createMatrix(int rowBlockSize, int colBlockSize, Mat *matrix) const
{
    // Create the matrix
#if BITPIT_ENABLE_MPI == 1
    MatCreate(m_communicator, matrix);
#else
    MatCreate(PETSC_COMM_SELF, matrix);
#endif
 
    // Set matrix type
    bool creatBlockMatrix = false;
    if (m_flatten) {
        creatBlockMatrix = false;
    } else {
        creatBlockMatrix = (rowBlockSize == colBlockSize) && (rowBlockSize != 1);
    }
 
#if BITPIT_ENABLE_MPI == 1
    if (m_partitioned) {
        if (creatBlockMatrix) {
            MatSetType(*matrix, MATMPIBAIJ);
        } else {
            MatSetType(*matrix, MATMPIAIJ);
        }
    } else
#endif
    {
        if (creatBlockMatrix) {
            MatSetType(*matrix, MATSEQBAIJ);
        } else {
            MatSetType(*matrix, MATSEQAIJ);
        }
    }
 
    // Set block size
    if (rowBlockSize == colBlockSize) {
        MatSetBlockSize(*matrix, rowBlockSize);
    } else if (rowBlockSize != 1 || colBlockSize != 1) {
        MatSetBlockSizes(*matrix, rowBlockSize, colBlockSize);
    }
}

void SystemSolver::createMatrix(int rowBlockSize, int colBlockSize, int nNestRows, int nNestCols,
                                Mat *subMatrices, Mat *matrix) const
{
    // Create the matrix
#if BITPIT_ENABLE_MPI == 1
    MatCreateNest(getCommunicator(), nNestRows, PETSC_NULLPTR, nNestCols, PETSC_NULLPTR, subMatrices, matrix);
#else
    MatCreateNest(PETSC_COMM_SELF, nNestRows, PETSC_NULLPTR, nNestCols, PETSC_NULLPTR, subMatrices, matrix);
#endif
 
    // Set block size
    if (rowBlockSize == colBlockSize) {
        MatSetBlockSize(*matrix, rowBlockSize);
    } else if (rowBlockSize != 1 || colBlockSize != 1) {
        MatSetBlockSizes(*matrix, rowBlockSize, colBlockSize);
    }
}

void SystemSolver::fillMatrix(Mat matrix, const std::string &filePath) const
{
    // Check if the file exists
    std::ifstream fileStream(filePath.c_str());
    if (!fileStream.good()) {
        throw std::runtime_error("The PETSc matrix file \"" + filePath + "\" doesn't exists.");
    }
    fileStream.close();
 
    // Fill matrix content
    PetscViewer viewer;
#if BITPIT_ENABLE_MPI==1
    PetscViewerCreate(m_communicator, &viewer);
#else
    PetscViewerCreate(PETSC_COMM_SELF, &viewer);
#endif
    PetscViewerSetType(viewer, PETSCVIEWERBINARY);
    PetscViewerFileSetMode(viewer, FILE_MODE_READ);
    PetscViewerFileSetName(viewer, filePath.c_str());
    MatLoad(matrix, viewer);
    PetscViewerDestroy(&viewer);
}

#if BITPIT_ENABLE_MPI==1
#endif
void SystemSolver::dumpMatrix(Mat matrix, const std::string &directory, const std::string &name) const
{
    int partitioningBlock = getBinaryArchiveBlock();
 
    // Store information needed to create the matrix
    OBinaryArchive infoArchive;
    if (partitioningBlock <= 0) {
        openBinaryArchive("", directory, name + ".info", -1, &infoArchive);
        std::ostream &infoStream = infoArchive.getStream();
 
        bool matrixExists = matrix;
        utils::binary::write(infoStream, matrixExists);
        if (!matrixExists) {
            closeBinaryArchive(&infoArchive);
            return;
        }
 
        PetscInt rowBlockSize;
        PetscInt colBlockSize;
        MatGetBlockSizes(matrix, &rowBlockSize, &colBlockSize);
        utils::binary::write(infoStream, static_cast<int>(rowBlockSize));
        utils::binary::write(infoStream, static_cast<int>(colBlockSize));
 
        closeBinaryArchive(&infoArchive);
    }
 
    if (!matrix) {
        return;
    }
 
#if BITPIT_ENABLE_MPI==1
    // Store partitioning information
    if (isPartitioned()) {
        OBinaryArchive partitioningArchive;
        openBinaryArchive("", directory, name + ".partitioning", partitioningBlock, &partitioningArchive);
        std::ostream &partitioningStream = partitioningArchive.getStream();
 
        PetscInt nLocalRows;
        PetscInt nLocalCols;
        MatGetLocalSize(matrix, &nLocalRows, &nLocalCols);
        utils::binary::write(partitioningStream, static_cast<std::size_t>(nLocalRows));
        utils::binary::write(partitioningStream, static_cast<std::size_t>(nLocalCols));
 
        PetscInt nGlobalRows;
        PetscInt nGlobalCols;
        MatGetSize(matrix, &nGlobalRows, &nGlobalCols);
        utils::binary::write(partitioningStream, static_cast<std::size_t>(nGlobalRows));
        utils::binary::write(partitioningStream, static_cast<std::size_t>(nGlobalCols));
 
        closeBinaryArchive(&partitioningArchive);
    }
#endif
 
    // Store matrix data
    std::string filePath = getDataFilePath(directory, name);
    exportMatrix(matrix, filePath, FILE_BINARY);
}

#if BITPIT_ENABLE_MPI==1
#endif
#if BITPIT_ENABLE_MPI==1
void SystemSolver::restoreMatrix(const std::string &directory, const std::string &name,
                                 bool redistribute, Mat *matrix) const
{
#else
void SystemSolver::restoreMatrix(const std::string &directory, const std::string &name,
                                 Mat *matrix) const
{
#endif
 
    // Create matrix
    IBinaryArchive infoArchive;
    openBinaryArchive(directory, name + ".info", -1, &infoArchive);
    std::istream &infoStream = infoArchive.getStream();
 
    bool matrixExists;
    utils::binary::read(infoStream, matrixExists);
    if (!matrixExists) {
        *matrix = PETSC_NULLPTR;
        closeBinaryArchive(&infoArchive);
        return;
    }
 
    int rowBlockSize;
    int colBlockSize;
    utils::binary::read(infoStream, rowBlockSize);
    utils::binary::read(infoStream, colBlockSize);
    createMatrix(rowBlockSize, colBlockSize, matrix);
 
    closeBinaryArchive(&infoArchive);
 
#if BITPIT_ENABLE_MPI==1
    // Set partitioning information
    if (isPartitioned() && !redistribute) {
        int partitioningBlock = getBinaryArchiveBlock();
 
        IBinaryArchive partitioningArchive;
        openBinaryArchive(directory, name + ".partitioning", partitioningBlock, &partitioningArchive);
        std::istream &partitioningStream = partitioningArchive.getStream();
 
        std::size_t nLocalRows;
        std::size_t nLocalCols;
        utils::binary::read(partitioningStream, nLocalRows);
        utils::binary::read(partitioningStream, nLocalCols);
        std::size_t nGlobalRows;
        std::size_t nGlobalCols;
        utils::binary::read(partitioningStream, nGlobalRows);
        utils::binary::read(partitioningStream, nGlobalCols);
        MatSetSizes(*matrix, nLocalRows, nLocalCols, nGlobalRows, nGlobalCols);
 
        closeBinaryArchive(&partitioningArchive);
    }
#endif
 
    // Fill matrix
    std::string filePath = getDataFilePath(directory, name);
    fillMatrix(*matrix, filePath);
}

void SystemSolver::exportMatrix(Mat matrix, const std::string &filePath, FileFormat fileFormat) const
{
    // Early return if the matrix doesn't exist
    bool matrixExists = matrix;
    if (!matrixExists) {
        std::ofstream dataFile(filePath);
        dataFile.close();
 
        std::ofstream infoFile(filePath + ".info");
        infoFile.close();
 
        return;
    }
 
    // Create the viewer
    PetscViewerType viewerType;
    PetscViewerFormat viewerFormat;
    if (fileFormat == FILE_BINARY) {
        viewerType   = PETSCVIEWERBINARY;
        viewerFormat = PETSC_VIEWER_DEFAULT;
    } else {
        viewerType   = PETSCVIEWERASCII;
        viewerFormat = PETSC_VIEWER_ASCII_MATLAB;
    }
 
    PetscViewer matViewer;
#if BITPIT_ENABLE_MPI==1
    PetscViewerCreate(m_communicator, &matViewer);
#else
    PetscViewerCreate(PETSC_COMM_SELF, &matViewer);
#endif
    PetscViewerSetType(matViewer, viewerType);
    PetscViewerFileSetMode(matViewer, FILE_MODE_WRITE);
    PetscViewerPushFormat(matViewer, viewerFormat);
    PetscViewerFileSetName(matViewer, filePath.c_str());
 
    // Export the matrix
    MatView(matrix, matViewer);
 
    // Destroy the viewer
    PetscViewerDestroy(&matViewer);
}

void SystemSolver::destroyMatrix(Mat *matrix) const
{
    if (matrix) {
        MatDestroy(matrix);
        *matrix = PETSC_NULLPTR;
    }
}

void SystemSolver::createVector(int blockSize, Vec *vector) const
{
    // Create the vector
#if BITPIT_ENABLE_MPI == 1
    VecCreate(m_communicator, vector);
#else
    VecCreate(PETSC_COMM_SELF, vector);
#endif
 
    // Set vector type
#if BITPIT_ENABLE_MPI == 1
    if (m_partitioned) {
        VecSetType(*vector, VECMPI);
    } else
#endif
    {
        VecSetType(*vector, VECSEQ);
    }
 
    // Set block size
    if (blockSize != 1) {
        VecSetBlockSize(*vector, blockSize);
    }
}

void SystemSolver::createVector(int blockSize, int nestSize, Vec *subVectors, Vec *vector) const
{
    // Create the vector
#if BITPIT_ENABLE_MPI == 1
    VecCreateNest(getCommunicator(), nestSize, PETSC_NULLPTR, subVectors, vector);
#else
    VecCreateNest(PETSC_COMM_SELF, nestSize, PETSC_NULLPTR, subVectors, vector);
#endif
 
    // Set block size
    if (blockSize != 1) {
        VecSetBlockSize(*vector, blockSize);
    }
}

void SystemSolver::fillVector(Vec vector, const std::string &filePath) const
{
    // Check if the file exists
    std::ifstream fileStream(filePath);
    if (!fileStream.good()) {
        throw std::runtime_error("The file \"" + filePath + "\" cannot be read.");
    }
    fileStream.close();
 
    // Fill the vector
    PetscViewer viewer;
#if BITPIT_ENABLE_MPI==1
    PetscViewerCreate(m_communicator, &viewer);
#else
    PetscViewerCreate(PETSC_COMM_SELF, &viewer);
#endif
    PetscViewerSetType(viewer, PETSCVIEWERBINARY);
    PetscViewerFileSetMode(viewer, FILE_MODE_READ);
    PetscViewerFileSetName(viewer, filePath.c_str());
    VecLoad(vector, viewer);
    PetscViewerDestroy(&viewer);
}

void SystemSolver::fillVector(Vec vector, const std::vector<double> &data) const
{
    int size;
    VecGetLocalSize(vector, &size);
 
    PetscScalar *petscData;
    VecGetArray(vector, &petscData);
    for (int i = 0; i < size; ++i) {
        petscData[i] = data[i];
    }
    VecRestoreArray(vector, &petscData);
}

#if BITPIT_ENABLE_MPI==1
#endif
void SystemSolver::dumpVector(Vec vector, const std::string &directory, const std::string &name) const
{
    int partitioningBlock = getBinaryArchiveBlock();
 
    // Store information needed to create the matrix
    if (partitioningBlock <= 0) {
        OBinaryArchive infoArchive;
        openBinaryArchive("", directory, name + ".info", -1, &infoArchive);
        std::ostream &infoStream = infoArchive.getStream();
 
        bool vectorExists = vector;
        utils::binary::write(infoStream, vectorExists);
        if (!vectorExists) {
            closeBinaryArchive(&infoArchive);
            return;
        }
 
        PetscInt blockSize;
        VecGetBlockSize(vector, &blockSize);
        utils::binary::write(infoStream, static_cast<int>(blockSize));
 
        closeBinaryArchive(&infoArchive);
    }
 
    if (!vector) {
        return;
    }
 
#if BITPIT_ENABLE_MPI==1
    // Store partitioning information
    if (isPartitioned()) {
        int partitioningBlock = getBinaryArchiveBlock();
 
        OBinaryArchive partitioningArchive;
        openBinaryArchive("", directory, name + ".partitioning", partitioningBlock, &partitioningArchive);
        std::ostream &partitioningStream = partitioningArchive.getStream();
 
        PetscInt localSize;
        VecGetLocalSize(vector, &localSize);
        utils::binary::write(partitioningStream, static_cast<std::size_t>(localSize));
 
        PetscInt globalSize;
        VecGetSize(vector, &globalSize);
        utils::binary::write(partitioningStream, static_cast<std::size_t>(globalSize));
 
        closeBinaryArchive(&partitioningArchive);
    }
#endif
 
    // Store vector data
    std::string filePath = getDataFilePath(directory, name);
    exportVector(vector, filePath, FILE_BINARY);
}

void SystemSolver::restoreVector(const std::string &directory, const std::string &name,
                                 Mat matrix, VectorSide side, Vec *vector) const
{
    // Get size information
    PetscInt nLocalRows;
    PetscInt nLocalCols;
    MatGetLocalSize(matrix, &nLocalRows, &nLocalCols);
 
    PetscInt nGlobalRows;
    PetscInt nGlobalCols;
    MatGetSize(matrix, &nGlobalRows, &nGlobalCols);
 
    std::size_t localSize  = std::numeric_limits<std::size_t>::max();
    std::size_t globalSize = std::numeric_limits<std::size_t>::max();
    if (!m_transpose) {
        if (side == VECTOR_SIDE_RIGHT) {
            localSize  = nLocalCols;
            globalSize = nGlobalCols;
        } else if (side == VECTOR_SIDE_LEFT) {
            localSize  = nLocalRows;
            globalSize = nGlobalRows;
        }
    } else {
        if (side == VECTOR_SIDE_RIGHT) {
            localSize  = nLocalRows;
            globalSize = nGlobalRows;
        } else if (side == VECTOR_SIDE_LEFT) {
            localSize  = nLocalCols;
            globalSize = nGlobalCols;
        }
    }
 
    // Restore the vector
    restoreVector(directory, name, localSize, globalSize, vector);
}

#if BITPIT_ENABLE_MPI==1
#endif
#if BITPIT_ENABLE_MPI==1
void SystemSolver::restoreVector(const std::string &directory, const std::string &name,
                                 bool redistribute, Vec *vector) const
#else
void SystemSolver::restoreVector(const std::string &directory, const std::string &name,
                                 Vec *vector) const
#endif
{
 
    // Get size information
    std::size_t localSize  = std::numeric_limits<std::size_t>::max();
    std::size_t globalSize = std::numeric_limits<std::size_t>::max();
#if BITPIT_ENABLE_MPI==1
    if (isPartitioned()) {
        if (!redistribute) {
            int partitioningBlock = getBinaryArchiveBlock();
 
            IBinaryArchive partitioningArchive;
            openBinaryArchive(directory, name + ".partitioning", partitioningBlock, &partitioningArchive);
            std::istream &partitioningStream = partitioningArchive.getStream();
 
            utils::binary::read(partitioningStream, localSize);
            utils::binary::read(partitioningStream, globalSize);
 
            closeBinaryArchive(&partitioningArchive);
        }
    }
#endif
 
    restoreVector(directory, name, localSize, globalSize, vector);
}

void SystemSolver::restoreVector(const std::string &directory, const std::string &name,
                                 std::size_t localSize, std::size_t globalSize, Vec *vector) const
{
    // Create vector
    IBinaryArchive infoArchive;
    openBinaryArchive(directory, name + ".info", -1, &infoArchive);
    std::istream &infoStream = infoArchive.getStream();
 
    bool vectorExists;;
    utils::binary::read(infoStream, vectorExists);
    if (!vectorExists) {
        *vector = PETSC_NULLPTR;
        closeBinaryArchive(&infoArchive);
        return;
    }
 
    int blockSize;
    utils::binary::read(infoStream, blockSize);
    createVector(blockSize, vector);
 
    closeBinaryArchive(&infoArchive);
 
    // Set size information
    PetscInt petscLocalSize;
    if (localSize != std::numeric_limits<std::size_t>::max()) {
        petscLocalSize = static_cast<PetscInt>(localSize);
    } else {
        petscLocalSize = PETSC_DECIDE;
    }
 
    PetscInt petscGlobalSize;
    if (globalSize != std::numeric_limits<std::size_t>::max()) {
        petscGlobalSize = static_cast<PetscInt>(globalSize);
    } else {
        petscGlobalSize = PETSC_DETERMINE;
    }
 
    if (petscLocalSize != PETSC_DECIDE || petscGlobalSize != PETSC_DETERMINE) {
        VecSetSizes(*vector, petscLocalSize, petscGlobalSize);
    }
 
    // Fill vector data
    std::string filePath = getDataFilePath(directory, name);
    fillVector(*vector, filePath);
}

void SystemSolver::reorderVector(Vec vector, IS permutations, bool invert) const
{
    if (!permutations) {
        return;
    }
 
    PetscBool petscInvert;
    if (invert) {
        petscInvert = PETSC_TRUE;
    } else {
        petscInvert = PETSC_FALSE;
    }
 
    VecPermute(vector, permutations, petscInvert);
}

void SystemSolver::exportVector(Vec vector, const std::string &filePath, FileFormat fileFormat) const
{
    // Early return if the matrix doesn't exist
    bool vectorExists = vector;
    if (!vectorExists) {
        std::ofstream dataFile(filePath);
        dataFile.close();
 
        std::ofstream infoFile(filePath + ".info");
        infoFile.close();
 
        return;
    }
 
    // Create the viewer
    PetscViewerType viewerType;
    PetscViewerFormat viewerFormat;
    if (fileFormat == FILE_BINARY) {
        viewerType   = PETSCVIEWERBINARY;
        viewerFormat = PETSC_VIEWER_DEFAULT;
    } else {
        viewerType   = PETSCVIEWERASCII;
        viewerFormat = PETSC_VIEWER_ASCII_MATLAB;
    }
 
    PetscViewer viewer;
#if BITPIT_ENABLE_MPI==1
    PetscViewerCreate(m_communicator, &viewer);
#else
    PetscViewerCreate(PETSC_COMM_SELF, &viewer);
#endif
    PetscViewerSetType(viewer, viewerType);
    PetscViewerFileSetMode(viewer, FILE_MODE_WRITE);
    PetscViewerPushFormat(viewer, viewerFormat);
    PetscViewerFileSetName(viewer, filePath.c_str());
 
    // Export the vector
    VecView(vector, viewer);
 
    // Destroy the viewer
    PetscViewerDestroy(&viewer);
}

void SystemSolver::exportVector(Vec vector, std::vector<double> *data) const
{
    int size;
    VecGetLocalSize(m_solution, &size);
 
    const PetscScalar *vectorData;
    VecGetArrayRead(m_solution, &vectorData);
    for (int i = 0; i < size; ++i) {
        (*data)[i] = vectorData[i];
    }
    VecRestoreArrayRead(vector, &vectorData);
}

void SystemSolver::destroyVector(Vec *vector) const
{
    if (vector) {
        VecDestroy(vector);
        *vector = PETSC_NULLPTR;
    }
}

std::string SystemSolver::getInfoFilePath(const std::string &directory, const std::string &name) const
{
    return getFilePath(directory, name + ".info");
}

std::string SystemSolver::getDataFilePath(const std::string &directory, const std::string &name) const
{
    return getFilePath(directory, name + ".dat");
}

std::string SystemSolver::getFilePath(const std::string &directory, const std::string &name) const
{
    std::string path = directory + "/" + name;
 
    return path;
}

void SystemSolver::setNullSpace()
{
    MatNullSpace nullspace;
#if BITPIT_ENABLE_MPI==1
    MatNullSpaceCreate(m_communicator, PETSC_TRUE, 0, NULL, &nullspace);
#else
    MatNullSpaceCreate(PETSC_COMM_SELF, PETSC_TRUE, 0, NULL, &nullspace);
#endif
    MatSetNullSpace(m_A, nullspace);
    MatNullSpaceDestroy(&nullspace);
}

void SystemSolver::unsetNullSpace()
{
    MatSetNullSpace(m_A, NULL);
}

void SystemSolver::setReordering(long nRows, long nCols, const SystemMatrixOrdering &reordering)
{
    // Clear existing reordering
    clearReordering();
 
    // Early return if natural ordering is used
    try {
        dynamic_cast<const NaturalSystemMatrixOrdering &>(reordering);
        return;
    } catch(const std::bad_cast &exception) {
        BITPIT_UNUSED(exception);
 
        // A reordering other than the natural one has been passed
    }
 
    // Set row reordering
    PetscInt *rowReorderingStorage;
    PetscMalloc(nRows * sizeof(PetscInt), &rowReorderingStorage);
    for (long i = 0; i < nRows; ++i) {
        rowReorderingStorage[reordering.getRowPermutationRank(i)] = i;
    }
 
#if BITPIT_ENABLE_MPI == 1
    ISCreateGeneral(m_communicator, nRows, rowReorderingStorage, PETSC_OWN_POINTER, &m_rowReordering);
#else
    ISCreateGeneral(PETSC_COMM_SELF, nRows, rowReorderingStorage, PETSC_OWN_POINTER, &m_rowReordering);
#endif
    ISSetPermutation(m_rowReordering);
 
    // Create new permutations
    PetscInt *colReorderingStorage;
    PetscMalloc(nCols * sizeof(PetscInt), &colReorderingStorage);
    for (long j = 0; j < nCols; ++j) {
        colReorderingStorage[reordering.getColPermutationRank(j)] = j;
    }
 
#if BITPIT_ENABLE_MPI == 1
    ISCreateGeneral(m_communicator, nCols, colReorderingStorage, PETSC_OWN_POINTER, &m_colReordering);
#else
    ISCreateGeneral(PETSC_COMM_SELF, nCols, colReorderingStorage, PETSC_OWN_POINTER, &m_colReordering);
#endif
    ISSetPermutation(m_colReordering);
}

void SystemSolver::clearReordering()
{
    if (m_rowReordering) {
        ISDestroy(&m_rowReordering);
    }
 
    if (m_colReordering) {
        ISDestroy(&m_colReordering);
    }
}

void SystemSolver::prepareKSP()
{
    // Check if the system is assembled
    if (!isAssembled()) {
        throw std::runtime_error("Unable to solve the system. The system is not yet assembled.");
    }
 
    // Early return if the KSP can be reused
    if (!m_KSPDirty) {
        return;
    }
 
    // Create KSP
    bool setupNeeded = false;
    if (!m_KSP) {
        createKSP();
        setupNeeded = true;
    }
 
    // Set the matrix associated with the linear system
    KSPSetOperators(m_KSP, m_A, m_A);
 
    // Set up
    if (setupNeeded) {
        // Perform actions before preconditioner set up
        prePreconditionerSetupActions();
 
        // Initialize preconditioner from options
        PC pc;
        KSPGetPC(m_KSP, &pc);
        PCSetFromOptions(pc);
 
        // Set up preconditioner
        setupPreconditioner();
 
        // Perform actions after preconditioner set up
        postPreconditionerSetupActions();
 
        // Perform actions before Krylov subspace method set up set up
        preKrylovSetupActions();
 
        // Initialize Krylov subspace from options
        KSPSetFromOptions(m_KSP);
 
        // Set up the Krylov subspace method
        setupKrylov();
 
        // Perform actions after Krylov subspace method set up
        postKrylovSetupActions();
    }
 
    // KSP is now ready
    m_KSPDirty = false;
 
    // Reset KSP status
    resetKSPStatus();
}

void SystemSolver::finalizeKSP()
{
    // Nothing to do
}

void SystemSolver::createKSP()
{
    // Create KSP object
#if BITPIT_ENABLE_MPI==1
    KSPCreate(m_communicator, &m_KSP);
#else
    KSPCreate(PETSC_COMM_SELF, &m_KSP);
#endif
 
    // Set options prefix
    if (!m_prefix.empty()) {
        KSPSetOptionsPrefix(m_KSP, m_prefix.c_str());
    }
}

void SystemSolver::destroyKSP()
{
    m_KSPDirty = true;
    KSPDestroy(&m_KSP);
    m_KSP = PETSC_NULLPTR;
}

void SystemSolver::setupPreconditioner()
{
    PC pc;
    KSPGetPC(m_KSP, &pc);
    setupPreconditioner(pc, getKSPOptions());
}

void SystemSolver::setupPreconditioner(PC pc, const KSPOptions &options) const
{
    // Set preconditioner type
    PCType pcType;
#if BITPIT_ENABLE_MPI == 1
    if (isPartitioned()) {
        pcType = PCASM;
    } else {
        pcType = PCILU;
    }
#else
    pcType = PCILU;
#endif
 
    PCSetType(pc, pcType);
 
    // Configure preconditioner
    if (strcmp(pcType, PCASM) == 0) {
        if (options.overlap != PETSC_DEFAULT) {
            PCASMSetOverlap(pc, options.overlap);
        }
    } else if (strcmp(pcType, PCILU) == 0) {
        if (options.levels != PETSC_DEFAULT) {
            PCFactorSetLevels(pc, options.levels);
        }
    }
 
    PCSetUp(pc);
 
    if (strcmp(pcType, PCASM) == 0) {
        KSP *subKSPs;
        PetscInt nSubKSPs;
        PCASMGetSubKSP(pc, &nSubKSPs, PETSC_NULLPTR, &subKSPs);
 
        for (PetscInt i = 0; i < nSubKSPs; ++i) {
            KSP subKSP = subKSPs[i];
            KSPSetType(subKSP, KSPPREONLY);
 
            PC subPC;
            KSPGetPC(subKSP, &subPC);
            PCSetType(subPC, PCILU);
            if (options.levels != PETSC_DEFAULT) {
                PCFactorSetLevels(subPC, options.levels);
            }
        }
    }
}

void SystemSolver::prePreconditionerSetupActions()
{
    // Nothing to do
}

void SystemSolver::postPreconditionerSetupActions()
{
    // Nothing to do
}

void SystemSolver::setupKrylov()
{
    setupKrylov(m_KSP, getKSPOptions());
}

void SystemSolver::setupKrylov(KSP ksp, const KSPOptions &options) const
{
    KSPSetType(ksp, KSPFGMRES);
    if (options.restart != PETSC_DEFAULT) {
        KSPGMRESSetRestart(ksp, options.restart);
    }
    if (options.rtol != PETSC_DEFAULT || options.atol != PETSC_DEFAULT || options.maxits != PETSC_DEFAULT) {
        KSPSetTolerances(ksp, options.rtol, options.atol, PETSC_DEFAULT, options.maxits);
    }
    KSPSetInitialGuessNonzero(ksp, options.initial_non_zero);
}

void SystemSolver::preKrylovSetupActions()
{
    // Enable convergence monitor
    if (m_convergenceMonitorEnabled) {
#if (PETSC_VERSION_MAJOR >= 3 && PETSC_VERSION_MINOR >= 7)
        PetscOptionsSetValue(PETSC_NULLPTR, ("-" + m_prefix + "ksp_monitor_true_residual").c_str(), "");
        PetscOptionsSetValue(PETSC_NULLPTR, ("-" + m_prefix + "ksp_monitor_singular_value").c_str(), "");
        PetscOptionsSetValue(PETSC_NULLPTR, ("-" + m_prefix + "ksp_converged_reason").c_str(), "");
#else
        PetscOptionsSetValue(("-" + m_prefix + "ksp_monitor_true_residual").c_str(), "");
        PetscOptionsSetValue(("-" + m_prefix + "ksp_monitor_singular_value").c_str(), "");
        PetscOptionsSetValue(("-" + m_prefix + "ksp_converged_reason").c_str(), "");
#endif
    }
}

void SystemSolver::postKrylovSetupActions()
{
}

KSPOptions & SystemSolver::getKSPOptions()
{
    if (!isAssembled()) {
        throw std::runtime_error("The options associated with the Krylov solver can only be accessed after assembling the system.");
    }
 
    return m_KSPOptions;
}

const KSPOptions & SystemSolver::getKSPOptions() const
{
    if (!isAssembled()) {
        throw std::runtime_error("The options associated with the Krylov solver can only be accessed after assembling the system.");
    }
 
    return m_KSPOptions;
}

void SystemSolver::initializeKSPOptions()
{
    resetKSPOptions(&m_KSPOptions);
}

void SystemSolver::resetKSPOptions(KSPOptions *options) const
{
    *options = KSPOptions();
}

void SystemSolver::destroyKSPOptions()
{
    resetKSPOptions(&m_KSPOptions);
}

const KSPStatus & SystemSolver::getKSPStatus() const
{
    if (!isAssembled()) {
        throw std::runtime_error("The status of the Krylov solver can only be accessed after assembling the system.");
    }
 
    return m_KSPStatus;
}

void SystemSolver::initializeKSPStatus()
{
    resetKSPStatus();
}

void SystemSolver::fillKSPStatus()
{
    fillKSPStatus(m_KSP, &m_KSPStatus);
}

void SystemSolver::fillKSPStatus(KSP ksp, KSPStatus *status) const
{
    status->error = 0;
    KSPGetIterationNumber(ksp, &(status->its));
    KSPGetConvergedReason(ksp, &(status->convergence));
}

void SystemSolver::resetKSPStatus()
{
    resetKSPStatus(&m_KSPStatus);
}

void SystemSolver::resetKSPStatus(KSPStatus *status) const
{
    status->error       = 0;
    status->its         = -1;
    status->convergence = KSP_CONVERGED_ITERATING;
}

void SystemSolver::destroyKSPStatus()
{
    resetKSPStatus();
}
 
#if BITPIT_ENABLE_MPI==1
const MPI_Comm & SystemSolver::getCommunicator() const
{
    return m_communicator;
}

void SystemSolver::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 SystemSolver::freeCommunicator()
{
    if (m_communicator != MPI_COMM_SELF) {
        int finalizedCalled;
        MPI_Finalized(&finalizedCalled);
        if (!finalizedCalled) {
            MPI_Comm_free(&m_communicator);
        }
    }
}
#endif

void SystemSolver::removeNullSpaceFromRHS()
{
    // Get the null space from the matrix
    MatNullSpace nullspace;
    MatGetNullSpace(m_A, &nullspace);
 
    // Remove null space components from right hand side
    MatNullSpaceRemove(nullspace, m_rhs);
}

bool SystemSolver::isForceConsistencyEnabled() const
{
    return m_forceConsistency;
}

void SystemSolver::enableForceConsistency(bool enable)
{
    m_forceConsistency = enable;
}

void SystemSolver::openBinaryArchive(const std::string &directory, const std::string &prefix,
                                     int block, IBinaryArchive *archive) const
{
    int expectedVersion = getDumpVersion();
 
    std::string path = getFilePath(directory, prefix);
    archive->open(path, "dat", block);
    bool versionsMatch = archive->checkVersion(expectedVersion);
    if (!versionsMatch) {
        std::string message = "Version stored in binary archive " + path + " does not match";
        message += " the expected version " + std::to_string(expectedVersion) + ".";
 
        throw std::runtime_error(message);
    }
}

void SystemSolver::openBinaryArchive(const std::string &header, const std::string &directory,
                                     const std::string &prefix, int block, OBinaryArchive *archive) const
{
    int version = getDumpVersion();
    std::string path = getFilePath(directory, prefix);
    archive->open(path, "dat", version, header, block);
}

void SystemSolver::closeBinaryArchive(BinaryArchive *archive) const
{
    archive->close();
}

int SystemSolver::getBinaryArchiveBlock() const
{
#if BITPIT_ENABLE_MPI == 1
    if (isPartitioned()) {
        int nProcesses;
        MPI_Comm_size(getCommunicator(), &nProcesses);
        if (nProcesses <= 1) {
            return -1;
        }
 
        int archiveBlock;
        MPI_Comm_rank(getCommunicator(), &archiveBlock);
 
        return archiveBlock;
    } else {
        return -1;
    }
#else
    return -1;
#endif
}
 
}