#+title: ATRIP: An MPI-asynchronous implementation of CCSD(T)
#+PROPERTY: header-args+ :noweb yes :comments noweb :mkdirp t

* Implementation

The algorithm uses two main data types, the =Slice= and the
=SliceUnion= as a container and resource manager of the =Slice=.

** The slice

#+begin_src c++ :tangle (atrip-slice-h)
#pragma once
#include <iostream>
#include <algorithm>
#include <vector>
#include <mpi.h>

#include <atrip/Tuples.hpp>
#include <atrip/Utils.hpp>
#include <atrip/Blas.hpp>

namespace atrip {

template <typename FF> FF maybeConjugate(const FF a) { return a; }
template <> Complex maybeConjugate(const Complex a) { return std::conj(a); }

namespace traits {
  template <typename FF> bool isComplex() { return false; };
  template <> bool isComplex<Complex>() { return true; };
namespace mpi {
  template <typename FF> MPI_Datatype datatypeOf(void);
  template <> MPI_Datatype datatypeOf<double>() { return MPI_DOUBLE; }
  template <> MPI_Datatype datatypeOf<Complex>() { return MPI_DOUBLE_COMPLEX; }
}
}


template <typename F=double>
struct Slice {

#+end_src

A slice is the concept of a subset of values of a given tensor.
As an example, for the doubles amplitudes \( T^{ab}_{ij} \), one need two kinds of objects:
- the object \( \mathsf{T}(a)^b_{ij} \) which for every \( a \) gets assigned the
  tensor \( T^{ab}_{ij} \) of size \( N_\mathrm{o}^2 \times N_\mathrm{v} \)

- the object \( \mathsf{T}(a,b)_{ij} \) which for every pair of \( a, b \)
  corresponds the \( N_\mathrm{o}^2 \)-sized tensor \( T^{ab}_{ij} \).



#+begin_src c++ :tangle (atrip-slice-h)
  // ASSOCIATED TYPES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  struct Location { size_t rank; size_t source; };

  enum Type
    { A = 10
    , B
    , C
    // Two-parameter slices
    , AB = 20
    , BC
    , AC
    // for abci and the doubles
    , CB
    , BA
    , CA
    // The non-typed slice
    , Blank = 404
    };

  enum State {
    // Fetch represents the state where a slice is to be fetched
    // and has a valid data pointer that can be written to
    Fetch = 0,
    // Dispatches represents the state that an MPI call has been
    // dispatched in order to get the data, but the data has not been
    // yet unwrapped, the data might be there or we might have to wait.
    Dispatched = 2,
    // Ready means that the data pointer can be read from
    Ready = 1,
    // Self sufficient is a slice when its contents are located
    // in the same rank that it lives, so that it does not have to
    // fetch from no one else.
    SelfSufficient = 911,
    // Recycled means that this slice gets its data pointer from another
    // slice, so it should not be written to
    Recycled = 123,
    // Acceptor means that the Slice can accept a new Slice, it is
    // the counterpart of the Blank type, but for states
    Acceptor = 405
  };

  struct Info {
    // which part of a,b,c the slice holds
    PartialTuple tuple;
    // The type of slice for the user to retrieve the correct one
    Type type;
    // What is the state of the slice
    State state;
    // Where the slice is to be retrieved
    // NOTE: this can actually be computed from tuple
    Location from;
    // If the data are actually to be found in this other slice
    Type recycling;

    Info() : tuple{0,0}
           , type{Blank}
           , state{Acceptor}
           , from{0,0}
           , recycling{Blank}
           {}
  };

  using Ty_x_Tu = std::pair< Type, PartialTuple >;

  // Names of the integrals that are considered in CCSD(T)
  enum Name
    { TA    = 100
    , VIJKA = 101
    , VABCI = 200
    , TABIJ = 201
    , VABIJ = 202
    };

  // DATABASE ==========================================================={{{1
  struct LocalDatabaseElement {
    Slice<F>::Name name;
    Slice<F>::Info info;
  };
  using LocalDatabase = std::vector<LocalDatabaseElement>;
  using Database = LocalDatabase;


    // STATIC METHODS ===========================================================
    //
    // They are useful to organize the structure of slices

    struct mpi {

      static MPI_Datatype vector(size_t n, MPI_Datatype const& DT) {
        MPI_Datatype dt;
        MPI_Type_vector(n, 1, 1, DT, &dt);
        MPI_Type_commit(&dt);
        return dt;
      }

      static MPI_Datatype sliceLocation () {
        constexpr int n = 2;
        // create a sliceLocation to measure in the current architecture
        // the packing of the struct
        Slice<F>::Location measure;
        MPI_Datatype dt;
        const std::vector<int> lengths(n, 1);
        const MPI_Datatype types[n] = {usizeDt(), usizeDt()};

        // measure the displacements in the struct
        size_t j = 0;
        MPI_Aint displacements[n];
        MPI_Get_address(&measure.rank,   &displacements[j++]);
        MPI_Get_address(&measure.source, &displacements[j++]);
        for (size_t i = 1; i < n; i++) displacements[i] -= displacements[0];
        displacements[0] = 0;

        MPI_Type_create_struct(n, lengths.data(), displacements, types, &dt);
        MPI_Type_commit(&dt);
        return dt;
      }

      static MPI_Datatype enumDt() { return MPI_INT; }
      static MPI_Datatype usizeDt() { return MPI_UINT64_T; }

      static MPI_Datatype sliceInfo () {
        constexpr int n = 5;
        MPI_Datatype dt;
        Slice<F>::Info measure;
        const std::vector<int> lengths(n, 1);
        const MPI_Datatype types[n]
          = { vector(2, usizeDt())
            , enumDt()
            , enumDt()
            , sliceLocation()
            , enumDt()
            };

        // create the displacements from the info measurement struct
        size_t j = 0;
        MPI_Aint displacements[n];
        MPI_Get_address(measure.tuple.data(), &displacements[j++]);
        MPI_Get_address(&measure.type,        &displacements[j++]);
        MPI_Get_address(&measure.state,       &displacements[j++]);
        MPI_Get_address(&measure.from,        &displacements[j++]);
        MPI_Get_address(&measure.recycling,   &displacements[j++]);
        for (size_t i = 1; i < n; i++) displacements[i] -= displacements[0];
        displacements[0] = 0;

        MPI_Type_create_struct(n, lengths.data(), displacements, types, &dt);
        MPI_Type_commit(&dt);
        return dt;
      }

      static MPI_Datatype localDatabaseElement () {
        constexpr int n = 2;
        MPI_Datatype dt;
        LocalDatabaseElement measure;
        const std::vector<int> lengths(n, 1);
        const MPI_Datatype types[n]
          = { enumDt()
            , sliceInfo()
            };

        // measure the displacements in the struct
        size_t j = 0;
        MPI_Aint displacements[n];
        MPI_Get_address(&measure.name, &displacements[j++]);
        MPI_Get_address(&measure.info, &displacements[j++]);
        for (size_t i = 1; i < n; i++) displacements[i] -= displacements[0];
        displacements[0] = 0;

        MPI_Type_create_struct(n, lengths.data(), displacements, types, &dt);
        MPI_Type_commit(&dt);
        return dt;
      }

    };

  static
  PartialTuple subtupleBySlice(ABCTuple abc, Type sliceType) {
    switch (sliceType) {
      case AB: return {abc[0], abc[1]};
      case BC: return {abc[1], abc[2]};
      case AC: return {abc[0], abc[2]};
      case CB: return {abc[2], abc[1]};
      case BA: return {abc[1], abc[0]};
      case CA: return {abc[2], abc[0]};
      case  A: return {abc[0], 0};
      case  B: return {abc[1], 0};
      case  C: return {abc[2], 0};
      default: throw "Switch statement not exhaustive!";
    }
  }


    /**
     ,* It is important here to return a reference to a Slice
     ,* not to accidentally copy the associated buffer of the slice.
     ,*/
    static Slice<F>& findOneByType(std::vector<Slice<F>> &slices, Slice<F>::Type type) {
        const auto sliceIt
          = std::find_if(slices.begin(), slices.end(),
                         [&type](Slice<F> const& s) {
                           return type == s.info.type;
                         });
        WITH_CRAZY_DEBUG
        WITH_RANK
          << "\t__ looking for " << type << "\n";
        if (sliceIt == slices.end())
          throw std::domain_error("Slice by type not found!");
        return *sliceIt;
    }

    /*
     ,* Check if an info has
     ,*
     ,*/
    static std::vector<Slice<F>*> hasRecycledReferencingToIt
      ( std::vector<Slice<F>> &slices
      , Info const& info
      ) {
      std::vector<Slice<F>*> result;

      for (auto& s: slices)
        if (  s.info.recycling == info.type
           && s.info.tuple == info.tuple
           && s.info.state == Recycled
           ) result.push_back(&s);

      return result;
    }

    static Slice<F>&
    findRecycledSource (std::vector<Slice<F>> &slices, Slice<F>::Info info) {
      const auto sliceIt
        = std::find_if(slices.begin(), slices.end(),
                       [&info](Slice<F> const& s) {
                         return info.recycling == s.info.type
                             && info.tuple == s.info.tuple
                             && State::Recycled != s.info.state
                             ;
                       });

      WITH_CRAZY_DEBUG
      WITH_RANK << "__slice__:find: recycling source of "
                << pretty_print(info) << "\n";
      if (sliceIt == slices.end())
        throw std::domain_error( "Slice not found: "
                               + pretty_print(info)
                               + " rank: "
                               + pretty_print(Atrip::rank)
                               );
      WITH_RANK << "__slice__:find: " << pretty_print(sliceIt->info) << "\n";
      return *sliceIt;
    }

    static Slice<F>& findByTypeAbc
      ( std::vector<Slice<F>> &slices
      , Slice<F>::Type type
      , ABCTuple const& abc
      ) {
        const auto tuple = Slice<F>::subtupleBySlice(abc, type);
        const auto sliceIt
          = std::find_if(slices.begin(), slices.end(),
                         [&type, &tuple](Slice<F> const& s) {
                           return type == s.info.type
                               && tuple == s.info.tuple
                               ;
                         });
        WITH_CRAZY_DEBUG
        WITH_RANK << "__slice__:find:" << type << " and tuple "
                  << pretty_print(tuple)
                  << "\n";
        if (sliceIt == slices.end())
          throw std::domain_error( "Slice not found: "
                                 + pretty_print(tuple)
                                 + ", "
                                 + pretty_print(type)
                                 + " rank: "
                                 + pretty_print(Atrip::rank)
                                 );
        return *sliceIt;
    }

    static Slice<F>& findByInfo(std::vector<Slice<F>> &slices,
                             Slice<F>::Info const& info) {
        const auto sliceIt
          = std::find_if(slices.begin(), slices.end(),
                         [&info](Slice<F> const& s) {
                           // TODO: maybe implement comparison in Info struct
                           return info.type == s.info.type
                               && info.state == s.info.state
                               && info.tuple == s.info.tuple
                               && info.from.rank == s.info.from.rank
                               && info.from.source == s.info.from.source
                                ;
                         });
        WITH_CRAZY_DEBUG
        WITH_RANK << "__slice__:find:looking for " << pretty_print(info) << "\n";
        if (sliceIt == slices.end())
          throw std::domain_error( "Slice by info not found: "
                                 + pretty_print(info));
        return *sliceIt;
    }

