add_subdirectory(CXX11)

set(Eigen_HEADERS Tensor TensorSymmetry Core MeasureCacheSizes

   )

install(FILES

  ${Eigen_HEADERS}

  DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/CXX11 COMPONENT Devel

  )

add_subdirectory(src)

// This file is part of Eigen, a lightweight C++ template library

// for linear algebra.

//

// Copyright (C) 2015 Benoit Jacob <benoitjacob@google.com>

//

// This Source Code Form is subject to the terms of the Mozilla

// Public License v. 2.0. If a copy of the MPL was not distributed

// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef EIGEN_MEASURE_CACHE_SIZES_MODULE

#define EIGEN_MEASURE_CACHE_SIZES_MODULE

/** \defgroup MeasureCacheSizes module

  *

  * Offers a simple measure_cache_sizes benchmark measuring CPU cache sizes

  * 

  * This is useful as a fallback in situations where this information

  * can't be queried directly. This information can then be passed to Eigen

  * to tune its algorithms.

  *

  * \code

  * #include <unsupported/Eigen/CXX11/MeasureCacheSizes>

  * \endcode

  */

#include "src/MeasureCacheSizes/MeasureCacheSizes.h"

#endif // EIGEN_MEASURE_CACHE_LEVELS_MODULE

#ADD_SUBDIRECTORY(Tensor)

#ADD_SUBDIRECTORY(TensorSymmetry)

#ADD_SUBDIRECTORY(Core)

ADD_SUBDIRECTORY(MeasureCacheSizes)

FILE(GLOB Eigen_MeasureCacheSizes_SRCS "*.h")

INSTALL(FILES

  ${Eigen_MeasureCacheSizes_SRCS}

  DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/CXX11/src/MeasureCacheSizes COMPONENT Devel

  )

// This file is part of Eigen, a lightweight C++ template library

// for linear algebra.

//

// Copyright (C) 2015 Benoit Jacob <benoitjacob@google.com>

//

// This Source Code Form is subject to the terms of the Mozilla

// Public License v. 2.0. If a copy of the MPL was not distributed

// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef EIGEN_MEASURE_CACHE_SIZES_H

#define EIGEN_MEASURE_CACHE_SIZES_H

#include <Eigen/Core>

#include <ctime>

#include <cassert>

#include <cstdint>

#include <cmath>

#include <cstdlib>

#include <cstring>

#include <limits>

#include <algorithm>

#include <vector>

#ifndef __linux__

#include <sys/time.h>

#endif

namespace Eigen {

namespace internal {

using std::size_t;

using std::ptrdiff_t;

using std::uint8_t;

using std::uint64_t;

const size_t megabyte = 1 << 20;

const size_t kilobyte = 1 << 10;

// We must repeat memory accesses enough times to give

// the CPU a good chance to have our stuff in cache for

// most of the measurements, so that the effect of caches

// can be measured. Since our accesses are by large memcpy's,

// the number of times that each cache line is touched is

// roughly equal to the number of times that each byte is touched.

const int minimum_times_accessed_each_byte = 8;

// Even with clock_gettime, time measurement accuracy is often

// worse than 1e-5 s. We avoid measuring times shorter than 1e-2 s

// so as to make this a non-issue. Increasing this minimum further

// significantly slows down the measurements.

const uint64_t minimum_measurement_time_in_ns = 1e7;

// We repeat all measurements a few times and keep the best throughput

// for each workload. This removes random noise and warm-up effects.

const int measurement_passes = 3;

// Minimum ratio of throughputs of different cache measurements.

// For example, 1.15 means that any difference smaller than 15%

// will be interpreted as measurement inaccuracy, not as

// indicating different cache measurements.

const float never_separate_levels_closer_than = 1.15;

// We need to be able to distinguish the highest level of cache (e.g. L3)

// from RAM. To do so, we need to know the throughput of RAM, so we need

// to know a working set size that cannot plausibly fit in cache. That is

// implausibly_large_cache_size. It should be set to a value larger than

// any CPU's cache on a given architecture.

//

// Conversely implausibly_small_cache_size should be set to a size

// smaller than any CPU's L1 cache, in order for us to correctly detect

// the L1 cache.

//

// Keep in mind that we require each cache level to be

// reflected in at least two consecutive power-of-two measurements,

// so for instance if the L1 cache size is 4k, implausibly_small_cache_size

// must be set to at most 2k in order for us to detect it.

#if defined __aarch64__ || defined __arm__

