p2p_precompute.hpp

Go to the documentation of this file.

 
#ifndef NIHU_P2P_PRECOMPUTE_HPP_INCLUDED
#define NIHU_P2P_PRECOMPUTE_HPP_INCLUDED
 
#include "cluster_tree.hpp"
#include "fmm_operator.hpp"
#include "lists.hpp"
#include "local_operator.hpp"
#include "util/eigen_utils.hpp"
#include "util/matrix_traits.hpp"
 
#include <Eigen/SparseCore>
 
#ifdef NIHU_FMM_PARALLEL
#include <omp.h>
#endif
 
#include <chrono>
#include <memory>
#include <type_traits>
#include <utility>
#include <vector>
 
namespace NiHu
{
namespace fmm
{
 
namespace internal
{
 
template <class Operator, class ClusterDerived, bool isLocal>
class sparse_index_computer;
 
template <class Operator, class ClusterDerived>
class sparse_index_computer<Operator, ClusterDerived, false>
{
    typedef typename std::decay<Operator>::type operator_t;
    typedef unsigned int index_t;
    typedef std::vector<std::pair<index_t, index_t> > indices_t;
    typedef ClusterDerived cluster_t;
    typedef cluster_tree<cluster_t> tree_t;
 
public:
 
    static indices_t eval(tree_t const &tree,
        interaction_lists::list_t const &p2p_list)
    {
        indices_t indices;
 
        // determine number of source- receiver pairs
        size_t s = 0;
        for (unsigned to = 0; to < p2p_list.size(); ++to)
            for (auto from : p2p_list[to])
                s += tree[to].get_rec_node_idx().size() *
                tree[from].get_src_node_idx().size();
 
        indices.reserve(s);
 
        for (unsigned to = 0; to < p2p_list.size(); ++to)
            for (auto from : p2p_list[to])
                for (auto i : tree[to].get_rec_node_idx())          // loop over receiver nodes
                    for (auto j : tree[from].get_src_node_idx())    // loop over source nodes
                        indices.push_back(std::make_pair(i, j));
 
        return indices;
    }
};
 
template <class Operator, class ClusterDerived>
class sparse_index_computer<Operator, ClusterDerived, true>
{
    typedef typename std::decay<Operator>::type operator_t;
    typedef std::vector<std::pair<unsigned int, unsigned int> > indices_t;
    typedef ClusterDerived cluster_t;
    typedef cluster_tree<cluster_t> tree_t;
 
public:
    static indices_t eval(tree_t const &tree,
        interaction_lists::list_t const &p2p_list)
    {
        indices_t indices;
 
        // determine number of src receiver pairs
        size_t s = 0;
        for (unsigned to = 0; to < p2p_list.size(); ++to)
            if (std::find(p2p_list[to].begin(), p2p_list[to].end(), to) != p2p_list[to].end())
                s += tree[to].get_rec_node_idx().size();
 
        indices.reserve(s);
 
        for (unsigned to = 0; to < p2p_list.size(); ++to)
            if (std::find(p2p_list[to].begin(), p2p_list[to].end(), to) != p2p_list[to].end())
                for (auto i : tree[to].get_rec_node_idx())
                    indices.push_back(std::make_pair(i, i));
 
        return indices;
    }
};
 
template <class Operator, class ClusterDerived, bool isResultEigen = is_eigen<
    typename std::decay<Operator>::type::result_t>::value>
class sparse_computer;
 
 
template <class Operator, class ClusterDerived>
class sparse_computer<Operator, ClusterDerived, true>
{
    typedef typename std::decay<Operator>::type operator_t;
    typedef typename operator_t::result_t result_t;
    typedef typename scalar<result_t>::type scalar_t;
    static size_t const rows = num_rows<result_t>::value;
    static size_t const cols = num_cols<result_t>::value;
    typedef std::ptrdiff_t storage_index_t;
    typedef Eigen::SparseMatrix<scalar_t, Eigen::RowMajor, storage_index_t> sparse_t;
    typedef ClusterDerived cluster_t;
    typedef cluster_tree<cluster_t> tree_t;
 