    // SLICE DEFINITION  =================================================={{{1

    // ATTRIBUTES ============================================================
    Info info;
    F  *data;
    MPI_Request request;
    const size_t size;

    void markReady() noexcept {
      info.state = Ready;
      info.recycling = Blank;
    }

    /*
     ,* This means that the data is there
     ,*/
    bool isUnwrapped() const noexcept {
      return info.state == Ready
          || info.state == SelfSufficient
          ;
    }

    bool isUnwrappable() const noexcept {
      return isUnwrapped()
          || info.state == Recycled
          || info.state == Dispatched
          ;
    }

    inline bool isDirectlyFetchable() const noexcept {
      return info.state == Ready || info.state == Dispatched;
    }

    void free() noexcept {
      info.tuple      = {0, 0};
      info.type       = Blank;
      info.state      = Acceptor;
      info.from       = {0, 0};
      info.recycling  = Blank;
      data            = nullptr;
    }

    inline bool isFree() const noexcept {
      return info.tuple       == PartialTuple{0, 0}
          && info.type        == Blank
          && info.state       == Acceptor
          && info.from.rank   == 0
          && info.from.source == 0
          && info.recycling   == Blank
          && data             == nullptr
           ;
    }


    /*
     ,* This function answers the question, which slices can be recycled.
     ,*
     ,* A slice can only be recycled if it is Fetch or Ready and has
     ,* a valid datapointer.
     ,*
     ,* In particular, SelfSufficient are not recyclable, since it is easier
     ,* just to create a SelfSufficient slice than deal with data dependencies.
     ,*
     ,* Furthermore, a recycled slice is not recyclable, if this is the case
     ,* then it is either bad design or a bug.
     ,*/
    inline bool isRecyclable() const noexcept {
      return (  info.state == Dispatched
             || info.state == Ready
             || info.state == Fetch
             )
          && hasValidDataPointer()
          ;
    }

    /*
     ,* This function describes if a slice has a valid data pointer.
     ,*
     ,* This is important to know if the slice has some data to it, also
     ,* some structural checks are done, so that it should not be Acceptor
     ,* or Blank, if this is the case then it is a bug.
     ,*/
    inline bool hasValidDataPointer() const noexcept {
      return data       != nullptr
          && info.state != Acceptor
          && info.type  != Blank
          ;
    }

    void unwrapAndMarkReady() {
      if (info.state == Ready) return;
      if (info.state != Dispatched)
        throw
          std::domain_error("Can't unwrap a non-ready, non-dispatched slice!");
      markReady();
      MPI_Status status;
#ifdef HAVE_OCD
        WITH_RANK << "__slice__:mpi: waiting " << "\n";
#endif
      const int errorCode = MPI_Wait(&request, &status);
      if (errorCode != MPI_SUCCESS)
        throw "MPI ERROR HAPPENED....";

#ifdef HAVE_OCD
      char errorString[MPI_MAX_ERROR_STRING];
      int errorSize;
      MPI_Error_string(errorCode, errorString, &errorSize);

      WITH_RANK << "__slice__:mpi: status "
                << "{ .source="    << status.MPI_SOURCE
                << ", .tag="       << status.MPI_TAG
                << ", .error="     << status.MPI_ERROR
                << ", .errCode="   << errorCode
                << ", .err="       << errorString
                << " }"
                << "\n";
#endif
    }

    Slice(size_t size_)
      : info({})
      , data(nullptr)
      , size(size_)
      {}


  }; // struct Slice


template <typename F=double>
std::ostream& operator<<(std::ostream& out, typename Slice<F>::Location const& v) {
  // TODO: remove me
  out << "{.r(" << v.rank << "), .s(" << v.source << ")};";
  return out;
}

template <typename F=double>
std::ostream& operator<<(std::ostream& out, typename Slice<F>::Info const& i) {
  out << "«t" << i.type << ", s" << i.state << "»"
      << " ⊙ {" << i.from.rank << ", " << i.from.source << "}"
      << " ∴ {" << i.tuple[0] << ", " << i.tuple[1] << "}"
      << " ♲t" << i.recycling
      ;
  return out;
}

} // namespace atrip
#+end_src

** Utils
#+begin_src c++ :tangle (atrip-utils-h)
#pragma once
#include <sstream>
#include <string>
#include <map>
#include <chrono>

#include <ctf.hpp>

namespace atrip {


  template <typename T>
  std::string pretty_print(T&& value) {
    std::stringstream stream;
#if ATRIP_DEBUG > 1
    dbg::pretty_print(stream, std::forward<T>(value));
#endif
    return stream.str();
  }

#define WITH_CHRONO(__chrono, ...) \
  __chrono.start(); __VA_ARGS__ __chrono.stop();

  struct Timer {
    using Clock = std::chrono::high_resolution_clock;
    using Event = std::chrono::time_point<Clock>;
    std::chrono::duration<double> duration;
    Event _start;
    inline void start() noexcept { _start = Clock::now(); }
    inline void stop() noexcept { duration += Clock::now() - _start; }
    inline void clear() noexcept { duration *= 0; }
    inline double count() const noexcept { return duration.count(); }
  };
  using Timings = std::map<std::string, Timer>;
}

#+end_src

** The rank mapping
#+begin_src c++ :tangle (atrip-rankmap-h)
#pragma once

#include <vector>
#include <algorithm>

#include <atrip/Slice.hpp>

namespace atrip {

  template <typename F=double>
  struct RankMap {

    std::vector<size_t> const lengths;
    size_t const np, size;

    RankMap(std::vector<size_t> lens, size_t np_)
      : lengths(lens)
      , np(np_)
      , size(std::accumulate(lengths.begin(), lengths.end(),
                            1UL, std::multiplies<size_t>()))
    { assert(lengths.size() <= 2); }

    size_t find(typename Slice<F>::Location const& p) const noexcept {
      return p.source * np + p.rank;
    }

    size_t nSources() const noexcept {
      return size / np + size_t(size % np != 0);
    }


    bool isPaddingRank(size_t rank) const noexcept {
      return size % np == 0
          ? false
          : rank > (size % np - 1)
          ;
    }

    bool isSourcePadding(size_t rank, size_t source) const noexcept {
      return source == nSources() && isPaddingRank(rank);
    }

    typename Slice<F>::Location
    find(ABCTuple const& abc, typename Slice<F>::Type sliceType) const noexcept {
      // tuple = {11, 8} when abc = {11, 8, 9} and sliceType = AB
      const auto tuple = Slice<F>::subtupleBySlice(abc, sliceType);

      const size_t index
        = tuple[0]
        + tuple[1] * (lengths.size() > 1 ? lengths[0] : 0)
        ;

      return
        { index % np
        , index / np
        };
    }

  };

}
#+end_src

** The slice union
#+begin_src c++ :tangle (atrip-slice-union-h)
#pragma once
#include <atrip/Debug.hpp>
#include <atrip/Slice.hpp>
#include <atrip/RankMap.hpp>

namespace atrip {

  template <typename F=double>
  struct SliceUnion {
    using Tensor = CTF::Tensor<F>;

    virtual void
    sliceIntoBuffer(size_t iteration, Tensor &to, Tensor const& from) = 0;

    /*
     * This function should enforce an important property of a SliceUnion.
     * Namely, there can be no two Slices of the same nature.
     *
     * This means that there can be at most one slice with a given Ty_x_Tu.
     */
    void checkForDuplicates() const {
      std::vector<typename Slice<F>::Ty_x_Tu> tytus;
      for (auto const& s: slices) {
        if (s.isFree()) continue;
        tytus.push_back({s.info.type, s.info.tuple});
      }

      for (auto const& tytu: tytus) {
        if (std::count(tytus.begin(), tytus.end(), tytu) > 1)
          throw "Invariance violated, more than one slice with same Ty_x_Tu";
      }

    }

    std::vector<typename Slice<F>::Ty_x_Tu> neededSlices(ABCTuple const& abc) {
      std::vector<typename Slice<F>::Ty_x_Tu> needed(sliceTypes.size());
      // build the needed vector
      std::transform(sliceTypes.begin(), sliceTypes.end(),
                     needed.begin(),
                     [&abc](typename Slice<F>::Type const type) {
                       auto tuple = Slice<F>::subtupleBySlice(abc, type);
                       return std::make_pair(type, tuple);
                     });
      return needed;
    }

    /* buildLocalDatabase
     *
     * It should build a database of slices so that we know what is needed
     * to fetch in the next iteration represented by the tuple 'abc'.
     *
     * 1. The algorithm works as follows, we build a database of the all
     * the slice types that we need together with their tuple.
     *
     * 2. Look in the SliceUnion if we already have this tuple,
     * if we already have it mark it (TODO)
     *
     * 3. If we don't have the tuple, look for a (state=acceptor, type=blank)
     * slice and mark this slice as type=Fetch with the corresponding type
     * and tuple.
     *
     * NOTE: The algorithm should certify that we always have enough blank
     * slices.
     *
     */
    typename
    Slice<F>::LocalDatabase buildLocalDatabase(ABCTuple const& abc) {
      typename Slice<F>::LocalDatabase result;

      auto const needed = neededSlices(abc);

      WITH_RANK << "__db__:needed:" << pretty_print(needed) << "\n";
      // BUILD THE DATABASE
      // we need to loop over all sliceTypes that this TensorUnion
      // is representing and find out how we will get the corresponding
      // slice for the abc we are considering right now.
      for (auto const& pair: needed) {
        auto const type = pair.first;
        auto const tuple = pair.second;
        auto const from  = rankMap.find(abc, type);

#ifdef HAVE_OCD
        WITH_RANK << "__db__:want:" << pretty_print(pair) << "\n";
        for (auto const& s: slices)
          WITH_RANK << "__db__:guts:ocd "
                    << s.info << " pt " << s.data
                    << "\n";
#endif

#ifdef HAVE_OCD
        WITH_RANK << "__db__: checking if exact match" << "\n";
#endif
        {
          // FIRST: look up if there is already a *Ready* slice matching what we
          // need
          auto const& it
            = std::find_if(slices.begin(), slices.end(),
                           [&tuple, &type](Slice<F> const& other) {
                             return other.info.tuple == tuple
                                 && other.info.type == type
                                    // we only want another slice when it
                                    // has already ready-to-use data
                                 && other.isUnwrappable()
                                 ;
                           });
          if (it != slices.end()) {
            // if we find this slice, it means that we don't have to do anything
            WITH_RANK << "__db__: EXACT: found EXACT in name=" << name
                      << " for tuple " << tuple[0] << ", " << tuple[1]
                      << " ptr " << it->data
                      << "\n";
            result.push_back({name, it->info});
            continue;
          }
        }

#ifdef HAVE_OCD
        WITH_RANK << "__db__: checking if recycle" << "\n";
#endif
        // Try to find a recyling possibility ie. find a slice with the same
        // tuple and that has a valid data pointer.
        auto const& recycleIt
          = std::find_if(slices.begin(), slices.end(),
                         [&tuple, &type](Slice<F> const& other) {
                           return other.info.tuple == tuple
                               && other.info.type != type
                               && other.isRecyclable()
                               ;
                         });

        // if we find this recylce, then we find a Blank slice
        // (which should exist by construction :THINK)
        //
        if (recycleIt != slices.end()) {
          auto& blank = Slice<F>::findOneByType(slices, Slice<F>::Blank);
          // TODO: formalize this through a method to copy information
          //       from another slice
          blank.data = recycleIt->data;
          blank.info.type = type;
          blank.info.tuple = tuple;
          blank.info.state = Slice<F>::Recycled;
          blank.info.from = from;
          blank.info.recycling = recycleIt->info.type;
          result.push_back({name, blank.info});
          WITH_RANK << "__db__: RECYCLING: n" << name
                    << " " << pretty_print(abc)
                    << " get " << pretty_print(blank.info)
                    << " from " << pretty_print(recycleIt->info)
                    << " ptr " << recycleIt->data
                    << "\n"
                    ;
          continue;
        }

        // in this case we have to create a new slice
        // this means that we should have a blank slice at our disposal
        // and also the freePointers should have some elements inside,
        // so we pop a data pointer from the freePointers container
#ifdef HAVE_OCD
        WITH_RANK << "__db__: none work, doing new" << "\n";
#endif
        {
          WITH_RANK << "__db__: NEW: finding blank in " << name
                    << " for type " << type
                    << " for tuple " << tuple[0] << ", " << tuple[1]
                    << "\n"
                    ;
          auto& blank = Slice<F>::findOneByType(slices, Slice<F>::Blank);
          blank.info.type = type;
          blank.info.tuple = tuple;
          blank.info.from = from;

          // Handle self sufficiency
          blank.info.state = Atrip::rank == from.rank
                           ? Slice<F>::SelfSufficient
                           : Slice<F>::Fetch
                           ;
          if (blank.info.state == Slice<F>::SelfSufficient) {
            blank.data = sources[from.source].data();
          } else {
            if (freePointers.size() == 0)
              throw std::domain_error("No more free pointers!");
            auto dataPointer = freePointers.begin();
            freePointers.erase(dataPointer);
            blank.data = *dataPointer;
          }

          result.push_back({name, blank.info});
          continue;
        }

      }

#ifdef HAVE_OCD
      for (auto const& s: slices)
        WITH_RANK << "__db__:guts:ocd:__end__ " << s.info << "\n";
#endif


      return result;

    }

