Commit 9f226fd4 authored by Utkarsh Ayachit's avatar Utkarsh Ayachit

Adding histogram.

parent 46c31e47
#------------------------------------------------------------------------------
# Since we need MPI left and right, this makes it easier to deal with MPI.
find_package(MPI REQUIRED)
function(mpi_link target)
target_include_directories(${target}
SYSTEM PRIVATE ${MPI_C_INCLUDE_PATH} ${MPI_CXX_INCLUDE_PATH})
target_compile_definitions(${target}
PRIVATE ${MPI_C_COMPILE_FLAGS} ${MPI_CXX_COMPILE_FLAGS})
target_link_libraries(${target}
PRIVATE ${MPI_C_LIBRARIES} ${MPI_CXX_LIBRARIES})
if(MPI_C_LINK_FLAGS OR MPI_CXX_LINK_FLAGS)
set_target_properties(${target}
PROPERTIES LINK_FLAGS ${MPI_C_LINK_FLAGS} ${MPI_CXX_LINK_FLAGS})
endif()
endfunction()
cmake_minimum_required(VERSION 2.8.12)
project(sensei)
list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/CMake")
include(CMakeDependentOption)
# Default to release.
......@@ -8,7 +9,6 @@ if(NOT CMAKE_BUILD_TYPE)
CACHE STRING "Choose the type of build, options are: Debug Release RelWithDebInfo MinSizeRel." FORCE)
set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS "Debug" "Release" "MinSizeRel" "RelWithDebInfo")
endif()
#------------------------------------------------------------------------------
# Set default output directories.
set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib)
......@@ -31,21 +31,6 @@ option(ENABLE_HISTOGRAM "Enable histogram analysis" ON)
# Options to enable various infrastructures.
#------------------------------------------------------------------------------
# Since we need MPI left and right, this makes it easier to deal with MPI.
find_package(MPI REQUIRED)
function(mpi_link target)
target_include_directories(${target}
SYSTEM PRIVATE ${MPI_C_INCLUDE_PATH} ${MPI_CXX_INCLUDE_PATH})
target_compile_definitions(${target}
PRIVATE ${MPI_C_COMPILE_FLAGS} ${MPI_CXX_COMPILE_FLAGS})
target_link_libraries(${target}
PRIVATE ${MPI_C_LIBRARIES} ${MPI_CXX_LIBRARIES})
if(MPI_C_LINK_FLAGS OR MPI_CXX_LINK_FLAGS)
set_target_properties(${target}
PROPERTIES LINK_FLAGS ${MPI_C_LINK_FLAGS} ${MPI_CXX_LINK_FLAGS})
endif()
endfunction()
#------------------------------------------------------------------------------
# Process subdirectories.
......
......@@ -3,7 +3,19 @@ set (sources
histogram.cpp
histogram.h
)
if(ENABLE_SENSEI)
list(APPEND sources
vtk_histogram.cxx
vtk_histogram.h)
endif()
#-----------------------------------------------------------------------
add_library(histogram STATIC ${sources})
target_compile_definitions(histogram INTERFACE ENABLE_HISTOGRAM)
include(mpi)
mpi_link(histogram)
if(ENABLE_SENSEI)
target_link_libraries(histogram PRIVATE core)
endif()
......@@ -5,6 +5,7 @@
extern "C" {
#endif
/// Simple parallel histogram implementation.
void histogram(MPI_Comm comm, double* data, long sz, int bins);
#ifdef __cplusplus
......
#include "vtk_histogram.h"
#include <vtkAOSDataArrayTemplate.h>
#include <vtkArrayDispatch.h>
#include <vtkSOADataArrayTemplate.h>
namespace {
// Private worker for histogram method. Computes the local histogram on
// array (passed to operator()). To be used with vtkArrayDispatch.
//
// Inputs:
// range: Global range of data
// bins: Number of histogram bins
// array: Local data.
//
// Outputs:
// Histogram: The histogram of the local data.
struct HistogramWorker
{
const double *Range;
int Bins;
std::vector<unsigned int> Histogram;
HistogramWorker(const double *range, int bins) : Range(range), Bins(bins) {}
template <typename ArrayT>
void operator()(ArrayT *array)
{
assert(array);
assert(array->GetNumberOfComponents() == 1);
typedef typename ArrayT::ValueType ValueType;
ValueType width = static_cast<ValueType>((this->Range[1] - this->Range[0]) / this->Bins);
ValueType min = static_cast<ValueType>(this->Range[0]);
vtkIdType numTuples = array->GetNumberOfTuples();
this->Histogram.clear();
// + 1 to store val == max. These will be moved to this last bin before
// returning (Avoids having to branch in the loop below);
this->Histogram.resize(this->Bins + 1, 0);
for (vtkIdType tIdx = 0; tIdx < numTuples; ++tIdx)
{
int bin = static_cast<int>((array->GetComponent(tIdx, 0) - min) / width);
++this->Histogram[bin];
}
// Merge the last two bins (the last is only when val == max)
this->Histogram[this->Bins] += this->Histogram[this->Bins + 1];
this->Histogram.resize(this->Bins);
}
};
} // end anon namespace
void histogram(MPI_Comm comm, vtkDataArray* array, int bins)
{
int rank;
MPI_Comm_rank(comm, &rank);
// find max and min.
double range[2];
array->GetRange(range);
double g_range[2];
// Find the global max/min
MPI_Allreduce(&range[0], &g_range[0], 1, MPI_DOUBLE, MPI_MIN, comm);
MPI_Allreduce(&range[1], &g_range[1], 1, MPI_DOUBLE, MPI_MAX, comm);
// Compute local histogram
HistogramWorker worker(g_range, bins);
// vtkArrayDispatch downcasts the array to a concrete typed implementation
// that provides faster access methods.
if (!vtkArrayDispatch::Dispatch::Execute(array, worker))
{
// This happens if vtkArrayDispatch doesn't know about the array subclass
// in use.
std::cerr << "HistogramWorker dispatch failed on rank " << rank << "!\n";
worker.Histogram.resize(bins, 0);
}
// Global histogram:
std::vector<unsigned int> gHist(bins, 0);
MPI_Reduce(&worker.Histogram[0], &gHist[0], bins, MPI_UNSIGNED, MPI_SUM, 0, comm);
if (rank == 0)
{
std::cout << "Histogram:\n";
double width = (g_range[1] - g_range[0]) / bins;
for (int i = 0; i < bins; ++i)
{
printf(" %f-%f: %d\n", g_range[0] + i*width, g_range[0] + (i+1)*width, gHist[i]);
}
}
}
#ifndef VTK_HISTOGRAM_H
#define VTK_HISTOGRAM_H
#include <mpi.h>
class vtkDataArray;
/// Parallel histogram implementation using array-layout independent
/// VTK's Generic Array infrastructure.
void vtk_histogram(MPI_Comm comm, vtkDataArray* array, int bins);
#endif
......@@ -6,6 +6,8 @@ set (sources
add_executable(3D_Grid ${sources})
target_link_libraries(3D_Grid PRIVATE histogram)
target_compile_definitions(3D_Grid PRIVATE _FILE_OFFSET_BITS=64 _LARGEFILE64_SOURCE)
include(mpi)
mpi_link(3D_Grid)
if(ENABLE_SENSEI)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment