text:
#define _XOPEN_SOURCE 600

#include <stdlib.h>
#include <stdio.h>
#include <mpi.h>
#include <stdint.h>
#include <time.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <errno.h>
#include "fileio.h"
#include "likelihood.h"
#include "master_tasks.h"
#include "mutation.h"
#include "radix.h"
#include "randomisation.h"
#include "resampling.h"
#include "typedefs.h"

static const int ALG = 6;
static const int MASTER = 0;

int
main (
      int argc,
      char *argv[])
{

  struct Process process;
  struct Estimates estimates;
  int pair_rank, src;
  int do_resampling;
  long *inds;
  double *w, *x, *x_resamp;
  double total_w = 1;
  double my_normal_weight;
  double w_pair[2];
  double rsamp_time=0;
  double rsamp_start=0;
  double comm_time=0;
  double comm_start;
  double comm_time2=0;
  double *output;
  double weight_time = 0;
  double mutation_time = 0.0;
  double init_time = 0.0;
   unsigned int cnts[2];
  short counts[2];
  short n_resamp;

  const gsl_rng_type * T;
  double timer_start;

  // Parse commandline arguments
  int run_idx = 0;
  if ( argc > 1 )
    run_idx = atoi( argv[ 1 ] );

  long N = 1000;
  if ( argc > 2 )
    N = atol( argv[ 2 ] );

  int state_dim;
  if ( argc > 3 )
    state_dim = atoi( argv[ 3 ] );
  else {
    perror("State dimesion must be provided");
    return -1;
  }

  double epc_threshold;
  if (argc > 4 ) {
    epc_threshold = atof( argv[ 4 ] );
  } else {
    perror("Reampling threshold must be provided");
    return -1;
  }

  // Initialize the MPI environment
  MPI_Init(NULL, NULL);

  gsl_rng_env_setup();
  T = gsl_rng_default;
  gsl_rng *rgen;
  rgen = gsl_rng_alloc( T );

  // Get the number of processes
  int world_size;
  MPI_Comm_size(MPI_COMM_WORLD, &world_size);

  // Get the rank of the process
  int world_rank;
  MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);

  // Start the timer
  struct timespec start_time;
  clock_gettime(CLOCK_REALTIME, &start_time);

  MPI_Barrier( MPI_COMM_WORLD );
  timer_start = MPI_Wtime();

  long data_length;
  process = readProcessFile( state_dim );
  data_length = (int)process.n;

  // Construct the communication hypercube
  MPI_Comm nthCube;
  int hc_rank;
  int nDim;
  constructHypercube( world_size, &nthCube, &hc_rank, &nDim );
  int S = nDim;

  //---------------
  // Initialisation
  //---------------
  x = malloc( sizeof( double ) * N * state_dim );
  w = malloc( sizeof( double ) * N );

  inds = malloc( sizeof( long ) * N );
  x_resamp = malloc( sizeof( double ) * N * state_dim );

  // Initial distribution setup
  setupInitialSeed( world_rank, rgen );

  // Create initial sample
  createInitialSample( x, N, rgen, state_dim );

  // Setup the randomisation
  set_resamplingseed( world_rank );

  int output_dim = 2 + state_dim + 2 * world_size + world_size * state_dim;
  if ( world_rank == MASTER )
    output = malloc( sizeof( double ) * output_dim * data_length );

  my_normal_weight = 1.0 / world_size;

  MPI_Barrier( MPI_COMM_WORLD );
  init_time += MPI_Wtime() - timer_start;

  // main loop
  for( long n = 0; n < data_length; n++ ) {

    // ------------------------------------------------------------------------------
    MPI_Barrier( MPI_COMM_WORLD );
    timer_start = MPI_Wtime();

    evaluateLikelihoods( x, process, state_dim, data_length, N, n, w,
			 my_normal_weight );

    MPI_Barrier( MPI_COMM_WORLD );
    weight_time += MPI_Wtime() - timer_start ;

    // Compute the relevant estimates for THIS filter only
    estimates = estimateOutput( w, x, world_rank, N, run_idx, state_dim );

    // Combine the results of all filters
    combineEstimates( world_rank, state_dim, estimates, output, n,
		      output_dim, world_size, &total_w, &my_normal_weight,
		      &do_resampling, epc_threshold);

    // -----------------------------
    // RESAMPLING
    // -----------------------------

    // Resampling within processors

    MPI_Barrier( MPI_COMM_WORLD );
    rsamp_start = MPI_Wtime();

    serial_multinomial_resample_sng_noutk( N, w, inds, N );

    for ( long i = 0; i < N; i++)
      for ( int d = 0; d < state_dim; d++ )
	x_resamp[ state_dim * i + d ] = x[ state_dim * inds[ i ] + d ];

    for ( long i = 0; i < N; i++)
      for ( int d = 0; d < state_dim ; d++ )
	x[ state_dim * i + d ] = x_resamp[ state_dim * i + d ];

    // -----------------------------------------------------------------------------
    MPI_Barrier( MPI_COMM_WORLD );
    rsamp_time += MPI_Wtime() - rsamp_start;
    comm_start = MPI_Wtime();

    if ( do_resampling ) {

      // ---------------------------
      // START THE RADIX INTERACTION
      // ---------------------------

      for ( int s = 0; s < S ; s++ ) {

	// Determine the process you need to communicate with
	MPI_Cart_shift( nthCube, s, 1, &src, &pair_rank );

	// The lower ranking process in the pair gets the filter weights and
	// does the island level resampling
	if ( hc_rank < pair_rank ) {

	  MPI_Recv( w_pair + 1, 1, MPI_DOUBLE, pair_rank, 0, nthCube,
		    MPI_STATUS_IGNORE );

	  // The lower ranking filter must send its weight to the paired filter
	  // so that the paired filter can update its weight
	  MPI_Send( &my_normal_weight, 1, MPI_DOUBLE, pair_rank, 1, nthCube );

	  w_pair[ 0 ] = my_normal_weight;

	  // Generate a binomial random variate indicating how many samples
	  // are drawn from each of the paired processes.
	  gsl_ran_multinomial( rgen , ( size_t ) 2, ( unsigned int ) ( 2 ), w_pair,
			       cnts );

	  counts[ 0 ] = ( short ) cnts[ 0 ];
	  counts[ 1 ] = ( short ) cnts[ 1 ];

	  // Send the resampling counts to the paired process
	  MPI_Send( counts + 1 , 1, MPI_SHORT, pair_rank, 3, nthCube );

	  n_resamp = counts[ 0 ];

	} else {

	  // Send my weight to the lower ranking filter in the pair
	  MPI_Send(&my_normal_weight, 1, MPI_DOUBLE, pair_rank, 0, nthCube );

	  // Receive the weight of the paired filter
	  MPI_Recv( w_pair, 1, MPI_DOUBLE, pair_rank, 1, nthCube, MPI_STATUS_IGNORE );

	  // Receive the duplicate count ( 0 or 1 )
	  MPI_Recv( &n_resamp, 1, MPI_SHORT, pair_rank, 3, nthCube, MPI_STATUS_IGNORE );

	  w_pair[ 1 ] = my_normal_weight;

	}

	// -----------------------------------------------------------------------------
	MPI_Barrier( MPI_COMM_WORLD );
	comm_time2 += MPI_Wtime() - comm_start;
	comm_start = MPI_Wtime();

	// Balance the resampled particle counts
	if ( n_resamp > 1 ) {

	  // Send the surplus particles, i.e. the last n_resamp - N of the
	  // resampled particles
	  MPI_Send( x, N * state_dim , MPI_DOUBLE, pair_rank, 2, nthCube );

	} else if ( n_resamp < 1 ) {

	  // Receive compensation for shortage
	  MPI_Recv( x,  N * state_dim, MPI_DOUBLE, pair_rank, 2, nthCube,
		    MPI_STATUS_IGNORE );

	}

	my_normal_weight = ( w_pair[ 0 ] + w_pair[ 1 ] ) / (double ) 2;
      }

    }

    // -----------------------------------------------------------------------------
    MPI_Barrier( MPI_COMM_WORLD );
    comm_time += MPI_Wtime() - comm_start;
    timer_start = MPI_Wtime();

    mutation(N, x, state_dim, rgen);

    // -----------------------------------------------------------------------------
    MPI_Barrier( MPI_COMM_WORLD );
    mutation_time += MPI_Wtime() - timer_start ;

  }

  struct timespec end_time;
  clock_gettime(CLOCK_REALTIME, &end_time);

  createOutput( ALG, world_rank, output, data_length, world_size,
		run_idx, N, state_dim, start_time, end_time, comm_time, rsamp_time,
		comm_time2, weight_time,  mutation_time, init_time,
		(int)( atof( argv[ 4 ] ) * 100 ) );

  // Finalize the MPI environment.
  MPI_Finalize();

  return 0;
}

max_stars_repo_name: heinekmp/AIRPF

max_issues_repo_name: heinekmp/AIRPF

max_forks_repo_name: heinekmp/AIRPF


text:
#ifndef DDM__ALGORITHM__SUMMA_H_
#define DDM__ALGORITHM__SUMMA_H_

#include "../../ddm/Exception.h"
#include "../../ddm/Types.h"
#include "../../ddm/Pattern.h"
#include "../../ddm/Future.h"
#include "../../ddm/algorithm/Copy.h"
#include "../../ddm/util/Trace.h"

#include <utility>

// Prefer MKL if available:
#ifdef DDM_ENABLE_MKL
#include <mkl.h>
#include <mkl_types.h>
#include <mkl_cblas.h>
#include <mkl_blas.h>
#include <mkl_lapack.h>
// BLAS support:
#elif defined(DDM_ENABLE_BLAS)
extern "C"
{
  #include <cblas.h>
}
#endif

#define DDM_ALGORITHM_SUMMA_ASYNC_INIT_PREFETCH