    /*
     * Garbage collect slices not needed for the next iteration.
     *
     * It will throw if it tries to gc a slice that has not been
     * previously unwrapped, as a safety mechanism.
     */
    void clearUnusedSlicesForNext(ABCTuple const& abc) {
      auto const needed = neededSlices(abc);

      // CLEAN UP SLICES, FREE THE ONES THAT ARE NOT NEEDED ANYMORE
      for (auto& slice: slices) {
        // if the slice is free, then it was not used anyways
        if (slice.isFree()) continue;


        // try to find the slice in the needed slices list
        auto const found
          = std::find_if(needed.begin(), needed.end(),
                         [&slice] (typename Slice<F>::Ty_x_Tu const& tytu) {
                           return slice.info.tuple == tytu.second
                               && slice.info.type == tytu.first
                               ;
                         });

        // if we did not find slice in needed, then erase it
        if (found == needed.end()) {

          // We have to be careful about the data pointer,
          // for SelfSufficient, the data pointer is a source pointer
          // of the slice, so we should just wipe it.
          //
          // For Ready slices, we have to be careful if there are some
          // recycled slices depending on it.
          bool freeSlicePointer = true;

          // allow to gc unwrapped and recycled, never Fetch,
          // if we have a Fetch slice then something has gone very wrong.
          if (!slice.isUnwrapped() && slice.info.state != Slice<F>::Recycled)
            throw
              std::domain_error("Trying to garbage collect "
                                " a non-unwrapped slice! "
                                + pretty_print(&slice)
                                + pretty_print(slice.info));

          // it can be that our slice is ready, but it has some hanging
          // references lying around in the form of a recycled slice.
          // Of course if we need the recycled slice the next iteration
          // this would be fatal, because we would then free the pointer
          // of the slice and at some point in the future we would
          // overwrite it. Therefore, we must check if slice has some
          // references in slices and if so then
          //
          //  - we should mark those references as the original (since the data
          //    pointer should be the same)
          //
          //  - we should make sure that the data pointer of slice
          //    does not get freed.
          //
          if (slice.info.state == Slice<F>::Ready) {
            WITH_OCD WITH_RANK
              << "__gc__:" << "checking for data recycled dependencies\n";
            auto recycled
              = Slice<F>::hasRecycledReferencingToIt(slices, slice.info);
            if (recycled.size()) {
              Slice<F>* newReady = recycled[0];
              WITH_OCD WITH_RANK
                << "__gc__:" << "swaping recycled "
                << pretty_print(newReady->info)
                << " and "
                << pretty_print(slice.info)
                << "\n";
              newReady->markReady();
              assert(newReady->data == slice.data);
              freeSlicePointer = false;

              for (size_t i = 1; i < recycled.size(); i++) {
                auto newRecyled = recycled[i];
                newRecyled->info.recycling = newReady->info.type;
                WITH_OCD WITH_RANK
                  << "__gc__:" << "updating recycled "
                  << pretty_print(newRecyled->info)
                  << "\n";
              }

            }
          }

          // if the slice is self sufficient, do not dare touching the
          // pointer, since it is a pointer to our sources in our rank.
          if (  slice.info.state == Slice<F>::SelfSufficient
             || slice.info.state == Slice<F>::Recycled
             ) {
            freeSlicePointer = false;
          }

          // make sure we get its data pointer to be used later
          // only for non-recycled, since it can be that we need
          // for next iteration the data of the slice that the recycled points
          // to
          if (freeSlicePointer) {
            freePointers.insert(slice.data);
            WITH_RANK << "~~~:cl(" << name << ")"
                      << " added to freePointer "
                      << pretty_print(freePointers)
                      << "\n";
          } else {
            WITH_OCD WITH_RANK << "__gc__:not touching the free Pointer\n";
          }

          // at this point, let us blank the slice
          WITH_RANK << "~~~:cl(" << name << ")"
                    << " freeing up slice "
                    // TODO: make this possible
                    // << " info " << slice.info
                    << "\n";
          slice.free();
        }

      }
    }

    // CONSTRUCTOR
    SliceUnion( Tensor const& sourceTensor
              , std::vector<typename Slice<F>::Type> sliceTypes_
              , std::vector<size_t> sliceLength_
              , std::vector<size_t> paramLength
              , size_t np
              , MPI_Comm child_world
              , MPI_Comm global_world
              , typename Slice<F>::Name name_
              , size_t nSliceBuffers = 4
              )
              : rankMap(paramLength, np)
              , world(child_world)
              , universe(global_world)
              , sliceLength(sliceLength_)
              , sources(rankMap.nSources(),
                        std::vector<F>
                          (std::accumulate(sliceLength.begin(),
                                           sliceLength.end(),
                                           1UL, std::multiplies<size_t>())))
              , name(name_)
              , sliceTypes(sliceTypes_)
              , sliceBuffers(nSliceBuffers, sources[0])
              //, slices(2 * sliceTypes.size(), Slice<F>{ sources[0].size() })
    { // constructor begin

      LOG(0,"Atrip") << "INIT SliceUnion: " << name << "\n";

      slices
        = std::vector<Slice<F>>(2 * sliceTypes.size(), { sources[0].size() });
      // TODO: think exactly ^------------------- about this number

      // initialize the freePointers with the pointers to the buffers
      std::transform(sliceBuffers.begin(), sliceBuffers.end(),
                     std::inserter(freePointers, freePointers.begin()),
                     [](std::vector<F> &vec) { return vec.data(); });



      LOG(1,"Atrip") << "rankMap.nSources "
                           << rankMap.nSources() << "\n";
      LOG(1,"Atrip") << "#slices "
                           << slices.size() << "\n";
      LOG(1,"Atrip") << "#slices[0] "
                           << slices[0].size << "\n";
      LOG(1,"Atrip") << "#sources "
                           << sources.size() << "\n";
      LOG(1,"Atrip") << "#sources[0] "
                           << sources[0].size() << "\n";
      LOG(1,"Atrip") << "#freePointers "
                           << freePointers.size() << "\n";
      LOG(1,"Atrip") << "#sliceBuffers "
                           << sliceBuffers.size() << "\n";
      LOG(1,"Atrip") << "#sliceBuffers[0] "
                           << sliceBuffers[0].size() << "\n";
      LOG(1,"Atrip") << "#sliceLength "
                           << sliceLength.size() << "\n";
      LOG(1,"Atrip") << "#paramLength "
                           << paramLength.size() << "\n";
      LOG(1,"Atrip") << "GB*" << np << " "
                           << double(sources.size() + sliceBuffers.size())
                            * sources[0].size()
                            * 8 * np
                            / 1073741824.0
                           << "\n";
    } // constructor ends

    void init(Tensor const& sourceTensor) {

      CTF::World w(world);
      const int rank = Atrip::rank
              , order = sliceLength.size()
              ;
      std::vector<int> const syms(order, NS);
      std::vector<int> __sliceLength(sliceLength.begin(), sliceLength.end());
      Tensor toSliceInto(order,
                         __sliceLength.data(),
                         syms.data(),
                         w);
      LOG(1,"Atrip") << "slicing... \n";

      // setUp sources
      for (size_t it(0); it < rankMap.nSources(); ++it) {
        const size_t
          source = rankMap.isSourcePadding(rank, source) ? 0 : it;
        WITH_OCD
        WITH_RANK
          << "Init:toSliceInto it-" << it
          << " :: source " << source << "\n";
        sliceIntoBuffer(source, toSliceInto, sourceTensor);
      }

    }

    /**
     * \brief Send asynchronously only if the state is Fetch
     */
    void send( size_t otherRank
             , typename Slice<F>::Info const& info
             , size_t tag) const noexcept {
      MPI_Request request;
      bool sendData_p = false;

      if (info.state == Slice<F>::Fetch) sendData_p = true;
      // TODO: remove this because I have SelfSufficient
      if (otherRank == info.from.rank)      sendData_p = false;
      if (!sendData_p) return;

      MPI_Isend( sources[info.from.source].data()
               , sources[info.from.source].size()
               , traits::mpi::datatypeOf<F>()
               , otherRank
               , tag
               , universe
               , &request
               );
      WITH_CRAZY_DEBUG
      WITH_RANK << "sent to " << otherRank << "\n";

    }

    /**
     * \brief Receive asynchronously only if the state is Fetch
     */
    void receive(typename Slice<F>::Info const& info, size_t tag) noexcept {
      auto& slice = Slice<F>::findByInfo(slices, info);

      if (Atrip::rank == info.from.rank) return;

      if (slice.info.state == Slice<F>::Fetch) {
        // TODO: do it through the slice class
        slice.info.state = Slice<F>::Dispatched;
        MPI_Request request;
        slice.request = request;
        MPI_Irecv( slice.data
                 , slice.size
                 , traits::mpi::datatypeOf<F>()
                 , info.from.rank
                 , tag
                 , universe
                 , &slice.request
                //, MPI_STATUS_IGNORE
                 );
      }
    }