    static bool const is_local = is_local_operator<operator_t>::value;
    typedef sparse_index_computer<Operator, ClusterDerived, is_local> index_computer_t;
 
public:
    static sparse_t *eval(Operator &&op, tree_t const &tree, interaction_lists::list_t const &list, size_t &assembly_time)
    {
        auto tstart = std::chrono::steady_clock::now();
 
        auto indices = index_computer_t::eval(tree, list);
        size_t s = indices.size();
 
        sparse_t *pmat = new sparse_t(rows * tree.get_n_rec_nodes(), cols * tree.get_n_src_nodes());
        Eigen::Matrix<unsigned, Eigen::Dynamic, 1> num_nonzeros_per_row(pmat->rows(),1);
        num_nonzeros_per_row.setZero();
 
        for (size_t isrcrec = 0; isrcrec < s; ++isrcrec)
        {
            auto i = indices[isrcrec].first;
            for (int ii = 0; ii < rows; ++ii)
                num_nonzeros_per_row(i * rows + ii) += cols;
        }
        pmat->reserve(num_nonzeros_per_row);
        
#ifdef NIHU_FMM_PARALLEL
// #pragma omp parallel for
#endif
        for (int isrcrec = 0; isrcrec < s; ++isrcrec)
        {
            auto i = indices[isrcrec].first;
            auto j = indices[isrcrec].second;
            result_t opmat = op(i, j);
            for (size_t ii = 0; ii < rows; ++ii)    // loop over matrix rows
                for (size_t jj = 0; jj < cols; ++jj) // loop over matrix cols 
                    pmat->insert(i * rows + ii, j * cols + jj) =  opmat(ii, jj);
        }
#ifdef NIHU_FMM_PARALLEL
// #pragma omp barrier
#endif
        pmat->makeCompressed();
 
        auto tend = std::chrono::steady_clock::now();
        assembly_time = std::chrono::duration_cast<std::chrono::microseconds>(tend - tstart).count();
 
        return pmat;
    }
};
 
template <class Operator, class ClusterDerived>
class sparse_computer<Operator, ClusterDerived, false>
{
    typedef typename std::decay<Operator>::type operator_t;
    typedef typename operator_t::result_t result_t;
    using scalar_t = typename scalar<result_t>::type;
    typedef std::ptrdiff_t storage_index_t;
    typedef Eigen::SparseMatrix<scalar_t, Eigen::RowMajor, storage_index_t> sparse_t;
    typedef ClusterDerived cluster_t;
    typedef cluster_tree<cluster_t> tree_t;
 
    static bool const is_local = is_local_operator<operator_t>::value;
    typedef sparse_index_computer<Operator, ClusterDerived, is_local> index_computer_t;
 
public:
 
    static sparse_t *eval(
        Operator &&op, 
        tree_t const &tree, 
        interaction_lists::list_t const &list, 
        size_t &assembly_time)
    {
        auto tstart = std::chrono::steady_clock::now();
 
        // i-j indices of final sparse matrix
        auto indices = index_computer_t::eval(tree, list);
        size_t s = indices.size();
 
        sparse_t *pmat = new sparse_t(tree.get_n_rec_nodes(), tree.get_n_src_nodes());
        Eigen::Matrix<unsigned, Eigen::Dynamic, 1> num_nonzeros_per_row(pmat->rows(),1);
        num_nonzeros_per_row.setZero();
 
        for (size_t isrcrec = 0; isrcrec < s; ++isrcrec)
            num_nonzeros_per_row(indices[isrcrec].first)++;
        pmat->reserve(num_nonzeros_per_row);
        
#ifdef NIHU_FMM_PARALLEL
// #pragma omp parallel for
#endif
        for (int isrcrec = 0; isrcrec < s; ++isrcrec)
        {
            size_t i = indices[isrcrec].first;
            size_t j = indices[isrcrec].second;
            pmat->insert(i, j) = op(i, j);
        }
#ifdef NIHU_FMM_PARALLEL
// #pragma omp barrier
#endif
        pmat->makeCompressed();
 
        auto tend = std::chrono::steady_clock::now();
        assembly_time = std::chrono::duration_cast<std::chrono::microseconds>(tend - tstart).count();
        
        return pmat;
    }
};
 
 
} // end of namespace internal
 
template <class Scalar, size_t NumDofPerSrc, size_t NumDofPerRec>
class p2p_precompute
    : public fmm_operator<p2p_tag>
{
private:
 
public:
    typedef Scalar scalar_t;
 
    // RowMajor is needed to allow Eigen's parallelization of sparse matrix * dense vector
    typedef Eigen::SparseMatrix<scalar_t, Eigen::RowMajor, std::ptrdiff_t> sparse_t;
    typedef Eigen::Matrix<scalar_t, Eigen::Dynamic, 1> response_t;
 
    static size_t const num_dof_per_src = NumDofPerSrc;
    static size_t const num_dof_per_rec = NumDofPerRec;
 
    template <class Operator, class ClusterDerived>
    p2p_precompute(Operator &&op, cluster_tree<ClusterDerived> const &tree, interaction_lists::list_t const &list)
        : m_pmat(internal::sparse_computer<Operator, ClusterDerived>::eval(std::forward<Operator>(op), tree, list, m_assembly_time))
    {
    }
 
    sparse_t const &get_sparse_matrix() const
    {
        return *m_pmat;
    }
    
    size_t get_assembly_time() const
    {
        return m_assembly_time;
    }
 
    template <class RHS>
    response_t operator*(RHS const &rhs) const
    {
        return *m_pmat * rhs;
    }
 
private:
    std::shared_ptr<sparse_t> m_pmat;
    size_t m_assembly_time;
};
 
template <class Operator, class ClusterDerived>
auto create_p2p_precompute(Operator &&op,
    cluster_tree<ClusterDerived> const &tree,
    interaction_lists::list_t const &list)
{
    typedef typename std::decay<Operator>::type operator_t;
    return p2p_precompute<
        typename scalar<typename operator_t::result_t>::type,
        operator_t::num_dof_per_src,
        operator_t::num_dof_per_rec
    >(std::forward<Operator>(op), tree, list);
}
 
 
} // end of namespace fmm
} // end of namespace NiHu
 
#endif /* NIHU_P2P_PRECOMPUTE_HPP_INCLUDED */