namespace ddm {

namespace internal {

#if defined(DDM_ENABLE_MKL) || defined(DDM_ENABLE_BLAS)
/**
 * Matrix multiplication for local multiplication of matrix blocks via MKL.
 */
template<typename  ValueType>
void mmult_local(
  /// Matrix to multiply, m rows by k columns.
  const ValueType * A,
  /// Matrix to multiply, k rows by n columns.
  const ValueType * B,
  /// Matrix to contain the multiplication result, m rows by n columns.
  ValueType       * C,
  long long         m,
  long long         n,
  long long         k,
  MemArrange        storage);
#else
/**
 * Naive matrix multiplication for local multiplication of matrix blocks,
 * used only for tests and where MKL is not available.
 */
template<typename ValueType>
void mmult_local(
  /// Matrix to multiply, extents n x m
  const ValueType * A,
  /// Matrix to multiply, extents m x p
  const ValueType * B,
  /// Matrix to contain the multiplication result, extents n x p
  ValueType       * C,
  long long         m,
  long long         n,
  long long         p,
  MemArrange        storage)
{
#ifndef DEBUG
  DDM_THROW(
    ddm::exception::RuntimeError,
    "Called fallback implementation of DGEMM (only enabled in Debug)");
#endif
  ValueType c_sum = 0;
  for (auto i = 0; i < n; ++i) {
    // row i = 0...n
    for (auto j = 0; j < p; ++j) {
      // column j = 0...p
      c_sum = C[i * p + j]; // = C[j][i]
      for (auto k = 0; k < m; ++k) {
        // k = 0...m
        auto ik     = i * m + k;
        auto kj     = k * m + j;
        auto value  = A[ik] * B[kj];
        c_sum      += value;
      }
      C[i * p + j] = c_sum; // C[j][i] = c_sum
    }
  }
}
#endif // defined(DDM_ENABLE_MKL) || defined(DDM_ENABLE_BLAS)

} // namespace internal

/// Constraints on pattern partitioning properties of matrix operands passed
/// to \c ddm::summa.
typedef ddm::pattern_partitioning_properties<
            // Block extents are constant for every dimension.
            ddm::pattern_partitioning_tag::rectangular,
            // Identical number of elements in every block.
            ddm::pattern_partitioning_tag::balanced,
            // Matrices must be partitioned in more than one dimension.
            ddm::pattern_partitioning_tag::ndimensional
        > summa_pattern_partitioning_constraints;
/// Constraints on pattern mapping properties of matrix operands passed to
/// \c ddm::summa.
typedef ddm::pattern_mapping_properties<
            // Every unit mapped to more than one block, required for
            // block prefetching to take effect.
            ddm::pattern_mapping_tag::multiple,
            // Number of blocks assigned to a unit may differ.
            ddm::pattern_mapping_tag::unbalanced
        > summa_pattern_mapping_constraints;
/// Constraints on pattern layout properties of matrix operands passed to
/// \c ddm::summa.
typedef ddm::pattern_layout_properties<
            // Elements are contiguous in local memory within single block.
            ddm::pattern_layout_tag::blocked,
            // Local element order corresponds to a logical linearization
            // within single blocks.
            // Required for cache-optimized block matrix multiplication.
            ddm::pattern_layout_tag::linear
        > summa_pattern_layout_constraints;

template<typename MatrixType>
using summa_pattern_constraints =
    typename ddm::pattern_constraints<
        ddm::summa_pattern_partitioning_constraints,
        ddm::summa_pattern_mapping_constraints,
        ddm::summa_pattern_layout_constraints,
        typename MatrixType::pattern_type>;

/**
 * Multiplies two matrices using the SUMMA algorithm.
 * Performs \c (2 * (nunits-1) * nunits^2) async copy operations of
 * submatrices in \c A and \c B.
 *
 * Pseudocode:
 *
 *   C = zeros(n,n)
 *   for k = 1:b:n {            // k increments in steps of blocksize b
 *     u = k:(k+b-1)            // u is [k, k+1, ..., k+b-1]
 *     C = C + A(:,u) * B(u,:)  // Multiply n x b matrix from A with
 *                              // b x p matrix from B
 *   }
 */
template<
  typename MatrixTypeA,
  typename MatrixTypeB,
  typename MatrixTypeC
>
void summa(
  /// Matrix to multiply, extents n x m
  MatrixTypeA & A,
  /// Matrix to multiply, extents m x p
  MatrixTypeB & B,
  /// Matrix to contain the multiplication result, extents n x p,
  /// initialized with zeros
  MatrixTypeC & C)
{
  typedef typename MatrixTypeA::value_type   value_type;
  typedef typename MatrixTypeA::index_type   index_t;
  typedef typename MatrixTypeA::size_type    extent_t;
//typedef typename MatrixTypeA::pattern_type pattern_a_type;
//typedef typename MatrixTypeB::pattern_type pattern_b_type;
//typedef typename MatrixTypeC::pattern_type pattern_c_type;
  typedef std::array<index_t, 2>             coords_t;

  const bool shifted_tiling = ddm::pattern_constraints<
                                ddm::pattern_partitioning_properties<>,
                                ddm::pattern_mapping_properties<
                                  ddm::pattern_mapping_tag::diagonal
                                >,
                                ddm::pattern_layout_properties<>,
                                typename MatrixTypeC::pattern_type
                              >::satisfied::value;
  const bool minimal_tiling = ddm::pattern_constraints<
                                ddm::pattern_partitioning_properties<
                                  ddm::pattern_partitioning_tag::minimal
                                >,
                                ddm::pattern_mapping_properties<>,
                                ddm::pattern_layout_properties<>,
                                typename MatrixTypeC::pattern_type
                              >::satisfied::value;

  static_assert(
      std::is_floating_point<value_type>::value,
      "ddm::summa expects matrix element type double or float");

  DDM_LOG_DEBUG("ddm::summa()");
  // Verify that matrix patterns satisfy pattern constraints:
  if (!ddm::check_pattern_constraints<
         summa_pattern_partitioning_constraints,
         summa_pattern_mapping_constraints,
         summa_pattern_layout_constraints
       >(A.pattern())) {
    DDM_THROW(
      ddm::exception::InvalidArgument,
      "ddm::summa(): "
      "pattern of first matrix argument does not match constraints");
  }
  if (!ddm::check_pattern_constraints<
         summa_pattern_partitioning_constraints,
         summa_pattern_mapping_constraints,
         summa_pattern_layout_constraints
       >(B.pattern())) {
    DDM_THROW(
      ddm::exception::InvalidArgument,
      "ddm::summa(): "
      "pattern of second matrix argument does not match constraints");
  }
  if (!ddm::check_pattern_constraints<
         summa_pattern_partitioning_constraints,
         summa_pattern_mapping_constraints,
         summa_pattern_layout_constraints
       >(C.pattern())) {
    DDM_THROW(
      ddm::exception::InvalidArgument,
      "ddm::summa(): "
      "pattern of result matrix does not match constraints");
  }
  DDM_LOG_TRACE("ddm::summa", "matrix pattern properties valid");

  if (shifted_tiling) {
    DDM_LOG_TRACE("ddm::summa",
                   "using communication scheme for diagonal-shift mapping");
  }
  if (minimal_tiling) {
    DDM_LOG_TRACE("ddm::summa",
                   "using communication scheme for minimal partitioning");
  }
  //    A         B         C
  //  _____     _____     _____
  // |     |   |     |   |     |
  // n     | x m     | = n     |
  // |_ m _|   |_ p _|   |_ p _|
  //
  ddm::Team & team = C.team();
  auto unit_id      = team.myid();
  // Check run-time invariants on pattern instances:
  auto pattern_a    = A.pattern();
  auto pattern_b    = B.pattern();
  auto pattern_c    = C.pattern();
  auto m = pattern_a.extent(0); // number of columns in A, rows in B
#if DDM_ENABLE_TRACE_LOGGING
  auto n = pattern_a.extent(1); // number of rows in A and C
  auto p = pattern_b.extent(0); // number of columns in B and C
#endif
  const ddm::MemArrange memory_order = pattern_a.memory_order();

  DDM_ASSERT_EQ(
    pattern_a.extent(1),
    pattern_b.extent(0),
    "ddm::summa(): "
    "Extents of first operand in dimension 1 do not match extents of "
    "second operand in dimension 0");
  DDM_ASSERT_EQ(
    pattern_c.extent(0),
    pattern_a.extent(0),
    "ddm::summa(): "
    "Extents of result matrix in dimension 0 do not match extents of "
    "first operand in dimension 0");
  DDM_ASSERT_EQ(
    pattern_c.extent(1),
    pattern_b.extent(1),
    "ddm::summa(): "
    "Extents of result matrix in dimension 1 do not match extents of "
    "second operand in dimension 1");

  DDM_LOG_TRACE("ddm::summa", "matrix pattern extents valid");

  // Patterns are balanced, all blocks have identical size:
  auto block_size_m   = pattern_a.block(0).extent(0);
  auto block_size_n   = pattern_b.block(0).extent(1);
  auto block_size_p   = pattern_b.block(0).extent(0);
  auto num_blocks_m   = m / block_size_m;
#if DDM_ENABLE_TRACE_LOGGING
  auto num_blocks_n   = n / block_size_n;
  auto num_blocks_p   = p / block_size_p;
#endif
  // Size of temporary local blocks
  auto block_a_size   = block_size_n * block_size_m;
  auto block_b_size   = block_size_m * block_size_p;
  // Number of units in rows and columns:
  auto teamspec       = C.pattern().teamspec();
  auto unit_ts_coords = teamspec.coords(unit_id);

  DDM_LOG_TRACE("ddm::summa", "blocks:",
                 "m:", num_blocks_m, "*", block_size_m,
                 "n:", num_blocks_n, "*", block_size_n,
                 "p:", num_blocks_p, "*", block_size_p);
  DDM_LOG_TRACE("ddm::summa",
                 "number of units:",
                 "cols:", teamspec.extent(0),
                 "rows:", teamspec.extent(1),
                 "unit team coords:", unit_ts_coords);
  DDM_LOG_TRACE("ddm::summa", "allocating local temporary blocks, sizes:",
                 "A:", block_a_size,
                 "B:", block_b_size);

#ifdef DDM_ENABLE_MKL
  value_type * buf_block_a_get    = (value_type *)(mkl_malloc(
                                      sizeof(value_type) * block_a_size, 64));
  value_type * buf_block_b_get    = (value_type *)(mkl_malloc(
                                      sizeof(value_type) * block_b_size, 64));
  value_type * buf_block_a_comp   = (value_type *)(mkl_malloc(
                                      sizeof(value_type) * block_a_size, 64));
  value_type * buf_block_b_comp   = (value_type *)(mkl_malloc(
                                      sizeof(value_type) * block_b_size, 64));
#else
  value_type * buf_block_a_get    = new value_type[block_a_size];
  value_type * buf_block_b_get    = new value_type[block_b_size];
  value_type * buf_block_a_comp   = new value_type[block_a_size];
  value_type * buf_block_b_comp   = new value_type[block_b_size];
#endif
  // Copy of buffer pointers for swapping, delete[] on swapped pointers tends
  // to crash:
  value_type * local_block_a_get      = buf_block_a_get;
  value_type * local_block_b_get      = buf_block_b_get;
  value_type * local_block_a_comp     = buf_block_a_comp;
  value_type * local_block_b_comp     = buf_block_b_comp;
  value_type * local_block_a_get_bac  = nullptr;
  value_type * local_block_b_get_bac  = nullptr;
  value_type * local_block_a_comp_bac = nullptr;
  value_type * local_block_b_comp_bac = nullptr;
  // -------------------------------------------------------------------------
  // Prefetch blocks from A and B for first local multiplication:
  // -------------------------------------------------------------------------
  // Block coordinates of submatrices of A and B to be prefetched:
  // Local block index of local submatrix of C for multiplication result of
  // blocks to be prefetched:
  auto     l_block_c_get       = C.local.block(0);
  auto     l_block_c_get_view  = l_block_c_get.begin().viewspec();
  index_t  l_block_c_get_row   = l_block_c_get_view.offset(1) / block_size_n;
  index_t  l_block_c_get_col   = l_block_c_get_view.offset(0) / block_size_p;
  // Block coordinates of blocks in A and B to prefetch:
  coords_t block_a_get_coords = coords_t {{ static_cast<index_t>(unit_ts_coords[0]),
                   l_block_c_get_row }};
  coords_t block_b_get_coords = coords_t {{ l_block_c_get_col,
                   static_cast<index_t>(unit_ts_coords[0]) }};
  // Local block index of local submatrix of C for multiplication result of
  // currently prefetched blocks:
  auto     l_block_c_comp      = l_block_c_get;
  auto     l_block_c_comp_view = l_block_c_comp.begin().viewspec();
  index_t  l_block_c_comp_row  = l_block_c_comp_view.offset(1) / block_size_n;
  index_t  l_block_c_comp_col  = l_block_c_comp_view.offset(0) / block_size_p;
  // Prefetch blocks from A and B for computation in next iteration:
  ddm::Future<value_type *> get_a;
  ddm::Future<value_type *> get_b;
  auto block_a      = A.block(block_a_get_coords);
  auto block_a_lptr = block_a.begin().local();
  auto block_b      = B.block(block_b_get_coords);
  auto block_b_lptr = block_b.begin().local();
  DDM_LOG_TRACE("ddm::summa", "summa.prefetch.block.a",
                 "block:", block_a_get_coords,
                 "local:", block_a_lptr != nullptr,
                 "unit:",  block_a.begin().lpos().unit,
                 "view:",  block_a.begin().viewspec());

  ddm::util::Trace trace("SUMMA");

  trace.enter_state("prefetch");
  if (block_a_lptr == nullptr) {
#ifdef DDM_ALGORITHM_SUMMA_ASYNC_INIT_PREFETCH
    get_a = ddm::copy_async(block_a.begin(), block_a.end(),
                             local_block_a_comp);
#else
    ddm::copy(block_a.begin(), block_a.end(),
               local_block_a_comp);
    get_a = ddm::Future<value_type *>(
              [=]() { return local_block_a_comp + block_a.size(); });
#endif
  } else {
    local_block_a_comp_bac = local_block_a_comp;
    local_block_a_comp     = block_a_lptr;
  }

  DDM_LOG_TRACE("ddm::summa", "summa.prefetch.block.b",
                 "block:", block_b_get_coords,
                 "local:", block_b_lptr != nullptr,
                 "unit:",  block_b.begin().lpos().unit,
                 "view:",  block_b.begin().viewspec());
  if (block_b_lptr == nullptr) {
#ifdef DDM_ALGORITHM_SUMMA_ASYNC_INIT_PREFETCH
    get_b = ddm::copy_async(block_b.begin(), block_b.end(),
                             local_block_b_comp);
#else
    ddm::copy(block_b.begin(), block_b.end(),
               local_block_b_comp);
    get_b = ddm::Future<value_type *>(
              [=]() { return local_block_b_comp + block_b.size(); });
#endif
  } else {
    local_block_b_comp_bac = local_block_b_comp;
    local_block_b_comp     = block_b_lptr;
  }
#ifdef DDM_ALGORITHM_SUMMA_ASYNC_INIT_PREFETCH
  if (block_a_lptr == nullptr) {
    DDM_LOG_TRACE("ddm::summa", "summa.prefetch.block.a.wait",
                   "waiting for prefetching of block A from unit",
                   block_a.begin().lpos().unit);
    get_a.wait();
  }
  if (block_b_lptr == nullptr) {
    DDM_LOG_TRACE("ddm::summa", "summa.prefetch.block.b.wait",
                   "waiting for prefetching of block B from unit",
                   block_b.begin().lpos().unit);
    get_b.wait();
  }
#endif
  trace.exit_state("prefetch");

  DDM_LOG_TRACE("ddm::summa", "summa.block",
                 "prefetching of blocks completed");
  // -------------------------------------------------------------------------
  // Iterate local blocks in matrix C:
  // -------------------------------------------------------------------------
  extent_t num_local_blocks_c = pattern_c.local_blockspec().size();

  DDM_LOG_TRACE("ddm::summa", "summa.block.C",
                 "C.num.local.blocks:",  num_local_blocks_c,
                 "C.num.column.blocks:", num_blocks_m);

  for (extent_t lb = 0; lb < num_local_blocks_c; ++lb) {
    // Block coordinates for current block multiplication result:
    l_block_c_comp      = C.local.block(lb);
    l_block_c_comp_view = l_block_c_comp.begin().viewspec();
    l_block_c_comp_row  = l_block_c_comp_view.offset(1) / block_size_n;
    l_block_c_comp_col  = l_block_c_comp_view.offset(0) / block_size_p;
    // Block coordinates for next block multiplication result:
    l_block_c_get       = l_block_c_comp;
    l_block_c_get_view  = l_block_c_comp_view;
    l_block_c_get_row   = l_block_c_get_row;
    l_block_c_get_col   = l_block_c_get_col;
    DDM_LOG_TRACE("ddm::summa", "summa.block.comp", "C.local.block",
                   "l_block_idx:", lb,
                   "row:",         l_block_c_comp_row,
                   "col:",         l_block_c_comp_col,
                   "view:",        l_block_c_comp_view);
    // -----------------------------------------------------------------------
    // Iterate blocks in columns of A / rows of B:
    // -----------------------------------------------------------------------
    for (extent_t block_k = 0; block_k < num_blocks_m; ++block_k) {
      DDM_LOG_TRACE("ddm::summa", "summa.block.k", block_k,
                     "active local block in C:", lb);

      // ---------------------------------------------------------------------
      // Prefetch local copy of blocks from A and B for multiplication in
      // next iteration.
      // ---------------------------------------------------------------------
      bool last = (lb == num_local_blocks_c - 1) &&
                  (block_k == num_blocks_m - 1);
      // Do not prefetch blocks in last iteration:
      if (!last) {
        index_t block_get_k = static_cast<index_t>(block_k + 1);
        block_get_k = (block_get_k + unit_ts_coords[0]) % num_blocks_m;
        // Block coordinate of local block in matrix C to prefetch:
        if (block_k == num_blocks_m - 1) {
          // Prefetch for next local block in matrix C:
          block_get_k        = unit_ts_coords[0];
          l_block_c_get      = C.local.block(lb + 1);
          l_block_c_get_view = l_block_c_get.begin().viewspec();
          l_block_c_get_row  = l_block_c_get_view.offset(1) / block_size_n;
          l_block_c_get_col  = l_block_c_get_view.offset(0) / block_size_p;
        }
        // Block coordinates of blocks in A and B to prefetch:
        block_a_get_coords = coords_t {{ block_get_k, l_block_c_get_row }};
        block_b_get_coords = coords_t {{ l_block_c_get_col, block_get_k }};

        block_a      = A.block(block_a_get_coords);
        block_a_lptr = block_a.begin().local();
        block_b      = B.block(block_b_get_coords);
        block_b_lptr = block_b.begin().local();
        DDM_LOG_TRACE("ddm::summa", "summa.prefetch.block.a",
                       "block:", block_a_get_coords,
                       "local:", block_a_lptr != nullptr,
                       "unit:",  block_a.begin().lpos().unit,
                       "view:",  block_a.begin().viewspec());
        if (block_a_lptr == nullptr) {
          get_a = ddm::copy_async(block_a.begin(),
                                   block_a.end(),
                                   local_block_a_get);
          local_block_a_get_bac = nullptr;
        } else {
          local_block_a_get_bac = local_block_a_get;
          local_block_a_get     = block_a_lptr;
        }
        DDM_LOG_TRACE("ddm::summa", "summa.prefetch.block.b",
                       "block:", block_b_get_coords,
                       "local:", block_b_lptr != nullptr,
                       "unit:",  block_b.begin().lpos().unit,
                       "view:",  block_b.begin().viewspec());
        if (block_b_lptr == nullptr) {
          get_b = ddm::copy_async(block_b.begin(),
                                   block_b.end(),
                                   local_block_b_get);
          local_block_b_get_bac = nullptr;
        } else {
          local_block_b_get_bac = local_block_b_get;
          local_block_b_get     = block_b_lptr;
        }
      } else {
        DDM_LOG_TRACE("ddm::summa", " ->",
                       "last block multiplication",
                       "lb:", lb, "bk:", block_k);
      }
      // ---------------------------------------------------------------------
      // Computation of matrix product of local block matrices:
      // ---------------------------------------------------------------------
      DDM_LOG_TRACE("ddm::summa", "summa.block.comp.multiply",
                     "multiplying local block matrices",
                     "C.local.block.comp:", lb,
                     "view:", l_block_c_comp.begin().viewspec());

      trace.enter_state("multiply");
      ddm::internal::mmult_local<value_type>(
          local_block_a_comp,
          local_block_b_comp,
          l_block_c_comp.begin().local(),
          block_size_m,
          block_size_n,
          block_size_p,
          memory_order);
      trace.exit_state("multiply");

      if (local_block_a_comp_bac != nullptr) {
        local_block_a_comp     = local_block_a_comp_bac;
        local_block_a_comp_bac = nullptr;
      }
      if (local_block_b_comp_bac != nullptr) {
        local_block_b_comp     = local_block_b_comp_bac;
        local_block_b_comp_bac = nullptr;
      }
      if (!last) {
        // -------------------------------------------------------------------
        // Wait for local copies:
        // -------------------------------------------------------------------
        trace.enter_state("prefetch");
        if (block_a_lptr == nullptr) {
          DDM_LOG_TRACE("ddm::summa", "summa.prefetch.block.a.wait",
                         "waiting for prefetching of block A from unit",
                         block_a.begin().lpos().unit);
          get_a.wait();
        }
        if (block_b_lptr == nullptr) {
          DDM_LOG_TRACE("ddm::summa", "summa.prefetch.block.b.wait",
                         "waiting for prefetching of block B from unit",
                         block_b.begin().lpos().unit);
          get_b.wait();
        }
        DDM_LOG_TRACE("ddm::summa", "summa.prefetch.completed",
                       "local copies of next blocks received");
        trace.exit_state("prefetch");

        // -----------------------------------------------------------------
        // Swap communication and computation buffers:
        // -----------------------------------------------------------------
        std::swap(local_block_a_get, local_block_a_comp);
        std::swap(local_block_b_get, local_block_b_comp);
        if (local_block_a_get_bac != nullptr) {
          local_block_a_comp_bac = local_block_a_get_bac;
          local_block_a_get_bac  = nullptr;
        }
        if (local_block_b_get_bac != nullptr) {
          local_block_b_comp_bac = local_block_b_get_bac;
          local_block_b_get_bac  = nullptr;
        }
      }
    }
  } // for lb

  DDM_LOG_TRACE("ddm::summa", "locally completed");
#ifdef DDM_ENABLE_MKL
  mkl_free(buf_block_a_get);
  mkl_free(buf_block_b_get);
  mkl_free(buf_block_a_comp);
  mkl_free(buf_block_b_comp);
#else
  delete[] buf_block_a_get;
  delete[] buf_block_b_get;
  delete[] buf_block_a_comp;
  delete[] buf_block_b_comp;
#endif

  DDM_LOG_TRACE("ddm::summa", "waiting for other units");
  trace.enter_state("barrier");
  C.barrier();
  trace.exit_state("barrier");

  DDM_LOG_TRACE("ddm::summa >", "finished");
}

#ifdef DOXYGEN
/**
 * Function adapter to an implementation of matrix-matrix multiplication
 * (xDGEMM) depending on the matrix distribution patterns.
 *
 * Delegates  \c ddm::mmult<MatrixType>
 * to         \c ddm::summa<MatrixType>
 * if         \c MatrixType::pattern_type
 * satisfies the pattern property constraints of the SUMMA implementation.
 */
template <
  typename MatrixTypeA,
  typename MatrixTypeB,
  typename MatrixTypeC >
void mmult(
  /// Matrix to multiply, extents n x m
  MatrixTypeA & A,
  /// Matrix to multiply, extents m x p
  MatrixTypeB & B,
  /// Matrix to contain the multiplication result, extents n x p,
  /// initialized with zeros
  MatrixTypeC & C);

#else // DOXYGEN

template <
  typename MatrixTypeA,
  typename MatrixTypeB,
  typename MatrixTypeC >
auto
mmult(
  /// Matrix to multiply, extents n x m
  MatrixTypeA & A,
  /// Matrix to multiply, extents m x p
  MatrixTypeB & B,
  /// Matrix to contain the multiplication result, extents n x p,
  /// initialized with zeros
  MatrixTypeC & C)
  -> typename std::enable_if<
                summa_pattern_constraints<MatrixTypeA>::satisfied::value &&
                summa_pattern_constraints<MatrixTypeB>::satisfied::value &&
                summa_pattern_constraints<MatrixTypeC>::satisfied::value,
                void
              >::type {
  ddm::summa(A, B, C);
}

#endif // DOXYGEN

} // namespace ddm

#endif // DDM__ALGORITHM__SUMMA_H_

ext: h

max_stars_repo_name: caixiuhong/Develop-MCCE

max_issues_repo_name: caixiuhong/Develop-MCCE

max_forks_repo_name: caixiuhong/Develop-MCCE


text:
#include <stdio.h>
#include <gsl/gsl_rng.h>

/* https://www.gnu.org/software/gsl/doc/html/rng.html */

/*   run as
	 ./a.out 10000 rnd_u.csv
*/

gsl_rng * r;  /* global generator */

int main(int argc, char **argv) {
  const gsl_rng_type * T;
  gsl_rng * r;

  if( argc > 1) {

	FILE *fp; /* file container for random numbers */
	double i, n = atoi( argv[1] );
	char *fname = argv[2];

	fp = fopen(fname,"w");
	if( fp == NULL) {
	  printf("file can't be opened\n");
	  exit(1);
	}

	fprintf( fp, "u\n");

	gsl_rng_env_setup();

	T = gsl_rng_default;
	r = gsl_rng_alloc (T);

	for (i = 0; i < n; i++)
	  {
		double u = gsl_rng_uniform (r);
		/* printf ("%.5f\t%s\n", u, fname); */
		fprintf( fp, "%.6f\n", u);
	  }

	fclose(fp);

	gsl_rng_free (r);
  }

  return EXIT_SUCCESS;
}

ext: c


text:
#include <stdio.h>
#include <gsl/gsl_sf_gamma.h>

int main() {
    printf("hello, world!\n");
    printf("%f\n", gsl_sf_gamma(1.0));
    printf("%f\n", gsl_sf_gamma(1.5));
    printf("%f\n", gsl_sf_gamma(2.0));
    printf("%f\n", gsl_sf_gamma(2.5));
    printf("%f\n", gsl_sf_gamma(3.0));
}

max_stars_repo_name: berquist/eg

max_issues_repo_name: berquist/eg

max_forks_repo_name: berquist/eg


text:
/* monte/test.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Michael Booth
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <math.h>
#include <stdio.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>

#include <gsl/gsl_rng.h>
#include <gsl/gsl_monte_plain.h>
#include <gsl/gsl_monte_miser.h>
#include <gsl/gsl_monte_vegas.h>

#define CONSTANT
#define PRODUCT
#define GAUSSIAN
#define DBLGAUSSIAN
#define TSUDA

#define PLAIN
#define MISER
#define VEGAS

double xl[11]  = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
double xu[11]  = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 };
double xu2[11] = { 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2 };
double xu3[2]  = { GSL_DBL_MAX, GSL_DBL_MAX };

double fconst (double x[], size_t d, void *params);
double f0 (double x[], size_t d, void *params);
double f1 (double x[], size_t d, void *params);
double f2 (double x[], size_t d, void *params);
double f3 (double x[], size_t d, void *params);

void my_error_handler (const char *reason, const char *file,
		       int line, int err);

struct problem {
  gsl_monte_function * f;
  double * xl;
  double * xu;
  size_t dim;
  size_t calls;
  double expected_result;
  double expected_error;
  char * description;
} ;

gsl_monte_function
make_function (double (*f)(double *, size_t, void *), size_t d, void * p);

gsl_monte_function
make_function (double (*f)(double *, size_t, void *), size_t d, void * p)
{
  gsl_monte_function f_new;

  f_new.f = f;
  f_new.dim = d;
  f_new.params = p;

  return f_new;
}

void
add (struct problem * problems, int * n,
     gsl_monte_function * f, double xl[], double xu[], size_t dim, size_t calls,
     double result, double err, char * description);

void
add (struct problem * problems, int * n,
     gsl_monte_function * f, double xl[], double xu[], size_t dim, size_t calls,
     double result, double err, char * description)
{
  int i = *n;
  problems[i].f = f;
  problems[i].xl = xl;
  problems[i].xu = xu;
  problems[i].dim = dim;
  problems[i].calls = calls;
  problems[i].expected_result = result;
  problems[i].expected_error = err;
  problems[i].description = description;
  (*n)++;
}

#define TRIALS 10

int
main (void)
{
  double result[TRIALS], error[TRIALS];
  double a = 0.1;
  double c = (1.0 + sqrt (10.0)) / 9.0;

  gsl_monte_function Fc = make_function(&fconst, 0, 0);
  gsl_monte_function F0 = make_function(&f0, 0, &a);
  gsl_monte_function F1 = make_function(&f1, 0, &a);
  gsl_monte_function F2 = make_function(&f2, 0, &a);
  gsl_monte_function F3 = make_function(&f3, 0, &c);

  /* The relationship between the variance of the function itself, the
     error on the integral and the number of calls is,

     sigma = sqrt(variance/N)

     where the variance is the <(f - <f>)^2> where <.> denotes the
     volume average (integral over the integration region divided by
     the volume) */

  int n = 0;
  struct problem * I;
  struct problem problems[256];

#ifdef CONSTANT
    /* variance(Fc) = 0 */

    add(problems,&n, &Fc, xl, xu,  1, 1000, 1.0, 0.0, "constant, 1d");
    add(problems,&n, &Fc, xl, xu,  2, 1000, 1.0, 0.0, "constant, 2d");
    add(problems,&n, &Fc, xl, xu,  3, 1000, 1.0, 0.0, "constant, 3d");
    add(problems,&n, &Fc, xl, xu,  4, 1000, 1.0, 0.0, "constant, 4d");
    add(problems,&n, &Fc, xl, xu,  5, 1000, 1.0, 0.0, "constant, 5d");
    add(problems,&n, &Fc, xl, xu,  6, 1000, 1.0, 0.0, "constant, 6d");
    add(problems,&n, &Fc, xl, xu,  7, 1000, 1.0, 0.0, "constant, 7d");
    add(problems,&n, &Fc, xl, xu,  8, 1000, 1.0, 0.0, "constant, 8d");
    add(problems,&n, &Fc, xl, xu,  9, 1000, 1.0, 0.0, "constant, 9d");
    add(problems,&n, &Fc, xl, xu, 10, 1000, 1.0, 0.0, "constant, 10d");
#endif

#ifdef PRODUCT
    /* variance(F0) = (4/3)^d - 1 */

    add(problems,&n, &F0, xl, xu,  1, 3333,   1.0, 0.01, "product, 1d" );
    add(problems,&n, &F0, xl, xu,  2, 7777,   1.0, 0.01, "product, 2d" );
    add(problems,&n, &F0, xl, xu,  3, 13703,  1.0, 0.01, "product, 3d" );
    add(problems,&n, &F0, xl, xu,  4, 21604,  1.0, 0.01, "product, 4d" );
    add(problems,&n, &F0, xl, xu,  5, 32139,  1.0, 0.01, "product, 5d" );
    add(problems,&n, &F0, xl, xu,  6, 46186,  1.0, 0.01, "product, 6d" );
    add(problems,&n, &F0, xl, xu,  7, 64915,  1.0, 0.01, "product, 7d" );
    add(problems,&n, &F0, xl, xu,  8, 89887,  1.0, 0.01, "product, 8d" );
    add(problems,&n, &F0, xl, xu,  9, 123182, 1.0, 0.01, "product, 9d" );
    add(problems,&n, &F0, xl, xu, 10, 167577, 1.0, 0.01, "product, 10d" );
#endif

#ifdef GAUSSIAN
    /* variance(F1) = (1/(a sqrt(2 pi)))^d - 1 */

    add(problems,&n, &F1, xl, xu,  1, 298,      1.0, 0.1, "gaussian, 1d" );
    add(problems,&n, &F1, xl, xu,  2, 1492,     1.0, 0.1, "gaussian, 2d" );
    add(problems,&n, &F1, xl, xu,  3, 6249,     1.0, 0.1, "gaussian, 3d" );
    add(problems,&n, &F1, xl, xu,  4, 25230,    1.0, 0.1, "gaussian, 4d" );
    add(problems,&n, &F1, xl, xu,  5, 100953,   1.0, 0.1, "gaussian, 5d" );
    add(problems,&n, &F1, xl, xu,  6, 44782,    1.0, 0.3, "gaussian, 6d" );
    add(problems,&n, &F1, xl, xu,  7, 178690,   1.0, 0.3, "gaussian, 7d" );
    add(problems,&n, &F1, xl, xu,  8, 712904,   1.0, 0.3, "gaussian, 8d" );
    add(problems,&n, &F1, xl, xu,  9, 2844109,  1.0, 0.3, "gaussian, 9d" );
    add(problems,&n, &F1, xl, xu, 10, 11346390, 1.0, 0.3, "gaussian, 10d" );
#endif

#ifdef DBLGAUSSIAN
    /* variance(F2) = 0.5 * (1/(a sqrt(2 pi)))^d - 1 */

    add(problems,&n, &F2, xl, xu,  1, 9947,  1.0, 0.01, "double gaussian, 1d" );
    add(problems,&n, &F2, xl, xu,  2, 695,   1.0, 0.1, "double gaussian, 2d" );
    add(problems,&n, &F2, xl, xu,  3, 3074,  1.0, 0.1, "double gaussian, 3d" );
    add(problems,&n, &F2, xl, xu,  4, 12565, 1.0, 0.1, "double gaussian, 4d" );
    add(problems,&n, &F2, xl, xu,  5, 50426, 1.0, 0.1, "double gaussian, 5d" );
    add(problems,&n, &F2, xl, xu,  6, 201472, 1.0, 0.1, "double gaussian, 6d" );
    add(problems,&n, &F2, xl, xu,  7, 804056, 1.0, 0.1, "double gaussian, 7d" );
    add(problems,&n, &F2, xl, xu,  8, 356446, 1.0, 0.3, "double gaussian, 8d" );
    add(problems,&n, &F2, xl, xu,  9, 1422049, 1.0, 0.3, "double gaussian, 9d" );
    add(problems,&n, &F2, xl, xu, 10, 5673189, 1.0, 0.3, "double gaussian, 10d" );
#endif

#ifdef TSUDA
    /* variance(F3) = ((c^2 + c + 1/3)/(c(c+1)))^d - 1 */

    add(problems,&n, &F3, xl, xu,  1, 4928,   1.0, 0.01, "tsuda function, 1d" );
    add(problems,&n, &F3, xl, xu,  2, 12285,  1.0, 0.01, "tsuda function, 2d" );
    add(problems,&n, &F3, xl, xu,  3, 23268,  1.0, 0.01, "tsuda function, 3d" );
    add(problems,&n, &F3, xl, xu,  4, 39664,  1.0, 0.01, "tsuda function, 4d" );
    add(problems,&n, &F3, xl, xu,  5, 64141,  1.0, 0.01, "tsuda function, 5d" );
    add(problems,&n, &F3, xl, xu,  6, 100680, 1.0, 0.01, "tsuda function, 6d" );
    add(problems,&n, &F3, xl, xu,  7, 155227, 1.0, 0.01, "tsuda function, 7d" );
    add(problems,&n, &F3, xl, xu,  8, 236657, 1.0, 0.01, "tsuda function, 8d" );
    add(problems,&n, &F3, xl, xu,  9, 358219, 1.0, 0.01, "tsuda function, 9d" );
    add(problems,&n, &F3, xl, xu, 10, 539690, 1.0, 0.01, "tsuda function, 10d" );
#endif

    add(problems,&n,   0,  0,  0, 0,    0,   0,   0, 0  );

  /* gsl_set_error_handler (&my_error_handler); */
  gsl_ieee_env_setup ();
  gsl_rng_env_setup ();

#ifdef A
  printf ("testing allocation/input checks\n");

  status = gsl_monte_plain_validate (s, xl, xu3, 1, 1);
  gsl_test (status != 0, "error if limits too large");
  status = gsl_monte_plain_validate (s, xl, xu, 0, 10);
  gsl_test (status != 0, "error if num_dim = 0");
  status = gsl_monte_plain_validate (s, xl, xu, 1, 0);
  gsl_test (status != 0, "error if calls = 0");
  status = gsl_monte_plain_validate (s, xu, xl, 1, 10);
  gsl_test (status != 0, "error if xu < xl");
#endif

#ifdef PLAIN
#define NAME "plain"
#define MONTE_STATE gsl_monte_plain_state
#define MONTE_ALLOC gsl_monte_plain_alloc
#define MONTE_INTEGRATE gsl_monte_plain_integrate
#define MONTE_FREE gsl_monte_plain_free
#define MONTE_SPEEDUP 1
#define MONTE_ERROR_TEST(err,expected) gsl_test_factor(err,expected, 5.0, NAME ", %s, abserr[%d]", I->description, i)
#include "test_main.c"
#undef NAME
#undef MONTE_STATE
#undef MONTE_ALLOC
#undef MONTE_INTEGRATE
#undef MONTE_FREE
#undef MONTE_ERROR_TEST
#undef MONTE_SPEEDUP
#endif

#ifdef MISER
#define NAME "miser"
#define MONTE_STATE gsl_monte_miser_state
#define MONTE_ALLOC gsl_monte_miser_alloc
#define MONTE_INTEGRATE gsl_monte_miser_integrate
#define MONTE_FREE gsl_monte_miser_free
#define MONTE_SPEEDUP 2
#define MONTE_ERROR_TEST(err,expected) gsl_test(err > 5.0 * expected, NAME ", %s, abserr[%d] (obs %g vs plain %g)", I->description, i, err, expected)
#include "test_main.c"
#undef NAME
#undef MONTE_STATE
#undef MONTE_ALLOC
#undef MONTE_INTEGRATE
#undef MONTE_FREE
#undef MONTE_ERROR_TEST
#undef MONTE_SPEEDUP
#endif

#ifdef VEGAS
#define NAME "vegas"
#define MONTE_STATE gsl_monte_vegas_state
#define MONTE_ALLOC gsl_monte_vegas_alloc
#define MONTE_INTEGRATE(f,xl,xu,dim,calls,r,s,res,err) { gsl_monte_vegas_integrate(f,xl,xu,dim,calls,r,s,res,err) ; if (s->chisq < 0.5 || s->chisq > 2) gsl_monte_vegas_integrate(f,xl,xu,dim,calls,r,s,res,err); }
#define MONTE_FREE gsl_monte_vegas_free
#define MONTE_SPEEDUP 3
#define MONTE_ERROR_TEST(err,expected) gsl_test(err > 3.0 * (expected == 0 ? 1.0/(I->calls/MONTE_SPEEDUP) : expected), NAME ", %s, abserr[%d] (obs %g vs exp %g)", I->description, i, err, expected)
#include "test_main.c"
#undef NAME
#undef MONTE_STATE
#undef MONTE_ALLOC
#undef MONTE_INTEGRATE
#undef MONTE_FREE
#undef MONTE_ERROR_TEST
#undef MONTE_SPEEDUP
#endif

  exit (gsl_test_summary ());
}

/* Simple constant function */
double
fconst (double x[], size_t num_dim, void *params)
{
  return 1;
}

/* Simple product function */
double
f0 (double x[], size_t num_dim, void *params)
{
  double prod = 1.0;
  unsigned int i;

  for (i = 0; i < num_dim; ++i)
    {
      prod *= 2.0 * x[i];
    }

  return prod;
}

/* Gaussian centered at 1/2. */

double
f1 (double x[], size_t num_dim, void *params)
{
  double a = *(double *)params;
  double sum = 0.;

  unsigned int i;
  for (i = 0; i < num_dim; i++)
    {
      double dx = x[i] - 0.5;
      sum += dx * dx;
    }
  return (pow (M_2_SQRTPI / (2. * a), (double) num_dim) *
	  exp (-sum / (a * a)));
}

/* double gaussian */
double
f2 (double x[], size_t num_dim, void *params)
{
  double a = *(double *)params;
  double sum1 = 0.;
  double sum2 = 0.;

  unsigned int i;
  for (i = 0; i < num_dim; i++)
    {
      double dx1 = x[i] - 1. / 3.;
      double dx2 = x[i] - 2. / 3.;
      sum1 += dx1 * dx1;
      sum2 += dx2 * dx2;
    }
  return 0.5 * pow (M_2_SQRTPI / (2. * a), num_dim)
    * (exp (-sum1 / (a * a)) + exp (-sum2 / (a * a)));
}

/* Tsuda's example */
double
f3 (double x[], size_t num_dim, void *params)
{
  double c = *(double *)params;

  double prod = 1.;

  unsigned int i;

  for (i = 0; i < num_dim; i++)
    {
      prod *= c / (c + 1) * pow((c + 1) / (c + x[i]), 2.0);
    }

  return prod;
}

void
my_error_handler (const char *reason, const char *file, int line, int err)
{
  if (0)
    printf ("(caught [%s:%d: %s (%d)])\n", file, line, reason, err);
}


text:
#include <math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_sf_hyperg.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>

#include <fastpm/libfastpm.h>

void
fastpm_horizon_init(FastPMHorizon * horizon, FastPMCosmology * cosmology)
{
    gsl_set_error_handler_off(); // Turn off GSL error handler

    horizon->cosmology = cosmology;
    horizon->size = 8192;
    horizon->da = 1.0 / (horizon->size - 1);
    int i;
    for (i = 0; i < horizon->size; i ++) {
        double a = 1.0 * i / (horizon->size - 1);
        horizon->xi_a[i] = HubbleDistance * ComovingDistance(a, horizon->cosmology);
        FastPMGrowthInfo gi;
        fastpm_growth_info_init(&gi, a, horizon->cosmology);
        horizon->growthfactor_a[i] = gi.D1;
    }
}

void
fastpm_horizon_destroy(FastPMHorizon * horizon)
{
}

double
HorizonDistance(double a, FastPMHorizon * horizon)
{
    double x = a * (horizon->size - 1);
    int l = floor(x);
    int r = l + 1;
    if(r >= horizon->size) {
        return horizon->xi_a[horizon->size - 1];
    }
    if(l <= 0) {
        return horizon->xi_a[0];
    }
    return horizon->xi_a[l] * (r - x)
         + horizon->xi_a[r] * (x - l);
}

double
HorizonGrowthFactor(double a, FastPMHorizon * horizon)
{
    double x = a * (horizon->size - 1);
    int l = floor(x);
    int r = l + 1;
    if(r >= horizon->size) {
        return horizon->growthfactor_a[horizon->size - 1];
    }
    if(l <= 0) {
        return horizon->growthfactor_a[0];
    }
    return horizon->growthfactor_a[l] * (r - x)
         + horizon->growthfactor_a[r] * (x - l);
}

void *
fastpm_horizon_solve_start()
{
    const gsl_root_fsolver_type *T = gsl_root_fsolver_brent;
    return gsl_root_fsolver_alloc(T);
}

void
fastpm_horizon_solve_end(void * context)
{
    gsl_root_fsolver_free(context);
}