    void unwrapAll(ABCTuple const& abc) {
      for (auto type: sliceTypes) unwrapSlice(type, abc);
    }

    F* unwrapSlice(typename Slice<F>::Type type, ABCTuple const& abc) {
      WITH_CRAZY_DEBUG
      WITH_RANK << "__unwrap__:slice " << type << " w n "
                << name
                << " abc" << pretty_print(abc)
                << "\n";
      auto& slice = Slice<F>::findByTypeAbc(slices, type, abc);
      //WITH_RANK << "__unwrap__:info " << slice.info << "\n";
      switch  (slice.info.state) {
        case Slice<F>::Dispatched:
          WITH_RANK << "__unwrap__:Fetch: " << &slice
                    << " info " << pretty_print(slice.info)
                    << "\n";
          slice.unwrapAndMarkReady();
          return slice.data;
          break;
        case Slice<F>::SelfSufficient:
          WITH_RANK << "__unwrap__:SelfSufficient: " << &slice
                    << " info " << pretty_print(slice.info)
                    << "\n";
          return slice.data;
          break;
        case Slice<F>::Ready:
          WITH_RANK << "__unwrap__:READY: UNWRAPPED ALREADY" << &slice
                    << " info " << pretty_print(slice.info)
                    << "\n";
          return slice.data;
          break;
        case Slice<F>::Recycled:
          WITH_RANK << "__unwrap__:RECYCLED " << &slice
                    << " info " << pretty_print(slice.info)
                    << "\n";
          return unwrapSlice(slice.info.recycling, abc);
          break;
        case Slice<F>::Fetch:
        case Slice<F>::Acceptor:
          throw std::domain_error("Can't unwrap an acceptor or fetch slice!");
          break;
        default:
          throw std::domain_error("Unknown error unwrapping slice!");
      }
      return slice.data;
    }

    const RankMap<F> rankMap;
    const MPI_Comm world;
    const MPI_Comm universe;
    const std::vector<size_t> sliceLength;
    std::vector< std::vector<F> > sources;
    std::vector< Slice<F> > slices;
    typename Slice<F>::Name name;
    const std::vector<typename Slice<F>::Type> sliceTypes;
    std::vector< std::vector<F> > sliceBuffers;
    std::set<F*> freePointers;

  };

  template <typename F=double>
  SliceUnion<F>&
  unionByName(std::vector<SliceUnion<F>*> const& unions,
              typename Slice<F>::Name name) {
      const auto sliceUnionIt
        = std::find_if(unions.begin(), unions.end(),
                      [&name](SliceUnion<F> const* s) {
                        return name == s->name;
                      });
      if (sliceUnionIt == unions.end())
        throw std::domain_error("SliceUnion not found!");
      return **sliceUnionIt;
  }

}
#+end_src

** Tuples
#+begin_src c++ :tangle (atrip-tuples-h)
#pragma once

#include <vector>
#include <array>
#include <numeric>

#include <atrip/Utils.hpp>
#include <atrip/Debug.hpp>

namespace atrip {

  using ABCTuple = std::array<size_t, 3>;
  using PartialTuple = std::array<size_t, 2>;
  using ABCTuples = std::vector<ABCTuple>;

  ABCTuples getTuplesList(size_t Nv) {
    const size_t n = Nv * (Nv + 1) * (Nv + 2) / 6 - Nv;
    ABCTuples result(n);
    size_t u(0);

    for (size_t a(0); a < Nv; a++)
    for (size_t b(a); b < Nv; b++)
    for (size_t c(b); c < Nv; c++){
      if ( a == b && b == c ) continue;
      result[u++] = {a, b, c};
    }

    return result;

  }


  std::pair<size_t, size_t>
  getABCRange(size_t np, size_t rank, ABCTuples const& tuplesList) {

    std::vector<size_t> n_tuples_per_rank(np, tuplesList.size()/np);
    const size_t
        // how many valid tuples should we still verteilen to nodes
        // since the number of tuples is not divisible by the number of nodes
        nRoundRobin = tuplesList.size() % np
        // every node must have the sanme amount of tuples in order for the
        // other nodes to receive and send somewhere, therefore
        // some nodes will get extra tuples but that are dummy tuples
      , nExtraInvalid = (np - nRoundRobin) % np
      ;

    if (nRoundRobin) for (int i = 0; i < np; i++) n_tuples_per_rank[i]++;

  #if defined(TODO)
    assert( tuplesList.size()
            ==
            ( std::accumulate(n_tuples_per_rank.begin(),
                              n_tuples_per_rank.end(),
                              0UL,
                              std::plus<size_t>())
            + nExtraInvalid
            ));
  #endif

    WITH_RANK << "nRoundRobin = " << nRoundRobin << "\n";
    WITH_RANK << "nExtraInvalid = " << nExtraInvalid << "\n";
    WITH_RANK << "ntuples = " << n_tuples_per_rank[rank] << "\n";

    auto const& it = n_tuples_per_rank.begin();

    return
      { std::accumulate(it, it + rank    , 0)
      , std::accumulate(it, it + rank + 1, 0)
      };

  }

}
#+end_src

** Unions
Since every tensor slice in a different way, we can override the slicing procedure
and define subclasses of slice unions.

#+begin_src c++ :tangle (atrip-unions-h)
#pragma once
#include <atrip/SliceUnion.hpp>

namespace atrip {

  template <typename F=double>
  void sliceIntoVector
    ( std::vector<F> &v
    , CTF::Tensor<F> &toSlice
    , std::vector<int64_t> const low
    , std::vector<int64_t> const up
    , CTF::Tensor<F> const& origin
    , std::vector<int64_t> const originLow
    , std::vector<int64_t> const originUp
    ) {
    // Thank you CTF for forcing me to do this
    struct { std::vector<int> up, low; }
        toSlice_ = { {up.begin(), up.end()}
                   , {low.begin(), low.end()} }
      , origin_ = { {originUp.begin(), originUp.end()}
                  , {originLow.begin(), originLow.end()} }
      ;

    WITH_OCD
    WITH_RANK << "slicing into " << pretty_print(toSlice_.up)
                          << "," << pretty_print(toSlice_.low)
              << " from " << pretty_print(origin_.up)
                   << "," << pretty_print(origin_.low)
              << "\n";

#ifndef ATRIP_DONT_SLICE
    toSlice.slice( toSlice_.low.data()
                 , toSlice_.up.data()
                 , 0.0
                 , origin
                 , origin_.low.data()
                 , origin_.up.data()
                 , 1.0);
    memcpy(v.data(), toSlice.data, sizeof(F) * v.size());
#endif

  }


  template <typename F=double>
  struct TAPHH : public SliceUnion<F> {
    TAPHH( CTF::Tensor<F> const& sourceTensor
         , size_t No
         , size_t Nv
         , size_t np
         , MPI_Comm child_world
         , MPI_Comm global_world
         ) : SliceUnion<F>( sourceTensor
                          , {Slice<F>::A, Slice<F>::B, Slice<F>::C}
                          , {Nv, No, No} // size of the slices
                          , {Nv}
                          , np
                          , child_world
                          , global_world
                          , Slice<F>::TA
                          , 4) {
           init(sourceTensor);
         }

    void sliceIntoBuffer(size_t it, CTF::Tensor<F> &to, CTF::Tensor<F> const& from) override
    {
      const int Nv = this->sliceLength[0]
              , No = this->sliceLength[1]
              , a = this->rankMap.find({static_cast<size_t>(Atrip::rank), it});
              ;


      sliceIntoVector<F>( this->sources[it]
                        , to,   {0, 0, 0},    {Nv, No, No}
                        , from, {a, 0, 0, 0}, {a+1, Nv, No, No}
                        );

    }

  };


  template <typename F=double>
  struct HHHA : public SliceUnion<F> {
    HHHA( CTF::Tensor<F> const& sourceTensor
        , size_t No
        , size_t Nv
        , size_t np
        , MPI_Comm child_world
        , MPI_Comm global_world
        ) : SliceUnion<F>( sourceTensor
                         , {Slice<F>::A, Slice<F>::B, Slice<F>::C}
                         , {No, No, No} // size of the slices
                         , {Nv}         // size of the parametrization
                         , np
                         , child_world
                         , global_world
                         , Slice<F>::VIJKA
                         , 4) {
           init(sourceTensor);
         }

    void sliceIntoBuffer(size_t it, CTF::Tensor<F> &to, CTF::Tensor<F> const& from) override
    {

      const int No = this->sliceLength[0]
              , a = this->rankMap.find({static_cast<size_t>(Atrip::rank), it})
              ;

      sliceIntoVector<F>( this->sources[it]
                        , to,   {0, 0, 0},    {No, No, No}
                        , from, {0, 0, 0, a}, {No, No, No, a+1}
                        );

    }
  };

  template <typename F=double>
  struct ABPH : public SliceUnion<F> {
    ABPH( CTF::Tensor<F> const& sourceTensor
        , size_t No
        , size_t Nv
        , size_t np
        , MPI_Comm child_world
        , MPI_Comm global_world
        ) : SliceUnion<F>( sourceTensor
                         , { Slice<F>::AB, Slice<F>::BC, Slice<F>::AC
                           , Slice<F>::BA, Slice<F>::CB, Slice<F>::CA
                           }
                         , {Nv, No} // size of the slices
                         , {Nv, Nv} // size of the parametrization
                         , np
                         , child_world
                         , global_world
                         , Slice<F>::VABCI
                         , 2*6) {
           init(sourceTensor);
         }

    void sliceIntoBuffer(size_t it, CTF::Tensor<F> &to, CTF::Tensor<F> const& from) override {

      const int Nv = this->sliceLength[0]
              , No = this->sliceLength[1]
              , el = this->rankMap.find({static_cast<size_t>(Atrip::rank), it})
              , a = el % Nv
              , b = el / Nv
              ;


      sliceIntoVector<F>( this->sources[it]
                        , to,   {0, 0},       {Nv, No}
                        , from, {a, b, 0, 0}, {a+1, b+1, Nv, No}
                        );

    }

  };

  template <typename F=double>
  struct ABHH : public SliceUnion<F> {
    ABHH( CTF::Tensor<F> const& sourceTensor
        , size_t No
        , size_t Nv
        , size_t np
        , MPI_Comm child_world
        , MPI_Comm global_world
        ) : SliceUnion<F>( sourceTensor
                         , {Slice<F>::AB, Slice<F>::BC, Slice<F>::AC}
                         , {No, No} // size of the slices
                         , {Nv, Nv} // size of the parametrization
                         , np
                         , child_world
                         , global_world
                         , Slice<F>::VABIJ
                         , 6) {
           init(sourceTensor);
         }

    void sliceIntoBuffer(size_t it, CTF::Tensor<F> &to, CTF::Tensor<F> const& from) override {

      const int Nv = from.lens[0]
              , No = this->sliceLength[1]
              , el = this->rankMap.find({static_cast<size_t>(Atrip::rank), it})
              , a = el % Nv
              , b = el / Nv
              ;

      sliceIntoVector<F>( this->sources[it]
                        , to,   {0, 0},       {No, No}
                        , from, {a, b, 0, 0}, {a+1, b+1, No, No}
                        );


    }

  };


  template <typename F=double>
  struct TABHH : public SliceUnion<F> {
    TABHH( CTF::Tensor<F> const& sourceTensor
         , size_t No
         , size_t Nv
         , size_t np
         , MPI_Comm child_world
         , MPI_Comm global_world
         ) : SliceUnion<F>( sourceTensor
                          , {Slice<F>::AB, Slice<F>::BC, Slice<F>::AC}
                          , {No, No} // size of the slices
                          , {Nv, Nv} // size of the parametrization
                          , np
                          , child_world
                          , global_world
                          , Slice<F>::TABIJ
                          , 6) {
           init(sourceTensor);
         }

    void sliceIntoBuffer(size_t it, CTF::Tensor<F> &to, CTF::Tensor<F> const& from) override {
      // TODO: maybe generalize this with ABHH

      const int Nv = from.lens[0]
              , No = this->sliceLength[1]
              , el = this->rankMap.find({static_cast<size_t>(Atrip::rank), it})
              , a = el % Nv
              , b = el / Nv
              ;

      sliceIntoVector<F>( this->sources[it]
                        , to,   {0, 0},       {No, No}
                        , from, {a, b, 0, 0}, {a+1, b+1, No, No}
                        );


    }

  };

}
#+end_src


** Equations
#+begin_src c++ :tangle (atrip-equations-h)
#pragma once

#include<atrip/Slice.hpp>
#include<atrip/Blas.hpp>

namespace atrip {

  template <typename F=double>
  double getEnergyDistinct
    ( const F epsabc
    , std::vector<F> const& epsi
    , std::vector<F> const& Tijk_
    , std::vector<F> const& Zijk_
    ) {
    constexpr size_t blockSize=16;
    F energy(0.);
    const size_t No = epsi.size();
    for (size_t kk=0; kk<No; kk+=blockSize){
      const size_t kend( std::min(No, kk+blockSize) );
      for (size_t jj(kk); jj<No; jj+=blockSize){
        const size_t jend( std::min( No, jj+blockSize) );
        for (size_t ii(jj); ii<No; ii+=blockSize){
          const size_t iend( std::min( No, ii+blockSize) );
          for (size_t k(kk); k < kend; k++){
            const F ek(epsi[k]);
            const size_t jstart = jj > k ? jj : k;
            for (size_t j(jstart); j < jend; j++){
              F const ej(epsi[j]);
              F const facjk = j == k ? F(0.5) : F(1.0);
              size_t istart = ii > j ? ii : j;
              for (size_t i(istart); i < iend; i++){
                const F
                    ei(epsi[i])
                  , facij = i == j ? F(0.5) : F(1.0)
                  , denominator(epsabc - ei - ej - ek)
                  , U(Zijk_[i + No*j + No*No*k])
                  , V(Zijk_[i + No*k + No*No*j])
                  , W(Zijk_[j + No*i + No*No*k])
                  , X(Zijk_[j + No*k + No*No*i])
                  , Y(Zijk_[k + No*i + No*No*j])
                  , Z(Zijk_[k + No*j + No*No*i])
                  , A(maybeConjugate<F>(Tijk_[i + No*j + No*No*k]))
                  , B(maybeConjugate<F>(Tijk_[i + No*k + No*No*j]))
                  , C(maybeConjugate<F>(Tijk_[j + No*i + No*No*k]))
                  , D(maybeConjugate<F>(Tijk_[j + No*k + No*No*i]))
                  , E(maybeConjugate<F>(Tijk_[k + No*i + No*No*j]))
                  , F(maybeConjugate<F>(Tijk_[k + No*j + No*No*i]))
                  , value
                    = 3.0 * ( A * U
                              + B * V
                              + C * W
                              + D * X
                              + E * Y
                              + F * Z )
                   + ( ( U + X + Y )
                     - 2.0 * ( V + W + Z )
                     ) * ( A + D + E )
                   + ( ( V + W + Z )
                     - 2.0 * ( U + X + Y )
                     ) * ( B + C + F )
                  ;
                energy += 2.0 * value / denominator * facjk * facij;
              } // i
            } // j
          } // k
        } // ii
      } // jj
    } // kk
    return std::real(energy);
  }


  template <typename F=double>
  double getEnergySame
    ( const F epsabc
    , std::vector<F> const& epsi
    , std::vector<F> const& Tijk_
    , std::vector<F> const& Zijk_
    ) {
    constexpr size_t blockSize = 16;
    const size_t No = epsi.size();
    F energy = F(0.);
    for (size_t kk=0; kk<No; kk+=blockSize){
      const size_t kend( std::min( kk+blockSize, No) );
      for (size_t jj(kk); jj<No; jj+=blockSize){
        const size_t jend( std::min( jj+blockSize, No) );
        for (size_t ii(jj); ii<No; ii+=blockSize){
          const size_t iend( std::min( ii+blockSize, No) );
          for (size_t k(kk); k < kend; k++){
            const F ek(epsi[k]);
            const size_t jstart = jj > k ? jj : k;
            for(size_t j(jstart); j < jend; j++){
              const F facjk( j == k ? F(0.5) : F(1.0));
              const F ej(epsi[j]);
              const size_t istart = ii > j ? ii : j;
              for(size_t i(istart); i < iend; i++){
                const F
                  ei(epsi[i])
                , facij ( i==j ? F(0.5) : F(1.0))
                , denominator(epsabc - ei - ej - ek)
                , U(Zijk_[i + No*j + No*No*k])
                , V(Zijk_[j + No*k + No*No*i])
                , W(Zijk_[k + No*i + No*No*j])
                , A(maybeConjugate<F>(Tijk_[i + No*j + No*No*k]))
                , B(maybeConjugate<F>(Tijk_[j + No*k + No*No*i]))
                , C(maybeConjugate<F>(Tijk_[k + No*i + No*No*j]))
                , value
                  = F(3.0) * ( A * U
                             + B * V
                             + C * W
                             )
                  - ( A + B + C ) * ( U + V + W )
                ;
                energy += F(2.0) * value / denominator * facjk * facij;
              } // i
            } // j
          } // k
        } // ii
      } // jj
    } // kk
    return std::real(energy);
  }

  template <typename F=double>
  void singlesContribution
    ( size_t No
    , size_t Nv
    , const ABCTuple &abc
    , F const* Tph
    , F const* VABij
    , F const* VACij
    , F const* VBCij
    , F *Zijk
    ) {
    const size_t a(abc[0]), b(abc[1]), c(abc[2]);
    for (size_t k=0; k < No; k++)
    for (size_t i=0; i < No; i++)
    for (size_t j=0; j < No; j++) {
      const size_t ijk = i + j*No + k*No*No
                ,  jk = j + No * k
                ;
      Zijk[ijk] += Tph[ a + i * Nv ] * VBCij[ j + k * No ];
      Zijk[ijk] += Tph[ b + j * Nv ] * VACij[ i + k * No ];
      Zijk[ijk] += Tph[ c + k * Nv ] * VABij[ i + j * No ];
    }
  }

