csr_concat.cc

Go to the documentation of this file.

// build csr from initializer list (c++ 2011)
//
#include "rheolef/csr_concat.h"
 
namespace rheolef { namespace details {
 
// =========================================================================
// 1rst case : one-line matrix initializer
//  A = {a, u};                 // matrix & vector
//  A = {trans(u), 3.2};        // trans(vector) & scalar
// =========================================================================
// ------------------------
// sizes utilities
// ------------------------
template <class T, class M>
void
csr_concat_line<T,M>::set_sizes (
    std::pair<size_type,size_type>&         row_sizes,
    std::pair<size_type,size_type>&         col_sizes,
    const std::pair<size_type,size_type>&   new_row_sizes,
    const std::pair<size_type,size_type>&   new_col_sizes)
{
  if (row_sizes.first == undef) {
    row_sizes = new_row_sizes;
  } else if (new_row_sizes.first != undef) {
    check_macro (row_sizes.first == new_row_sizes.first,
            "matrix initializer list: matrix row argument [0:"
               << new_row_sizes.first << "|" << new_row_sizes.second << "["
               << " has incompatible size: expect [0:"
               << row_sizes.first  << "|" << row_sizes.second << "[");
  }
  if (col_sizes.first == undef) {
    col_sizes = new_col_sizes;
  } else if (new_col_sizes.first != undef) {
    check_macro (col_sizes.first == new_col_sizes.first,
            "matrix initializer list: matrix col argument [0:"
               << new_col_sizes.first << "|" << new_col_sizes.second << "["
               << " has incompatible size: expect [0:"
               << col_sizes.first  << "|" << col_sizes.second << "[");
  }
}
template <class T, class M>
void
csr_concat_line<T,M>::finalize_sizes (
  std::pair<size_type,size_type>& sizes,
  const communicator&             comm)
{
  // when "0" has still unknown block size, set it to 1
  if (sizes.first == undef) {
    size_type my_proc = comm.rank();
    size_type iproc0  = constraint_process_rank (comm);
    sizes.first  = (my_proc == iproc0) ? 1 : 0;
    sizes.second = 1;
  }
}
template <class T, class M>
void
csr_concat_line<T,M>::finalize_sizes (
  std::vector<std::pair<size_type,size_type> >&  sizes,
  const communicator&                            comm)
{
  for (size_type i_comp = 0; i_comp < sizes.size(); i_comp++) {
    finalize_sizes (sizes [i_comp], comm);
  }
}
// ------------------------
// pass 0 : compute sizes
// ------------------------
template <class T, class M>
void
csr_concat_line<T,M>::build_csr_pass0 (
    std::pair<size_type,size_type>&                row_sizes,
    std::vector<std::pair<size_type,size_type> >&  col_sizes,
    communicator&                                  comm) const
{
  size_type my_proc = comm.rank();
  size_type iproc0 = constraint_process_rank (comm);
  if (col_sizes.size() == 0) {
    col_sizes.resize (_l.size(), std::pair<size_type,size_type>(undef,undef));
  } else {
    check_macro (col_sizes.size() == _l.size(),
            "matrix initializer list: unexpected matrix line size "
               << _l.size() << ": size " << col_sizes.size()
               << " was expected");
  }
  size_type j_comp = 0;
  for (typename std::list<value_type>::const_iterator iter = _l.begin(); iter != _l.end(); ++iter, j_comp++) {
    const value_type& x = *iter;
    switch (x.variant) {
      case value_type::empty: {
        const sizes_type& nrows = x.e.first;
        const sizes_type& ncols = x.e.second;
        trace_macro ("block col  j_comp="<<j_comp<<": e: "<<nrows.first<<"x"<<ncols.first<<"...");
        set_sizes (row_sizes, col_sizes[j_comp], nrows, ncols);
        trace_macro ("block col  j_comp="<<j_comp<<": E: "<<row_sizes.first<<"x"<<col_sizes[j_comp].first);
        break;
      }
      case value_type::scalar: {
        size_type     nrow = undef;
        size_type dis_nrow = undef;
        size_type     ncol = undef;
        size_type dis_ncol = undef;
#ifdef TODO
        // TODO: non-zero value => force size to one ?
        if (x.s != 0) {
              nrow = (my_proc == iproc0) ? 1 : 0;
          dis_nrow = 1;
              ncol = (my_proc == iproc0) ? 1 : 0;
          dis_ncol = 1;
        }
#endif // TODO
        trace_macro ("block col  j_comp="<<j_comp<<": s: "<<nrow<<"x"<<ncol<<"...");
        set_sizes (row_sizes, col_sizes[j_comp], std::make_pair (nrow, dis_nrow), std::make_pair (ncol, dis_ncol));
        trace_macro ("block col  j_comp="<<j_comp<<": S: "<<row_sizes.first<<"x"<<col_sizes[j_comp].first);
        break;
      }
      case value_type::vector: {
        // vertical vec
        size_type     nrow = x.v.ownership().size();
        size_type dis_nrow = x.v.ownership().dis_size();
        size_type     ncol = (my_proc == iproc0) ? 1 : 0;
        size_type dis_ncol = 1;
        trace_macro ("block col  j_comp="<<j_comp<<": v: "<<nrow<<"x"<<ncol<<"...");
        set_sizes (row_sizes, col_sizes[j_comp], std::make_pair (nrow, dis_nrow), std::make_pair (ncol, dis_ncol));
        trace_macro ("block col  j_comp="<<j_comp<<": V: "<<row_sizes.first<<"x"<<col_sizes[j_comp].first);
        comm = x.v.comm(); // TODO: check if compatible comm
        break;
      }
      case value_type::vector_transpose: {
        // horizontal vec
        size_type     ncol = x.v.ownership().size();
        size_type dis_ncol = x.v.ownership().dis_size();
        size_type     nrow = (my_proc == iproc0) ? 1 : 0;
        size_type dis_nrow = 1;
        set_sizes (row_sizes, col_sizes[j_comp], std::make_pair (nrow, dis_nrow), std::make_pair (ncol, dis_ncol));
        trace_macro ("block col  j_comp="<<j_comp<<": VT: "<<row_sizes.first<<"x"<<col_sizes[j_comp].first);
        comm = x.v.comm(); // TODO: check if compatible comm
        comm = x.v.comm(); // TODO: check if compatible comm
        break;
      }
      case value_type::vector_vec: {
        // horizontal vector<vec>
        size_type     ncol = undef;
        size_type dis_ncol = undef;
        if (x.vv.size() > 0) {
              ncol = x.vv[0].ownership().size();
          dis_ncol = x.vv[0].ownership().dis_size();
              comm = x.vv[0].ownership().comm(); // TODO: check if compatible comm
        }
        size_type     nrow = (my_proc == iproc0) ? x.vv.size() : 0;
        size_type dis_nrow = x.vv.size();
        set_sizes (row_sizes, col_sizes[j_comp], std::make_pair (nrow, dis_nrow), std::make_pair (ncol, dis_ncol));
        trace_macro ("block col  j_comp="<<j_comp<<": VV: "<<row_sizes.first<<"x"<<col_sizes[j_comp].first);
        break;
      }
      case value_type::vector_vec_transpose: {
        // vertical vector<vec>
        size_type     nrow = undef;
        size_type dis_nrow = undef;
        if (x.vv.size() > 0) {
              nrow = x.vv[0].ownership().size();
          dis_nrow = x.vv[0].ownership().dis_size();
              comm = x.vv[0].ownership().comm(); // TODO: check if compatible comm
        }
        size_type     ncol = (my_proc == iproc0) ? x.vv.size() : 0;
        size_type dis_ncol = x.vv.size();
        set_sizes (row_sizes, col_sizes[j_comp], std::make_pair (nrow, dis_nrow), std::make_pair (ncol, dis_ncol));
        trace_macro ("block col  j_comp="<<j_comp<<": VVT: "<<row_sizes.first<<"x"<<col_sizes[j_comp].first);
        break;
      }
      case value_type::matrix: {
        size_type     nrow = x.m.row_ownership().size();
        size_type dis_nrow = x.m.row_ownership().dis_size();
        size_type     ncol = x.m.col_ownership().size();
        size_type dis_ncol = x.m.col_ownership().dis_size();
        set_sizes (row_sizes, col_sizes[j_comp], std::make_pair (nrow, dis_nrow), std::make_pair (ncol, dis_ncol));
        trace_macro ("block col  j_comp="<<j_comp<<": m: "<<nrow<<"x"<<ncol<<"...");
        comm = x.m.row_ownership().comm(); // TODO: check if compatible comm
        trace_macro ("block col  j_comp="<<j_comp<<": M: "<<row_sizes.first<<"x"<<col_sizes[j_comp].first);
        break;
        break;
      } 
    }
    trace_macro ("block col  j_comp="<<j_comp<<" done");
  }
}
// ---------------------------
// pass 1 : compute ownership
// ---------------------------
template <class T, class M>
void
csr_concat_line<T,M>::build_csr_pass1 (
    const std::pair<size_type,size_type>&                row_sizes,
    const std::vector<std::pair<size_type,size_type> >&  col_sizes,
    const communicator&                                  comm,
    distributor&                                         row_ownership,
    distributor&                                         col_ownership,
    std::vector<distributor>&                            col_start_by_component) const
{
  row_ownership = distributor (row_sizes.second, comm, row_sizes.first);
  size_type col_comp_start_j     = 0;
  size_type col_comp_start_dis_j = 0;
  col_start_by_component.resize (col_sizes.size());
  for (size_type j_comp = 0; j_comp < col_sizes.size(); j_comp++) {
    col_start_by_component [j_comp] = distributor (col_comp_start_dis_j, comm, col_comp_start_j);
    col_comp_start_j     += col_sizes [j_comp].first;
    col_comp_start_dis_j += col_sizes [j_comp].second;
  }
  col_ownership = distributor (col_comp_start_dis_j, comm, col_comp_start_j);
}
// -----------------------------------
// pass 2 : copy at the right location
// -----------------------------------
template <class T, class M>
void
csr_concat_line<T,M>::build_csr_pass2 (
    asr<T,M>&                                            A,
    const std::pair<size_type,size_type>&                row_sizes,
    const distributor&                                   row_start_by_component,
    const std::vector<std::pair<size_type,size_type> >&  col_sizes,
    const std::vector<distributor>&                      col_start_by_component) const
{
  communicator comm = A.row_ownership().comm();
  size_type my_proc = comm.rank();
  size_type iproc0 = constraint_process_rank (comm);
  size_type j_comp = 0;
  for (typename std::list<value_type>::const_iterator iter = _l.begin(); iter != _l.end(); ++iter, j_comp++) {
    const value_type& x = *iter;
    switch (x.variant) {
      case value_type::empty: {
        // no value to insert
        break;
      }
      case value_type::scalar: {
        // do not add eplicit zero values
        size_type dis_i = A.row_ownership().first_index(iproc0)
                        + row_start_by_component.size(iproc0);
        size_type dis_j = A.col_ownership().first_index(iproc0)
                        + col_start_by_component[j_comp].size(iproc0);
        if (x.s != 0 && my_proc == iproc0) A.dis_entry (dis_i, dis_j) += x.s;
        break;
      }
      case value_type::vector: {
        // vertical vec : moved column to iproc0, but still distributed by lines
        size_type dis_i = A.row_ownership().first_index()
                        + row_start_by_component.size();
        size_type dis_j = A.col_ownership().first_index(iproc0)
                        + col_start_by_component[j_comp].size(iproc0);
        for (typename vec<T,M>::const_iterator iter = x.v.begin(), last = x.v.end(); iter != last; iter++, dis_i++) {
          if (*iter != 0) A.dis_entry (dis_i, dis_j) += *iter;
        }
        break;
      }
      case value_type::vector_transpose: {
        // horizontal vec: row integrally hosted by iproc=0 ; col still distributed
        size_type dis_i = A.row_ownership().first_index(iproc0)
                        + row_start_by_component.size(iproc0);
        size_type dis_j = A.col_ownership().first_index()
                        + col_start_by_component[j_comp].size();
        for (typename vec<T,M>::const_iterator iter = x.v.begin(), last = x.v.end(); iter != last; iter++, dis_j++) {
          if (*iter != 0) A.dis_entry (dis_i, dis_j) += *iter;
        }
        break;
      }
      case value_type::vector_vec: {
        // horizontal vector<vec>: rows integrally hosted by iproc=0 ; col still distributed
        size_type n = x.vv.size(); 
        for (size_type k = 0; k < n; k++) {
          size_type dis_i = A.row_ownership().first_index(iproc0)
                          + row_start_by_component.size(iproc0)
                          + k;
          size_type dis_j = A.col_ownership().first_index()
                          + col_start_by_component[j_comp].size();
          for (typename vec<T,M>::const_iterator iter = x.vv[k].begin(), last = x.vv[k].end(); iter != last; iter++, dis_j++) {
            if (*iter != 0) A.dis_entry (dis_i, dis_j) += *iter;
          }
        }
        break;
      }
      case value_type::vector_vec_transpose: {
        // vertical vector<vec>: moved column to iproc0, but still distributed by lines
        size_type n = x.vv.size(); 
        for (size_type k = 0; k < n; k++) {
          size_type dis_i = A.row_ownership().first_index()
                          + row_start_by_component.size();
          size_type dis_j = A.col_ownership().first_index(iproc0)
                          + col_start_by_component[j_comp].size(iproc0)
                          + k;
          for (typename vec<T,M>::const_iterator iter = x.vv[k].begin(), last = x.vv[k].end(); iter != last; iter++, dis_i++) {
            if (*iter != 0) A.dis_entry (dis_i, dis_j) += *iter;
          }
        }
        break;
      }
      case value_type::matrix: {
        typename csr<T,M>::const_iterator dia_ia = x.m.begin();
        typename csr<T,M>::const_iterator ext_ia = x.m.ext_begin();
        for (size_type i = 0, n = x.m.nrow(); i < n; i++) {
          size_type dis_i = A.row_ownership().first_index()
                          + row_start_by_component.size()
                          + i;
          for (typename csr<T,M>::const_data_iterator p = dia_ia[i]; p < dia_ia[i+1]; p++) {
            size_type j = (*p).first;
            size_type dis_j = A.col_ownership().first_index()
                            + col_start_by_component[j_comp].size()
                            + j;
            A.dis_entry (dis_i, dis_j) += (*p).second;
          }
          if (is_distributed<M>::value && x.m.ext_nnz() != 0) {
            for (typename csr<T,M>::const_data_iterator p = ext_ia[i]; p < ext_ia[i+1]; p++) {
              size_type comp_dis_j = x.m.jext2dis_j ((*p).first);
              size_type iproc = x.m.col_ownership().find_owner (comp_dis_j);
              size_type first_comp_dis_j_iproc = x.m.col_ownership().first_index (iproc);
              assert_macro (comp_dis_j >= first_comp_dis_j_iproc, "invalid index");
              size_type j = comp_dis_j - first_comp_dis_j_iproc;
              size_type dis_j = A.col_ownership().first_index(iproc)
                              + col_start_by_component[j_comp].size(iproc)
                              + j;
              A.dis_entry (dis_i, dis_j) += (*p).second;
            }
          }
        }
        break;
      } 
    }
  }
}
// ------------------------------
// merge passes 1 & 2
// ------------------------------
template <class T, class M>
csr<T,M>
csr_concat_line<T,M>::build_csr() const
{
  communicator comm;
  std::pair<size_type,size_type>               row_sizes = std::pair<size_type,size_type>(undef,undef);
  distributor                                  row_start_by_component (0, comm, 0);
  std::vector<std::pair<size_type,size_type> > col_sizes;
  std::vector<distributor>                     col_start_by_component;
  build_csr_pass0 (row_sizes, col_sizes, comm);
  finalize_sizes  (row_sizes, comm);
  finalize_sizes  (col_sizes, comm);
  distributor              row_ownership;
  distributor              col_ownership;
  build_csr_pass1 (row_sizes, col_sizes, comm, row_ownership, col_ownership, col_start_by_component);
  asr<T,M> A (row_ownership, col_ownership);
  build_csr_pass2 (A, row_sizes, row_start_by_component, col_sizes, col_start_by_component);
  A.dis_entry_assembly();
  return csr<T,M>(A);
}
// =========================================================================
// 2nd case : multi-line matrix initializer
//  A = { {a,        u  },
//        {trans(u), 3.2} };
// =========================================================================
template <class T, class M>
csr<T,M>
csr_concat<T,M>::build_csr() const
{
  // ---------------------------
  // pass 0 : compute sizes
  // ---------------------------
  std::vector<std::pair<size_type,size_type> > row_sizes (_l.size(), std::pair<size_type,size_type>(undef,undef));
  std::vector<std::pair<size_type,size_type> > col_sizes;
  communicator comm;
  size_type i_comp = 0;
  for (typename std::list<line_type>::const_iterator iter = _l.begin(); iter != _l.end(); ++iter, i_comp++) {
    const line_type& line = *iter;
    trace_macro ("block line i_comp="<<i_comp<<"...");
    line.build_csr_pass0 (row_sizes[i_comp], col_sizes, comm);
    trace_macro ("block line i_comp="<<i_comp<<" done");
  }
  csr_concat_line<T,M>::finalize_sizes (row_sizes, comm);
  csr_concat_line<T,M>::finalize_sizes (col_sizes, comm);
  // ---------------------------
  // pass 1 : compute ownership
  // ---------------------------
  size_type row_comp_start_i     = 0;
  size_type row_comp_start_dis_i = 0;
  std::vector<distributor> row_start_by_component (row_sizes.size()+1);
  for (size_type i_comp = 0; i_comp <= row_sizes.size(); i_comp++) {
    row_start_by_component [i_comp] = distributor (row_comp_start_dis_i, comm, row_comp_start_i);
    if (i_comp == row_sizes.size()) break;
    row_comp_start_i     += row_sizes [i_comp].first;
    row_comp_start_dis_i += row_sizes [i_comp].second;
  }
  distributor row_ownership = row_start_by_component[row_sizes.size()];
 
  size_type col_comp_start_j     = 0;
  size_type col_comp_start_dis_j = 0;
  std::vector<distributor> col_start_by_component (col_sizes.size()+1);
  for (size_type j_comp = 0; j_comp <= col_sizes.size(); j_comp++) {
    col_start_by_component [j_comp] = distributor (col_comp_start_dis_j, comm, col_comp_start_j);
    if (j_comp == col_sizes.size()) break;
    col_comp_start_j     += col_sizes [j_comp].first;
    col_comp_start_dis_j += col_sizes [j_comp].second;
  }
  distributor col_ownership = col_start_by_component[col_sizes.size()];
  // ------------------------
  // pass 2 : copy
  // ------------------------
  asr<T,M> a (row_ownership, col_ownership);
  i_comp = 0;
  for (typename std::list<line_type>::const_iterator iter = _l.begin(); iter != _l.end(); ++iter, ++i_comp) {
    const line_type& line = *iter;
    line.build_csr_pass2 (a, row_sizes[i_comp], row_start_by_component[i_comp], col_sizes, col_start_by_component);
  }
  a.dis_entry_assembly();
  return csr<T,M>(a);
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T,M)                                     \
template class csr_concat_line<T,M>;                                    \
template class csr_concat<T,M>;
 
_RHEOLEF_instanciation(Float,sequential)
#ifdef _RHEOLEF_HAVE_MPI
_RHEOLEF_instanciation(Float,distributed)
#endif // _RHEOLEF_HAVE_MPI
 
}} // namespace rheolef::details