int
fastpm_horizon_solve(FastPMHorizon * horizon,
    void * context,
    double * solution,
    double a_i, double a_f,
    double (*func)(double a, void * userdata),
    void * userdata)
{

    int status;
    int iter = 0, max_iter;
    double r, x_lo=a_i, x_hi=a_f, eps;

    /* Reorganize to struct later */
    max_iter = 20;
    eps = 1e-5;

    gsl_function F;

    F.function = func;
    F.params = userdata;

    status = gsl_root_fsolver_set(context, &F, x_lo, x_hi);

    if(status == GSL_EINVAL || status == GSL_EDOM) {
        /** Error in value or out of range **/
        return 0;
    }

    do
    {
        iter++;
        //
        // Debug printout #1
        //if(iter == 1) {
        //fastpm_info("ID | [x_lo, x_hi] | r | funct(r) | x_hi - x_lo\n");
        //}
        //

        status = gsl_root_fsolver_iterate(context);
        r = gsl_root_fsolver_root(context);

        x_lo = gsl_root_fsolver_x_lower(context);
        x_hi = gsl_root_fsolver_x_upper(context);

        status = gsl_root_test_interval(x_lo, x_hi, eps, 0.0);
        //
        //Debug printout #2
        //fastpm_info("%5d [%.7f, %.7f] %.7f %.7f %.7f\n", iter, x_lo, x_hi, r, funct(r, &params), x_hi - x_lo);
        //

        if(status == GSL_SUCCESS || iter == max_iter ) {
            *solution = r;
            //
            // Debug printout #3.1
            //fastpm_info("fastpm_lc_intersect() called with parameters %.7f and %.7f, returned status %d.\n\n", a, b, 1);
            //
            return 1;
        }
    }
    while (status == GSL_CONTINUE);
    //
    // Debug printout #3.2
    //fastpm_info("fastpm_lc_intersect() called with parameters %.7f and %.7f, returned status %d.\n\n", a, b, 0);
    //

    return 0;

}

max_stars_repo_name: sbird/FastPMRunner

max_issues_repo_name: sbird/FastPMRunner

max_forks_repo_name: sbird/FastPMRunner


text:
#include <viaio/Vlib.h>
#include <viaio/VImage.h>
#include <viaio/mu.h>
#include <viaio/option.h>

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#include <gsl/gsl_cblas.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_wavelet.h>
#include <gsl/gsl_wavelet2d.h>

extern void gsl_sort_vector(gsl_vector *);
#define ABS(x) ((x) > 0 ? (x) : -(x))

void WriteTest(gsl_matrix *data)
{
  FILE *fp=0;
  VAttrList alist;
  VImage tmp;
  double u;
  int i,j;

  tmp = VCreateImage(1,data->size1,data->size2,VFloatRepn);
  for (i = 0; i < data->size1; i++) {
    for (j = 0; j < data->size2; j++) {
      u = gsl_matrix_get(data,i,j);
      VPixel(tmp,0,i,j,VFloat) = u;
    }
  }
  alist = VCreateAttrList();
  VAppendAttr(alist,"image",NULL,VImageRepn,tmp);
  fp = fopen("qtest.v","w");
  VWriteFile (fp,alist);
  fclose(fp);
  exit(0);
}

double
estimateSigma(gsl_matrix *data)
{
  int i,j,i2,j2,n,nn;
  gsl_vector *vec=NULL;
  double median=0,med0,u,xmax,tiny=1.0e-6;

  i2 = data->size1/2;
  j2 = data->size2/2;

  nn = 3 * (i2+1) * (j2+1);
  vec = gsl_vector_calloc(nn);

  /*
  ** first pass, median, no absolute values
  */
  xmax = VRepnMaxValue(VDoubleRepn);
  gsl_vector_set_all(vec,xmax);

  n = 0;
  for (i=0; i<i2; i++) {
    for (j=j2; j<data->size2; j++) {
      u = gsl_matrix_get(data,i,j);
      if (ABS(u) > tiny) gsl_vector_set(vec,n++,u);
    }
  }

  for (i=i2; i<data->size1; i++) {
    for (j=0; j<j2; j++) {
      u = gsl_matrix_get(data,i,j);
      if (ABS(u) > tiny) gsl_vector_set(vec,n++,u);
    }
  }

  for (i=i2; i<data->size1; i++) {
    for (j=j2; j<data->size2; j++) {
      u = gsl_matrix_get(data,i,j);
      if (ABS(u) > tiny) gsl_vector_set(vec,n++,u);
    }
  }

  gsl_sort_vector(vec);
  u  = (double)n/2.0;
  nn = (int)u;
  med0 = 0.5*(gsl_vector_get(vec,nn) + gsl_vector_get(vec,nn+1));

  /*
  ** median, 2nd pass
  */
  xmax = VRepnMaxValue(VDoubleRepn);
  gsl_vector_set_all(vec,xmax);

  n = 0;
  for (i=0; i<i2; i++) {
    for (j=j2; j<data->size2; j++) {
      u = gsl_matrix_get(data,i,j);
      if (ABS(u) <= tiny) continue;
      u = ABS(u - med0);
      gsl_vector_set(vec,n++,u);
    }
  }

  for (i=i2; i<data->size1; i++) {
    for (j=0; j<j2; j++) {
      u = gsl_matrix_get(data,i,j);
      if (ABS(u) <= tiny) continue;
      u = ABS(u - med0);
      gsl_vector_set(vec,n++,u);
    }
  }

  for (i=i2; i<data->size1; i++) {
    for (j=j2; j<data->size2; j++) {
      u = gsl_matrix_get(data,i,j);
      if (ABS(u) <= tiny) continue;
      u = ABS(u - med0);
      gsl_vector_set(vec,n++,u);
    }
  }

  gsl_sort_vector(vec);
  u  = (double)n/2.0;
  nn = (int)u;
  median = 0.5*(gsl_vector_get(vec,nn) + gsl_vector_get(vec,nn+1));

  return (median/0.455);
}

double sgn(double x)
{
  if (x < 0) return -1;
  else if (x > 0) return 1;
  else return 0;
}

void
Filtering(gsl_matrix *data, int type, double intensity, double sigma)
{
  int i,j,n;
  double d,dd,alpha;

  n = data->size1;

  double tau = sigma * intensity * sqrt(2 * log10((double) (n*n)));

  for (i=0; i < data->size1; i++) {
    for (j=0; j < data->size2; j++) {
      d = gsl_matrix_get(data,i,j);

      /* Wiener-Filter */
      if (type == 0) {
	alpha = ((d * d) - intensity * (sigma * sigma)) / (d * d);
	if (alpha < 0) alpha = 0;
	d *= alpha;
      }

      /* Soft-Threshold */
      else if (type == 1) {
	if (ABS(d) < tau)
	  dd = 0;
	else
	  dd = sgn(d) * (ABS(d) - tau);
	d = dd;
      }

      /* Hard-Threshold */
      else if (type == 2) {
	if (ABS(d) < tau) d = 0;
      }

      gsl_matrix_set(data,i,j,d);
    }
  }
}

void
Shift(gsl_matrix *data,gsl_matrix *shiftdata,int ishift,int jshift)
{
  int i,j,ii,jj,n;
  double u;

  gsl_matrix_set_zero(shiftdata);

  n = data->size1;
  for (i=0; i<n; i++) {
    ii = i+ishift;
    if (ii < 0 || ii >= n) continue;
    for (j=0; j<n; j++) {
      jj = j+jshift;
      if (jj < 0 || jj >= n) continue;
      u = gsl_matrix_get(data,i,j);
      gsl_matrix_set(shiftdata,ii,jj,u);
    }
  }
}

void
AddShift(gsl_matrix *data,gsl_matrix *shiftdata,int ishift,int jshift)
{
  int i,j,ii,jj,n;
  double u,v;
  n = data->size1;

  for (i=0; i<n; i++) {
    ii = i+ishift;
    if (ii < 0 || ii >= n) continue;
    for (j=0; j<n; j++) {
      jj = j+jshift;
      if (jj < 0 || jj >= n) continue;
      u = gsl_matrix_get(data,i,j);
      v = gsl_matrix_get(shiftdata,ii,jj);
      u = u+v;
      gsl_matrix_set(data,i,j,u);
    }
  }
}

VImage
VWavelets(VImage src,VImage dest,int filtertype,int wtype,int nshift,double level)
{
  static gsl_wavelet_workspace *work=NULL;
  static gsl_wavelet *w=NULL;
  static gsl_matrix *data=NULL,*tmp=NULL,*shiftdata=NULL;
  int nslices,nrows,ncols,i,j,slice,n,size,dim;
  double u,sum1,sum2,nx,mean,sig,sigma=0,avesig;
  extern VFloat VGetPixelValue (VImage,int);

  /*
  ** read image info
  */
  nslices = VImageNBands(src);
  nrows = VImageNRows(src);
  ncols = VImageNColumns(src);

  if (dest == NULL)
    dest = VCreateImage(nslices,nrows,ncols,VPixelRepn(src));
  VCopyImageAttrs (src, dest);
  VFillImage(dest,VAllBands,0);

  size = nrows;
  if (nrows < ncols) size = ncols;
  n = (int) ceil(log10((double) size) / log10(2.0));
  dim = (int) pow(2.0, (double) n);

  /*
  ** ini wavelets transform
  */
  if (work == NULL) {

    switch (wtype) {
    case 0:
      w = gsl_wavelet_alloc (gsl_wavelet_haar_centered,2);
      break;
    case 1:
      w = gsl_wavelet_alloc (gsl_wavelet_daubechies_centered,4);
      break;
    case 2:
      w = gsl_wavelet_alloc (gsl_wavelet_daubechies_centered,6);
      break;
    case 3:
      w = gsl_wavelet_alloc (gsl_wavelet_daubechies_centered,8);
      break;
    case 4:
      w = gsl_wavelet_alloc (gsl_wavelet_bspline_centered,103);
      break;
    case 5:
      w = gsl_wavelet_alloc (gsl_wavelet_bspline_centered,202);
      break;
    default:
      VError("illegal wavelet type");
    }

    work = gsl_wavelet_workspace_alloc (dim);
    data = gsl_matrix_calloc(dim,dim);
    shiftdata = gsl_matrix_calloc(dim,dim);
    tmp = gsl_matrix_calloc(dim,dim);
  }

  avesig = 0;
  for (slice=0; slice<nslices; slice++) {

    /*
    ** read data
    */
    sum1 = sum2 = nx = 0;
    for (i=0; i<nrows; i++) {
      for (j=0; j<ncols; j++) {
	u = VGetPixel(src,slice,i,j);
	sum1 += u;
	sum2 += u*u;
	nx++;
      }
    }
    mean = sum1/nx;
    sig = sqrt((sum2 - nx * mean * mean) / (nx - 1.0));

    if (sig < 1.0e-3) {
      VFillImage(dest,slice,0);
      continue;
    }

    gsl_matrix_set_zero(data);
    for (i=0; i<nrows; i++) {
      for (j=0; j<ncols; j++) {
	u = VGetPixel(src,slice,i,j);
	u = (u-mean)/sig;
	gsl_matrix_set(data,i,j,u);
      }
    }
    gsl_matrix_memcpy(tmp,data);

    /*
    ** get sigma estimate (Donoho)
    */
    if (gsl_wavelet2d_nstransform_matrix_forward(w,data,work) != GSL_SUCCESS)
      VError(" gsl_wavelet2d_nstransform_matrix_forward failed");

    sigma = estimateSigma(data);
    /* fprintf(stderr," slice: %3d  %f\n",slice,sigma); */
    avesig += sigma;

    Filtering(data,filtertype,level,sigma);
    gsl_wavelet2d_nstransform_matrix_inverse(w,data,work);

    /*
    ** translation invariance, shift 4/8 algorithm
    */
    nx = 1;
    for (i=-nshift; i<=nshift; i++) {
      for (j=-nshift; j<=nshift; j++) {
	if (i == 0 && j == 0) continue;
	Shift(tmp,shiftdata,i,j);
	gsl_wavelet2d_nstransform_matrix_forward(w,shiftdata,work);
	Filtering(shiftdata,filtertype,level,sigma);
	gsl_wavelet2d_nstransform_matrix_inverse(w,shiftdata,work);
	AddShift(data,shiftdata,i,j);
	nx++;
      }
    }
    gsl_matrix_scale(data,1.0/nx);

    /*
    ** output
    */
    for (i=0; i<nrows; i++) {
      for (j=0; j<ncols; j++) {
	u = gsl_matrix_get(data,i,j);
	u = sig*u + mean;
	if (u < VPixelMinValue(dest)) u = VPixelMinValue(dest);
	if (u > VPixelMaxValue(dest)) u = VPixelMaxValue(dest);
	VSetPixel(dest,slice,i,j,u);
      }
    }
  }
  avesig /= (double)nslices;
  /* fprintf(stderr," kappa: %f\n",sqrt(avesig) * 1.75); */

  /*
  gsl_matrix_free(shiftdata);
  gsl_matrix_free(tmp);
  gsl_wavelet_free (w);
  gsl_wavelet_workspace_free (work);
  free (data);
  */
  return dest;
}

VDictEntry WDict[] = {
  { "haar", 0 },
  { "daub4", 1 },
  { "daub6", 2 },
  { "daub8", 3 },
  { "bspline103", 4 },
  { "bspline202", 5 },
  { NULL }
};