  template <typename F=double>
  void doublesContribution
    ( const ABCTuple &abc
    , size_t const No
    , size_t const Nv
    // -- VABCI
    , F const* VABph
    , F const* VACph
    , F const* VBCph
    , F const* VBAph
    , F const* VCAph
    , F const* VCBph
    // -- VHHHA
    , F const* VhhhA
    , F const* VhhhB
    , F const* VhhhC
    // -- TA
    , F const* TAphh
    , F const* TBphh
    , F const* TCphh
    // -- TABIJ
    , F const* TABhh
    , F const* TAChh
    , F const* TBChh
    // -- TIJK
    , F *Tijk
    , atrip::Timings& chrono
    ) {

    auto& t_reorder = chrono["doubles:reorder"];
    const size_t a = abc[0], b = abc[1], c = abc[2]
              , NoNo = No*No, NoNv = No*Nv
              ;

  #if defined(ATRIP_USE_DGEMM)
  #define _IJK_(i, j, k) i + j*No + k*NoNo
  #define REORDER(__II, __JJ, __KK)                                 \
    t_reorder.start();                                              \
    for (size_t k = 0; k < No; k++)                                 \
    for (size_t j = 0; j < No; j++)                                 \
    for (size_t i = 0; i < No; i++) {                               \
      Tijk[_IJK_(i, j, k)] += _t_buffer[_IJK_(__II, __JJ, __KK)];   \
    }                                                               \
    t_reorder.stop();
  #define DGEMM_PARTICLES(__A, __B)      \
    atrip::xgemm<F>( "T"                 \
                   , "N"                 \
                   , (int const*)&NoNo   \
                   , (int const*)&No     \
                   , (int const*)&Nv     \
                   , &one                \
                   , __A                 \
                   , (int const*)&Nv     \
                   , __B                 \
                   , (int const*)&Nv     \
                   , &zero               \
                   , _t_buffer.data()    \
                   , (int const*)&NoNo   \
                   );
  #define DGEMM_HOLES(__A, __B, __TRANSB)    \
    atrip::xgemm<F>( "N"                     \
                   , __TRANSB                \
                   , (int const*)&NoNo       \
                   , (int const*)&No         \
                   , (int const*)&No         \
                   , &m_one                  \
                   , __A                     \
                   , (int const*)&NoNo       \
                   , __B                     \
                   , (int const*)&No         \
                   , &zero                   \
                   , _t_buffer.data()        \
                   , (int const*)&NoNo       \
                   );
  #define MAYBE_CONJ(_conj, _buffer)                 \
    for (size_t __i = 0; __i < NoNoNo; ++__i)        \
      _conj[__i] = maybeConjugate<F>(_buffer[__i]);  \

    const size_t NoNoNo = No*NoNo;
    std::vector<F> _t_buffer;
    _t_buffer.reserve(NoNoNo);
    F one{1.0}, m_one{-1.0}, zero{0.0};

    t_reorder.start();
    for (size_t k = 0; k < NoNoNo; k++) {
      // zero the Tijk
      Tijk[k] = 0.0;
    }
    t_reorder.stop();

    chrono["doubles:holes"].start();
    { // Holes part ============================================================

      std::vector<F> _vhhh(NoNoNo);

      // VhhhC[i + k*No + L*NoNo] * TABhh[L + j*No]; H1
      MAYBE_CONJ(_vhhh, VhhhC)
      chrono["doubles:holes:1"].start();
      DGEMM_HOLES(_vhhh.data(), TABhh, "N")
      REORDER(i, k, j)
      chrono["doubles:holes:1"].stop();
      // VhhhC[j + k*No + L*NoNo] * TABhh[i + L*No]; H0
      chrono["doubles:holes:2"].start();
      DGEMM_HOLES(_vhhh.data(), TABhh, "T")
      REORDER(j, k, i)
      chrono["doubles:holes:2"].stop();

      // VhhhB[i + j*No + L*NoNo] * TAChh[L + k*No]; H5
      MAYBE_CONJ(_vhhh, VhhhB)
      chrono["doubles:holes:3"].start();
      DGEMM_HOLES(_vhhh.data(), TAChh, "N")
      REORDER(i, j, k)
      chrono["doubles:holes:3"].stop();
      // VhhhB[k + j*No + L*NoNo] * TAChh[i + L*No]; H3
      chrono["doubles:holes:4"].start();
      DGEMM_HOLES(_vhhh.data(), TAChh, "T")
      REORDER(k, j, i)
      chrono["doubles:holes:4"].stop();

      // VhhhA[j + i*No + L*NoNo] * TBChh[L + k*No]; H1
      MAYBE_CONJ(_vhhh, VhhhA)
      chrono["doubles:holes:5"].start();
      DGEMM_HOLES(_vhhh.data(), TBChh, "N")
      REORDER(j, i, k)
      chrono["doubles:holes:5"].stop();
      // VhhhA[k + i*No + L*NoNo] * TBChh[j + L*No]; H4
      chrono["doubles:holes:6"].start();
      DGEMM_HOLES(_vhhh.data(), TBChh, "T")
      REORDER(k, i, j)
      chrono["doubles:holes:6"].stop();

    }
    chrono["doubles:holes"].stop();
  #undef MAYBE_CONJ

    chrono["doubles:particles"].start();
    { // Particle part =========================================================
      // TAphh[E + i*Nv + j*NoNv] * VBCph[E + k*Nv]; P0
      chrono["doubles:particles:1"].start();
      DGEMM_PARTICLES(TAphh, VBCph)
      REORDER(i, j, k)
      chrono["doubles:particles:1"].stop();
      // TAphh[E + i*Nv + k*NoNv] * VCBph[E + j*Nv]; P3
      chrono["doubles:particles:2"].start();
      DGEMM_PARTICLES(TAphh, VCBph)
      REORDER(i, k, j)
      chrono["doubles:particles:2"].stop();
      // TCphh[E + k*Nv + i*NoNv] * VABph[E + j*Nv]; P5
      chrono["doubles:particles:3"].start();
      DGEMM_PARTICLES(TCphh, VABph)
      REORDER(k, i, j)
      chrono["doubles:particles:3"].stop();
      // TCphh[E + k*Nv + j*NoNv] * VBAph[E + i*Nv]; P2
      chrono["doubles:particles:4"].start();
      DGEMM_PARTICLES(TCphh, VBAph)
      REORDER(k, j, i)
      chrono["doubles:particles:4"].stop();
      // TBphh[E + j*Nv + i*NoNv] * VACph[E + k*Nv]; P1
      chrono["doubles:particles:5"].start();
      DGEMM_PARTICLES(TBphh, VACph)
      REORDER(j, i, k)
      chrono["doubles:particles:5"].stop();
      // TBphh[E + j*Nv + k*NoNv] * VCAph[E + i*Nv]; P4
      chrono["doubles:particles:6"].start();
      DGEMM_PARTICLES(TBphh, VCAph)
      REORDER(j, k, i)
      chrono["doubles:particles:6"].stop();
    }
    chrono["doubles:particles"].stop();

  #undef REORDER
  #undef DGEMM_HOLES
  #undef DGEMM_PARTICLES
  #undef _IJK_
  #else
    for (size_t k = 0; k < No; k++)
    for (size_t j = 0; j < No; j++)
    for (size_t i = 0; i < No; i++){
      const size_t ijk = i + j*No + k*NoNo
                ,  jk = j + k*No
                ;
      Tijk[ijk] = 0.0; // :important
      // HOLE DIAGRAMS: TABHH and VHHHA
      for (size_t L = 0; L < No; L++){
        // t[abLj] * V[Lcik]        H1
        // t[baLi] * V[Lcjk]        H0      TODO: conjugate T for complex
        Tijk[ijk] -= TABhh[L + j*No] * VhhhC[i + k*No + L*NoNo];
        Tijk[ijk] -= TABhh[i + L*No] * VhhhC[j + k*No + L*NoNo];

        // t[acLk] * V[Lbij]        H5
        // t[caLi] * V[Lbkj]        H3
        Tijk[ijk] -= TAChh[L + k*No] * VhhhB[i + j*No + L*NoNo];
        Tijk[ijk] -= TAChh[i + L*No] * VhhhB[k + j*No + L*NoNo];

        // t[bcLk] * V[Laji]        H2
        // t[cbLj] * V[Laki]        H4
        Tijk[ijk] -= TBChh[L + k*No] * VhhhA[j + i*No + L*NoNo];
        Tijk[ijk] -= TBChh[j + L*No] * VhhhA[k + i*No + L*NoNo];
      }
      // PARTILCE DIAGRAMS: TAPHH and VABPH
      for (size_t E = 0; E < Nv; E++) {
        // t[aEij] * V[bcEk]        P0
        // t[aEik] * V[cbEj]        P3 // TODO: CHECK THIS ONE, I DONT KNOW
        Tijk[ijk] += TAphh[E + i*Nv + j*NoNv] * VBCph[E + k*Nv];
        Tijk[ijk] += TAphh[E + i*Nv + k*NoNv] * VCBph[E + j*Nv];

        // t[cEki] * V[abEj]        P5
        // t[cEkj] * V[baEi]        P2
        Tijk[ijk] += TCphh[E + k*Nv + i*NoNv] * VABph[E + j*Nv];
        Tijk[ijk] += TCphh[E + k*Nv + j*NoNv] * VBAph[E + i*Nv];

        // t[bEji] * V[acEk]        P1
        // t[bEjk] * V[caEi]        P4
        Tijk[ijk] += TBphh[E + j*Nv + i*NoNv] * VACph[E + k*Nv];
        Tijk[ijk] += TBphh[E + j*Nv + k*NoNv] * VCAph[E + i*Nv];
      }

    }
  #endif
  }

}
#+end_src

** Blas
The main matrix-matrix multiplication method used in this algorithm
is mainly using the =DGEMM= function, which we declare as
=extern= since it should be known only at link-time.
#+begin_src c++ :tangle (atrip-blas-h)
#pragma once
namespace atrip {

  using Complex = std::complex<double>;

  extern "C" {
    void dgemm_(
      const char *transa,
      const char *transb,
      const int *m,
      const int *n,
      const int *k,
      double *alpha,
      const double *a,
      const int *lda,
      const double *b,
      const int *ldb,
      double *beta,
      double *c,
      const int *ldc
    );

    void zgemm_(
      const char *transa,
      const char *transb,
      const int *m,
      const int *n,
      const int *k,
      Complex *alpha,
      const Complex *A,
      const int *lda,
      const Complex *B,
      const int *ldb,
      Complex *beta,
      Complex *C,
      const int *ldc
    );
  }


  template <typename F=double>
  void xgemm(const char *transa,
             const char *transb,
             const int *m,
             const int *n,
             const int *k,
             F *alpha,
             const F *A,
             const int *lda,
             const F *B,
             const int *ldb,
             F *beta,
             F *C,
             const int *ldc) {
    dgemm_(transa, transb,
           m, n, k,
           alpha, A, lda,
           B, ldb, beta,
           C, ldc);
  }

  template <>
  void xgemm(const char *transa,
             const char *transb,
             const int *m,
             const int *n,
             const int *k,
             Complex *alpha,
             const Complex *A,
             const int *lda,
             const Complex *B,
             const int *ldb,
             Complex *beta,
             Complex *C,
             const int *ldc) {
    zgemm_(transa, transb,
           m, n, k,
           alpha, A, lda,
           B, ldb, beta,
           C, ldc);
  }
}
#+end_src

** Atrip
*** Header
#+begin_src c++ :tangle (atrip-atrip-h)
#pragma once
#include <sstream>
#include <string>
#include <map>
#include <chrono>

#include <ctf.hpp>

namespace atrip {

  struct Atrip {

    static int rank;
    static int np;
    static void init();

    template <typename F=double>
    struct Input {
      CTF::Tensor<F> *ei = nullptr
                        , *ea = nullptr
                        , *Tph = nullptr
                        , *Tpphh = nullptr
                        , *Vpphh = nullptr
                        , *Vhhhp = nullptr
                        , *Vppph = nullptr
                        ;
      int maxIterations = 0, iterationMod = -1, percentageMod = -1;
      bool barrier = false;
      bool chrono = false;
      Input& with_epsilon_i(CTF::Tensor<F> * t) { ei = t; return *this; }
      Input& with_epsilon_a(CTF::Tensor<F> * t) { ea = t; return *this; }
      Input& with_Tai(CTF::Tensor<F> * t) { Tph = t; return *this; }
      Input& with_Tabij(CTF::Tensor<F> * t) { Tpphh = t; return *this; }
      Input& with_Vabij(CTF::Tensor<F> * t) { Vpphh = t; return *this; }
      Input& with_Vijka(CTF::Tensor<F> * t) { Vhhhp = t; return *this; }
      Input& with_Vabci(CTF::Tensor<F> * t) { Vppph = t; return *this; }
      Input& with_maxIterations(int i) { maxIterations = i; return *this; }
      Input& with_iterationMod(int i) { iterationMod = i; return *this; }
      Input& with_percentageMod(int i) { percentageMod = i; return *this; }
      Input& with_barrier(bool i) { barrier = i; return *this; }
      Input& with_chrono(bool i) { chrono = i; return *this; }
    };

    struct Output {
      double energy;
    };
    template <typename F=double>
    static Output run(Input<F> const& in);
  };

}
#+end_src

*** Main

#+begin_src c++ :tangle (atrip-atrip-cxx)
#include <iomanip>

#include <atrip/Atrip.hpp>
#include <atrip/Utils.hpp>
#include <atrip/Equations.hpp>
#include <atrip/SliceUnion.hpp>
#include <atrip/Unions.hpp>

using namespace atrip;

int Atrip::rank;
int Atrip::np;

// user printing block
IterationDescriptor IterationDescription::descriptor;
void atrip::registerIterationDescriptor(IterationDescriptor d) {
  IterationDescription::descriptor = d;
}

void Atrip::init()  {
  MPI_Comm_rank(MPI_COMM_WORLD, &Atrip::rank);
  MPI_Comm_size(MPI_COMM_WORLD, &Atrip::np);
}

