packages/Alea/mc_solvers/MC_Data.cc
//---------------------------------------------------------------------------//
// Build probability and weight matrices.
//---------------------------------------------------------------------------//
Build probability and weight matrices.
Build probability and weight matrices .
void MC_Data::buildMonteCarloMatrices()
{
    REQUIRE( d_H != Teuchos::null );
    Teuchos::RCP<Teuchos::ParameterList> mc_pl =
        Teuchos::sublist(d_pl,"Monte Carlo");
    SCALAR abs_prob = mc_pl->get("absorption_probability",0.0);
    SCALAR trans_factor = mc_pl->get("transition_factor",1.0);
    VALIDATE(mc_pl->isType<std::string>("mc_type"),"Must specify mc_type.");
    std::string type = mc_pl->get<std::string>("mc_type");
    INSIST( type == "forward" || type == "adjoint",
            "mc_type must be forward or adjoint." );
    if( type == "adjoint" )
    {
        typedef Kokkos::Compat::KokkosDeviceWrapperNode<HOST> HostNode;
        auto host_node = KokkosClassic::Details::getNode<HostNode>();
        auto dev_node = d_H->getNode();
        auto H_host = d_H->clone(host_node);
        Tpetra::RowMatrixTransposer<SCALAR,LO,GO,HostNode> transposer(H_host);
        H_host = transposer.createTranspose();
        d_H = H_host->clone(dev_node);
    }
    CHECK( d_H != Teuchos::null );
    LO N = d_H->getNodeNumRows();
    LO max_nnz = d_H->getNodeMaxNumRowEntries();
    d_P = Teuchos::rcp( new CRS_MATRIX(d_H->getDomainMap(),max_nnz) );
    d_W = Teuchos::rcp( new CRS_MATRIX(d_H->getDomainMap(),max_nnz) );
    Teuchos::ArrayRCP<SCALAR> H_vals(max_nnz), P_vals(max_nnz), W_vals(max_nnz);
    Teuchos::ArrayRCP<GO>     H_inds(max_nnz);
    size_t num_entries;
    for( LO irow=0; irow<N; ++irow )
    {
        d_H->getGlobalRowCopy(irow,H_inds(),H_vals(),num_entries);
        MAGNITUDE row_sum=0.0;
        for( size_t icol=0; icol<num_entries; ++icol )
        {
            if( SCALAR_TRAITS::magnitude(H_vals[icol]) > 0.0 )
            {
                row_sum += SCALAR_TRAITS::pow(SCALAR_TRAITS::magnitude(
                            H_vals[icol]),trans_factor);
            }
        }
        if( row_sum > 0.0 )
        {
            SCALAR cdf = 0.0;
            for( size_t icol=0; icol<num_entries; ++icol )
            {
                MAGNITUDE this_prob = 0.0;
                if( SCALAR_TRAITS::magnitude(H_vals[icol]) > 0.0 )
                {
                    this_prob = SCALAR_TRAITS::pow(
                        SCALAR_TRAITS::magnitude( H_vals[icol] ), trans_factor )
                        / row_sum;
                }
                this_prob /= (1.0 - abs_prob);
                cdf += this_prob;
                P_vals[icol] = cdf;
                if( this_prob > 0.0 )
                {
                    W_vals[icol] = H_vals[icol] / this_prob;
                }
                else
                {
                    W_vals[icol] = 0.0;
                }
            }
        }
        else
        {
            std::fill( P_vals.begin(), P_vals.end(), 0.0 );
            std::fill( W_vals.begin(), W_vals.end(), 0.0 );
        }
        d_P->insertGlobalValues(irow,H_inds(0,num_entries),P_vals(0,num_entries));
        d_W->insertGlobalValues(irow,H_inds(0,num_entries),W_vals(0,num_entries));
    }
    d_P->fillComplete();
    d_W->fillComplete();
}
void MC_Data :: buildMonteCarloMatrices ( ) { REQUIRE ( d_H != Teuchos :: null ) ; Teuchos :: RCP < Teuchos :: ParameterList > mc_pl = Teuchos :: sublist ( d_pl , " " ) ; SCALAR abs_prob = mc_pl -> get ( " " , 0.0 ) ; SCALAR trans_factor = mc_pl -> get ( " " , 1.0 ) ; VALIDATE ( mc_pl -> isType < std :: string > ( " " ) , " " ) ; std :: string type = mc_pl -> get < std :: string > ( " " ) ; INSIST ( type == " " || type == " " , " " ) ; if ( type == " " ) { typedef Kokkos :: Compat :: KokkosDeviceWrapperNode < HOST > HostNode ; auto host_node = KokkosClassic :: Details :: getNode < HostNode > ( ) ; auto dev_node = d_H -> getNode ( ) ; auto H_host = d_H -> clone ( host_node ) ; Tpetra :: RowMatrixTransposer < SCALAR , LO , GO , HostNode > transposer ( H_host ) ; H_host = transposer . createTranspose ( ) ; d_H = H_host -> clone ( dev_node ) ; } CHECK ( d_H != Teuchos :: null ) ; LO N = d_H -> getNodeNumRows ( ) ; LO max_nnz = d_H -> getNodeMaxNumRowEntries ( ) ; d_P = Teuchos :: rcp ( new CRS_MATRIX ( d_H -> getDomainMap ( ) , max_nnz ) ) ; d_W = Teuchos :: rcp ( new CRS_MATRIX ( d_H -> getDomainMap ( ) , max_nnz ) ) ; Teuchos :: ArrayRCP < SCALAR > H_vals ( max_nnz ) , P_vals ( max_nnz ) , W_vals ( max_nnz ) ; Teuchos :: ArrayRCP < GO > H_inds ( max_nnz ) ; size_t num_entries ; for ( LO irow = 0 ; irow < N ; ++ irow ) { d_H -> getGlobalRowCopy ( irow , H_inds ( ) , H_vals ( ) , num_entries ) ; MAGNITUDE row_sum = 0.0 ; for ( size_t icol = 0 ; icol < num_entries ; ++ icol ) { if ( SCALAR_TRAITS :: magnitude ( H_vals [ icol ] ) > 0.0 ) { row_sum += SCALAR_TRAITS :: pow ( SCALAR_TRAITS :: magnitude ( H_vals [ icol ] ) , trans_factor ) ; } } if ( row_sum > 0.0 ) { SCALAR cdf = 0.0 ; for ( size_t icol = 0 ; icol < num_entries ; ++ icol ) { MAGNITUDE this_prob = 0.0 ; if ( SCALAR_TRAITS :: magnitude ( H_vals [ icol ] ) > 0.0 ) { this_prob = SCALAR_TRAITS :: pow ( SCALAR_TRAITS :: magnitude ( H_vals [ icol ] ) , trans_factor ) / row_sum ; } this_prob /= ( 1.0 - abs_prob ) ; cdf += this_prob ; P_vals [ icol ] = cdf ; if ( this_prob > 0.0 ) { W_vals [ icol ] = H_vals [ icol ] / this_prob ; } else { W_vals [ icol ] = 0.0 ; } } } else { std :: fill ( P_vals . begin ( ) , P_vals . end ( ) , 0.0 ) ; std :: fill ( W_vals . begin ( ) , W_vals . end ( ) , 0.0 ) ; } d_P -> insertGlobalValues ( irow , H_inds ( 0 , num_entries ) , P_vals ( 0 , num_entries ) ) ; d_W -> insertGlobalValues ( irow , H_inds ( 0 , num_entries ) , W_vals ( 0 , num_entries ) ) ; } d_P -> fillComplete ( ) ; d_W -> fillComplete ( ) ; }
Build probability and weight matrices.
Build probability and weight matrices .
// Get absorption probability // Determine if we want forward or adjoint MC // Determine "Base" matrix B, such that B = P \circ W // Forward MC -> B=H // Adjoint MC -> B=H^T // Transpose H /*
        Tpetra::RowMatrixTransposer<SCALAR,LO,GO,NODE> transposer(d_H);
        d_H = transposer.createTranspose();
        */ // Transposer is broken for Cuda, this will be fun // Now loop over rows in B to build probability/weight matrices // Get row from H^T // Compute row sum // Never include zeros