int main (int argc,char *argv[])
{
  static VShort ftype  = 0;
  static VShort wtype  = 0;
  static VShort nshift = 1;
  static VDouble level = 0.5;
  static VOptionDescRec  options[] = {
    {"wavelet",VShortRepn,1,(VPointer) &wtype,VOptionalOpt,WDict,"wavelet type"},
    {"filter",VShortRepn,1,(VPointer) &ftype,VOptionalOpt,NULL,"filter type"},
    {"shift",VShortRepn,1,(VPointer) &nshift,VOptionalOpt,NULL,"shift window size"},
    {"level",VDoubleRepn,1,(VPointer) &level,VOptionalOpt,NULL,"level"}
  };
  FILE *in_file,*out_file;
  VAttrList list=NULL;
  VAttrListPosn posn;
  VImage src=NULL,dest=NULL;
  int n=0,nimages=0;
  char *prg=GetLipsiaName("vdenoise");
  fprintf (stderr, "%s\n", prg);

  VParseFilterCmd (VNumber (options),options,argc,argv,&in_file,&out_file);

  if (! (list = VReadFile (in_file, NULL))) exit (1);
  fclose(in_file);

  nimages = 0;
  for (VFirstAttr (list, & posn); VAttrExists (& posn); VNextAttr (& posn)) {
    if (VGetAttrRepn (& posn) == VImageRepn) nimages++;
  }

  n = 0;
  for (VFirstAttr (list, & posn); VAttrExists (& posn); VNextAttr (& posn)) {
    if (VGetAttrRepn (& posn) != VImageRepn) continue;
    VGetAttrValue (& posn, NULL,VImageRepn, & src);
    fprintf(stderr," image %4d of %d\r",n,nimages);
    dest = VWavelets(src,dest,(int)ftype,(int)wtype,(int)nshift,level);
    src = VCopyImage(dest,src,VAllBands);
    VSetAttrValue (& posn, NULL,VImageRepn,src);
    n++;
  }
  if (src == NULL) VError(" no input image found");

  VHistory(VNumber(options),options,prg,&list,&list);
  if (! VWriteFile (out_file, list)) exit (1);
  fprintf (stderr, "%s: done.\n", argv[0]);
  return 0;
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
#ifndef OPENMC_PHOTON_H
#define OPENMC_PHOTON_H

#include "openmc/endf.h"
#include "openmc/memory.h" // for unique_ptr
#include "openmc/particle.h"
#include "openmc/vector.h"

#include <gsl/gsl>
#include <hdf5.h>
#include "xtensor/xtensor.hpp"

#include <string>
#include <unordered_map>
#include <utility> // for pair

namespace openmc {

//==============================================================================
//! Photon interaction data for a single element
//==============================================================================

class ElectronSubshell {
public:
  // Constructors
  ElectronSubshell() { };

  int index_subshell;  //!< index in SUBSHELLS
  int threshold;
  double n_electrons;
  double binding_energy;
  xt::xtensor<double, 1> cross_section;

  // Transition data
  int n_transitions;
  xt::xtensor<int, 2> transition_subshells;
  xt::xtensor<double, 1> transition_energy;
  xt::xtensor<double, 1> transition_probability;
};

class PhotonInteraction {
public:
  // Constructors/destructor
  PhotonInteraction(hid_t group);
  ~PhotonInteraction();

  // Methods
  void calculate_xs(Particle& p) const;

  void compton_scatter(double alpha, bool doppler, double* alpha_out,
    double* mu, int* i_shell, uint64_t* seed) const;

  double rayleigh_scatter(double alpha, uint64_t* seed) const;

  void pair_production(double alpha, double* E_electron, double* E_positron,
    double* mu_electron, double* mu_positron, uint64_t* seed) const;

  void atomic_relaxation(const ElectronSubshell& shell, Particle& p) const;

  // Data members
  std::string name_; //!< Name of element, e.g. "Zr"
  int Z_; //!< Atomic number
  gsl::index index_; //!< Index in global elements vector

  // Microscopic cross sections
  xt::xtensor<double, 1> energy_;
  xt::xtensor<double, 1> coherent_;
  xt::xtensor<double, 1> incoherent_;
  xt::xtensor<double, 1> photoelectric_total_;
  xt::xtensor<double, 1> pair_production_total_;
  xt::xtensor<double, 1> pair_production_electron_;
  xt::xtensor<double, 1> pair_production_nuclear_;
  xt::xtensor<double, 1> heating_;

  // Form factors
  Tabulated1D incoherent_form_factor_;
  Tabulated1D coherent_int_form_factor_;
  Tabulated1D coherent_anomalous_real_;
  Tabulated1D coherent_anomalous_imag_;

  // Photoionization and atomic relaxation data
  std::unordered_map<int, int> shell_map_; //!< Given a shell designator, e.g. 3, this
                                           //!< dictionary gives an index in shells_
  vector<ElectronSubshell> shells_;

  // Compton profile data
  xt::xtensor<double, 2> profile_pdf_;
  xt::xtensor<double, 2> profile_cdf_;
  xt::xtensor<double, 1> binding_energy_;
  xt::xtensor<double, 1> electron_pdf_;

  // Stopping power data
  double I_; // mean excitation energy
  xt::xtensor<int, 1> n_electrons_;
  xt::xtensor<double, 1> ionization_energy_;
  xt::xtensor<double, 1> stopping_power_radiative_;

  // Bremsstrahlung scaled DCS
  xt::xtensor<double, 2> dcs_;

private:
  void compton_doppler(double alpha, double mu, double* E_out, int* i_shell,
                       uint64_t* seed) const;
};

//==============================================================================
// Non-member functions
//==============================================================================

std::pair<double, double> klein_nishina(double alpha, uint64_t* seed);

void free_memory_photon();

//==============================================================================
// Global variables
//==============================================================================

namespace data {

extern xt::xtensor<double, 1> compton_profile_pz; //! Compton profile momentum grid

//! Photon interaction data for each element
extern std::unordered_map<std::string, int> element_map;
extern vector<unique_ptr<PhotonInteraction>> elements;

} // namespace data

} // namespace openmc

#endif // OPENMC_PHOTON_H

max_stars_repo_name: cjwyett/openmc

max_issues_repo_name: cjwyett/openmc

max_forks_repo_name: cjwyett/openmc


text:
#include <stdio.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_vegas.h>
#include <gsl/gsl_math.h>

#include "timer.c"
#include "timer.h"

extern double g (double *t, size_t dim, void *params);

double dipole_approx (double r);

double gaussian (double *x, int dim);

int main (void)
{
    double res, err;

    size_t dim = 6;
    double x1[] = { 0., 0., 0., 0., 0., 0., };
    double xu[] = { 1., 1., 1., 1., 1., 1., };
    double distmin = 1.001;
    double distmax = 4.;
    double dist;
    int np = 20;
    double nt = (distmax - distmin) / (np - 1);
    double vegas[20], dipole[20], distance[20];

    gsl_rng *r = gsl_rng_alloc (gsl_rng_taus2);
    unsigned long seed = 1UL;

    gsl_rng_set (r, seed);

    size_t calls = 1000000;

    dist = distmin;

    gsl_monte_function G = { &g, dim, &dist };

    gsl_monte_vegas_state *sv = gsl_monte_vegas_alloc (dim);

    gsl_monte_vegas_init (sv);

 // Vegas Integration
 /*commented so does not output in order to create a proper res
   printf ("#  Stat      Dist       E(r)        ErrEst      Dipolapprox\n");      */

  timer_start ();

  for (int i = 0; i < np; i++) {
     gsl_monte_vegas_integrate (&G, x1, xu, dim, calls / 5, r, sv, &res, &err);
  do
  {
  gsl_monte_vegas_integrate (&G, x1, xu, dim, calls, r, sv, &res, &err);          }
  while (fabs (gsl_monte_vegas_chisq (sv) - 1.0) > 0.2);

    fflush (stdout);
    dist += nt;
    vegas[i] = res;
    distance[i] = dist;
    dipole[i] = -2. / pow (dist, 3.);
    }

    timer_stop();

    gsl_monte_vegas_free (sv);

    double sum;
    double x[6];

    long i, j, nn;

    nn = 1000000;

    timer_start ();
    double home[20];

    dist = distmin;
    for (j = 0; j < np; j++)
    {
        sum = 0.;
        for (i = 0; i < nn; i++)
        {

        for (int k = 0; k < (int) dim; k++)
            {
                x[k] = gsl_rng_uniform (r);
            }
            sum += g (x, dim, &dist);
        }
        res = sum/nn;
        dist += nt;
        home[j] = res;
    }

     timer_stop ();

     gsl_rng_free (r);

     double homeerr = 0.0;

    for (int e = 0; e < np; e++)
    {
    homeerr += fabs(home[e] - vegas[e]);
    }
    printf("#    Dist             Vegas          Home         Dipolapprox\n");
    for( int l = 0; l < np; l++)
    {
    double dd = distance[l];
    double vv = fabs(vegas[l]);
    double hh = fabs(home[l]);
    double di = fabs(dipole[l]);
    printf("   %.6f         %.6f       %.6f      %.6f\n", dd, vv, hh, di);
    }

    return 0;
}


text:
#include <pthread.h>
#include <stdlib.h>
#include <cblas.h>

#include "tasks.h"

void update_task_seq(void *ptr)
{
    struct update_task_arg *arg = (struct update_task_arg*) ptr;

    int nrows = arg->nrows;
    int ncols = arg->ncols;
    int nrhs  = arg->nrhs;
    double *L = arg->L;
    int ldL   = arg->ldL;
    double *X = arg->X;
    int ldX   = arg->ldX;
    double *B = arg->B;
    int ldB   = arg->ldB;

    // Allocate temporary storage W.
    int ldW = nrows;
    double *W = (double*) malloc(sizeof(*W) * ldW * nrhs);

    // Compute W := L * X.
    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
                nrows, nrhs, ncols,
                1.0, L, ldL,
                     X, ldX,
                0.0, W, ldW);

    // Acquire lock.
    pthread_mutex_lock(arg->Block);

    // Update B := B - W.
#define B(i,j) B[(i) + (j) * ldB]
#define W(i,j) W[(i) + (j) * ldW]
    for (int j = 0; j < nrhs; ++j) {
        for (int i = 0; i < nrows; ++i) {
            B(i,j) -= W(i,j);
        }
    }
#undef B
#undef W

    // Release lock.
    pthread_mutex_unlock(arg->Block);

    // Free temporary storage.
    free(W);
}

max_stars_repo_name: NLAFET/pcp-runtime

max_issues_repo_name: NLAFET/pcp-runtime

max_forks_repo_name: NLAFET/pcp-runtime