template <typename F>
Atrip::Output Atrip::run(Atrip::Input<F> const& in) {

  const int np = Atrip::np;
  const int rank = Atrip::rank;
  MPI_Comm universe = in.ei->wrld->comm;

  // Timings in seconds ================================================{{{1
  Timings chrono{};

  const size_t No = in.ei->lens[0];
  const size_t Nv = in.ea->lens[0];
  LOG(0,"Atrip") << "No: " << No << "\n";
  LOG(0,"Atrip") << "Nv: " << Nv << "\n";

  // allocate the three scratches, see piecuch
  std::vector<F>   Tijk(No*No*No) // doubles only (see piecuch)
                 , Zijk(No*No*No) // singles + doubles (see piecuch)
                 // we need local copies of the following tensors on every
                 // rank
                 , epsi(No)
                 , epsa(Nv)
                 , Tai(No * Nv)
                 ;

  in.ei->read_all(epsi.data());
  in.ea->read_all(epsa.data());
  in.Tph->read_all(Tai.data());

  // COMMUNICATOR CONSTRUCTION ========================================={{{1
  //
  // Construct a new communicator living only on a single rank
  int child_size = 1
    , child_rank
    ;
  const
  int color = rank / child_size
    , crank = rank % child_size
    ;
  MPI_Comm child_comm;
  if (np == 1) {
    child_comm = universe;
  } else {
    MPI_Comm_split(universe, color, crank, &child_comm);
    MPI_Comm_rank(child_comm, &child_rank);
    MPI_Comm_size(child_comm, &child_size);
  }


  chrono["nv-slices"].start();
  // BUILD SLICES PARAMETRIZED BY NV ==================================={{{1
  LOG(0,"Atrip") << "BUILD NV-SLICES\n";
  TAPHH<F> taphh(*in.Tpphh, (size_t)No, (size_t)Nv, (size_t)np, child_comm, universe);
  HHHA<F>  hhha(*in.Vhhhp, (size_t)No, (size_t)Nv, (size_t)np, child_comm, universe);
  chrono["nv-slices"].stop();

  chrono["nv-nv-slices"].start();
  // BUILD SLICES PARAMETRIZED BY NV x NV =============================={{{1
  LOG(0,"Atrip") << "BUILD NV x NV-SLICES\n";
  ABPH<F> abph(*in.Vppph, (size_t)No, (size_t)Nv, (size_t)np, child_comm, universe);
  ABHH<F> abhh(*in.Vpphh, (size_t)No, (size_t)Nv, (size_t)np, child_comm, universe);
  TABHH<F> tabhh(*in.Tpphh, (size_t)No, (size_t)Nv, (size_t)np, child_comm, universe);
  chrono["nv-nv-slices"].stop();

  // all tensors
  std::vector< SliceUnion<F>* > unions = {&taphh, &hhha, &abph, &abhh, &tabhh};

  //CONSTRUCT TUPLE LIST ==============================================={{{1
  LOG(0,"Atrip") << "BUILD TUPLE LIST\n";
  const auto tuplesList = std::move(getTuplesList(Nv));
  WITH_RANK << "tupList.size() = " << tuplesList.size() << "\n";

  // GET ABC INDEX RANGE FOR RANK ======================================{{{1
  auto abcIndex = getABCRange(np, rank, tuplesList);
  size_t nIterations = abcIndex.second - abcIndex.first;

  WITH_RANK << "abcIndex = " << pretty_print(abcIndex) << "\n";
  LOG(0,"Atrip") << "#iterations: " << nIterations << "\n";

  // first abc
  const ABCTuple firstAbc = tuplesList[abcIndex.first];


  double energy(0.);

  const size_t
      iterationMod = (in.percentageMod > 0)
                  ? nIterations * in.percentageMod / 100
                  : in.iterationMod

    , iteration1Percent = nIterations * 0.01
    ;



  auto const isFakeTuple
    = [&tuplesList](size_t const i) { return i >= tuplesList.size(); };


  using Database = typename Slice<F>::Database;
  using LocalDatabase = typename Slice<F>::LocalDatabase;
  auto communicateDatabase
    = [ &unions
      , np
      , &chrono
      ] (ABCTuple const& abc, MPI_Comm const& c) -> Database {

        chrono["db:comm:type:do"].start();
        auto MPI_LDB_ELEMENT = Slice<F>::mpi::localDatabaseElement();
        chrono["db:comm:type:do"].stop();

        chrono["db:comm:ldb"].start();
        LocalDatabase ldb;

        for (auto const& tensor: unions) {
          auto const& tensorDb = tensor->buildLocalDatabase(abc);
          ldb.insert(ldb.end(), tensorDb.begin(), tensorDb.end());
        }
        chrono["db:comm:ldb"].stop();

        Database db(np * ldb.size(), ldb[0]);

        chrono["oneshot-db:comm:allgather"].start();
        chrono["db:comm:allgather"].start();
        MPI_Allgather( ldb.data()
                     , ldb.size()
                     , MPI_LDB_ELEMENT
                     , db.data()
                     , ldb.size()
                     , MPI_LDB_ELEMENT
                     , c);
        chrono["db:comm:allgather"].stop();
        chrono["oneshot-db:comm:allgather"].stop();

        chrono["db:comm:type:free"].start();
        MPI_Type_free(&MPI_LDB_ELEMENT);
        chrono["db:comm:type:free"].stop();

        return db;
      };

  auto doIOPhase
    = [&unions, &rank, &np, &universe, &chrono] (Database const& db) {

    const size_t localDBLength = db.size() / np;

    size_t sendTag = 0
         , recvTag = rank * localDBLength
         ;

    // RECIEVE PHASE ======================================================
    {
      // At this point, we have already send to everyone that fits
      auto const& begin = &db[rank * localDBLength]
                , end   = begin + localDBLength
                ;
      for (auto it = begin; it != end; ++it) {
        recvTag++;
        auto const& el = *it;
        auto& u = unionByName(unions, el.name);

        WITH_DBG std::cout
          << rank << ":r"
          << "♯" << recvTag << " =>"
          << " «n" << el.name
          << ", t" << el.info.type
          << ", s" << el.info.state
          << "»"
          << " ⊙ {" << rank << "⇐" << el.info.from.rank
                    << ", "
                    << el.info.from.source << "}"
          << " ∴ {" << el.info.tuple[0]
                    << ", "
                    << el.info.tuple[1]
                    << "}"
          << "\n"
          ;

        chrono["db:io:recv"].start();
        u.receive(el.info, recvTag);
        chrono["db:io:recv"].stop();

      } // recv
    }

    // SEND PHASE =========================================================
    for (size_t otherRank = 0; otherRank<np; otherRank++) {
      auto const& begin = &db[otherRank * localDBLength]
                , end = begin + localDBLength
                ;
      for (auto it = begin; it != end; ++it) {
        sendTag++;
        typename Slice<F>::LocalDatabaseElement const& el = *it;

        if (el.info.from.rank != rank) continue;

        auto& u = unionByName(unions, el.name);
        WITH_DBG std::cout
          << rank << ":s"
          << "♯" << sendTag << " =>"
          << " «n" << el.name
          << ", t" << el.info.type
          << ", s" << el.info.state
          << "»"
          << " ⊙ {" << el.info.from.rank << "⇒" << otherRank
                    << ", "
                    << el.info.from.source << "}"
          << " ∴ {" << el.info.tuple[0]
                    << ", "
                    << el.info.tuple[1]
                    << "}"
          << "\n"
          ;

        chrono["db:io:send"].start();
        u.send(otherRank, el.info, sendTag);
        chrono["db:io:send"].stop();

      } // send phase

    } // otherRank


  };

#if defined(HAVE_OCD) || defined(ATRIP_PRINT_TUPLES)
  std::map<ABCTuple, double> tupleEnergies;
#endif

  const double doublesFlops
    = double(No)
    * double(No)
    * double(No)
    * (double(No) + double(Nv))
    * 1
    * (traits::isComplex<F>() ? 2.0 : 1.0)
    * 6
    / 1e9
    ;

  // START MAIN LOOP ======================================================{{{1

  for ( size_t i = abcIndex.first, iteration = 1
      ; i < abcIndex.second
      ; i++, iteration++
      ) {
    chrono["iterations"].start();


    // check overhead from chrono over all iterations
    chrono["start:stop"].start(); chrono["start:stop"].stop();

    // check overhead of doing a barrier at the beginning
    chrono["oneshot-mpi:barrier"].start();
    chrono["mpi:barrier"].start();
    // TODO: REMOVE
    if (in.barrier == 1)
    MPI_Barrier(universe);
    chrono["mpi:barrier"].stop();
    chrono["oneshot-mpi:barrier"].stop();

    if (iteration % iterationMod == 0 || iteration == iteration1Percent) {

      if (IterationDescription::descriptor) {
        IterationDescription::descriptor({
          iteration,
          nIterations,
          chrono["iterations"].count()
        });
      }

      LOG(0,"Atrip")
        << "iteration " << iteration
        << " [" << 100 * iteration / nIterations << "%]"
        << " (" << doublesFlops * iteration / chrono["doubles"].count()
        << "GF)"
        << " (" << doublesFlops * iteration / chrono["iterations"].count()
        << "GF)"
        << " ===========================\n";

      // PRINT TIMINGS
      if (in.chrono)
      for (auto const& pair: chrono)
        LOG(1, " ") << pair.first << " :: "
                    << pair.second.count()
                    << std::endl;

    }

    const ABCTuple abc = isFakeTuple(i)
                       ? tuplesList[tuplesList.size() - 1]
                       : tuplesList[i]
                 , *abcNext = i == (abcIndex.second - 1)
                            ? nullptr
                            : isFakeTuple(i + 1)
                            ? &tuplesList[tuplesList.size() - 1]
                            : &tuplesList[i + 1]
                 ;

    chrono["with_rank"].start();
    WITH_RANK << " :it " << iteration
              << " :abc " << pretty_print(abc)
              << " :abcN "
              << (abcNext ? pretty_print(*abcNext) : "None")
              << "\n";
    chrono["with_rank"].stop();


    // COMM FIRST DATABASE ================================================{{{1
    if (i == abcIndex.first) {
      WITH_RANK << "__first__:first database ............ \n";
      const auto __db = communicateDatabase(abc, universe);
      WITH_RANK << "__first__:first database communicated \n";
      WITH_RANK << "__first__:first database io phase \n";
      doIOPhase(__db);
      WITH_RANK << "__first__:first database io phase DONE\n";
      WITH_RANK << "__first__::::Unwrapping all slices for first database\n";
      for (auto& u: unions) u->unwrapAll(abc);
      WITH_RANK << "__first__::::Unwrapping all slices for first database DONE\n";
      MPI_Barrier(universe);
    }

    // COMM NEXT DATABASE ================================================={{{1
    if (abcNext) {
      WITH_RANK << "__comm__:" << iteration << "th communicating database\n";
      chrono["db:comm"].start();
      //const auto db = communicateDatabase(*abcNext, universe);
      Database db = communicateDatabase(*abcNext, universe);
      chrono["db:comm"].stop();
      chrono["db:io"].start();
      doIOPhase(db);
      chrono["db:io"].stop();
      WITH_RANK << "__comm__:" <<  iteration << "th database io phase DONE\n";
    }

    // COMPUTE DOUBLES ===================================================={{{1
    OCD_Barrier(universe);
    if (!isFakeTuple(i)) {
      WITH_RANK << iteration << "-th doubles\n";
      WITH_CHRONO(chrono["oneshot-unwrap"],
      WITH_CHRONO(chrono["unwrap"],
      WITH_CHRONO(chrono["unwrap:doubles"],
        for (auto& u: decltype(unions){&abph, &hhha, &taphh, &tabhh}) {
          u->unwrapAll(abc);
        }
      )))
      chrono["oneshot-doubles"].start();
      chrono["doubles"].start();
      doublesContribution<F>( abc, (size_t)No, (size_t)Nv
                            // -- VABCI
                            , abph.unwrapSlice(Slice<F>::AB, abc)
                            , abph.unwrapSlice(Slice<F>::AC, abc)
                            , abph.unwrapSlice(Slice<F>::BC, abc)
                            , abph.unwrapSlice(Slice<F>::BA, abc)
                            , abph.unwrapSlice(Slice<F>::CA, abc)
                            , abph.unwrapSlice(Slice<F>::CB, abc)
                            // -- VHHHA
                            , hhha.unwrapSlice(Slice<F>::A, abc)
                            , hhha.unwrapSlice(Slice<F>::B, abc)
                            , hhha.unwrapSlice(Slice<F>::C, abc)
                            // -- TA
                            , taphh.unwrapSlice(Slice<F>::A, abc)
                            , taphh.unwrapSlice(Slice<F>::B, abc)
                            , taphh.unwrapSlice(Slice<F>::C, abc)
                            // -- TABIJ
                            , tabhh.unwrapSlice(Slice<F>::AB, abc)
                            , tabhh.unwrapSlice(Slice<F>::AC, abc)
                            , tabhh.unwrapSlice(Slice<F>::BC, abc)
                            // -- TIJK
                            , Tijk.data()
                            , chrono
                            );
      WITH_RANK << iteration << "-th doubles done\n";
      chrono["doubles"].stop();
      chrono["oneshot-doubles"].stop();
    }

    // COMPUTE SINGLES =================================================== {{{1
    OCD_Barrier(universe);
    if (!isFakeTuple(i)) {
      WITH_CHRONO(chrono["oneshot-unwrap"],
      WITH_CHRONO(chrono["unwrap"],
      WITH_CHRONO(chrono["unwrap:singles"],
        abhh.unwrapAll(abc);
      )))
      chrono["reorder"].start();
      for (size_t I(0); I < Zijk.size(); I++) Zijk[I] = Tijk[I];
      chrono["reorder"].stop();
      chrono["singles"].start();
      singlesContribution<F>( No, Nv, abc
                            , Tai.data()
                            , abhh.unwrapSlice(Slice<F>::AB, abc)
                            , abhh.unwrapSlice(Slice<F>::AC, abc)
                            , abhh.unwrapSlice(Slice<F>::BC, abc)
                            , Zijk.data());
      chrono["singles"].stop();
    }


    // COMPUTE ENERGY ==================================================== {{{1
    if (!isFakeTuple(i)) {
      double tupleEnergy(0.);

      int distinct(0);
      if (abc[0] == abc[1]) distinct++;
      if (abc[1] == abc[2]) distinct--;
      const F epsabc(epsa[abc[0]] + epsa[abc[1]] + epsa[abc[2]]);

      chrono["energy"].start();
      if ( distinct == 0)
        tupleEnergy = getEnergyDistinct<F>(epsabc, epsi, Tijk, Zijk);
      else
        tupleEnergy = getEnergySame<F>(epsabc, epsi, Tijk, Zijk);
      chrono["energy"].stop();

#if defined(HAVE_OCD) || defined(ATRIP_PRINT_TUPLES)
      tupleEnergies[abc] = tupleEnergy;
#endif

      energy += tupleEnergy;

    }

    if (isFakeTuple(i)) {
      // fake iterations should also unwrap whatever they got
      WITH_RANK << iteration
                << "th unwrapping because of fake in "
                << i << "\n";
      for (auto& u: unions) u->unwrapAll(abc);
    }

#ifdef HAVE_OCD
    for (auto const& u: unions) {
      WITH_RANK << "__dups__:"
                << iteration
                << "-th n" << u->name << " checking duplicates\n";
      u->checkForDuplicates();
    }
#endif


    // CLEANUP UNIONS ===================================================={{{1
    OCD_Barrier(universe);
    if (abcNext) {
      chrono["gc"].start();
      WITH_RANK << "__gc__:" << iteration << "-th cleaning up.......\n";
      for (auto& u: unions) {

        u->unwrapAll(abc);
        WITH_RANK << "__gc__:n" << u->name  << " :it " << iteration
                  << " :abc " << pretty_print(abc)
                  << " :abcN " << pretty_print(*abcNext)
                  << "\n";
        // for (auto const& slice: u->slices)
        //   WITH_RANK << "__gc__:guts:" << slice.info << "\n";
        u->clearUnusedSlicesForNext(*abcNext);

        WITH_RANK << "__gc__: checking validity\n";

#ifdef HAVE_OCD
        // check for validity of the slices
        for (auto type: u->sliceTypes) {
          auto tuple = Slice<F>::subtupleBySlice(abc, type);
        for (auto& slice: u->slices) {
          if ( slice.info.type == type
             && slice.info.tuple == tuple
             && slice.isDirectlyFetchable()
             ) {
            if (slice.info.state == Slice<F>::Dispatched)
              throw std::domain_error( "This slice should not be undispatched! "
                                     + pretty_print(slice.info));
          }
        }
        }
#endif


      }
      chrono["gc"].stop();
    }

      WITH_RANK << iteration << "-th cleaning up....... DONE\n";

    chrono["iterations"].stop();
    // ITERATION END ====================================================={{{1

  }
    // END OF MAIN LOOP

  MPI_Barrier(universe);

  // PRINT TUPLES ========================================================={{{1
#if defined(HAVE_OCD) || defined(ATRIP_PRINT_TUPLES)
  LOG(0,"Atrip") << "tuple energies" << "\n";
  for (size_t i = 0; i < np; i++) {
    MPI_Barrier(universe);
    for (auto const& pair: tupleEnergies) {
      if (i == rank)
        std::cout << pair.first[0]
                  << " " << pair.first[1]
                  << " " << pair.first[2]
                  << std::setprecision(15) << std::setw(23)
                  << " tupleEnergy: " << pair.second
                  << "\n"
                  ;
    }
  }
#endif

  // COMMUNICATE THE ENERGIES ============================================={{{1
  LOG(0,"Atrip") << "COMMUNICATING ENERGIES \n";
  double globalEnergy = 0;
  MPI_Reduce(&energy, &globalEnergy, 1, MPI_DOUBLE, MPI_SUM, 0, universe);

  WITH_RANK << "local energy " << energy << "\n";
  LOG(0, "Atrip") << "Energy: "
    << std::setprecision(15) << std::setw(23)
    << (- globalEnergy) << std::endl;

  // PRINT TIMINGS {{{1
  if (in.chrono)
  for (auto const& pair: chrono)
    LOG(0,"atrip:chrono") << pair.first << " "
                          << pair.second.count() << std::endl;


  LOG(0, "atrip:flops(doubles)")
    << nIterations * doublesFlops / chrono["doubles"].count() << "\n";
  LOG(0, "atrip:flops(iterations)")
    << nIterations * doublesFlops / chrono["iterations"].count() << "\n";

  // TODO: change the sign in  the getEnergy routines
  return { - globalEnergy };

}
// instantiate
template Atrip::Output Atrip::run(Atrip::Input<double> const& in);
template Atrip::Output Atrip::run(Atrip::Input<Complex> const& in);