  const size_t implausibly_large_cache_size = 8 * megabyte;

  const size_t implausibly_small_cache_size = 1 * kilobyte;

#elif defined __x86_64__ || defined __i386__

  const size_t implausibly_large_cache_size = 64 * megabyte;

  const size_t implausibly_small_cache_size = 4 * kilobyte;

#else

#error "What's your architecture?"

#endif

// To avoid alignment-related noise, especially with small working set sizes,

// we keep our buffers on an alignment much higher than any plausible cache

// line size. It doesn't hurt, then, to make this the typical value of

// page size on most systems.

const size_t alignment = 0x1000;

class cache_measurement_context_t

{

  size_t allocated_buffer_size;

  uint8_t* buffer1;

  uint8_t* buffer2;

  __attribute__((noinline))

  static void consume(uint8_t&)

  {}

  static uint64_t get_raw_time_in_ns()

  {

  #ifdef __linux__

    timespec t;

    clock_gettime(CLOCK_THREAD_CPUTIME_ID, &t);

    return t.tv_nsec + 1e9 * t.tv_sec;

  #else

    timeval t;

    gettimeofday(&t, nullptr);

    return 1e3 * t.tv_usec + 1e9 * t.tv_sec;

  #endif

  }

public:

  cache_measurement_context_t(size_t working_set_size)

  {

    allocated_buffer_size = working_set_size / 2;

    bool alloc_failure =

      posix_memalign(reinterpret_cast<void**>(&buffer1), alignment, allocated_buffer_size) ||

      posix_memalign(reinterpret_cast<void**>(&buffer2), alignment, allocated_buffer_size);

    eigen_assert(!alloc_failure);

    EIGEN_ONLY_USED_FOR_DEBUG(alloc_failure);

    // Initial warm up and ensure that the compiler can't optimize all

    // the workload away as depending on uninitialized values

    for (size_t i = 0; i < allocated_buffer_size; i++) {

      buffer1[i] = buffer2[i] = uint8_t(i);

    }

  }

  ~cache_measurement_context_t()

  {

    free(buffer1);

    free(buffer2);

  }

  float measure_throughput(size_t working_set_size)

  {

    const size_t buffer_size = working_set_size / 2;

    eigen_assert(buffer_size <= allocated_buffer_size);

    size_t iters_at_a_time = 1;

    uint64_t bytes_copied = 0;

    const uint64_t minimum_bytes_copied = buffer_size * minimum_times_accessed_each_byte;

    uint64_t t1, t2;

    t1 = get_raw_time_in_ns();

    while (true) {

      t2 = get_raw_time_in_ns();

      uint64_t elapsed = t2 - t1;

      if (elapsed >= minimum_measurement_time_in_ns && bytes_copied >= minimum_bytes_copied) {

        consume(buffer1[buffer_size / 2]);

        consume(buffer2[buffer_size / 2]);

        return float(bytes_copied) / elapsed; // bytes per ns, equivalently GB/s

      }

      for (size_t i = 0; i < iters_at_a_time; i++) {

        // Inner benchmark loop. We use memcpy as it is probably well optimized

        // for optimal throughput.

        // We copy back and forth. This could be updated in the future

        // if we wanted to test more specific cache behavior.

        uint8_t* bufferA = (i & 1) ? buffer1 : buffer2;

        uint8_t* bufferB = (i & 1) ? buffer2 : buffer1;

        // We apply a variable offset to make it near-impossible for the compiler

        // to outsmart us and optimize away our memcpy calls.

        const size_t offset = 16 * (i & 15);

        const size_t copy_size = buffer_size - offset;

        memcpy(bufferB, bufferA + offset, copy_size);

        bytes_copied += copy_size;

      }

      eigen_assert(iters_at_a_time <= std::numeric_limits<size_t>::max() / 2);

      iters_at_a_time *= 2;

    }

  }

};

/*

 * A measurement_t is the result of a single experimental measurement of

 * sustained memory throughput for a buffer of a given size.

 */

struct measurement_t

{