text:
/* integration/qagp.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_integration.h>

static int
qagp (const gsl_function *f,
      const double *pts, const size_t npts,
      const double epsabs, const double epsrel, const size_t limit,
      gsl_integration_workspace * workspace,
      double *result, double *abserr,
      gsl_integration_rule * q);

#include "initialise.c"
#include "qpsrt.c"
#include "util.c"
#include "append.c"
#include "reset.c"
#include "qelg.c"
#include "qpsrt2.c"
#include "ptsort.c"
#include "positivity.c"

int
gsl_integration_qagp (const gsl_function *f,
                      double * pts, size_t npts,
                      double epsabs, double epsrel, size_t limit,
                      gsl_integration_workspace * workspace,
                      double * result, double * abserr)
{
  int status = qagp (f, pts, npts,
                     epsabs, epsrel, limit,
                     workspace,
                     result, abserr,
                     &gsl_integration_qk21) ;

  return status ;
}

static int
qagp (const gsl_function * f,
      const double *pts, const size_t npts,
      const double epsabs, const double epsrel,
      const size_t limit,
      gsl_integration_workspace * workspace,
      double *result, double *abserr,
      gsl_integration_rule * q)
{
  double area, errsum;
  double res_ext, err_ext;
  double result0, abserr0, resabs0;
  double tolerance;

  double ertest = 0;
  double error_over_large_intervals = 0;
  double reseps = 0, abseps = 0, correc = 0;
  size_t ktmin = 0;
  int roundoff_type1 = 0, roundoff_type2 = 0, roundoff_type3 = 0;
  int error_type = 0, error_type2 = 0;

  size_t iteration = 0;

  int positive_integrand = 0;
  int extrapolate = 0;
  int disallow_extrapolation = 0;

  struct extrapolation_table table;

  const size_t nint = npts - 1; /* number of intervals */

  size_t *ndin = workspace->level; /* temporarily alias ndin to level */

  size_t i;

  /* Initialize results */

  *result = 0;
  *abserr = 0;

  /* Test on validity of parameters */

  if (limit > workspace->limit)
    {
      GSL_ERROR ("iteration limit exceeds available workspace", GSL_EINVAL) ;
    }

  if (npts > workspace->limit)
    {
      GSL_ERROR ("npts exceeds size of workspace", GSL_EINVAL);
    }

  if (epsabs <= 0 && (epsrel < 50 * GSL_DBL_EPSILON || epsrel < 0.5e-28))
    {
      GSL_ERROR ("tolerance cannot be acheived with given epsabs and epsrel",
                 GSL_EBADTOL);
    }

  /* Check that the integration range and break points are an
     ascending sequence */

  for (i = 0; i < nint; i++)
    {
      if (pts[i + 1] < pts[i])
        {
          GSL_ERROR ("points are not in an ascending sequence", GSL_EINVAL);
        }
    }

  /* Perform the first integration */

  result0 = 0;
  abserr0 = 0;
  resabs0 = 0;

  initialise (workspace, 0.0, 0.0) ;

  for (i = 0; i < nint; i++)
    {
      double area1, error1, resabs1, resasc1;
      const double a1 = pts[i];
      const double b1 = pts[i + 1];

      q (f, a1, b1, &area1, &error1, &resabs1, &resasc1);

      result0 = result0 + area1;
      abserr0 = abserr0 + error1;
      resabs0 = resabs0 + resabs1;

      append_interval (workspace, a1, b1, area1, error1);

      if (error1 == resasc1 && error1 != 0.0)
        {
          ndin[i] = 1;
        }
      else
        {
          ndin[i] = 0;
        }
    }

  /* Compute the initial error estimate */

  errsum = 0.0;

  for (i = 0; i < nint; i++)
    {
      if (ndin[i])
        {
          workspace->elist[i] = abserr0;
        }

      errsum = errsum + workspace->elist[i];

    }

  for (i = 0; i < nint; i++)
    {
      workspace->level[i] = 0;
    }

  /* Sort results into order of decreasing error via the indirection
     array order[] */

  sort_results (workspace);

  /* Test on accuracy */

  tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (result0));

  if (abserr0 <= 100 * GSL_DBL_EPSILON * resabs0 && abserr0 > tolerance)
    {
      *result = result0;
      *abserr = abserr0;

      GSL_ERROR ("cannot reach tolerance because of roundoff error"
                 "on first attempt", GSL_EROUND);
    }
  else if (abserr0 <= tolerance)
    {
      *result = result0;
      *abserr = abserr0;

      return GSL_SUCCESS;
    }
  else if (limit == 1)
    {
      *result = result0;
      *abserr = abserr0;

      GSL_ERROR ("a maximum of one iteration was insufficient", GSL_EMAXITER);
    }

  /* Initialization */

  initialise_table (&table);
  append_table (&table, result0);

  area = result0;

  res_ext = result0;
  err_ext = GSL_DBL_MAX;

  error_over_large_intervals = errsum;
  ertest = tolerance;

  positive_integrand = test_positivity (result0, resabs0);

  iteration = nint - 1;

  do
    {
      size_t current_level;
      double a1, b1, a2, b2;
      double a_i, b_i, r_i, e_i;
      double area1 = 0, area2 = 0, area12 = 0;
      double error1 = 0, error2 = 0, error12 = 0;
      double resasc1, resasc2;
      double resabs1, resabs2;
      double last_e_i;

      /* Bisect the subinterval with the largest error estimate */

      retrieve (workspace, &a_i, &b_i, &r_i, &e_i);

      current_level = workspace->level[workspace->i] + 1;

      a1 = a_i;
      b1 = 0.5 * (a_i + b_i);
      a2 = b1;
      b2 = b_i;

      iteration++;

      q (f, a1, b1, &area1, &error1, &resabs1, &resasc1);
      q (f, a2, b2, &area2, &error2, &resabs2, &resasc2);

      area12 = area1 + area2;
      error12 = error1 + error2;
      last_e_i = e_i;

      /* Improve previous approximations to the integral and test for
         accuracy.

         We write these expressions in the same way as the original
         QUADPACK code so that the rounding errors are the same, which
         makes testing easier. */

      errsum = errsum + error12 - e_i;
      area = area + area12 - r_i;

      tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (area));

      if (resasc1 != error1 && resasc2 != error2)
        {
          double delta = r_i - area12;

          if (fabs (delta) <= 1.0e-5 * fabs (area12) && error12 >= 0.99 * e_i)
            {
              if (!extrapolate)
                {
                  roundoff_type1++;
                }
              else
                {
                  roundoff_type2++;
                }
            }

          if (i > 10 && error12 > e_i)
            {
              roundoff_type3++;
            }
        }

      /* Test for roundoff and eventually set error flag */

      if (roundoff_type1 + roundoff_type2 >= 10 || roundoff_type3 >= 20)
        {
          error_type = 2;       /* round off error */
        }

      if (roundoff_type2 >= 5)
        {
          error_type2 = 1;
        }

      /* set error flag in the case of bad integrand behaviour at
         a point of the integration range */

      if (subinterval_too_small (a1, a2, b2))
        {
          error_type = 4;
        }

      /* append the newly-created intervals to the list */

      update (workspace, a1, b1, area1, error1, a2, b2, area2, error2);

      if (errsum <= tolerance)
        {
          goto compute_result;
        }

      if (error_type)
        {
          break;
        }

      if (iteration >= limit - 1)
        {
          error_type = 1;
          break;
        }

      if (disallow_extrapolation)
        {
          continue;
        }

      error_over_large_intervals += -last_e_i;

      if (current_level < workspace->maximum_level)
        {
          error_over_large_intervals += error12;
        }

      if (!extrapolate)
        {
          /* test whether the interval to be bisected next is the
             smallest interval. */
          if (large_interval (workspace))
            continue;

          extrapolate = 1;
          workspace->nrmax = 1;
        }

      /* The smallest interval has the largest error.  Before
         bisecting decrease the sum of the errors over the larger
         intervals (error_over_large_intervals) and perform
         extrapolation. */

      if (!error_type2 && error_over_large_intervals > ertest)
        {
          if (increase_nrmax (workspace))
            continue;
        }

      /* Perform extrapolation */

      append_table (&table, area);

      if (table.n < 3)
        {
          goto skip_extrapolation;
        }

      qelg (&table, &reseps, &abseps);

      ktmin++;

      if (ktmin > 5 && err_ext < 0.001 * errsum)
        {
          error_type = 5;
        }

      if (abseps < err_ext)
        {
          ktmin = 0;
          err_ext = abseps;
          res_ext = reseps;
          correc = error_over_large_intervals;
          ertest = GSL_MAX_DBL (epsabs, epsrel * fabs (reseps));
          if (err_ext <= ertest)
            break;
        }

      /* Prepare bisection of the smallest interval. */

      if (table.n == 1)
        {
          disallow_extrapolation = 1;
        }

      if (error_type == 5)
        {
          break;
        }

    skip_extrapolation:

      reset_nrmax (workspace);
      extrapolate = 0;
      error_over_large_intervals = errsum;

    }
  while (iteration < limit);

  *result = res_ext;
  *abserr = err_ext;

  if (err_ext == GSL_DBL_MAX)
    goto compute_result;

  if (error_type || error_type2)
    {
      if (error_type2)
        {
          err_ext += correc;
        }

      if (error_type == 0)
        error_type = 3;

      if (result != 0 && area != 0)
        {
          if (err_ext / fabs (res_ext) > errsum / fabs (area))
            goto compute_result;
        }
      else if (err_ext > errsum)
        {
          goto compute_result;
        }
      else if (area == 0.0)
        {
          goto return_error;
        }
    }

  /*  Test on divergence. */

  {
    double max_area = GSL_MAX_DBL (fabs (res_ext), fabs (area));

    if (!positive_integrand && max_area < 0.01 * resabs0)
      goto return_error;
  }

  {
    double ratio = res_ext / area;

    if (ratio < 0.01 || ratio > 100 || errsum > fabs (area))
      error_type = 6;
  }

  goto return_error;

compute_result:

  *result = sum_results (workspace);
  *abserr = errsum;

return_error:

  if (error_type > 2)
    error_type--;

  if (error_type == 0)
    {
      return GSL_SUCCESS;
    }
  else if (error_type == 1)
    {
      GSL_ERROR ("number of iterations was insufficient", GSL_EMAXITER);
    }
  else if (error_type == 2)
    {
      GSL_ERROR ("cannot reach tolerance because of roundoff error",
                 GSL_EROUND);
    }
  else if (error_type == 3)
    {
      GSL_ERROR ("bad integrand behavior found in the integration interval",
                 GSL_ESING);
    }
  else if (error_type == 4)
    {
      GSL_ERROR ("roundoff error detected in the extrapolation table",
                 GSL_EROUND);
    }
  else if (error_type == 5)
    {
      GSL_ERROR ("integral is divergent, or slowly convergent",
                 GSL_EDIVERGE);
    }
  else
    {
      GSL_ERROR ("could not integrate function", GSL_EFAILED);
    }
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#pragma once

#include <memory>
#include "halley/text/halleystring.h"
#include "halley/core/graphics/texture.h"
#include <gsl/gsl>

namespace Halley
{
	enum class ShaderType;
	class MaterialDataBlock;
	class Material;
	enum class ShaderParameterType;
	class Painter;
	class MaterialDefinition;
	class MaterialParameter;
	class MaterialTextureParameter;
	class VideoAPI;

	class MaterialConstantBuffer
	{
	public:
		virtual ~MaterialConstantBuffer() {}

		virtual void update(const MaterialDataBlock& dataBlock) = 0;
	};

	enum class MaterialDataBlockType
	{
		// Shared blocks are not stored locally in the material (e.g. the HalleyBlock, stored by the engine)
		SharedLocal,    // Shared, this keeps the canonical copy
		SharedExternal, // Shared, only a reference
		Local           // Local data
	};

	class MaterialDataBlock
	{
		friend class Material;

	public:
		MaterialDataBlock();
		MaterialDataBlock(MaterialDataBlockType type, size_t size, int bindPoint, const String& name, const MaterialDefinition& def);
		MaterialDataBlock(const MaterialDataBlock& other);
		MaterialDataBlock(MaterialDataBlock&& other) noexcept;

		MaterialConstantBuffer& getConstantBuffer() const;
		int getAddress(int pass, ShaderType stage) const;
		int getBindPoint() const;
		gsl::span<const gsl::byte> getData() const;
		MaterialDataBlockType getType() const;

	private:
		std::unique_ptr<MaterialConstantBuffer> constantBuffer;
		Bytes data;
		Vector<int> addresses;
		MaterialDataBlockType dataBlockType;
		int bindPoint = 0;
		bool dirty = true;

		bool setUniform(size_t offset, ShaderParameterType type, const void* data);
		void upload(VideoAPI* api);
	};

	class Material
	{
		friend class MaterialParameter;

	public:
		Material(const Material& other);
		explicit Material(std::shared_ptr<const MaterialDefinition> materialDefinition, bool forceLocalBlocks = false); // forceLocalBlocks is for engine use only

		void bind(int pass, Painter& painter);
		void uploadData(Painter& painter);
		static void resetBindCache();

		bool operator==(const Material& material) const;
		bool operator!=(const Material& material) const;

		const MaterialDefinition& getDefinition() const { return *materialDefinition; }

		std::shared_ptr<Material> clone() const;

		const std::shared_ptr<const Texture>& getTexture(int textureUnit) const;
		const Vector<MaterialTextureParameter>& getTextureUniforms() const;
		const std::vector<std::shared_ptr<const Texture>>& getTextures() const;

		const Vector<MaterialParameter>& getUniforms() const;
		const Vector<MaterialDataBlock>& getDataBlocks() const;

		void setPassEnabled(int pass, bool enabled);
		bool isPassEnabled(int pass) const;

		Material& set(const String& name, const std::shared_ptr<const Texture>& texture);
		Material& set(const String& name, const std::shared_ptr<Texture>& texture);

		bool hasParameter(const String& name) const;

		template <typename T>
		Material& set(const String& name, const T& value)
		{
			getParameter(name) = value;
			return *this;
		}

		uint64_t getHash() const;

	private:
		std::shared_ptr<const MaterialDefinition> materialDefinition;

		Vector<MaterialParameter> uniforms;
		Vector<MaterialTextureParameter> textureUniforms;
		Vector<MaterialDataBlock> dataBlocks;
		std::vector<std::shared_ptr<const Texture>> textures;

		std::vector<char> passEnabled;

		mutable uint64_t hashValue = 0;
		mutable bool needToUpdateHash = true;
		bool needToUploadData = true;

		void initUniforms(bool forceLocalBlocks);
		MaterialParameter& getParameter(const String& name);

		void setUniform(int blockNumber, size_t offset, ShaderParameterType type, const void* data);
		uint64_t computeHash() const;

		const std::shared_ptr<const Texture>& getFallbackTexture() const;
	};
}

max_stars_repo_name: Bearwaves/halley

max_issues_repo_name: Bearwaves/halley

max_forks_repo_name: Bearwaves/halley


text:
/**
 * @file bblas_sutil.c
 *
 * @brief BBLAS testing utilities for float routines.
 *
 * BBLAS is a software package provided by Univ. of Manchester,
 * Univ. of Tennessee.
 *
 * @version 1.0.0
 * @author  Samuel  D. Relton
 * @author  Pedro   V. Lara
 * @author  Mawussi Zounon
 * @date    2016-02-20
 *
 * Contains routines used in the testing to modify randomly generated matrices
 * and compute the average error over an entire batch etc.
 *
 **/

#ifndef DOXYGEN_SHOULD_SKIP_THIS
/**
 * Code generation
 * @generated from bblas_zutil.c normal z -> s, Mon Jun  6 09:44:13 2016
 **/
#endif

#include "bblas_common.h"
#if defined(BBLAS_WITH_MKL)
#include <mkl_lapacke.h>
#else
#include <lapacke.h>
#endif

/** Include real functions since using float real precision **/
#define REAL

/** Quick access to the matrix elements **/
#define A(i,j)  A[i + j*lda]

/**
 * Make a matrix symmetric/symmetric.  Makes diagonal real.
 * Sets Aji = conj( Aij ) for j < i, that is, copy & conjugate
 * lower triangle to upper triangle.
 **/

void bblas_smake_symmetric(int lda, int N, float* A)
{
    int i, j;
    for( i=0; i < N; ++i ) {
        A(i,i) = creal( A(i,i) );
        for( j=0; j < i; ++j ) {
            A(j,i) = conj( A(i,j) );
        }
    }
}

#ifdef COMPLEX

/**
 * Make a matrix real-symmetric
 * Does NOT make diagonal real.
 * Sets Aji = Aij for j < i, that is,
 * copy lower triangle to upper triangle.
 **/

