packages/MC/mc/Fission_Matrix_Processor.cc
//---------------------------------------------------------------------------//
// PUBLIC FUNCTIONS
//---------------------------------------------------------------------------//
/*!
 * \brief Globally reduce and build the fission matrix.
 *
 * This builds a flattened fission matrix.  The order is stored in the graph.
 */
PUBLIC FUNCTIONS
\brief Globally reduce and build the fission matrix.

This builds a flattened fission matrix.  The order is stored in the graph.
PUBLIC FUNCTIONS \ brief Globally reduce and build the fission matrix . This builds a flattened fission matrix . The order is stored in the graph .
void Fission_Matrix_Processor::build_matrix(
    const Sparse_Matrix &local_matrix,
    const Denominator   &local_denominator)
{
    REQUIRE(!local_matrix.empty());
    REQUIRE(local_matrix.size() <= local_denominator.size() *
            local_denominator.size());
    reset();
    d_N = local_denominator.size();
    Ordered_Graph local_graph(local_matrix.size());
    auto itr = local_graph.begin();
    for (const auto &element : local_matrix)
    {
        *itr = element.first;
        ++itr;
    }
    std::sort(local_graph.begin(), local_graph.end());
    d_graph.insert(d_graph.end(), local_graph.begin(), local_graph.end());
    CHECK(d_graph.size() == local_graph.size());
    reduce();
    int size = d_graph.size();
    profugus::broadcast(&size, 1, 0);
    if (d_node > 0)
    {
        Ordered_Graph g(size);
        std::swap(d_graph, g);
    }
    CHECK(d_graph.size() == size);
    profugus::broadcast(&d_graph[0].first, size * 2, 0);
    d_matrix.resize(size);
    {
        auto g = d_graph.begin();
        auto l = local_graph.begin();
        auto m = d_matrix.begin();
        while (l != local_graph.end())
        {
            CHECK(m != d_matrix.end());
            CHECK(g != d_graph.end());
            if (*l == *g)
            {
                CHECK(local_matrix.find(*l) != local_matrix.end());
                CHECK(local_matrix.find(*l)->first == *l);
                *m = local_matrix.find(*l)->second;
                ++l; ++g; ++m;
            }
            else
            {
                ++g; ++m;
            }
        }
    }
    Denominator denominator(local_denominator.begin(), local_denominator.end());
    profugus::global_sum(&d_matrix[0], size);
    profugus::global_sum(&denominator[0], d_N);
    for (int n = 0, j = 0; n < size; ++n)
    {
        j = d_graph[n].second;
        CHECK(j < denominator.size());
        CHECK(denominator[j] > 0.0);
        d_matrix[n] /= denominator[j];
    }
}
void Fission_Matrix_Processor :: build_matrix ( const Sparse_Matrix & local_matrix , const Denominator & local_denominator ) { REQUIRE ( ! local_matrix . empty ( ) ) ; REQUIRE ( local_matrix . size ( ) <= local_denominator . size ( ) * local_denominator . size ( ) ) ; reset ( ) ; d_N = local_denominator . size ( ) ; Ordered_Graph local_graph ( local_matrix . size ( ) ) ; auto itr = local_graph . begin ( ) ; for ( const auto & element : local_matrix ) { * itr = element . first ; ++ itr ; } std :: sort ( local_graph . begin ( ) , local_graph . end ( ) ) ; d_graph . insert ( d_graph . end ( ) , local_graph . begin ( ) , local_graph . end ( ) ) ; CHECK ( d_graph . size ( ) == local_graph . size ( ) ) ; reduce ( ) ; int size = d_graph . size ( ) ; profugus :: broadcast ( & size , 1 , 0 ) ; if ( d_node > 0 ) { Ordered_Graph g ( size ) ; std :: swap ( d_graph , g ) ; } CHECK ( d_graph . size ( ) == size ) ; profugus :: broadcast ( & d_graph [ 0 ] . first , size * 2 , 0 ) ; d_matrix . resize ( size ) ; { auto g = d_graph . begin ( ) ; auto l = local_graph . begin ( ) ; auto m = d_matrix . begin ( ) ; while ( l != local_graph . end ( ) ) { CHECK ( m != d_matrix . end ( ) ) ; CHECK ( g != d_graph . end ( ) ) ; if ( * l == * g ) { CHECK ( local_matrix . find ( * l ) != local_matrix . end ( ) ) ; CHECK ( local_matrix . find ( * l ) -> first == * l ) ; * m = local_matrix . find ( * l ) -> second ; ++ l ; ++ g ; ++ m ; } else { ++ g ; ++ m ; } } } Denominator denominator ( local_denominator . begin ( ) , local_denominator . end ( ) ) ; profugus :: global_sum ( & d_matrix [ 0 ] , size ) ; profugus :: global_sum ( & denominator [ 0 ] , d_N ) ; for ( int n = 0 , j = 0 ; n < size ; ++ n ) { j = d_graph [ n ] . second ; CHECK ( j < denominator . size ( ) ) ; CHECK ( denominator [ j ] > 0.0 ) ; d_matrix [ n ] /= denominator [ j ] ; } }
PUBLIC FUNCTIONS
 \brief Globally reduce and build the fission matrix.
PUBLIC FUNCTIONS \ brief Globally reduce and build the fission matrix .
// reset the internal storage // the NxN fission matrix size is determined by the size of the // Denominator, which equals N // on every domain, write the local matrix into an ordered graph // sort the local graph // initialize the global graph on this domain with the local graph // do a parallel merge/sort on the global graph // broadcast the graph // resize the graph on the work nodes // broadcast the graph // walk the local graph and write local data into the global matrix // get iterators to each data field // since all of the graphs (local and global) are ordered, we can // traverse the local graph and simply advance the global graph until // we hit the correct index in the local graph; this algorithm is O(N) // where N is the number of entries in the global graph // write the local denominator into the global denominator // now do a global reduction on the matrix and denominator // finally normalize the fission matrix // get the j column index (source cell index in the denominator) // normalize this element of the fission matrix (i,j) with the // starting source weight in j