mpi_projection.h

Go to the documentation of this file.

#pragma once
 
#include "dg/backend/typedefs.h"
#include "dg/backend/index.h"
#include "dg/backend/mpi_matrix.h"
#include "dg/backend/mpi_gather.h"
#include "mpi_grid.h"
#include "projection.h"
 
namespace dg
{
 
template<class ConversionPolicy, class real_type>
dg::MIHMatrix_t<real_type> make_mpi_matrix(
    const dg::IHMatrix_t<real_type>& global_cols, const ConversionPolicy& col_policy)
{
    // Our approach here is:
    // 1. Convert to coo format
    // 2. mark rows in global_cols that are communicating
    // we need to grab the entire row to ensure reproducibility (which is guaranteed by order of computations)
    // 3. Convert inner rows to local and move outer rows to vectors (invalidating them in inner)
    // 4. Use global columns vector to construct MPIGather
    // 5. Use buffer index to construct local CooMat
    int size;
    MPI_Comm_size( col_policy.communicator(), &size);
    if( size == 1)
    {
        return { global_cols,
                 dg::SparseMatrix<int, real_type, thrust::host_vector>{},
                 MPIGather<thrust::host_vector>{},
                 thrust::host_vector<int>{} };
    }
    int rank;
    MPI_Comm_rank( col_policy.communicator(), &rank);
 
    const dg::IHMatrix_t<real_type>& A = global_cols;
    std::vector<int> global_row( A.num_rows(), 0);
    thrust::host_vector<int> pids(A.num_nnz()), lIdx(A.num_nnz());
    // 1st pass determine local rows
    for(int i = 0; i < (int)A.num_rows(); i++)
    for (int jj = A.row_offsets()[i]; jj < A.row_offsets()[i+1]; jj++)
    {
        assert( col_policy.global2localIdx( A.column_indices()[jj], lIdx[jj], pids[jj]));
        if( pids[jj] != rank)
        {
            global_row[i] = 1;
        }
    }
 
    // 2nd pass: distribute entries
    thrust::host_vector<std::array<int,2>> outer_col;
    thrust::host_vector<int> inner_row, inner_col, outer_row;
    thrust::host_vector<real_type> inner_val, outer_val;
    thrust::host_vector<int> row_scatter;
    inner_row.push_back(0);
    outer_row.push_back(0);
    for(int i = 0; i < (int)A.num_rows(); i++)
    {
        for (int jj = A.row_offsets()[i]; jj < A.row_offsets()[i+1]; jj++)
        {
            if( global_row[i] == 1)
            {
                outer_col.push_back( {pids[jj],lIdx[jj]});
                outer_val.push_back( A.values()[jj]);
            }
            else
            {
                inner_col.push_back( lIdx[jj]);
                inner_val.push_back( A.values()[jj]);
            }
        }
        if( global_row[i] == 1)
        {
            row_scatter.push_back(i);
            int old_end = outer_row.back();
            outer_row.push_back( old_end + A.row_offsets()[i+1] - A.row_offsets()[i]);
            inner_row.push_back( inner_row[i]);
        }
        else
        {
            inner_row.push_back( inner_row[i] + A.row_offsets()[i+1] - A.row_offsets()[i]);
        }
    }
    // 3. Now make MPI Gather object
    dg::SparseMatrix<int, real_type, thrust::host_vector> inner( A.num_rows(),
            col_policy.local_size(), inner_row, inner_col, inner_val);
 
    thrust::host_vector<int> lColIdx;
    auto gather_map = dg::gIdx2unique_idx( outer_col, lColIdx);
    MPIGather<thrust::host_vector> mpi_gather( gather_map,
            col_policy.communicator());
 
    dg::SparseMatrix<int, real_type, thrust::host_vector> outer( outer_row.size()-1,
            mpi_gather.buffer_size(), outer_row, lColIdx, outer_val);
 
    return { inner, outer, mpi_gather, row_scatter};
}
 
 
template<class ConversionPolicy, class real_type>
dg::IHMatrix_t<real_type> convertGlobal2LocalRows( const
    dg::IHMatrix_t<real_type>& global, const ConversionPolicy& row_policy)
{
    int size;
    MPI_Comm_size( row_policy.communicator(), &size);
    if( size == 1)
    {
        return global; // no conversion necessary
    }
    // 0. Convert to coo matrix
    thrust::host_vector<int> global_row_indices = dg::detail::csr2coo( global.row_offsets());
 
    // 1. For all rows determine pid to which it belongs
    thrust::host_vector<int> pids = dg::gIdx2gIdx( global_row_indices, row_policy);
    std::map<int, thrust::host_vector<int>> rows, cols;
    std::map<int, thrust::host_vector<real_type>> vals;
    for( unsigned u=0; u<global_row_indices.size(); u++)
    {
        rows[pids[u]].push_back( global_row_indices[u]);
        cols[pids[u]].push_back( global.column_indices()[u]);
        vals[pids[u]].push_back( global.values()[u]);
    }
    // 2. Now send those rows to where they belong
    auto row_buf = dg::mpi_permute( rows, row_policy.communicator());
    auto col_buf = dg::mpi_permute( cols, row_policy.communicator());
    auto val_buf = dg::mpi_permute( vals, row_policy.communicator());
 
    dg::SparseMatrix<int, real_type, thrust::host_vector> B;
    // indices come in any order ( so we need to sort)
    B.setFromCoo( row_policy.local_size(), global.num_cols(),
        dg::detail::flatten_map( row_buf),
        dg::detail::flatten_map(col_buf),
        dg::detail::flatten_map(val_buf), true);
    return B;
    // 4. Reduce on identical rows/cols
    // .... Maybe later
}
 
template<class ConversionPolicy, class real_type>
void convertLocal2GlobalCols( dg::IHMatrix_t<real_type>& local, const ConversionPolicy& policy)
{
    // 1. For all columns determine pid to which it belongs
    int rank=0;
    MPI_Comm_rank( policy.communicator(), &rank);
    thrust::host_vector<int>& local_cols = local.column_indices();
    local.num_cols() = policy.size();
 
    for(unsigned i=0; i<local.column_indices().size(); i++)
        assert( policy.local2globalIdx(local_cols[i], rank, local_cols[i]) );
    //local.set( local.num_rows(), policy.size(), local.row_offsets(), local_cols, local.values());
}
 
namespace create
{
 
//
template<class MPITopology, typename = std::enable_if_t<dg::is_vector_v<
    typename MPITopology::host_vector, MPIVectorTag>>>
dg::MIHMatrix_t<typename MPITopology::value_type> interpolation( const MPITopology&
        g_new, const MPITopology& g_old, std::string method = "dg")
{
    dg::IHMatrix_t<typename MPITopology::value_type> mat = dg::create::interpolation(
        g_new.local(), g_old.global(), method);
    return make_mpi_matrix(  mat, g_old);
}
// This is now dg::create::prolongation
//template<class real_type>
//dg::MIHMatrix_t<real_type> interpolation( const aRealMPITopology3d<real_type>&
//        g_new, const aRealMPITopology2d<real_type>& g_old,std::string method = "dg")
//{
//    dg::IHMatrix_t<real_type> mat = dg::create::interpolation(
//        g_new.local(), g_old.global(), method);
//    return make_mpi_matrix(  mat, g_old);
//}
 
template<class MPITopology, typename = std::enable_if_t<dg::is_vector_v<
    typename MPITopology::host_vector, MPIVectorTag>>>
dg::MIHMatrix_t<typename MPITopology::value_type> projection( const MPITopology&
        g_new, const MPITopology& g_old, std::string method = "dg")
{
    dg::IHMatrix_t<typename MPITopology::value_type> mat = dg::create::projection(
        g_new.global(), g_old.local(), method);
    convertLocal2GlobalCols( mat, g_old);
    auto mat_loc = convertGlobal2LocalRows( mat, g_new);
    return make_mpi_matrix(  mat_loc, g_old);
}
 
template<class RecursiveHostVector, class real_type, size_t Nd>
dg::MIHMatrix_t<real_type> interpolation(
        const RecursiveHostVector& x,
        const aRealMPITopology<real_type, Nd>& g,
        std::array<dg::bc, Nd> bcx,
        std::string method = "dg")
{
    dg::IHMatrix_t<real_type> mat = dg::create::interpolation( x, g.global(),
            bcx, method);
    return make_mpi_matrix(  mat, g);
}
 
template<class host_vector, class real_type>
dg::MIHMatrix_t<real_type> interpolation(
        const host_vector& x,
        const RealMPIGrid1d<real_type>& g,
        dg::bc bcx = dg::NEU,
        std::string method = "dg")
{
    dg::IHMatrix_t<real_type> mat = dg::create::interpolation( x, g.global(),
            bcx, method);
    return make_mpi_matrix(  mat, g);
}
template<class host_vector, class real_type>
dg::MIHMatrix_t<real_type> interpolation(
        const host_vector& x,
        const host_vector& y,
        const aRealMPITopology2d<real_type>& g,
        dg::bc bcx = dg::NEU, dg::bc bcy = dg::NEU,
        std::string method = "dg")
{
    dg::IHMatrix_t<real_type> mat = dg::create::interpolation( x,y, g.global(),
            bcx, bcy, method);
    return make_mpi_matrix(  mat, g);
}
 
template<class host_vector, class real_type>
dg::MIHMatrix_t<real_type> interpolation(
        const host_vector& x,
        const host_vector& y,
        const host_vector& z,
        const aRealMPITopology3d<real_type>& g,
        dg::bc bcx = dg::NEU, dg::bc bcy = dg::NEU, dg::bc bcz = dg::PER,
        std::string method = "dg")
{
    dg::IHMatrix_t<real_type> mat = dg::create::interpolation( x,y,z,
            g.global(), bcx, bcy, bcz, method);
    return make_mpi_matrix(  mat, g);
}
 
 
 
 
}//namespace create
}//namespace dg