  // The size of the memory buffer that this measurement used

  size_t size;

  // The sustained memory throughput that we measured, in GB/s

  float throughput;

  // Identifier of the cache level that this buffer size fits in.

  // The initial value 0 means that this has not yet been determined.

  // Identifiers are not ordered, and in particular are not

  // numerically equal to cache levels. Also, buffers too large to

  // fit in any cache still get an identifier, which corresponds to

  // RAM.

  int cache_level_id;

  measurement_t(size_t s = 0)

    : size(s)

    , throughput(0.0f)

    , cache_level_id(0)

  {}

};

/* Measures memory throughput for various buffer sizes.

 * This gathers all the experimental data that will be needed to

 * determine the cache structure.

 */

inline void measure_throughput_for_various_sizes(

              std::vector<measurement_t>& measurements)

{

  eigen_assert(measurements.empty());

  cache_measurement_context_t cache_measurement_context(implausibly_large_cache_size);

  for (size_t size = implausibly_small_cache_size;

       size <= implausibly_large_cache_size;

       size *= 2)

  {

    measurements.push_back(measurement_t(size));

  }

  for (size_t pass = 1; pass <= measurement_passes; pass++) {

    for (ptrdiff_t i = measurements.size() - 1; i >= 0; i--) {

      measurements[i].throughput = std::max(

        measurements[i].throughput,

        cache_measurement_context.measure_throughput(measurements[i].size));

      if (size_t(i) < measurements.size() - 1) {

        measurements[i].throughput = std::max(

          measurements[i].throughput,

          measurements[i+1].throughput);

      }

    }

  }

}

/* Partitions a vector of measurement levels into "cache levels", i.e.

 * groups having approximately the same throughput. The partitioning is

 * written into the cache_level_id fields of the measurement levels.

 */

inline void partition_measurements_into_cache_levels(

              std::vector<measurement_t>& measurements)

{

  eigen_assert(!measurements.empty());

  int current_cache_level_id = 0;

  int next_cache_level_id = 1;

  measurements.front().cache_level_id = next_cache_level_id++;

  measurements.back().cache_level_id = next_cache_level_id++;

  float threshold = std::log(never_separate_levels_closer_than);

  ptrdiff_t step = 1;

  // Step 1: we traverse the measurements back and forth, in the direction

  // of 'step'. In each traversal, we try to identify groups of

  // at least two consecutive measurements with approximately the same

  // throughput, as determined by the current value of 'threshold'. Such groups

  // are called cache levels, and are recorded in the 'cache_level' field

  // of each measurement. We continue until no two consecutive measurements

  // are left without a measurement.

  while(true) {

    ptrdiff_t start_level = step == 1 ? 0 : measurements.size() - 1;

    ptrdiff_t end_level = step == -1 ? 0 : measurements.size() - 1;

    for (ptrdiff_t i = start_level; i != end_level; i += step) {

      if (measurements[i].cache_level_id) {

        current_cache_level_id = measurements[i].cache_level_id;

        continue;

      }

      if (current_cache_level_id) {

        // Test if the current measurement can join the current cache level

        float ratio = measurements[i].throughput / measurements[i - step].throughput;

        if (std::abs(std::log(ratio)) < threshold) {

          measurements[i].cache_level_id = current_cache_level_id; 

        } else {

          current_cache_level_id = 0;

        }

      }

      if (!current_cache_level_id) {

        // Test if the current and next measurements are close enough to

        // create a new cache level for them

        if (!measurements[i + step].cache_level_id) {

          float ratio = measurements[i + step].throughput / measurements[i].throughput;

          if (std::abs(std::log(ratio)) < threshold)

          {

            current_cache_level_id = next_cache_level_id++;

            measurements[i].cache_level_id = current_cache_level_id;

          }

        }

      }

    }

    bool finished = true;

    for (size_t i = 0; i < measurements.size() - 1; i++) {

      // Test if two consecutive measurements are still left without a cache_level

      if (!measurements[i].cache_level_id && !measurements[i+1].cache_level_id) {

        finished = false;

      }

    }

    if (finished) {

      break;

    }

    step = -step;

    threshold *= 1.2f; // The value of this multiplier doesn't matter much.

                       // We just increase the threshold little by little

  }                    // until only isolated 1-wide gaps remain.

  // Step 2: there may still remain isolated measurements either without a cache level,

  // or with a cache level different from either of their neighbors'.

  // We resolve those isolated measurements by letting them join the cache level

  // immediately above or below, whichever is closer.

  bool might_have_isolated_level = true;

  while (might_have_isolated_level) {

    for (size_t i = 0; i < measurements.size(); i++) {

      bool isolated_below = i == 0 ||

        measurements[i].cache_level_id != measurements[i-1].cache_level_id;

      bool isolated_above = i == measurements.size() - 1 ||

        measurements[i].cache_level_id != measurements[i+1].cache_level_id;

      if (i == 0 && isolated_above) {

        measurements[0].cache_level_id = measurements[1].cache_level_id;

        continue;

      }

      if (i == measurements.size() - 1 && isolated_below) {

        measurements[measurements.size() - 1].cache_level_id =

          measurements[measurements.size() - 2].cache_level_id;

        continue;

      }

      if (isolated_below && isolated_above) {

        float distance_to_prev =

          std::abs(std::log(measurements[i-1].throughput / measurements[i].throughput));

        float distance_to_next =

          std::abs(std::log(measurements[i+1].throughput / measurements[i].throughput));

        if (distance_to_prev < distance_to_next) {

          measurements[i].cache_level_id = measurements[i-1].cache_level_id;

        } else {

          measurements[i].cache_level_id = measurements[i+1].cache_level_id;

        }

        continue;

      }

      might_have_isolated_level = false;

    }

  }

}

/* Takes a vector of measurement levels already partitioned into cache levels,

 * and outputs a plain vector of integers giving the size in bytes of each cache level.

 */

inline void get_cache_sizes_from_partitioned_measurements(

              std::vector<size_t>& cache_sizes,

              const std::vector<measurement_t>& measurements)

{

  eigen_assert(cache_sizes.empty());

  for (size_t i = 1; i < measurements.size() - 1; i++) {

    // should have generated cache_level_id's

    eigen_assert(measurements[i].cache_level_id);

    // pick biggest size in each cache level, and omit the last one, as it's

    // just RAM, not a real cache level.

    if (measurements[i].cache_level_id != measurements[i+1].cache_level_id &&

        measurements[i].cache_level_id != measurements.back().cache_level_id)

    {

      cache_sizes.push_back(measurements[i].size);

    }

  }

}

} // namespace internal

/**

 * Empirically determines the CPU's cache structure (number of levels and size of

 * each level) by running a benchmark. Takes roughly 1 second to complete.

 *

 * \param cache_sizes Output-parameter where to store the description of the

 *                     cache structure. This is a vector of size_t's, one

 *                     per cache level, giving its size in bytes.

 *

 * For example, on a Nexus 4 phone, the resulting cache_sizes vector is

 *   [16384, 1048576]  (size 2)

 * which means that there are 2 levels of cache: L1=16k, L2=1M.

 *

 * The sizes are approximate and may vary between runs. In particular,

 * only power-of-two sizes are ever returned.

 *

 * This function is only useful on systems where the cache structure cannot

 * be queried directly from the CPU. For example, userspace programs

 * on Android on ARM CPUs as of 2015.

 */

inline void measure_cache_sizes(std::vector<size_t>& cache_sizes)

{

  using namespace internal;

  std::vector<measurement_t> measurements;

  measure_throughput_for_various_sizes(measurements);

  partition_measurements_into_cache_levels(measurements);

  get_cache_sizes_from_partitioned_measurements(cache_sizes, measurements);

}

} // namespace Eigen

#endif // EIGEN_MEASURE_CACHE_SIZES_H