#+end_src


** Debug and Logging
*** Macros

#+begin_src c++ :tangle (atrip-debug-h)
#pragma once
#include <functional>
#define ATRIP_BENCHMARK
//#define ATRIP_DONT_SLICE
#ifndef ATRIP_DEBUG
#  define ATRIP_DEBUG 1
#endif
//#define ATRIP_WORKLOAD_DUMP
#define ATRIP_USE_DGEMM
//#define ATRIP_PRINT_TUPLES

#ifndef ATRIP_DEBUG
#define ATRIP_DEBUG 1
#endif

#if ATRIP_DEBUG == 4
#  pragma message("WARNING: You have OCD debugging ABC triples "    \
                  "expect GB of output and consult your therapist")
#  include <dbg.h>
#  define HAVE_OCD
#  define OCD_Barrier(com) MPI_Barrier(com)
#  define WITH_OCD
#  define WITH_ROOT if (atrip::Atrip::rank == 0)
#  define WITH_SPECIAL(r) if (atrip::Atrip::rank == r)
#  define WITH_RANK std::cout << atrip::Atrip::rank << ": "
#  define WITH_CRAZY_DEBUG
#  define WITH_DBG
#  define DBG(...) dbg(__VA_ARGS__)
#elif ATRIP_DEBUG == 3
#  pragma message("WARNING: You have crazy debugging ABC triples,"  \
                  " expect GB of output")
#  include <dbg.h>
#  define OCD_Barrier(com)
#  define WITH_OCD if (false)
#  define WITH_ROOT if (atrip::Atrip::rank == 0)
#  define WITH_SPECIAL(r) if (atrip::Atrip::rank == r)
#  define WITH_RANK std::cout << atrip::Atrip::rank << ": "
#  define WITH_CRAZY_DEBUG
#  define WITH_DBG
#  define DBG(...) dbg(__VA_ARGS__)
#elif ATRIP_DEBUG == 2
#  pragma message("WARNING: You have some debugging info for ABC triples")
#  include <dbg.h>
#  define OCD_Barrier(com)
#  define WITH_OCD if (false)
#  define WITH_ROOT if (atrip::Atrip::rank == 0)
#  define WITH_SPECIAL(r) if (atrip::Atrip::rank == r)
#  define WITH_RANK std::cout << atrip::Atrip::rank << ": "
#  define WITH_CRAZY_DEBUG if (false)
#  define WITH_DBG
#  define DBG(...) dbg(__VA_ARGS__)
#else
#  define OCD_Barrier(com)
#  define WITH_OCD if (false)
#  define WITH_ROOT if (false)
#  define WITH_SPECIAL(r) if (false)
#  define WITH_RANK if (false) std::cout << atrip::Atrip::rank << ": "
#  define WITH_DBG if (false)
#  define WITH_CRAZY_DEBUG if (false)
#  define DBG(...)
#endif
#+end_src

And users of the library can redefine the =LOG= macro
which in case of not being defined is defined as follows:

#+begin_src c++ :tangle (atrip-debug-h)
#ifndef LOG
#define LOG(level, name) if (Atrip::rank == 0) std::cout << name << ": "
#endif
#+end_src

Furthermore, if you do not wish to see any output from ATRIP, simply
define =ATRIP_NO_OUTPUT=


#+begin_src c++ :tangle (atrip-debug-h)
#ifdef ATRIP_NO_OUTPUT
#  undef LOG
#  define LOG(level, name) if (false) std::cout << name << ": "
#endif
#+end_src

*** Iteration informer

In general a code writer will want to write some messages in every iteration.
A developer then can register a function to be used in this sense.
The input of the function is an [[IterationDescriptor]] structure and the output
should be nothing.

#+name: IterationDescriptor
#+begin_src c++ :tangle (atrip-debug-h)
namespace atrip {

  struct IterationDescription;
  using IterationDescriptor = std::function<void(IterationDescription const&)>;
  struct IterationDescription {
    static IterationDescriptor descriptor;
    size_t currentIteration;
    size_t totalIterations;
    double currentElapsedTime;
  };

  void registerIterationDescriptor(IterationDescriptor);

}
#+end_src

** Include header

#+begin_src c++ :tangle (atrip-main-h)
#pragma once

#include <atrip/Atrip.hpp>

#+end_src