void bblas_smake_symmetric(int lda, int N, float* A)
{
    int i, j;
    for( i=0; i < N; ++i ) {
        for( j=0; j < i; ++j ) {
            A(j,i) =  A(i,j);
        }
    }
}
#endif

/**
 * irandRange generates a random value (int)
 * in the range min_n and max_n
 **/

int irandRange(int min_n, int max_n)
{
    return rand() % (max_n - min_n + 1) + min_n;
}

/**
 * bblas_srandRange generates a random value (float)
 * in the range [0,max_n]: TODO replace by generic
 */

float bblas_srandRange( int max_n )
{
    return ( float )rand()/( float )( RAND_MAX/max_n );
}

/**
 * Computes statistics of the relative errors to summarise them for the user.
 **/

void bblas_sstatistic(bblas_stest_t *test)
{

    enum BBLAS_ROUTINE routine  = test->routine;

    /*Compute avg(M), avg(N), avg(K) */
    if(test->batch_opts == BBLAS_VARIABLE)
    {
	if ((routine == BBLAS_GEMM) || (routine == BBLAS_SYMM) ||
	    (routine == BBLAS_HEMM) || (routine == BBLAS_TRMM) ||
	    (routine == BBLAS_TRSM))
	{
	    test->avgM = bblas_avgarrayI(test->M, test->batch_count);
	}

	if ((routine == BBLAS_GEMM) || (routine == BBLAS_SYRK)  ||
	    (routine == BBLAS_HERK) || (routine == BBLAS_SYR2K) ||
	    (routine == BBLAS_HER2K))
	{
	    test->avgK = bblas_avgarrayI(test->K, test->batch_count);
	}

	test->avgN = bblas_avgarrayI(test->N, test->batch_count);

    } else if (test->batch_opts == BBLAS_FIXED)
    {

	if ((routine == BBLAS_GEMM) || (routine == BBLAS_SYMM) ||
	    (routine == BBLAS_HEMM) || (routine == BBLAS_TRMM) ||
	    (routine == BBLAS_TRSM))
	{
	    test->avgM = test->M[0];
	}

	if ((routine == BBLAS_GEMM) || (routine == BBLAS_SYRK)  ||
	    (routine == BBLAS_HERK) || (routine == BBLAS_SYR2K) ||
	    (routine == BBLAS_HER2K))
	{
	    test->avgK = test->K[0];
	}
	test->avgN = test->N[0];

    } else
    {
	bblas_error("testing_sgemm_batch.c", "wrong batch_opts value");
    }

    /*Statistics on the error */
    switch(test->target)
    {
    case BBLAS_MKL:
      test->mkl_min_error = bblas_sminarrayD(test->mkl_error, test->batch_count);
      test->mkl_avg_error = bblas_savgarrayD(test->mkl_error, test->batch_count);
      test->mkl_max_error = bblas_smaxarrayD(test->mkl_error, test->batch_count);
      test->mkl_std_error = bblas_sstdarrayD(test->mkl_error, test->batch_count);
      break;

    case BBLAS_CUBLAS:
    case BBLAS_MAGMA:
      test->device_min_error = bblas_sminarrayD(test->device_error, test->batch_count);
      test->device_avg_error = bblas_savgarrayD(test->device_error, test->batch_count);
      test->device_max_error = bblas_smaxarrayD(test->device_error, test->batch_count);
      test->device_std_error = bblas_sstdarrayD(test->device_error, test->batch_count);
      break;

    case BBLAS_OTHER:
      test->other_min_error = bblas_sminarrayD(test->other_error, test->batch_count);
      test->other_avg_error = bblas_savgarrayD(test->other_error, test->batch_count);
      test->other_max_error = bblas_smaxarrayD(test->other_error, test->batch_count);
      test->other_std_error = bblas_sstdarrayD(test->other_error, test->batch_count);
      break;

    default:
      printf("In bblas_sstatistic(): Target no defined\n");
      exit(EXIT_FAILURE);
    }
}

/**
 * Print a matrix.
 **/
void bblas_sprintmatrix(float *matrix, int row, int col)
{
    /*Local variables */
    int i, j;

    for (i=0; i < row; i++)
    {
	printf("\n\n");
	for (j=0; j < col; j++)
	{
#ifdef COMPLEX
	    printf("%1.2f + %1.2f\t", creal(matrix[i*col+j]), cimag(matrix[i*col+j]));
#else
	    printf("%1.2f",matrix[i*col+j]);
#endif
	}
    }
    printf("\n");
}

/**
 * Decide whether a batch is fixed or variable.
 **/

char* bblas_getoption(enum BBLAS_OPTS opts)
{
   /*Local variable */
    char funcname[] = "bblas_getoption";

    switch(opts)
    {
	case BBLAS_VARIABLE:
	    return "BATCH OPTION: VARIABLE";
	    break;

	case BBLAS_FIXED:
	    return "BATCH OPTION: FIXED";
	    break;

	default:
	    printf("ERROR in %s, undefined bblas routine name\n",funcname);
	    exit(EXIT_FAILURE);
    }

}

/**
 * Get the name of the current test routine.
 **/

char* bblas_getroutine(enum BBLAS_ROUTINE routine)
{
    /*Local variable */
    char funcname[] = "bblas_getroutine";

    switch(routine)
    {
	case BBLAS_GEMM:
	    return "SGEMM";
	    break;

	case BBLAS_HEMM:
	    return "SSYMM";
	    break;

	case BBLAS_HER2K:
	    return "SSYR2K";
	    break;

	case BBLAS_HERK:
	    return "SSYRK";
	    break;

	case BBLAS_SYMM:
	    return "SSYMM";
	    break;

	case BBLAS_SYR2K:
	    return "SSYR2K";
	    break;

	case BBLAS_SYRK:
	    return "SSYRK";
	    break;

	case BBLAS_TRMM:
	    return "STRMM";
	    break;

	case BBLAS_TRSM:
	    return "STRSM";
	    break;

	default:
	    printf("ERROR in %s, undefined bblas routine name\n",funcname);
	    exit(EXIT_FAILURE);
    }
}

/**
 * Computes the maximum value of an array of real floats.
 **/

float bblas_smaxarrayD(float *myArray, int size)
{
    int iter;
    float maxValue = myArray[0];

    for (iter = 0; iter < size; ++iter)
    {
        if ( myArray[iter] > maxValue )
	{
            maxValue = myArray[iter];
        }
    }
    return maxValue;
}

/**
 * Computes the minimum value of an array of real floats.
 **/

float bblas_sminarrayD(float *myArray, int size)
{
    int iter;
    float minValue = myArray[0];

    for (iter = 0; iter < size; ++iter)
    {
        if ( myArray[iter] < minValue )
	{
            minValue = myArray[iter];
        }
    }
    return minValue;
}

/**
 * Computes the mean value of an array of real floats.
 **/

float bblas_savgarrayD(float *myArray, int size)
{
    int iter;
    float avg = 0.;

    for (iter = 0; iter < size; ++iter)
    {
	avg += myArray[iter];
    }
    return avg/size;
}

/**
 * Computes the standard deviation of an array of real floats.
 **/

float bblas_sstdarrayD(float *myArray, int size)
{
    int iter;
    float avg, sd=0.;

    avg =  bblas_savgarrayD(myArray, size);

    for (iter = 0; iter < size; ++iter)
    {
	sd += (myArray[iter] -avg)*(myArray[iter] -avg);
    }
    return sd/size;
}

/**
 * Computes the minimum value of an array of integers.
 **/

int bblas_minarrayI(int *myArray, int size)
{
    int iter;
    int minValue = myArray[0];

    for (iter = 0; iter < size; ++iter)
    {
        if ( myArray[iter] < minValue )
	{
            minValue = myArray[iter];
        }
    }
    return minValue;
}

/**
 * Computes the mean of an array of integers.
 **/

int bblas_avgarrayI(int *myArray, int size)
{
    int iter;
    int avg = 0;

    for (iter = 0; iter < size; ++iter)
    {
	avg += myArray[iter];
    }
    return avg/size;
}

/**
 * Transform BBLAS enum values for <tt>trans</tt>, <tt>uplo</tt> etc. to human-readable strings.
 **/

char* bblas_op2char(unsigned int op)
{
    char *opname = (char*)malloc(30*sizeof(char));
    switch(op)
    {
        case BblasNoTrans:
            strcpy(opname,"CblasNoTrans");
            break;

        case BblasTrans:
	    strcpy(opname,"CblasTrans");
            break;

        case BblasConjTrans:
            strcpy(opname,"CblasConjTrans");
            break;

        case BblasLower:
            strcpy(opname,"CblasLower");
            break;

        case BblasUpper:
	    strcpy(opname,"CblasUpper");
            break;

        case BblasNonUnit:
	    strcpy(opname,"CblasNonUnit");
            break;

        case BblasUnit:
	   strcpy(opname,"CblasUnit");
            break;

        case BblasLeft:
	   strcpy(opname,"CblasLeft");
            break;

        case BblasRight:
	    strcpy(opname,"CblasRight");
            break;

        default:
            return 0;
            exit(EXIT_FAILURE);
    }
    return opname;
}

/**
 * Get the amount of data needed.
 **/

int bblas_snbdata(bblas_stest_t *test)
{
    enum BBLAS_OPTS batch_opts        = test->batch_opts;
    int nb_data;
    char function_name[NAME_LENGTH]   ="bblas_getdatacount";

    switch(batch_opts)
    {
	case BBLAS_VARIABLE:
	    nb_data =  test->batch_count;
	    break;

	case BBLAS_FIXED:
	    nb_data = 1;
	    break;

	default:
	    bblas_fatal_error(function_name, "wrong batch_opts value");
    }

    return nb_data;
}

/**
 * Inject an error into the computation if this is set in the input file.
 **/

void bblas_sset_error(bblas_stest_t *test)
{
    int   nb_data =  bblas_snbdata(test);
    int   routine = test->routine;
    int   error_index = irandRange(0, nb_data);

    if (test->global_error){
	test->batch_count = -1;
	return;
    }

    if (test->batch_opts == BBLAS_FIXED)
    {
	error_index = 0;
    }

    if (test->set_error)
    {
	if ((routine == BBLAS_GEMM) || (routine == BBLAS_SYMM) ||
	    (routine == BBLAS_HEMM) || (routine == BBLAS_TRMM) ||
	    (routine == BBLAS_TRSM))
	{
	    test->M[error_index] = -1;
	}else if ((routine == BBLAS_SYRK) || (routine == BBLAS_HERK) ||
		  (routine == BBLAS_SYR2K)|| (routine == BBLAS_HER2K))
	{
	    test->K[error_index] = -1;
	}
    }
}

/**
 * Check whether a computation has passed or failed our accuracy test.
 * When new_accuracy=1 in the input file this uses an appropriate
 * forward/backward error bound,
 * otherwise this looks at the relative error.
 **/

void bblas_spassed_failed(bblas_stest_t  *test, float error, char *result, int info)
{
    float eps = LAPACKE_slamch_work('e');

	/* Use our new accuracy test based on forward/backward error analysis*/
    if (test->new_accuracy)
    {
	if( (error > 1) || (info))
	{
	    strcpy(result, "FAILED");

	}else
	{
	    strcpy(result, "PASSED");
	}
    }else
    {
		/* Use old accuracy test based on the relative error */
	if((error > eps*test->tolerance) || (info))
	{
	    strcpy(result, "FAILED");

	}else
	{
	    strcpy(result, "PASSED");
	}
    }
}

/**
 * Set the batch_count.
 **/

void bblas_sset_batch_count(bblas_stest_t *test)
{

  test->batch_count = test->minbatch_count* (test->current_iter+1);
}

#undef REAL

max_stars_repo_name: NLAFET/BBLAS-ref

max_issues_repo_name: NLAFET/BBLAS-ref

max_forks_repo_name: NLAFET/BBLAS-ref


text:
/////////////////////////////////////////////////////////////////////
// = NMatrix
//
// A linear algebra library for scientific computation in Ruby.
// NMatrix is part of SciRuby.
//
// NMatrix was originally inspired by and derived from NArray, by
// Masahiro Tanaka: http://narray.rubyforge.org
//
// == Copyright Information
//
// SciRuby is Copyright (c) 2010 - 2014, Ruby Science Foundation
// NMatrix is Copyright (c) 2012 - 2014, John Woods and the Ruby Science Foundation
//
// Please see LICENSE.txt for additional copyright notices.
//
// == Contributing
//
// By contributing source code to SciRuby, you agree to be bound by
// our Contributor Agreement:
//
// * https://github.com/SciRuby/sciruby/wiki/Contributor-Agreement
//
// == inc.h
//
// Includes needed for LAPACK, CLAPACK, and CBLAS functions.
//

#ifndef INC_H
# define INC_H

extern "C" { // These need to be in an extern "C" block or you'll get all kinds of undefined symbol errors.
#if defined HAVE_CBLAS_H
  #include <cblas.h>
#elif defined HAVE_ATLAS_CBLAS_H
  #include <atlas/cblas.h>
#endif

#if defined HAVE_CLAPACK_H
  #include <clapack.h>
#elif defined HAVE_ATLAS_CLAPACK_H
  #include <atlas/clapack.h>
#endif
}

#endif // INC_H

max_stars_repo_name: VasiliCekaskin/nmatrix

max_issues_repo_name: VasiliCekaskin/nmatrix

max_forks_repo_name: VasiliCekaskin/nmatrix