text: /**
   * \file
   * Compute the tensor transposition
   */
  
  #include <tuple>
  #include <string>
  #include <vector>
  #include <list>
  #include <algorithm>
  #include <numeric>
  #include <iostream>
  #include <cmath>
  #include <chrono>
  
  #include <float.h>
  #include <stdio.h>
  #include <assert.h>
  
  #ifdef _OPENMP
  #include <omp.h>
  #endif
  
  #include "../include/utils.h"
  #include "../include/macros.h"
  #include "../include/compute_node.h"
  #include "../include/utils.h"
  #include "../include/hptt_types.h"
  #include "../include/plan.h"
  #include "../include/transpose.h"
  
  namespace hptt {
  
  template <typename floatType, int betaIsZero, bool conjA>
  struct micro_kernel
  {
      static void execute(const floatType* __restrict__ A, const size_t lda, floatType* __restrict__ B, const size_t ldb, const floatType alpha, const floatType beta)
      {
         constexpr int n = (REGISTER_BITS/8) / sizeof(floatType);
  
         if( betaIsZero )
            for(int j=0; j < n; ++j)
               for(int i=0; i < n; ++i)
                  if( conjA )
                     B[i + j * ldb] = alpha * conj(A[i * lda + j]);
                  else
                     B[i + j * ldb] = alpha * A[i * lda + j];
         else
            for(int j=0; j < n; ++j)
               for(int i=0; i < n; ++i)
                  if( conjA )
                     B[i + j * ldb] = alpha * conj(A[i * lda + j]) + beta * B[i + j * ldb];
                  else
                     B[i + j * ldb] = alpha * A[i * lda + j] + beta * B[i + j * ldb];
      }
  };
  
  template<typename floatType>
  static void streamingStore( floatType* out, const floatType *in )
  {
     constexpr int n = REGISTER_BITS/8/sizeof(floatType);
     for(int i=0; i < n; ++i)
        out[i] = in[i];
  }
  
  
  #ifdef HPTT_ARCH_AVX
  #include <immintrin.h>
  
  template<typename floatType>
  static INLINE void prefetch(const floatType* A, const int lda)
  {
     constexpr int blocking_micro_ = REGISTER_BITS/8 / sizeof(floatType);
     for(int i=0;i < blocking_micro_; ++i)
        _mm_prefetch((char*)(A + i * lda), _MM_HINT_T2);
  }
  template <int betaIsZero, bool conjA>
  struct micro_kernel<double, betaIsZero, conjA>
  {
      static void execute(const double* __restrict__ A, const size_t lda, double* __restrict__ B, const size_t ldb, const double alpha ,const double beta)
      {
         __m256d reg_alpha = _mm256_set1_pd(alpha); // do not alter the content of B
         __m256d reg_beta = _mm256_set1_pd(beta); // do not alter the content of B
         //Load A
         __m256d rowA0 = _mm256_loadu_pd((A +0*lda));
         __m256d rowA1 = _mm256_loadu_pd((A +1*lda));
         __m256d rowA2 = _mm256_loadu_pd((A +2*lda));
         __m256d rowA3 = _mm256_loadu_pd((A +3*lda));
  
         //4x4 transpose micro kernel
         __m256d r4, r34, r3, r33;
         r33 = _mm256_shuffle_pd( rowA2, rowA3, 0x3 );
         r3 = _mm256_shuffle_pd( rowA0, rowA1, 0x3 );
         r34 = _mm256_shuffle_pd( rowA2, rowA3, 0xc );
         r4 = _mm256_shuffle_pd( rowA0, rowA1, 0xc );
         rowA0 = _mm256_permute2f128_pd( r34, r4, 0x2 );
         rowA1 = _mm256_permute2f128_pd( r33, r3, 0x2 );
         rowA2 = _mm256_permute2f128_pd( r33, r3, 0x13 );
         rowA3 = _mm256_permute2f128_pd( r34, r4, 0x13 );
  
         //Scale A
         rowA0 = _mm256_mul_pd(rowA0, reg_alpha);
         rowA1 = _mm256_mul_pd(rowA1, reg_alpha);
         rowA2 = _mm256_mul_pd(rowA2, reg_alpha);
         rowA3 = _mm256_mul_pd(rowA3, reg_alpha);
  
         //Load B
         if( !betaIsZero )
         {
            __m256d rowB0 = _mm256_loadu_pd((B +0*ldb));
            __m256d rowB1 = _mm256_loadu_pd((B +1*ldb));
            __m256d rowB2 = _mm256_loadu_pd((B +2*ldb));
            __m256d rowB3 = _mm256_loadu_pd((B +3*ldb));
  
            rowB0 = _mm256_add_pd( _mm256_mul_pd(rowB0, reg_beta), rowA0);
            rowB1 = _mm256_add_pd( _mm256_mul_pd(rowB1, reg_beta), rowA1);
            rowB2 = _mm256_add_pd( _mm256_mul_pd(rowB2, reg_beta), rowA2);
            rowB3 = _mm256_add_pd( _mm256_mul_pd(rowB3, reg_beta), rowA3);
            //Store B
            _mm256_storeu_pd((B + 0 * ldb), rowB0);
            _mm256_storeu_pd((B + 1 * ldb), rowB1);
            _mm256_storeu_pd((B + 2 * ldb), rowB2);
            _mm256_storeu_pd((B + 3 * ldb), rowB3);
         } else {
            //Store B
            _mm256_storeu_pd((B + 0 * ldb), rowA0);
            _mm256_storeu_pd((B + 1 * ldb), rowA1);
            _mm256_storeu_pd((B + 2 * ldb), rowA2);
            _mm256_storeu_pd((B + 3 * ldb), rowA3);
         }
      }
  };
  
  template <int betaIsZero, bool conjA>
  struct micro_kernel<float, betaIsZero, conjA>
  {
      static void execute(const float* __restrict__ A, const size_t lda, float* __restrict__ B, const size_t ldb, const float alpha ,const float beta)
      {
         __m256 reg_alpha = _mm256_set1_ps(alpha); // do not alter the content of B
         __m256 reg_beta = _mm256_set1_ps(beta); // do not alter the content of B
         //Load A
         __m256 rowA0 = _mm256_loadu_ps((A +0*lda));
         __m256 rowA1 = _mm256_loadu_ps((A +1*lda));
         __m256 rowA2 = _mm256_loadu_ps((A +2*lda));
         __m256 rowA3 = _mm256_loadu_ps((A +3*lda));
         __m256 rowA4 = _mm256_loadu_ps((A +4*lda));
         __m256 rowA5 = _mm256_loadu_ps((A +5*lda));
         __m256 rowA6 = _mm256_loadu_ps((A +6*lda));
         __m256 rowA7 = _mm256_loadu_ps((A +7*lda));
  
         //8x8 transpose micro kernel
         __m256 r121, r139, r120, r138, r71, r89, r70, r88, r11, r1, r55, r29, r10, r0, r54, r28;
         r28 = _mm256_unpacklo_ps( rowA4, rowA5 );
         r54 = _mm256_unpacklo_ps( rowA6, rowA7 );
         r0 = _mm256_unpacklo_ps( rowA0, rowA1 );
         r10 = _mm256_unpacklo_ps( rowA2, rowA3 );
         r29 = _mm256_unpackhi_ps( rowA4, rowA5 );
         r55 = _mm256_unpackhi_ps( rowA6, rowA7 );
         r1 = _mm256_unpackhi_ps( rowA0, rowA1 );
         r11 = _mm256_unpackhi_ps( rowA2, rowA3 );
         r88 = _mm256_shuffle_ps( r28, r54, 0x44 );
         r70 = _mm256_shuffle_ps( r0, r10, 0x44 );
         r89 = _mm256_shuffle_ps( r28, r54, 0xee );
         r71 = _mm256_shuffle_ps( r0, r10, 0xee );
         r138 = _mm256_shuffle_ps( r29, r55, 0x44 );
         r120 = _mm256_shuffle_ps( r1, r11, 0x44 );
         r139 = _mm256_shuffle_ps( r29, r55, 0xee );
         r121 = _mm256_shuffle_ps( r1, r11, 0xee );
         rowA0 = _mm256_permute2f128_ps( r88, r70, 0x2 );
         rowA1 = _mm256_permute2f128_ps( r89, r71, 0x2 );
         rowA2 = _mm256_permute2f128_ps( r138, r120, 0x2 );
         rowA3 = _mm256_permute2f128_ps( r139, r121, 0x2 );
         rowA4 = _mm256_permute2f128_ps( r88, r70, 0x13 );
         rowA5 = _mm256_permute2f128_ps( r89, r71, 0x13 );
         rowA6 = _mm256_permute2f128_ps( r138, r120, 0x13 );
         rowA7 = _mm256_permute2f128_ps( r139, r121, 0x13 );
  
         //Scale A
         rowA0 = _mm256_mul_ps(rowA0, reg_alpha);
         rowA1 = _mm256_mul_ps(rowA1, reg_alpha);
         rowA2 = _mm256_mul_ps(rowA2, reg_alpha);
         rowA3 = _mm256_mul_ps(rowA3, reg_alpha);
         rowA4 = _mm256_mul_ps(rowA4, reg_alpha);
         rowA5 = _mm256_mul_ps(rowA5, reg_alpha);
         rowA6 = _mm256_mul_ps(rowA6, reg_alpha);
         rowA7 = _mm256_mul_ps(rowA7, reg_alpha);
  
         //Load B
         if( !betaIsZero )
         {
            __m256 rowB0 = _mm256_loadu_ps((B +0*ldb));
            __m256 rowB1 = _mm256_loadu_ps((B +1*ldb));
            __m256 rowB2 = _mm256_loadu_ps((B +2*ldb));
            __m256 rowB3 = _mm256_loadu_ps((B +3*ldb));
            __m256 rowB4 = _mm256_loadu_ps((B +4*ldb));
            __m256 rowB5 = _mm256_loadu_ps((B +5*ldb));
            __m256 rowB6 = _mm256_loadu_ps((B +6*ldb));
            __m256 rowB7 = _mm256_loadu_ps((B +7*ldb));
  
            rowB0 = _mm256_add_ps( _mm256_mul_ps(rowB0, reg_beta), rowA0);
            rowB1 = _mm256_add_ps( _mm256_mul_ps(rowB1, reg_beta), rowA1);
            rowB2 = _mm256_add_ps( _mm256_mul_ps(rowB2, reg_beta), rowA2);
            rowB3 = _mm256_add_ps( _mm256_mul_ps(rowB3, reg_beta), rowA3);
            rowB4 = _mm256_add_ps( _mm256_mul_ps(rowB4, reg_beta), rowA4);
            rowB5 = _mm256_add_ps( _mm256_mul_ps(rowB5, reg_beta), rowA5);
            rowB6 = _mm256_add_ps( _mm256_mul_ps(rowB6, reg_beta), rowA6);
            rowB7 = _mm256_add_ps( _mm256_mul_ps(rowB7, reg_beta), rowA7);
            //Store B
            _mm256_storeu_ps((B + 0 * ldb), rowB0);
            _mm256_storeu_ps((B + 1 * ldb), rowB1);
            _mm256_storeu_ps((B + 2 * ldb), rowB2);
            _mm256_storeu_ps((B + 3 * ldb), rowB3);
            _mm256_storeu_ps((B + 4 * ldb), rowB4);
            _mm256_storeu_ps((B + 5 * ldb), rowB5);
            _mm256_storeu_ps((B + 6 * ldb), rowB6);
            _mm256_storeu_ps((B + 7 * ldb), rowB7);
         } else {
            _mm256_storeu_ps((B + 0 * ldb), rowA0);
            _mm256_storeu_ps((B + 1 * ldb), rowA1);
            _mm256_storeu_ps((B + 2 * ldb), rowA2);
            _mm256_storeu_ps((B + 3 * ldb), rowA3);
            _mm256_storeu_ps((B + 4 * ldb), rowA4);
            _mm256_storeu_ps((B + 5 * ldb), rowA5);
            _mm256_storeu_ps((B + 6 * ldb), rowA6);
            _mm256_storeu_ps((B + 7 * ldb), rowA7);
         }
      }
  };
  
  template<>
  void streamingStore<float>( float* out, const float*in ){
     _mm256_stream_ps(out, _mm256_loadu_ps(in));
  }
  template<>
  void streamingStore<double>( double* out, const double*in ){
     _mm256_stream_pd(out, _mm256_loadu_pd(in));
  }
  #else
  template<typename floatType>
  static INLINE void prefetch(const floatType* A, const int lda) { }
  #endif
  
  #ifdef HPTT_ARCH_ARM
  #include <arm_neon.h>
  
  template <int betaIsZero, bool conjA>
  struct micro_kernel<float, betaIsZero>
  {
      static void execute(const float* __restrict__ A, const size_t lda, float* __restrict__ B, const size_t ldb, const float alpha ,const float beta)
      {
         float32x4_t reg_alpha = vdupq_n_f32(alpha);
         float32x4_t reg_beta = vdupq_n_f32(beta);
  
         //Load A
         float32x4_t rowA0 = vld1q_f32((A +0*lda));
         float32x4_t rowA1 = vld1q_f32((A +1*lda));
         float32x4_t rowA2 = vld1q_f32((A +2*lda));
         float32x4_t rowA3 = vld1q_f32((A +3*lda));
  
         //4x4 transpose micro kernel
         float32x4x2_t t0,t1,t2,t3;
         t0 = vuzpq_f32(rowA0, rowA2);
         t1 = vuzpq_f32(rowA1, rowA3);
         t2 = vtrnq_f32(t0.val[0], t1.val[0]);
         t3 = vtrnq_f32(t0.val[1], t1.val[1]);
  
         //Scale A
         rowA0 = vmulq_f32(t2.val[0], reg_alpha);
         rowA1 = vmulq_f32(t3.val[0], reg_alpha);
         rowA2 = vmulq_f32(t2.val[1], reg_alpha);
         rowA3 = vmulq_f32(t3.val[1], reg_alpha);
  
         //Load B
         if( !betaIsZero )
         {
            float32x4_t rowB0 = vld1q_f32((B +0*ldb));
            float32x4_t rowB1 = vld1q_f32((B +1*ldb));
            float32x4_t rowB2 = vld1q_f32((B +2*ldb));
            float32x4_t rowB3 = vld1q_f32((B +3*ldb));
  
            rowB0 = vaddq_f32( vmulq_f32(rowB0, reg_beta), rowA0);
            rowB1 = vaddq_f32( vmulq_f32(rowB1, reg_beta), rowA1);
            rowB2 = vaddq_f32( vmulq_f32(rowB2, reg_beta), rowA2);
            rowB3 = vaddq_f32( vmulq_f32(rowB3, reg_beta), rowA3);
            //Store B
            vst1q_f32((B + 0 * ldb), rowB0);
            vst1q_f32((B + 1 * ldb), rowB1);
            vst1q_f32((B + 2 * ldb), rowB2);
            vst1q_f32((B + 3 * ldb), rowB3);
         } else {
            //Store B
            vst1q_f32((B + 0 * ldb), rowA0);
            vst1q_f32((B + 1 * ldb), rowA1);
            vst1q_f32((B + 2 * ldb), rowA2);
            vst1q_f32((B + 3 * ldb), rowA3);
         }
      }
  };
  #endif
  
  #ifdef HPTT_ARCH_IBM
  //#include <altivec.h> //vector conflicts with std::vector (TODO)
  //
  //template <int betaIsZero>
  //struct micro_kernel<float, betaIsZero>
  //{
  //    static void execute(const float* __restrict__ A, const size_t lda, float* __restrict__ B, const size_t ldb, const float alpha ,const float beta)
  //    {
  //       vector float reg_alpha = vec_splats(alpha);
  //
  //       //Load A
  //       vector float rowA0 = vec_ld(0,const_cast<float*>(A+0*lda));
  //       vector float rowA1 = vec_ld(0,const_cast<float*>(A+1*lda));
  //       vector float rowA2 = vec_ld(0,const_cast<float*>(A+2*lda));
  //       vector float rowA3 = vec_ld(0,const_cast<float*>(A+3*lda));
  //
  //       //4x4 transpose micro kernel
  //       vector float aa = (vector float) {2, 3, 2.5, 3.5};
  //       vector float bb = (vector float) {2.25, 3.25, 2.75, 3.75};
  //       vector float cc = (vector float) {2, 2.25, 3, 3.25};
  //       vector float dd = (vector float) {2.5, 2.75, 3.5, 3.75};
  //
  //       vector float r010 = vec_perm(rowA0,rowA1, aa); //0,4,2,6
  //       vector float r011 = vec_perm(rowA0,rowA1, bb); //1,5,3,7
  //       vector float r230 = vec_perm(rowA2,rowA3, aa); //8,12,10,14
  //       vector float r231 = vec_perm(rowA2,rowA3, bb); //9,13,11,15
  //
  //       rowA0 = vec_perm(r010, r230, cc); //0,4,8,12
  //       rowA1 = vec_perm(r011, r231, cc); //1,5,9,13
  //       rowA2 = vec_perm(r010, r230, dd); //2,6,10,14
  //       rowA3 = vec_perm(r011, r231, dd); //3,7,11,15
  //
  //       //Scale A
  //       rowA0 = vec_mul(rowA0, reg_alpha);
  //       rowA1 = vec_mul(rowA1, reg_alpha);
  //       rowA2 = vec_mul(rowA2, reg_alpha);
  //       rowA3 = vec_mul(rowA3, reg_alpha);
  //
  //       if( !betaIsZero )
  //       {
  //          vector float reg_beta = vec_splats(beta);
  //          //Load B
  //          vector float rowB0 = vec_ld(0,const_cast<float*>(B+0*ldb));
  //          vector float rowB1 = vec_ld(0,const_cast<float*>(B+1*ldb));
  //          vector float rowB2 = vec_ld(0,const_cast<float*>(B+2*ldb));
  //          vector float rowB3 = vec_ld(0,const_cast<float*>(B+3*ldb));
  //
  //          rowB0 = vec_madd( rowB0, reg_beta, rowA0);
  //          rowB1 = vec_madd( rowB1, reg_beta, rowA1);
  //          rowB2 = vec_madd( rowB2, reg_beta, rowA2);
  //          rowB3 = vec_madd( rowB3, reg_beta, rowA3);
  //
  //          //Store B
  //          vec_st(rowB0, 0, B + 0 * ldb);
  //          vec_st(rowB1, 0, B + 1 * ldb);
  //          vec_st(rowB2, 0, B + 2 * ldb);
  //          vec_st(rowB3, 0, B + 3 * ldb);
  //       } else {             
  //          //Store B         
  //          vec_st(rowA0, 0, B + 0 * ldb);
  //          vec_st(rowA1, 0, B + 1 * ldb);
  //          vec_st(rowA2, 0, B + 2 * ldb);
  //          vec_st(rowA3, 0, B + 3 * ldb);
  //       }
  //    }
  //};
  #endif
  
  
  template<int betaIsZero, typename floatType, bool conjA>
  static INLINE void macro_kernel_scalar(const floatType* __restrict__ A, const size_t lda, int blockingA,  
                                               floatType* __restrict__ B, const size_t ldb, int blockingB,
                                               const floatType alpha ,const floatType beta)
  {
  #ifdef DEBUG
     assert( blockingA > 0 && blockingB > 0);
  #endif
  
     if( betaIsZero )
        for(int j=0; j < blockingA; ++j)
           for(int i=0; i < blockingB; ++i)
              if( conjA )
                 B[i + j * ldb] = alpha * conj(A[i * lda + j]);
              else
                 B[i + j * ldb] = alpha * A[i * lda + j];
     else
        for(int j=0; j < blockingA; ++j)
           for(int i=0; i < blockingB; ++i)
              if( conjA )
                 B[i + j * ldb] = alpha * conj(A[i * lda + j]) + beta * B[i + j * ldb];
              else
                 B[i + j * ldb] = alpha * A[i * lda + j] + beta * B[i + j * ldb];
  }
  
  template<int blockingA, int blockingB, int betaIsZero, typename floatType, bool useStreamingStores_, bool conjA>
  static INLINE void macro_kernel(const floatType* __restrict__ A, const floatType* __restrict__ Anext, const size_t lda, 
                                     floatType* __restrict__ B, const floatType* __restrict__ Bnext, const size_t ldb,
                                     const floatType alpha ,const floatType beta)
  {
     constexpr int blocking_micro_ = REGISTER_BITS/8 / sizeof(floatType);
     constexpr int blocking_ = blocking_micro_ * 4;
  
     const bool useStreamingStores = useStreamingStores_ && betaIsZero && (blockingB*sizeof(floatType))%64 == 0 && ((uint64_t)B)%32 == 0 && (ldb*sizeof(floatType))%32 == 0;
  
     floatType *Btmp = B;
     size_t ldb_tmp = ldb;
     floatType buffer[blockingA * blockingB];// __attribute__((aligned(64)));
     if( (useStreamingStores_ && useStreamingStores) ){
        Btmp = buffer;
        ldb_tmp = blockingB;
     }
  
     if( blockingA == blocking_ && blockingB == blocking_ )
     {
        if( !(useStreamingStores_ && useStreamingStores) )
           prefetch<floatType>(Bnext + (0 * ldb_tmp + 0), ldb_tmp);
        prefetch<floatType>(Anext + (0 * lda + 0), lda);
        micro_kernel<floatType,betaIsZero, conjA>::execute(A + (0 * lda + 0), lda, Btmp + (0 * ldb_tmp + 0), ldb_tmp  , alpha , beta);
        prefetch<floatType>(Anext + (blocking_micro_ * lda + 0), lda);
        micro_kernel<floatType,betaIsZero, conjA>::execute(A + (blocking_micro_ * lda + 0), lda, Btmp + (0 * ldb_tmp + blocking_micro_), ldb_tmp  , alpha , beta);
        if( !(useStreamingStores_ && useStreamingStores) )
           prefetch<floatType>(Bnext + (0 * ldb_tmp + 2*blocking_micro_), ldb_tmp);
        prefetch<floatType>(Anext + (2*blocking_micro_ * lda + 0), lda);
        micro_kernel<floatType,betaIsZero, conjA>::execute(A + (2*blocking_micro_ * lda + 0), lda, Btmp + (0 * ldb_tmp + 2*blocking_micro_), ldb_tmp  , alpha , beta);
        prefetch<floatType>(Anext + (3*blocking_micro_ * lda + 0), lda);
        micro_kernel<floatType,betaIsZero, conjA>::execute(A + (3*blocking_micro_ * lda + 0), lda, Btmp + (0 * ldb_tmp + 3*blocking_micro_), ldb_tmp  , alpha , beta);
        if( !(useStreamingStores_ && useStreamingStores) )
           prefetch<floatType>(Bnext + (blocking_micro_ * ldb_tmp + 0), ldb_tmp);
        micro_kernel<floatType,betaIsZero, conjA>::execute(A + (0 * lda + blocking_micro_), lda, Btmp + (blocking_micro_ * ldb_tmp + 0), ldb_tmp  , alpha , beta);
        micro_kernel<floatType,betaIsZero, conjA>::execute(A + (blocking_micro_ * lda + blocking_micro_), lda, Btmp + (blocking_micro_ * ldb_tmp + blocking_micro_), ldb_tmp  , alpha , beta);
        if( !(useStreamingStores_ && useStreamingStores) )
           prefetch<floatType>(Bnext + (blocking_micro_ * ldb_tmp + 2*blocking_micro_), ldb_tmp);
        micro_kernel<floatType,betaIsZero, conjA>::execute(A + (2*blocking_micro_ * lda + blocking_micro_), lda, Btmp + (blocking_micro_ * ldb_tmp + 2*blocking_micro_), ldb_tmp  , alpha , beta);
        micro_kernel<floatType,betaIsZero, conjA>::execute(A + (3*blocking_micro_ * lda + blocking_micro_), lda, Btmp + (blocking_micro_ * ldb_tmp + 3*blocking_micro_), ldb_tmp  , alpha , beta);
        if( !(useStreamingStores_ && useStreamingStores) )
           prefetch<floatType>(Bnext + (2*blocking_micro_ * ldb_tmp + 0), ldb_tmp);
        prefetch<floatType>(Anext + (0 * lda + 2*blocking_micro_), lda);
        micro_kernel<floatType,betaIsZero, conjA>::execute(A + (0 * lda + 2*blocking_micro_), lda, Btmp + (2*blocking_micro_ * ldb_tmp + 0), ldb_tmp  , alpha , beta);
        prefetch<floatType>(Anext + (blocking_micro_ * lda + 2*blocking_micro_), lda);
        micro_kernel<floatType,betaIsZero, conjA>::execute(A + (blocking_micro_ * lda + 2*blocking_micro_), lda, Btmp + (2*blocking_micro_ * ldb_tmp + blocking_micro_), ldb_tmp  , alpha , beta);
        if( !(useStreamingStores_ && useStreamingStores) )
           prefetch<floatType>(Bnext + (2*blocking_micro_ * ldb_tmp + 2*blocking_micro_), ldb_tmp);
        prefetch<floatType>(Anext + (2*blocking_micro_ * lda + 2*blocking_micro_), lda);
        micro_kernel<floatType,betaIsZero, conjA>::execute(A + (2*blocking_micro_ * lda + 2*blocking_micro_), lda, Btmp + (2*blocking_micro_ * ldb_tmp + 2*blocking_micro_), ldb_tmp  , alpha , beta);
        prefetch<floatType>(Anext + (3*blocking_micro_ * lda + 2*blocking_micro_), lda);
        micro_kernel<floatType,betaIsZero, conjA>::execute(A + (3*blocking_micro_ * lda + 2*blocking_micro_), lda, Btmp + (2*blocking_micro_ * ldb_tmp + 3*blocking_micro_), ldb_tmp  , alpha , beta);
        if( !(useStreamingStores_ && useStreamingStores) )
           prefetch<floatType>(Bnext + (3*blocking_micro_ * ldb_tmp + 0), ldb_tmp);
        micro_kernel<floatType,betaIsZero, conjA>::execute(A + (0 * lda + 3*blocking_micro_), lda, Btmp + (3*blocking_micro_ * ldb_tmp + 0), ldb_tmp  , alpha , beta);
        micro_kernel<floatType,betaIsZero, conjA>::execute(A + (blocking_micro_ * lda + 3*blocking_micro_), lda, Btmp + (3*blocking_micro_ * ldb_tmp + blocking_micro_), ldb_tmp  , alpha , beta);
        if( !(useStreamingStores_ && useStreamingStores) )
           prefetch<floatType>(Bnext + (3*blocking_micro_ * ldb_tmp + 2*blocking_micro_), ldb_tmp);
        micro_kernel<floatType,betaIsZero, conjA>::execute(A + (2*blocking_micro_ * lda + 3*blocking_micro_), lda, Btmp + (3*blocking_micro_ * ldb_tmp + 2*blocking_micro_), ldb_tmp  , alpha , beta);
        micro_kernel<floatType,betaIsZero, conjA>::execute(A + (3*blocking_micro_ * lda + 3*blocking_micro_), lda, Btmp + (3*blocking_micro_ * ldb_tmp + 3*blocking_micro_), ldb_tmp  , alpha , beta);
     }else if( blockingA == 2*blocking_micro_ && blockingB == blocking_ ) {
        if( !(useStreamingStores_ && useStreamingStores) )
           prefetch<floatType>(Bnext + (0 * ldb_tmp + 0), ldb_tmp);
        prefetch<floatType>(Anext + (0 * lda + 0), lda);
        micro_kernel<floatType,betaIsZero, conjA>::execute(A + (0 * lda + 0), lda, Btmp + (0 * ldb_tmp + 0), ldb_tmp  , alpha , beta);
        prefetch<floatType>(Anext + (blocking_micro_ * lda + 0), lda);
        micro_kernel<floatType,betaIsZero, conjA>::execute(A + (blocking_micro_ * lda + 0), lda, Btmp + (0 * ldb_tmp + blocking_micro_), ldb_tmp  , alpha , beta);
        if( !(useStreamingStores_ && useStreamingStores) )
           prefetch<floatType>(Bnext + (0 * ldb_tmp + 2*blocking_micro_), ldb_tmp);
        prefetch<floatType>(Anext + (2*blocking_micro_ * lda + 0), lda);
        micro_kernel<floatType,betaIsZero, conjA>::execute(A + (2*blocking_micro_ * lda + 0), lda, Btmp + (0 * ldb_tmp + 2*blocking_micro_), ldb_tmp  , alpha , beta);
        prefetch<floatType>(Anext + (3*blocking_micro_ * lda + 0), lda);
        micro_kernel<floatType,betaIsZero, conjA>::execute(A + (3*blocking_micro_ * lda + 0), lda, Btmp + (0 * ldb_tmp + 3*blocking_micro_), ldb_tmp  , alpha , beta);
        if( !(useStreamingStores_ && useStreamingStores) )
           prefetch<floatType>(Bnext + (blocking_micro_ * ldb_tmp + 0), ldb_tmp);
        micro_kernel<floatType,betaIsZero, conjA>::execute(A + (0 * lda + blocking_micro_), lda, Btmp + (blocking_micro_ * ldb_tmp + 0), ldb_tmp  , alpha , beta);
        micro_kernel<floatType,betaIsZero, conjA>::execute(A + (blocking_micro_ * lda + blocking_micro_), lda, Btmp + (blocking_micro_ * ldb_tmp + blocking_micro_), ldb_tmp  , alpha , beta);
        if( !(useStreamingStores_ && useStreamingStores) )
           prefetch<floatType>(Bnext + (blocking_micro_ * ldb_tmp + 2*blocking_micro_), ldb_tmp);
        micro_kernel<floatType,betaIsZero, conjA>::execute(A + (2*blocking_micro_ * lda + blocking_micro_), lda, Btmp + (blocking_micro_ * ldb_tmp + 2*blocking_micro_), ldb_tmp  , alpha , beta);
        micro_kernel<floatType,betaIsZero, conjA>::execute(A + (3*blocking_micro_ * lda + blocking_micro_), lda, Btmp + (blocking_micro_ * ldb_tmp + 3*blocking_micro_), ldb_tmp  , alpha , beta);
     }else if( blockingA == blocking_ && blockingB == 2*blocking_micro_) {
        if( !(useStreamingStores_ && useStreamingStores) )
           prefetch<floatType>(Bnext + (0 * ldb_tmp + 0), ldb_tmp);
        prefetch<floatType>(Anext + (0 * lda + 0), lda);
        micro_kernel<floatType,betaIsZero, conjA>::execute(A + (0 * lda + 0), lda, Btmp + (0 * ldb_tmp + 0), ldb_tmp  , alpha , beta);
        prefetch<floatType>(Anext + (blocking_micro_ * lda + 0), lda);
        micro_kernel<floatType,betaIsZero, conjA>::execute(A + (blocking_micro_ * lda + 0), lda, Btmp + (0 * ldb_tmp + blocking_micro_), ldb_tmp  , alpha , beta);
        if( !(useStreamingStores_ && useStreamingStores) )
           prefetch<floatType>(Bnext + (blocking_micro_ * ldb_tmp + 0), ldb_tmp);
        micro_kernel<floatType,betaIsZero, conjA>::execute(A + (0 * lda + blocking_micro_), lda, Btmp + (blocking_micro_ * ldb_tmp + 0), ldb_tmp  , alpha , beta);
        micro_kernel<floatType,betaIsZero, conjA>::execute(A + (blocking_micro_ * lda + blocking_micro_), lda, Btmp + (blocking_micro_ * ldb_tmp + blocking_micro_), ldb_tmp  , alpha , beta);
        if( !(useStreamingStores_ && useStreamingStores) )
           prefetch<floatType>(Bnext + (2*blocking_micro_ * ldb_tmp + 0), ldb_tmp);
        prefetch<floatType>(Anext + (0 * lda + 2*blocking_micro_), lda);
        micro_kernel<floatType,betaIsZero, conjA>::execute(A + (0 * lda + 2*blocking_micro_), lda, Btmp + (2*blocking_micro_ * ldb_tmp + 0), ldb_tmp  , alpha , beta);
        prefetch<floatType>(Anext + (blocking_micro_ * lda + 2*blocking_micro_), lda);
        micro_kernel<floatType,betaIsZero, conjA>::execute(A + (blocking_micro_ * lda + 2*blocking_micro_), lda, Btmp + (2*blocking_micro_ * ldb_tmp + blocking_micro_), ldb_tmp  , alpha , beta);
        if( !(useStreamingStores_ && useStreamingStores) )
           prefetch<floatType>(Bnext + (3*blocking_micro_ * ldb_tmp + 0), ldb_tmp);
        micro_kernel<floatType,betaIsZero, conjA>::execute(A + (0 * lda + 3*blocking_micro_), lda, Btmp + (3*blocking_micro_ * ldb_tmp + 0), ldb_tmp  , alpha , beta);
        micro_kernel<floatType,betaIsZero, conjA>::execute(A + (blocking_micro_ * lda + 3*blocking_micro_), lda, Btmp + (3*blocking_micro_ * ldb_tmp + blocking_micro_), ldb_tmp  , alpha , beta);
     } else {
        //invoke micro-transpose
        if(blockingA > 0 && blockingB > 0 )
           micro_kernel<floatType,betaIsZero, conjA>::execute(A, lda, Btmp, ldb_tmp  , alpha , beta);
  
        //invoke micro-transpose
        if(blockingA > 0 && blockingB > blocking_micro_ )
           micro_kernel<floatType,betaIsZero, conjA>::execute(A + blocking_micro_ * lda, lda, Btmp + blocking_micro_, ldb_tmp  , alpha , beta);
  
        //invoke micro-transpose
        if(blockingA > 0 && blockingB > 2*blocking_micro_ )
           micro_kernel<floatType,betaIsZero, conjA>::execute(A + 2*blocking_micro_ * lda, lda, Btmp + 2*blocking_micro_, ldb_tmp  , alpha , beta);
  
        //invoke micro-transpose
        if(blockingA > 0 && blockingB > 3*blocking_micro_ )
           micro_kernel<floatType,betaIsZero, conjA>::execute(A + 3*blocking_micro_ * lda, lda, Btmp + 3*blocking_micro_, ldb_tmp  , alpha , beta);
  
        //invoke micro-transpose
        if(blockingA > blocking_micro_ && blockingB > 0 )
           micro_kernel<floatType,betaIsZero, conjA>::execute(A + blocking_micro_, lda, Btmp + blocking_micro_ * ldb_tmp, ldb_tmp  , alpha , beta);
  
        //invoke micro-transpose
        if(blockingA > blocking_micro_ && blockingB > blocking_micro_ )
           micro_kernel<floatType,betaIsZero, conjA>::execute(A + blocking_micro_ + blocking_micro_ * lda, lda, Btmp + blocking_micro_ + blocking_micro_ * ldb_tmp, ldb_tmp  , alpha , beta);
  
        //invoke micro-transpose
        if(blockingA > blocking_micro_ && blockingB > 2*blocking_micro_ )
           micro_kernel<floatType,betaIsZero, conjA>::execute(A + blocking_micro_ + 2*blocking_micro_ * lda, lda, Btmp + 2*blocking_micro_ + blocking_micro_ * ldb_tmp, ldb_tmp  , alpha , beta);
  
        //invoke micro-transpose
        if(blockingA > blocking_micro_ && blockingB > 3*blocking_micro_ )
           micro_kernel<floatType,betaIsZero, conjA>::execute(A + blocking_micro_ + 3*blocking_micro_ * lda, lda, Btmp + 3*blocking_micro_ + blocking_micro_ * ldb_tmp, ldb_tmp  , alpha , beta);
  
        //invoke micro-transpose
        if(blockingA > 2*blocking_micro_ && blockingB > 0 )
           micro_kernel<floatType,betaIsZero, conjA>::execute(A + 2*blocking_micro_, lda, Btmp + 2*blocking_micro_ * ldb_tmp, ldb_tmp  , alpha , beta);
  
        //invoke micro-transpose
        if(blockingA > 2*blocking_micro_ && blockingB > blocking_micro_ )
           micro_kernel<floatType,betaIsZero, conjA>::execute(A + 2*blocking_micro_ + blocking_micro_ * lda, lda, Btmp + blocking_micro_ + 2*blocking_micro_ * ldb_tmp, ldb_tmp  , alpha , beta);
  
        //invoke micro-transpose
        if(blockingA > 2*blocking_micro_ && blockingB > 2*blocking_micro_ )
           micro_kernel<floatType,betaIsZero, conjA>::execute(A + 2*blocking_micro_ + 2*blocking_micro_ * lda, lda, Btmp + 2*blocking_micro_ + 2*blocking_micro_ * ldb_tmp, ldb_tmp  , alpha , beta);
  
        //invoke micro-transpose
        if(blockingA > 2*blocking_micro_ && blockingB > 3*blocking_micro_ )
           micro_kernel<floatType,betaIsZero, conjA>::execute(A + 2*blocking_micro_ + 3*blocking_micro_ * lda, lda, Btmp + 3*blocking_micro_ + 2*blocking_micro_ * ldb_tmp, ldb_tmp  , alpha , beta);
  
        //invoke micro-transpose
        if(blockingA > 3*blocking_micro_ && blockingB > 0 )
           micro_kernel<floatType,betaIsZero, conjA>::execute(A + 3*blocking_micro_, lda, Btmp + 3*blocking_micro_ * ldb_tmp, ldb_tmp  , alpha , beta);
  
        //invoke micro-transpose
        if(blockingA > 3*blocking_micro_ && blockingB > blocking_micro_ )
           micro_kernel<floatType,betaIsZero, conjA>::execute(A + 3*blocking_micro_ + blocking_micro_ * lda, lda, Btmp + blocking_micro_ + 3*blocking_micro_ * ldb_tmp, ldb_tmp  , alpha , beta);
  
        //invoke micro-transpose
        if(blockingA > 3*blocking_micro_ && blockingB > 2*blocking_micro_ )
           micro_kernel<floatType,betaIsZero, conjA>::execute(A + 3*blocking_micro_ + 2*blocking_micro_ * lda, lda, Btmp + 2*blocking_micro_ + 3*blocking_micro_ * ldb_tmp, ldb_tmp  , alpha , beta);
  
        //invoke micro-transpose
        if(blockingA > 3*blocking_micro_ && blockingB > 3*blocking_micro_ )
           micro_kernel<floatType,betaIsZero, conjA>::execute(A + 3*blocking_micro_ + 3*blocking_micro_ * lda, lda, Btmp + 3*blocking_micro_ + 3*blocking_micro_ * ldb_tmp, ldb_tmp  , alpha , beta);
     }
  
     // write buffer to main-memory via non-temporal stores
     if( (useStreamingStores_ && useStreamingStores) )
        for( int i = 0; i < blockingA; i++){
           for( int j = 0; j < blockingB; j+=blocking_micro_)
              streamingStore<floatType>(B + i * ldb + j, buffer + i * ldb_tmp + j);
        }
  }
  
  template<int betaIsZero, typename floatType, bool conjA>
  void transpose_int_scalar( const floatType* __restrict__ A, int sizeStride1A, 
                                    floatType* __restrict__ B, int sizeStride1B, const floatType alpha, const floatType beta, const ComputeNode* plan)
  {
     const int32_t end = plan->end;
     const size_t lda = plan->lda;
     const size_t ldb = plan->ldb;
     if( plan->next->next != nullptr ){
        // recurse
        int i = plan->start;
        if( lda == 1)
           transpose_int_scalar<betaIsZero, floatType, conjA>( &A[i*lda], end - plan->start, &B[i*ldb], sizeStride1B, alpha, beta, plan->next);
        else if( ldb == 1)
           transpose_int_scalar<betaIsZero, floatType, conjA>( &A[i*lda], sizeStride1A, &B[i*ldb], end - plan->start, alpha, beta, plan->next);
        else
           for(; i < end; i++)
              transpose_int_scalar<betaIsZero, floatType, conjA>( &A[i*lda], sizeStride1A, &B[i*ldb], sizeStride1B, alpha, beta, plan->next);
     }else{
        // macro-kernel
        const size_t lda_macro = plan->next->lda;
        const size_t ldb_macro = plan->next->ldb;
        int i = plan->start;
        const size_t scalarRemainder = plan->end - plan->start;
        if( scalarRemainder > 0 ){
           if( lda == 1)
              macro_kernel_scalar<betaIsZero,floatType, conjA>(&A[i*lda], lda_macro, scalarRemainder, &B[i*ldb], ldb_macro, sizeStride1B, alpha, beta);
           else if( ldb == 1)
              macro_kernel_scalar<betaIsZero,floatType, conjA>(&A[i*lda], lda_macro, sizeStride1A, &B[i*ldb], ldb_macro, scalarRemainder, alpha, beta);
        }
     }
  }
  template<int blockingA, int blockingB, int betaIsZero, typename floatType, bool useStreamingStores, bool conjA>
  void transpose_int( const floatType* __restrict__ A, const floatType* __restrict__ Anext, 
                       floatType* __restrict__ B, const floatType* __restrict__ Bnext, const floatType alpha, const floatType beta, 
                       const ComputeNode* plan)
  {
     const int32_t end = plan->end - (plan->inc - 1);
     const int32_t inc = plan->inc;
     const size_t lda = plan->lda;
     const size_t ldb = plan->ldb;
  
     constexpr int blocking_micro_ = REGISTER_BITS/8 / sizeof(floatType);
     constexpr int blocking_ = blocking_micro_ * 4;
  
     if( plan->next->next != nullptr ){
        // recurse
        int32_t i;
        for(i = plan->start; i < end; i+= inc)
        {
           if( i + inc < end )
              transpose_int<blockingA, blockingB, betaIsZero, floatType, useStreamingStores, conjA>( &A[i*lda], &A[(i+1)*lda], &B[i*ldb], &B[(i+1)*ldb], alpha, beta, plan->next);
           else
              transpose_int<blockingA, blockingB, betaIsZero, floatType, useStreamingStores, conjA>( &A[i*lda], Anext, &B[i*ldb], Bnext, alpha, beta, plan->next);
        }
        // remainder
        if( blocking_/2 >= blocking_micro_ && (i + blocking_/2) <= plan->end ){
           if( lda == 1)
              transpose_int<blocking_/2, blockingB, betaIsZero, floatType, useStreamingStores, conjA>( &A[i*lda], Anext, &B[i*ldb], Bnext, alpha, beta, plan->next);
           else if( ldb == 1)
              transpose_int<blockingA, blocking_/2, betaIsZero, floatType, useStreamingStores, conjA>( &A[i*lda], Anext, &B[i*ldb], Bnext, alpha, beta, plan->next);
           i+=blocking_/2;
        }
        if( blocking_/4 >= blocking_micro_ && (i + blocking_/4) <= plan->end ){
           if( lda == 1)
              transpose_int<blocking_/4, blockingB, betaIsZero, floatType, useStreamingStores, conjA>( &A[i*lda], Anext, &B[i*ldb], Bnext, alpha, beta, plan->next);
           else if( ldb == 1)
              transpose_int<blockingA, blocking_/4, betaIsZero, floatType, useStreamingStores, conjA>( &A[i*lda], Anext, &B[i*ldb], Bnext, alpha, beta, plan->next);
           i+=blocking_/4;
        }
        const size_t scalarRemainder = plan->end - i;
        if( scalarRemainder > 0 ){
           if( lda == 1)
              transpose_int_scalar<betaIsZero, floatType, conjA>( &A[i*lda], scalarRemainder, &B[i*ldb], -1, alpha, beta, plan->next);
           else
              transpose_int_scalar<betaIsZero, floatType, conjA>( &A[i*lda], -1, &B[i*ldb], scalarRemainder, alpha, beta, plan->next);
        }
     } else {
        const size_t lda_macro = plan->next->lda;
        const size_t ldb_macro = plan->next->ldb;
        // invoke macro-kernel
        
        int32_t i;
        for(i = plan->start; i < end; i+= inc)
           if( i + inc < end )
              macro_kernel<blockingA, blockingB, betaIsZero,floatType, useStreamingStores, conjA>(&A[i*lda], &A[(i+1)*lda], lda_macro, &B[i*ldb], &B[(i+1)*ldb], ldb_macro, alpha, beta);
           else
              macro_kernel<blockingA, blockingB, betaIsZero,floatType, useStreamingStores, conjA>(&A[i*lda], Anext, lda_macro, &B[i*ldb], Bnext, ldb_macro, alpha, beta);
        // remainder
        if( blocking_/2 >= blocking_micro_ && (i + blocking_/2) <= plan->end ){
           if( lda == 1)
              macro_kernel<blocking_/2, blockingB, betaIsZero,floatType, useStreamingStores, conjA>(&A[i*lda], Anext, lda_macro, &B[i*ldb], Bnext, ldb_macro, alpha, beta);
           else if( ldb == 1)
              macro_kernel<blockingA, blocking_/2, betaIsZero,floatType, useStreamingStores, conjA>(&A[i*lda], Anext, lda_macro, &B[i*ldb], Bnext, ldb_macro, alpha, beta);
           i+=blocking_/2;
        }
        if( blocking_/4 >= blocking_micro_ && (i + blocking_/4) <= plan->end ){
           if( lda == 1)
              macro_kernel<blocking_/4, blockingB, betaIsZero,floatType, useStreamingStores, conjA>(&A[i*lda], Anext, lda_macro, &B[i*ldb], Bnext, ldb_macro, alpha, beta);
           else if( ldb == 1)
              macro_kernel<blockingA, blocking_/4, betaIsZero,floatType, useStreamingStores, conjA>(&A[i*lda], Anext, lda_macro, &B[i*ldb], Bnext, ldb_macro, alpha, beta);
           i+=blocking_/4;
        }
        const size_t scalarRemainder = plan->end - i;
        if( scalarRemainder > 0 ){
           if( lda == 1)
              macro_kernel_scalar<betaIsZero,floatType, conjA>(&A[i*lda], lda_macro, scalarRemainder, &B[i*ldb], ldb_macro, blockingB, alpha, beta);
           else if( ldb == 1)
              macro_kernel_scalar<betaIsZero,floatType, conjA>(&A[i*lda], lda_macro, blockingA, &B[i*ldb], ldb_macro, scalarRemainder, alpha, beta);
        }
     }
  }
  
  template<int betaIsZero, typename floatType, bool useStreamingStores, bool conjA>
  void transpose_int_constStride1( const floatType* __restrict__ A, floatType* __restrict__ B, const floatType alpha, const floatType beta, 
        const ComputeNode* plan)
  {
     const int32_t end = plan->end - (plan->inc - 1);
     constexpr int32_t inc = 1; // TODO
     const size_t lda = plan->lda;
     const size_t ldb = plan->ldb;
  
     if( plan->next != nullptr )
        for(int i = plan->start; i < end; i+= inc)
           // recurse
           transpose_int_constStride1<betaIsZero, floatType, useStreamingStores, conjA>( &A[i*lda], &B[i*ldb], alpha, beta, plan->next);
     else 
        if( !betaIsZero )
        {
           for(int32_t i = plan->start; i < end; i+= inc)
              if( conjA )
                 B[i] = alpha * conj(A[i]) + beta * B[i];
              else
                 B[i] = alpha * A[i] + beta * B[i];
        } else {
           if( useStreamingStores)
              if( conjA )
  #pragma vector nontemporal
                 for(int32_t i = plan->start; i < end; i+= inc)
                    B[i] = alpha * conj(A[i]);
              else
  #pragma vector nontemporal
                 for(int32_t i = plan->start; i < end; i+= inc)
                    B[i] = alpha * A[i];
           else
              if( conjA )
                 for(int32_t i = plan->start; i < end; i+= inc)
                    B[i] = alpha * conj(A[i]);
              else
                 for(int32_t i = plan->start; i < end; i+= inc)
                    B[i] = alpha * A[i];
        }
  }
  
  
     template<typename floatType>
     Transpose<floatType>::Transpose( const int *sizeA, 
                   const int *perm, 
                   const int *outerSizeA, 
                   const int *outerSizeB, 
                   const int dim,
                   const floatType *A,
                   const floatType alpha,
                   floatType *B,
                   const floatType beta,
                   const SelectionMethod selectionMethod,
                   const int numThreads, 
                   const int *threadIds, 
                   const bool useRowMajor) : 
           A_(A),
           B_(B),
           alpha_(alpha),
           beta_(beta),
           dim_(-1),
           numThreads_(numThreads), 
           masterPlan_(nullptr),
           selectionMethod_(selectionMethod),
           maxAutotuningCandidates_(-1),
           selectedParallelStrategyId_(-1),
           selectedLoopOrderId_(-1),
           conjA_(false)
        {
  #ifdef _OPENMP
           omp_init_lock(&writelock);
  #endif
           int tmpPerm[dim];
           int tmpSizeA[dim];
           int tmpOuterSizeA[dim];
           int tmpOuterSizeB[dim];
           accountForRowMajor(sizeA, outerSizeA, outerSizeB, perm, 
                 tmpSizeA, tmpOuterSizeA, tmpOuterSizeB, tmpPerm, dim, useRowMajor); 
  
           sizeA_.resize(dim);
           perm_.resize(dim);
           outerSizeA_.resize(dim);
           outerSizeB_.resize(dim);
           lda_.resize(dim);
           ldb_.resize(dim);
           if(threadIds){
              // compact threadIds. E.g., 1, 7, 5 -> local_id(1) = 0, local_id(7) = 2, local_id(5) = 1
              for(int i=0; i < numThreads; ++i)
                 threadIds_.push_back(threadIds[i]);
              std::sort(threadIds_.begin(), threadIds_.end());
           }else{
              for(int i=0; i < numThreads; ++i)
                 threadIds_.push_back(i);
           }
  
           verifyParameter(tmpSizeA, tmpPerm, tmpOuterSizeA, tmpOuterSizeB, dim);
  
           // initializes dim_, outerSizeA, outerSizeB, sizeA and perm 
           skipIndices(tmpSizeA, tmpPerm, tmpOuterSizeA, tmpOuterSizeB, dim);
           fuseIndices();
  
           // initializes lda_ and ldb_
           computeLeadingDimensions();
  
           // create plan
           this->createPlan();
        }
  
        template<typename floatType>
        Transpose<floatType>::Transpose(const Transpose<floatType> &other) : A_(other.A_), B_(other.B_), 
                                            alpha_(other.alpha_),
                                            beta_(other.beta_),
                                            dim_(other.dim_),
                                            numThreads_(other.numThreads_),
                                            masterPlan_(other.masterPlan_),
                                            selectionMethod_(other.selectionMethod_),
                                            selectedParallelStrategyId_(other.selectedParallelStrategyId_),
                                            selectedLoopOrderId_(other.selectedLoopOrderId_),
                                            maxAutotuningCandidates_(other.maxAutotuningCandidates_),
                                            sizeA_(other.sizeA_),
                                            perm_(other.perm_),
                                            outerSizeA_(other.outerSizeA_),
                                            outerSizeB_(other.outerSizeB_),
                                            lda_(other.lda_),
                                            ldb_(other.ldb_),
                                            threadIds_(other.threadIds_),
                                            conjA_(other.conjA_)
        { 
  #ifdef _OPENMP
           omp_init_lock(&writelock);
  #endif
        }
  
        template<typename floatType>
        Transpose<floatType>::~Transpose() { 
  #ifdef _OPENMP
           omp_destroy_lock(&writelock); 
  #endif
        }
  
  template<typename floatType>
  void Transpose<floatType>::executeEstimate(const Plan *plan) noexcept
  {
     if( plan == nullptr ) {
        fprintf(stderr,"[HPTT] ERROR: plan has not yet been created.\n");
        exit(-1);
     }
     
     constexpr bool useStreamingStores = false;
  
     const int numTasks = plan->getNumTasks();
  #ifdef _OPENMP
  #pragma omp parallel for num_threads(numThreads_)  if(numThreads_ > 1)
  #endif
     for( int taskId = 0; taskId < numTasks; taskId++)
        if ( perm_[0] != 0 ) {
           auto rootNode = plan->getRootNode_const( taskId );
           if( std::abs(beta_) < getZeroThreshold<floatType>() ) {
              if( conjA_ )
                 transpose_int<blocking_,blocking_,1,floatType, useStreamingStores, true>( A_,A_, B_, B_, 0.0, 1.0, rootNode);
              else
                 transpose_int<blocking_,blocking_,1,floatType, useStreamingStores, false>( A_,A_, B_, B_, 0.0, 1.0, rootNode);
           } else {
              if( conjA_ )
                 transpose_int<blocking_,blocking_,0,floatType, useStreamingStores, true>( A_,A_, B_, B_, 0.0, 1.0, rootNode);
              else
                 transpose_int<blocking_,blocking_,0,floatType, useStreamingStores, false>( A_,A_, B_, B_, 0.0, 1.0, rootNode);
           }
        } else {
           auto rootNode = plan->getRootNode_const( taskId );
           if( std::abs(beta_) < getZeroThreshold<floatType>() ) {
              if( conjA_ )
                 transpose_int_constStride1<1,floatType, useStreamingStores, true>( A_, B_, 0.0, 1.0, rootNode);
              else
                 transpose_int_constStride1<1,floatType, useStreamingStores, false>( A_, B_, 0.0, 1.0, rootNode);
           }else{
              if( conjA_ )
                 transpose_int_constStride1<0,floatType, useStreamingStores, true>( A_, B_, 0.0, 1.0, rootNode);
              else
                 transpose_int_constStride1<0,floatType, useStreamingStores, false>( A_, B_, 0.0, 1.0, rootNode);
           }
        }
  }
  
  
  
  template<int betaIsZero, typename floatType, bool useStreamingStores, bool spawnThreads, bool conjA>
  static void axpy_1D( const floatType* __restrict__ A, floatType* __restrict__ B, const int myStart, const int myEnd, const floatType alpha, const floatType beta, int numThreads)
  {
     if( !betaIsZero )
     {
        HPTT_DUPLICATE(spawnThreads,
           for(int32_t i = myStart; i < myEnd; i++)
              if( conjA )
                 B[i] = alpha * conj(A[i]) + beta * B[i];
              else
                 B[i] = alpha * A[i] + beta * B[i];
        )
     } else {
        if( useStreamingStores)
  #pragma vector nontemporal
           HPTT_DUPLICATE(spawnThreads,
              for(int32_t i = myStart; i < myEnd; i++)
                 if( conjA )
                    B[i] = alpha * conj(A[i]);
                 else
                    B[i] = alpha * A[i];
           )
        else
           HPTT_DUPLICATE(spawnThreads,
              for(int32_t i = myStart; i < myEnd; i++)
                 if( conjA )
                    B[i] = alpha * conj(A[i]);
                 else
                    B[i] = alpha * A[i];
           )
     }
  }
  
  template<int betaIsZero, typename floatType, bool useStreamingStores, bool spawnThreads, bool conjA>
  static void axpy_2D( const floatType* __restrict__ A, const int lda, 
                          floatType* __restrict__ B, const int ldb, 
                          const int n0, const int myStart, const int myEnd, const floatType alpha, const floatType beta, int numThreads)
  {
     if( !betaIsZero )
     {
        HPTT_DUPLICATE(spawnThreads,
           for(int32_t j = myStart; j < myEnd; j++)
              for(int32_t i = 0; i < n0; i++)
                 if( conjA )
                    B[i + j * ldb] = alpha * conj(A[i + j * lda]) + beta * B[i + j * ldb];
                 else
                    B[i + j * ldb] = alpha * A[i + j * lda] + beta * B[i + j * ldb];
        )
     } else {
        if( useStreamingStores)
           HPTT_DUPLICATE(spawnThreads,
              for(int32_t j = myStart; j < myEnd; j++)
  _Pragma("vector nontemporal")
              for(int32_t i = 0; i < n0; i++)
                 if( conjA )
                    B[i + j * ldb] = alpha * conj(A[i + j * lda]);
                 else
                    B[i + j * ldb] = alpha * A[i + j * lda];
           )
        else
           HPTT_DUPLICATE(spawnThreads,
              for(int32_t j = myStart; j < myEnd; j++)
                 for(int32_t i = 0; i < n0; i++)
                    if( conjA )
                       B[i + j * ldb] = alpha * conj(A[i + j * lda]);
                    else
                       B[i + j * ldb] = alpha * A[i + j * lda];
           )
     }
  }
  
  template<typename floatType> 
  template<bool spawnThreads>
  void Transpose<floatType>::getStartEnd(int n, int &myStart, int &myEnd) const
  {
  #ifdef _OPENMP
     int myLocalThreadId = getLocalThreadId(omp_get_thread_num());
  #else
     int myLocalThreadId = 0;
  #endif
  
     if(myLocalThreadId == -1 ) // skip those threads which do not participate in this plan
     {
        myStart = n;
        myEnd = n;
        return;
     }
     if(spawnThreads){ // worksharing will be handled by the OpenMP runtime
        myStart = 0;
        myEnd = n;
        return;
     }
  
     const int workPerThread = (n + numThreads_-1)/numThreads_;
     myStart = std::min(n, myLocalThreadId * workPerThread);
     myEnd = std::min(n, (myLocalThreadId+1) * workPerThread);
  }
  
  
  template<typename floatType> 
  int Transpose<floatType>::getLocalThreadId(int myThreadId) const
  {
     int myLocalId = -1;
     for(int i=0; i < numThreads_; ++i)
        if(myThreadId == threadIds_[i])
           myLocalId = i;
     return myLocalId;
  }
  
  template<typename floatType> 
  template<bool useStreamingStores, bool spawnThreads, bool betaIsZero>
  void Transpose<floatType>::execute_expert() noexcept
  {
     if( masterPlan_ == nullptr ) {
        fprintf(stderr,"[HPTT] ERROR: master plan has not yet been created.\n");
        exit(-1);
     }
  
     int myStart = 0;
     int myEnd = 0;
  
     if( dim_ == 1)
     {
        getStartEnd<spawnThreads>(sizeA_[0], myStart, myEnd);
        if( conjA_ )
           axpy_1D<betaIsZero, floatType, useStreamingStores, spawnThreads, true>( A_, B_, myStart, myEnd, alpha_, beta_, numThreads_ );
        else
           axpy_1D<betaIsZero, floatType, useStreamingStores, spawnThreads, false>( A_, B_, myStart, myEnd, alpha_, beta_, numThreads_ );
        return;
     }
     else if( dim_ == 2 && perm_[0] == 0)
     {
        getStartEnd<spawnThreads>(sizeA_[1], myStart, myEnd);
        if( conjA_ )
           axpy_2D<betaIsZero, floatType, useStreamingStores, spawnThreads, true>( A_, lda_[1], B_, ldb_[1], sizeA_[0], myStart, myEnd, alpha_, beta_, numThreads_ );
        else
           axpy_2D<betaIsZero, floatType, useStreamingStores, spawnThreads, false>( A_, lda_[1], B_, ldb_[1], sizeA_[0], myStart, myEnd, alpha_, beta_, numThreads_ );
        return;
     }
  
     const int numTasks = masterPlan_->getNumTasks();
     const int numThreads = numThreads_;
     getStartEnd<spawnThreads>(numTasks, myStart, myEnd);
  
     HPTT_DUPLICATE(spawnThreads,
        for( int taskId = myStart; taskId < myEnd; taskId++)
           if ( perm_[0] != 0 ) {
              auto rootNode = masterPlan_->getRootNode_const( taskId );
              if( conjA_ )
                 transpose_int<blocking_,blocking_,betaIsZero,floatType, useStreamingStores, true>( A_, A_, B_, B_, alpha_, beta_, rootNode);
              else
                 transpose_int<blocking_,blocking_,betaIsZero,floatType, useStreamingStores, false>( A_, A_, B_, B_, alpha_, beta_, rootNode);
           } else {
              auto rootNode = masterPlan_->getRootNode_const( taskId );
              if( conjA_ )
                 transpose_int_constStride1<betaIsZero,floatType, useStreamingStores, true>( A_, B_, alpha_, beta_, rootNode);
              else
                 transpose_int_constStride1<betaIsZero,floatType, useStreamingStores, false>( A_, B_, alpha_, beta_, rootNode);
           }
     )
  }
  template<typename floatType> 
  void Transpose<floatType>::execute() noexcept
  {
     if( masterPlan_ == nullptr ) {
        fprintf(stderr,"[HPTT] ERROR: master plan has not yet been created.\n");
        exit(-1);
     }
  
     bool spawnThreads = numThreads_ > 1;
     bool betaIsZero = (beta_ == (floatType) 0.0);
     constexpr bool useStreamingStores = true;
     if( spawnThreads )
     {
        if( betaIsZero) {
           this->execute_expert<useStreamingStores, true, true>();
        } else {
           this->execute_expert<useStreamingStores, true, false>();
        }
     } else {
        if( betaIsZero) {
           this->execute_expert<useStreamingStores, false, true>();
        } else {
           this->execute_expert<useStreamingStores, false, false>();
        }
     }
  }
  
  template<typename floatType> 
  void Transpose<floatType>::print() noexcept 
  { 
     masterPlan_->print(); 
  }
  
  template<typename floatType>
  int Transpose<floatType>::getIncrement( int loopIdx ) const
  {
     int inc = 1;
     if( perm_[0] != 0 ) {
        if( loopIdx == 0 || loopIdx == perm_[0] )
           inc = blocking_;
     }
     return inc;
  }
  
  template<typename floatType>
  void Transpose<floatType>::getAvailableParallelism( std::vector<int> &numTasksPerLoop ) const
  {
     numTasksPerLoop.resize(dim_);
     for(int loopIdx=0; loopIdx < dim_; ++loopIdx){
        int inc = this->getIncrement(loopIdx);
        numTasksPerLoop[loopIdx] = (sizeA_[loopIdx] + inc - 1) / inc;
     }
  }
  
  template<typename floatType>
  void Transpose<floatType>::getAllParallelismStrategies( std::list<int> &primeFactorsToMatch, 
                                               std::vector<int> &availableParallelismAtLoop, 
                                               std::vector<int> &achievedParallelismAtLoop, 
                                               std::vector<std::vector<int> > &parallelismStrategies) const
  {
     if( primeFactorsToMatch.size() > 0 ){
        // match every primefactor ...
        for( auto p : primeFactorsToMatch )
        {
           // ... with every loop
           for( int i = 0; i < dim_; i++ )
           {
              std::list<int> primeFactorsToMatch_(primeFactorsToMatch); 
              std::vector<int> availableParallelismAtLoop_(availableParallelismAtLoop); 
              std::vector<int> achievedParallelismAtLoop_(achievedParallelismAtLoop);
  
              primeFactorsToMatch_.erase( std::find(primeFactorsToMatch_.begin(), primeFactorsToMatch_.end(), p) );
              availableParallelismAtLoop_[i] = (availableParallelismAtLoop_[i] + p - 1) / p;
              achievedParallelismAtLoop_[i] *= p;
  
              this->getAllParallelismStrategies( primeFactorsToMatch_, 
                    availableParallelismAtLoop_, 
                    achievedParallelismAtLoop_, 
                    parallelismStrategies);
           }
        }
     } else {
        // avoid duplicates
        if( parallelismStrategies.end() == std::find(parallelismStrategies.begin(), parallelismStrategies.end(), achievedParallelismAtLoop) )
           parallelismStrategies.push_back(achievedParallelismAtLoop);
     }
  }
  
  // balancing if one tries to parallelize avail many tasks with req many threads
  // e.g., balancing(3,4) = 0.75
  static float getBalancing(int avail, int req)
  {
     return ((float)(avail))/(((avail+req-1)/req)*req);
  }
  
  template<typename floatType>
  float Transpose<floatType>::getLoadBalance( const std::vector<int> &parallelismStrategy ) const
  {
     float load_balance = 1.0;
     int totalTasks = 1;
     for(int i=0; i < dim_; ++i){
  
        int inc = this->getIncrement(i);
        while(sizeA_[i] < inc)
           inc /= 2;
        int availableParallelism = (sizeA_[i] + inc - 1) / inc;
  
        if( i == 0 || perm_[i] == 0) 
           // account for the load-imbalancing due to blocking
           load_balance *= getBalancing(sizeA_[i], inc); 
        load_balance *= getBalancing(availableParallelism, parallelismStrategy[i]); 
        totalTasks *= parallelismStrategy[i];
     }
  
     //how well can these tasks be distributed among numThreads_?
     // e.g., totalTasks = 3, numThreads = 8 => 3./8
     // e.g., totalTasks = 5, numThreads = 8 => 5./8
     // e.g., totalTasks = 15, numThreads = 8 => 15./16
     // e.g., totalTasks = 17, numThreads = 8 => 17./24
     float workDistribution = ((float)totalTasks) / (((totalTasks + numThreads_ - 1)/numThreads_)*numThreads_);
  
     load_balance *= workDistribution ;
     return load_balance; 
  }
  
  template<typename floatType>
  void Transpose<floatType>::getBestParallelismStrategy ( std::vector<int> &bestParallelismStrategy ) const
  {
     std::vector<int> availableParallelismAtLoop;
     this->getAvailableParallelism( availableParallelismAtLoop);
     int totalAvailableParallelism = std::accumulate(availableParallelismAtLoop.begin(), 
                                   availableParallelismAtLoop.end(), 1, std::multiplies<int>());
  
     // reduce the probability of parallelizing the stride-1 index
     // if this loops would be parallelized, these two statements ensure that each
     // thread would have at least two macro-kernels of work at this loop-level
     //
     // However, if the total available parallelism is too small, then we do not
     // artificially limit the available parallelism further
     int reduceParallelismB = 4; //avoid parallelization in stride-1 B more strongly
     int reduceParallelismA = 2;
     if( totalAvailableParallelism < 2 * numThreads_ )
        reduceParallelismB = 1;
     else if( totalAvailableParallelism < 4 * numThreads_ )
        reduceParallelismB = 2;
     totalAvailableParallelism = (totalAvailableParallelism / availableParallelismAtLoop[perm_[0]]) * 
                                  (availableParallelismAtLoop[perm_[0]]/reduceParallelismB);
     if( totalAvailableParallelism < 2 * numThreads_ )
        reduceParallelismA = 1;
     availableParallelismAtLoop[perm_[0]] = std::max(1, availableParallelismAtLoop[perm_[0]] 
                                                     / reduceParallelismB ); 
     availableParallelismAtLoop[0] = std::max(1, availableParallelismAtLoop[0] / reduceParallelismA );
  
     // Objectives: 1) load-balancing 
     //             2) avoid parallelizing stride-1 loops (rational: less consecutive memory accesses)
     //             3) avoid false sharing
  
     std::vector<int> loopsAllowed;
     for(int i=dim_-1; i >= 1; i--)
        if(perm_[i] != 0)
           loopsAllowed.push_back(perm_[i]);
     std::vector<int> loopsAllowedStride1 {0,perm_[0]};
  
     int totalTasks = 1; // goal: totalTasks should be a close multiple of numTasks_
     std::list<int> primeFactors;
     getPrimeFactors( numThreads_, primeFactors );
  
     // 1. parallelize using 100% load balancing
     parallelize( bestParallelismStrategy, availableParallelismAtLoop, totalTasks, primeFactors, 1.0, loopsAllowed);
   
     if( totalTasks != numThreads_ )
     { // no perfect match has been found
  
        // Option 1: keep parallelizing non-stride-1 loops only, but allowing load-imbalance
        std::vector<int> strat1(bestParallelismStrategy);
        std::vector<int> avail1(availableParallelismAtLoop);
        std::list<int> primes1(primeFactors);
        int totalTasks1 = totalTasks;
        parallelize( strat1, avail1, totalTasks1, primes1, 0.92, loopsAllowed);
        if( getLoadBalance(strat1) > 0.90 )
        {
           std::copy(strat1.begin(), strat1.end(), bestParallelismStrategy.begin());
           return;
        }
  
        if(perm_[0] != 0)
        {
           // Option 2: also parallelize stride-1 loops, enforcing perfect loop balancing
           std::vector<int> strat2(bestParallelismStrategy);
           std::vector<int> avail2(availableParallelismAtLoop);
           std::list<int> primes2(primeFactors);
           int totalTasks2 = totalTasks;
           parallelize( strat2, avail2, totalTasks2, primes2, 1.0, loopsAllowedStride1 );
           if( getLoadBalance(strat2) > 0.92 )
           {
              std::copy(strat2.begin(), strat2.end(), bestParallelismStrategy.begin());
              return;
           }
  
           //keep on going based on strat1
           parallelize( strat1, avail1, totalTasks1, primes1, 1.0, loopsAllowedStride1 );
           if( getLoadBalance(strat1) > 0.90 )
           {
              std::copy(strat1.begin(), strat1.end(), bestParallelismStrategy.begin());
              return;
           }
  
           //keep on going based on strat2
           parallelize( strat2, avail2, totalTasks2, primes2, 0.92, loopsAllowed);
           if( getLoadBalance(strat2) > 0.92 )
           {
              std::copy(strat2.begin(), strat2.end(), bestParallelismStrategy.begin());
              return;
           }
  
           if( getLoadBalance(strat1) > 0.80 ) //reduced threshold
           {
              std::copy(strat1.begin(), strat1.end(), bestParallelismStrategy.begin());
              return;
           }
           if( getLoadBalance(strat2) > 0.82 ) //reduced threshold
           {
              std::copy(strat2.begin(), strat2.end(), bestParallelismStrategy.begin());
              return;
           }
  
           parallelize( strat1, avail1, totalTasks1, primes1, 0.9, loopsAllowedStride1 );
           parallelize( strat2, avail2, totalTasks2, primes2, 0.8, loopsAllowed );
           float lb1 = getLoadBalance(strat1);
           float lb2 = getLoadBalance(strat2);
  //         printVector(strat2,"strat2");
  //         printf("strat2: %f\n",getLoadBalance(strat2));
           if( lb1 > 0.8 && lb2 < 0.85 || lb1 >lb2 && lb1 > 0.75 )
           {
              std::copy(strat1.begin(), strat1.end(), bestParallelismStrategy.begin());
              return;
           }
           if(lb2 >= 0.85)
           {
              std::copy(strat2.begin(), strat2.end(), bestParallelismStrategy.begin());
              return;
           }
  
           //fallback
           std::vector<int> allLoops;
           for(int i=dim_-1; i >= 1; i--)
              allLoops.push_back(perm_[i]);
           allLoops.push_back(0);
           allLoops.push_back(perm_[0]);
           parallelize( strat1, avail1, totalTasks1, primes1, 0., allLoops);
           std::copy(strat1.begin(), strat1.end(), bestParallelismStrategy.begin());
  
        }else{
           parallelize( strat1, avail1, totalTasks1, primes1, 0.0, loopsAllowed);
           std::copy(strat1.begin(), strat1.end(), bestParallelismStrategy.begin());
        }
     }
  }
  
  template<typename floatType>
  void Transpose<floatType>::parallelize( std::vector<int> &parallelismStrategy,
                                          std::vector<int> &availableParallelismAtLoop, 
                                          int &totalTasks,
                                          std::list<int> &primeFactors, 
                                          const float minBalancing, const std::vector<int> &loopsAllowed ) const
     
  {
     bool suboptimalParallelizationUsed = false;
     // find loop which minimizes load imbalance for the given prime factor
     for(auto it = primeFactors.begin(); it != primeFactors.end(); it++)
     {
        int suitedLoop = -1;
        float bestBalancing = 0;
  
        for(auto idx : loopsAllowed )
        {
           float balancing = getBalancing(availableParallelismAtLoop[idx], *it);
           if(balancing > bestBalancing){
              bestBalancing = balancing;
              suitedLoop = idx;
           }
        }
        //allow up to one slightly less optimal splitting to prefer parallelizing
        //idx=0 over idx=perm[0]
        if( suboptimalParallelizationUsed == false &&
              suitedLoop == perm_[0] && getBalancing(availableParallelismAtLoop[0], *it) >= 0.949 ){
           suitedLoop = 0;
           suboptimalParallelizationUsed = true;
        }
        if( suitedLoop != -1 && bestBalancing >= minBalancing ){
           availableParallelismAtLoop[suitedLoop] /= *it;
           parallelismStrategy[suitedLoop] *= *it;
           totalTasks *= *it;
           it = primeFactors.erase(it);
           it--;
        }
     }
  }
  
  template<typename floatType>
  double Transpose<floatType>::parallelismCostHeuristic( const std::vector<int> &achievedParallelismAtLoop ) const
  {
     std::vector<int> availableParallelismAtLoop;
     this->getAvailableParallelism( availableParallelismAtLoop);
  
     double cost = 1;
     // penalize load-imbalance
     for(int loopIdx=0; loopIdx < dim_ ; ++loopIdx){
        if ( achievedParallelismAtLoop[loopIdx] <= 1 ) 
           continue;
  
        const int blocksPerThread = (availableParallelismAtLoop[loopIdx] + achievedParallelismAtLoop[loopIdx] -1) / achievedParallelismAtLoop[loopIdx];
        int inc = this->getIncrement( loopIdx );
        const int effectiveSize = blocksPerThread * inc * achievedParallelismAtLoop[loopIdx];
        cost *= ((double)(effectiveSize ) / sizeA_[loopIdx]);
     } 
  
     // penalize parallelization of stride-1 loops
     if( perm_[0] == 0 )
        cost *= std::pow(1.01, achievedParallelismAtLoop[0] - 1); //strongly penalize this case
  
     cost *= std::pow(1.00010, std::min(16,achievedParallelismAtLoop[0] - 1));        // if at all, prefer ...
     cost *= std::pow(1.00015, std::min(16,achievedParallelismAtLoop[perm_[0]] - 1)); // parallelization in stride-1 of A
     
  
     const int workPerThread = (availableParallelismAtLoop[perm_[0]] + achievedParallelismAtLoop[perm_[0]] -1) / achievedParallelismAtLoop[perm_[0]];
     if( workPerThread * sizeof(floatType) % 64 != 0 && achievedParallelismAtLoop[perm_[0]] > 1 ){ //avoid false-sharing
        cost *= std::pow(1.00015, std::min(16,achievedParallelismAtLoop[perm_[0]] - 1)); // penalize this parallelization again
     }
     return cost;
  }
  
  template<typename floatType>
  void Transpose<floatType>::getParallelismStrategies(std::vector<std::vector<int> > &parallelismStrategies) const
  {
     parallelismStrategies.clear();
     if( numThreads_ == 1 ){
        parallelismStrategies.emplace_back(std::vector<int>(dim_, 1));
        return;
     }
     std::vector<int> bestParallelismStrategy(dim_, 1);
     getBestParallelismStrategy(bestParallelismStrategy);
     if( this->infoLevel_ > 0 )
        printf("Loadbalancing: %f\n", getLoadBalance(bestParallelismStrategy));
  
     if( selectionMethod_ == ESTIMATE )
     {
        parallelismStrategies.push_back(bestParallelismStrategy);
        return;
     }
  
     // ATTENTION: we don't care about the case where numThreads_ is a large prime number...
     // (sorry, KNC)
     //
     // we factorize numThreads into its prime factors because we have to match
     // every one to a certain loop. In principle every loop could be used to
     // match every primefactor, but some choices are preferable over others.
     // E.g., we want to achieve good load-balancing _and_ try to avoid the
     // stride-1 index of B (due to false sharing)
     std::list<int> primeFactors;
     getPrimeFactors( numThreads_, primeFactors );
     if( this->infoLevel_ > 0 )
        printVector(primeFactors,"primes");
  
     std::vector<int> availableParallelismAtLoop;
     this->getAvailableParallelism( availableParallelismAtLoop);
     if( this->infoLevel_ > 0 )
        printVector(availableParallelismAtLoop,"available Parallelism");
  
     std::vector<int> achievedParallelismAtLoop (dim_, 1);
  
     this->getAllParallelismStrategies( primeFactors, 
           availableParallelismAtLoop, 
           achievedParallelismAtLoop, 
           parallelismStrategies);
     
     // sort according to loop heuristic
     std::sort(parallelismStrategies.begin(), parallelismStrategies.end(), 
           [this](const std::vector<int> loopOrder1, const std::vector<int> loopOrder2)
           { 
              return this->parallelismCostHeuristic(loopOrder1) < this->parallelismCostHeuristic(loopOrder2); 
           });
  
     parallelismStrategies.insert(parallelismStrategies.begin(), bestParallelismStrategy);
  
     if( this->infoLevel_ > 1 )
        for( auto strat : parallelismStrategies ){
           printVector(strat,"parallelization");
           printf("cost: %f\n", this->parallelismCostHeuristic( strat ));
        }
  }
  
  template<typename floatType>
  void Transpose<floatType>::verifyParameter(const int *size, const int* perm, const int* outerSizeA, const int* outerSizeB, const int dim) const
  {
     if ( dim < 1 ) {
        fprintf(stderr,"[HPTT] ERROR: dimensionality too low.\n");
        exit(-1);
     }
  
     std::vector<int> found(dim, 0);
  
     for(int i=0;i < dim ; ++i)
     {
        if( size[i] <= 0 ) {
           fprintf(stderr,"[HPTT] ERROR: size at position %d is invalid\n", i);
           exit(-1);
        }
        found[ perm[i] ] = 1;
     }
  
     for(int i=0;i < dim ; ++i)
        if( found[i] <= 0 ) {
           fprintf(stderr,"[HPTT] ERROR: permutation invalid\n");
           exit(-1);
        }
  
     if ( outerSizeA != NULL )
        for(int i=0;i < dim ; ++i)
           if ( outerSizeA[i] < size[i] ) {
              fprintf(stderr,"[HPTT] ERROR: outerSizeA invalid\n");
              exit(-1);
           }
  
     if ( outerSizeB != NULL )
        for(int i=0;i < dim ; ++i)
           if ( outerSizeB[i] < size[perm[i]] ) {
              fprintf(stderr,"[HPTT] ERROR: outerSizeB invalid\n");
              exit(-1);
           }
  }
  
  template<typename floatType>
  void Transpose<floatType>::computeLeadingDimensions()
  {
     lda_[0] = 1;
     if( outerSizeA_[0] == -1 )
        for(int i=1;i < dim_ ; ++i)
           lda_[i] = lda_[i-1] * sizeA_[i-1];
     else
        for(int i=1;i < dim_ ; ++i)
           lda_[i] = outerSizeA_[i-1] * lda_[i-1];
  
     ldb_[0] = 1;
     if( outerSizeB_[0] == -1 )
        for(int i=1;i < dim_ ; ++i)
           ldb_[i] = ldb_[i-1] * sizeA_[perm_[i-1]];
     else
        for(int i=1;i < dim_ ; ++i)
           ldb_[i] = outerSizeB_[i-1] * ldb_[i-1];
  }
  
  template<typename floatType>
  void Transpose<floatType>::skipIndices(const int *sizeA, const int* perm, const int *outerSizeA, const int *outerSizeB, const int dim)
  {
     for(int i=0;i < dim ; ++i)
     {
        perm_[i] = perm[i];
        sizeA_[i] = sizeA[i];
        if(outerSizeA)
           outerSizeA_[i] = outerSizeA[i];
        else
           outerSizeA_[i] = sizeA[i];
        if(outerSizeB)
           outerSizeB_[i] = outerSizeB[i];
        else
           outerSizeB_[i] = sizeA[perm[i]];
     }
  
     int skipped = 0;
     for(int i=0;i < dim ; ++i)
     {
        int idxB = 0;
        for(;idxB < dim ; ++idxB)
           if(perm[idxB] == i)
              break;
        if( sizeA[i] == 1 && (!outerSizeA || outerSizeA[i] == 1) && (!outerSizeB || outerSizeB[idxB] == 1) )
        {
           sizeA_[i] = -1;
           outerSizeA_[i] = -1;
           outerSizeB_[idxB] = -1;
           perm_[idxB] = -1;
           skipped++;
        }
     }
     // compact arrays (remove -1)
     for(int i=0;i < dim ; ++i)
        if( sizeA_[i] == -1 )
        {
           int j=i+1;
           for(;j < dim ; ++j)
              if( sizeA_[j] != -1 )
                 break;
           if( j < dim )
              std::swap(sizeA_[i], sizeA_[j]);
        }
     for(int i=0;i < dim ; ++i)
        if( outerSizeA_[i] == -1 )
        {
           int j=i+1;
           for(;j < dim ; ++j)
              if( outerSizeA_[j] != -1 )
                 break;
           if( j < dim )
              std::swap(outerSizeA_[i], outerSizeA_[j]);
        }
     for(int i=0;i < dim ; ++i)
        if( outerSizeB_[i] == -1 )
        {
           int j=i+1;
           for(;j < dim ; ++j)
              if( outerSizeB_[j] != -1 )
                 break;
           if( j < dim )
              std::swap(outerSizeB_[i], outerSizeB_[j]);
        }
     for(int i=0;i < dim ; ++i)
        if( perm_[i] == -1 )
        {
           int j=i+1;
           for(;j < dim ; ++j)
              if( perm_[j] != -1 )
                 break;
           if( j < dim )
              std::swap(perm_[i], perm_[j]);
        }
  
     dim_ = dim - skipped;
     if( dim_ == 0 ){
        dim_ = 1;
        perm_.resize(dim_);
        sizeA_.resize(dim_);
        outerSizeA_.resize(dim_);
        outerSizeB_.resize(dim_);
        perm_[0] = 0;
        sizeA_[0] = 1;
        outerSizeA_[0] = 1;
        outerSizeB_[0] = 1;
     } else {
        perm_.resize(dim_);
        sizeA_.resize(dim_);
        outerSizeA_.resize(dim_);
        outerSizeB_.resize(dim_);
  
        // remove gaps in the perm, if requried (e.g., perm=3,1,0 -> 2,1,0)
        int currentValue = 0;
        for(int i=0;i < dim_; ++i){
           //find smallest element in perm_ and rename it to currentValue
           int minValue = 1000000;
           int minPos = -1;
           for(int pos=0; pos < dim_; ++pos){
              if ( perm_[pos] >= currentValue && perm_[pos] < minValue) {
                 minValue = perm_[pos];
                 minPos = pos;
              }
           }
           perm_[minPos] = currentValue; // minValue renamed to currentValue
           currentValue++;
        }
     }
  
  #ifdef DEBUG
     printVector(perm_,"perm");
     printVector(sizeA_,"sizeA");
     printVector(outerSizeA_,"outerSizeA");
     printVector(outerSizeB_,"outerSizeB");
  #endif
  }
  
  /**
    * \brief fuses indices whenever possible 
    * \detailed For instance:
    *           perm=3,1,2,0 & size=10,11,12,13  becomes: perm=2,1,0 & size=10,11*12,13
    * \return This function will initialize sizeA_, perm_, outerSizeA_, outersize_ and dim_
  */
  template<typename floatType>
  void Transpose<floatType>::fuseIndices()
  {
     std::list< std::tuple<int, int> > fusedIndices;
  
     std::vector<int> perm;
     // correct perm
     for(int i=0;i < dim_; ++i){
        // merge indices if the two consecutive entries are identical
        int toMerge = i;
        perm.push_back(perm_[i]);
        while(i+1 < dim_ && perm_[i] + 1 == perm_[i+1] 
              && (sizeA_[perm_[i]] == outerSizeA_[perm_[i]]) 
              && (sizeA_[perm_[i]] == outerSizeB_[i]) ){ 
  #ifdef DEBUG
           fprintf(stderr,"[HPTT] MERGING indices %d and %d\n",perm_[i], perm_[i+1]); 
  #endif
           fusedIndices.emplace_back( std::make_tuple(perm_[toMerge],perm_[i+1]) );
           i++;
        }
     }
  
     // correct sizes and outer-sizes
     for( auto tup : fusedIndices )
     {
        sizeA_[std::get<0>(tup)] *= sizeA_[std::get<1>(tup)];
        outerSizeA_[std::get<0>(tup)] *= outerSizeA_[std::get<1>(tup)];
        outerSizeA_[std::get<1>(tup)] = -1;
  
        auto pos1 = std::find(perm_.begin(), perm_.end(), std::get<0>(tup)) - perm_.begin();
        auto pos2 = std::find(perm_.begin(), perm_.end(), std::get<1>(tup)) - perm_.begin();
        outerSizeB_[pos1] *= outerSizeB_[pos2];
        outerSizeB_[pos2] = -1;
     }
  
     if ( fusedIndices.size() > 0 ) {
        perm_ = perm;
        // remove gaps in the perm, if requried (e.g., perm=3,1,0 -> 2,1,0)
        int currentValue = 0;
        for(int i=0;i < perm_.size(); ++i){
           //find smallest element in perm_ and rename it to currentValue
           int minValue = 1000000;
           int minPos = -1;
           for(int pos=0; pos < perm_.size(); ++pos){
              if ( perm_[pos] >= currentValue && perm_[pos] < minValue) {
                 minValue = perm_[pos];
                 minPos = pos;
              }
           }
  #ifdef DEBUG
           printf("perm[%d]: %d -> %d\n",minPos, perm_[minPos], currentValue);
  #endif
           perm_[minPos] = currentValue; // minValue renamed to currentValue
           sizeA_[currentValue] = sizeA_[minValue];
           currentValue++;
        }
  
  
        // compact outer size (e.g.: outerSizeA_[] = {24,-1,5,-1,13} -> {24,5,13,-1,-1} -> {24,5,13}
        for(int i=0;i < dim_ ; ++i)
           if( outerSizeA_[i] == -1 )
           {
              int j=i+1;
              for(;j < dim_ ; ++j)
                 if( outerSizeA_[j] != -1 )
                    break;
              if( j < dim_ )
                 std::swap(outerSizeA_[i], outerSizeA_[j]);
           }
        for(int i=0;i < dim_ ; ++i)
           if( outerSizeB_[i] == -1 )
           {
              int j=i+1;
              for(;j < dim_ ; ++j)
                 if( outerSizeB_[j] != -1 )
                    break;
              if( j < dim_ )
                 std::swap(outerSizeB_[i], outerSizeB_[j]);
           }
        dim_ -= fusedIndices.size();
        outerSizeA_.resize(dim_);
        outerSizeB_.resize(dim_);
        sizeA_.resize(dim_);
        perm_.resize(dim_);
  
  #ifdef DEBUG
        printf("\nperm_new: ");
        for(int i=0;i < dim_ ; ++i)
           printf("%d ",perm_[i]);
        printf("\nsizes_new: ");
        for(int i=0;i < dim_ ; ++i)
           printf("%d ",sizeA_[i]);
        printf("\nouterSizeA_new: ");
        for(int i=0;i < dim_ ; ++i)
           printf("%d ",outerSizeA_[i]);
        printf("\nouterSizeB_new: ");
        for(int i=0;i < dim_ ; ++i)
           printf("%d ",outerSizeB_[i]);
        printf("\n");
  #endif
     } 
  }
  
  // returns the best loop order (same as the best one with exhaustive search)
  template<typename floatType>
  void Transpose<floatType>::getBestLoopOrder( std::vector<int> &loopOrder ) const
  {
     auto totalOuterSizeA = std::accumulate(outerSizeA_.begin(), outerSizeA_.end(), 1, std::multiplies<size_t>()) * sizeof(floatType);
     auto totalOuterSizeB = std::accumulate(outerSizeB_.begin(), outerSizeB_.end(), 1, std::multiplies<size_t>()) * sizeof(floatType);
     if( totalOuterSizeA > totalOuterSizeB && totalOuterSizeB <= 22 * 1024. * 1024. ) // B is likely to fit into L3 cache
     {
        // prefer accesses to A over those to B (Rationale: reduce TLB misses)
        for(int i=0;i < dim_ ; ++i) 
           loopOrder[dim_-1-i] = i; //innermost loop idx is stored at dim_-1
        return;
     }else if( totalOuterSizeB > totalOuterSizeA && totalOuterSizeA <= 22 * 1024. * 1024. ) // B is likely to fit into L3 cache
     {
        // prefer accesses to B over those to A (Rationale: reduce TLB misses)
        for(int i=0;i < dim_ ; ++i) 
           loopOrder[dim_-1-i] = dim_-1-i; //innermost loop idx is stored at dim_-1
        return;
     }
  
     // create cost matrix; cost[i,idx] === cost for idx being at loop-level i
     double costs[dim_*dim_];
     for(int i=0;i < dim_ ; ++i){
        for(int idx=0;idx < dim_ ; ++idx){ //idx is at loop i
           double cost = 0;
           if( i != 0 ){
              const int posB = findPos(idx, perm_);
              const int importanceA = (1<<(dim_ - idx)); // stride-1 has the most importance ...
              const int importanceB = (1<<(dim_ - posB));                // subsequent indices are half as important
              const int penalty = 10 * (1<<(i-1));
              constexpr double bias = 1.01;
              cost = (importanceA + importanceB * bias) * penalty;
           }
           costs[i + idx * dim_] = cost;
        }
     }
     std::list<int> availLoopLevels; //available rows
     std::list<int> availIndices;
     for(int i=0;i < dim_ ; ++i){
        availLoopLevels.push_back(i);
        availIndices.push_back(i);
     }
  
     // create best loop order constructively without generating all
     for(int i=0;i < dim_ ; ++i){
        //find column with maximum cost
        int selectedIdx= 0;
        double maxValueAll = 0;
        for(auto c : availIndices)
        {
           double maxValue = 0;
           for(auto r : availLoopLevels){
              const double val = costs[c * dim_ + r];
              maxValue = (val > maxValue) ? val : maxValue;
           }
  
           if(maxValue > maxValueAll){
              maxValueAll = maxValue;
              selectedIdx= c;
           }
        }
        //find minimum in that column
        int selectedLoopLevel = 0;
        double minValue = 1e100;
        for(auto r : availLoopLevels){
           const double val = costs[selectedIdx * dim_ + r];
           if(val < minValue){
              minValue = val;
              selectedLoopLevel = r;
           }
        }
        // update loop order
        loopOrder[dim_-1-i] = selectedIdx; //innermost loop idx is stored at dim_-1
        // remove selected row
        for( auto it = availLoopLevels.begin(); it != availLoopLevels.end(); it++ )
           if(*it == selectedLoopLevel){
              availLoopLevels.erase(it);
              break;
           }
        // remove selected col
        for( auto it = availIndices.begin(); it != availIndices.end(); it++ )
           if(*it == selectedIdx){
              availIndices.erase(it);
              break;
           }
     }
  }
  
  template<typename floatType>
  double Transpose<floatType>::loopCostHeuristic( const std::vector<int> &loopOrder ) const
  {
     double loopCost= 0.0;
     for(int i=1;i < dim_ ; ++i){
        const int idx = loopOrder[dim_-1-i];
        const int posB = findPos(idx, perm_);
        const int importanceA = (1<<(dim_ - idx)); // stride-1 has the most importance ...
        const int importanceB = (1<<(dim_ - posB));                // subsequent indices are half as important
        const int penalty = 10 * (1<<(i-1));
        double bias = 1.01;
        loopCost += (importanceA + importanceB * bias) * penalty;
     } 
  
     return loopCost;
  }
  
  template<typename floatType>
  void Transpose<floatType>::getLoopOrders(std::vector<std::vector<int> > &loopOrders) const
  {
     loopOrders.clear();
     if( selectionMethod_ == ESTIMATE ){
        loopOrders.emplace_back(std::vector<int>(dim_));
        getBestLoopOrder( loopOrders[0] );
        return;
     }
  
     std::vector<int> loopOrder;
     for(int i = 0; i < dim_; i++)
        loopOrder.push_back(i);
  
     // create all loopOrders
     do {
        if ( perm_[0] == 0 && loopOrder[dim_-1] != 0 )
           continue; // ATTENTION: we skip all loop-orders where the stride-1 index is not the inner-most loop iff perm[0] == 0 (both for perf & correctness)
  
        loopOrders.push_back(loopOrder);
     } while(std::next_permutation(loopOrder.begin(), loopOrder.end()));
  
  
     // sort according to loop heuristic
     std::sort(loopOrders.begin(), loopOrders.end(), 
           [this](const std::vector<int> loopOrder1, const std::vector<int> loopOrder2)
           { 
              return this->loopCostHeuristic(loopOrder1) < this->loopCostHeuristic(loopOrder2); 
           });
  
     if( this->infoLevel_ > 1 )
        for(auto loopOrder : loopOrders) {
           printVector(loopOrder,"loop");
           printf("penalty: %f\n",loopCostHeuristic(loopOrder));
        }
  }
  
  template<typename floatType>
  void Transpose<floatType>::createPlan()
  {
  //   printf("entering createPlan()\n");
  #ifdef HPTT_TIMERS
     double timeStart = omp_get_wtime();
  #endif
     
     std::vector<std::shared_ptr<Plan> > allPlans;
     createPlans(allPlans);
  
  #ifdef HPTT_TIMERS
     printf("createPlans() took %f ms\n",(omp_get_wtime()-timeStart)*1000);
     timeStart = omp_get_wtime();
  #endif
     masterPlan_ = selectPlan( allPlans );
     if( this->infoLevel_ > 0 ){
        printf("Configuration of best plan:\n");
        masterPlan_->print();
     }
  #ifdef HPTT_TIMERS
     printf("SelectPlan() took %f ms\n",(omp_get_wtime()-timeStart)*1000);
  #endif
  }
  
  template<typename floatType>
  void Transpose<floatType>::createPlans( std::vector<std::shared_ptr<Plan> > &plans ) const
  {
     if( dim_ == 1 || (dim_ == 2 && perm_[0] == 0)) {
        plans.emplace_back(new Plan); //create dummy plan
        return; // handled within execute()
     }
  #ifdef HPTT_TIMERS
     double parallelStrategiesTime = omp_get_wtime();
  #endif
     std::vector<std::vector<int> > parallelismStrategies;
     this->getParallelismStrategies(parallelismStrategies);
  #ifdef HPTT_TIMERS
     printf("There exists %d parallel strategies. Time: %f ms\n",parallelismStrategies.size(), (omp_get_wtime()-parallelStrategiesTime)*1000);
  
     double loopOrdersTime = omp_get_wtime();
  #endif
     std::vector<std::vector<int> > loopOrders;
     this->getLoopOrders(loopOrders);
  #ifdef HPTT_TIMERS
     printf("There exists %d loop orders. Time: %f ms\n",loopOrders.size(), (omp_get_wtime()-loopOrdersTime)*1000);
  #endif
  
     if( selectedParallelStrategyId_ != -1 ){
        int selectedParallelStrategyId = std::min((int)parallelismStrategies.size()-1, selectedParallelStrategyId_);
        std::vector<int>parStrategy(parallelismStrategies[selectedParallelStrategyId]);
        printVector(parStrategy,"selected parallel: ");
        parallelismStrategies.clear();
        parallelismStrategies.push_back(parStrategy);
     }
     if( selectedLoopOrderId_ != -1 ){
        int selectedLoopOrderId = std::min((int)loopOrders.size()-1, selectedLoopOrderId_);
        std::vector<int>loopOrder(loopOrders[selectedLoopOrderId]);
        printVector(loopOrder,"selected loopOrder: ");
        loopOrders.clear();
        loopOrders.push_back(loopOrder);
     }
  
     const int posStride1A_inB = findPos(0, perm_);
     const int posStride1B_inA = perm_[0];
  
     // combine the loopOrder and parallelismStrategies according to their
     // heuristics, search the space with a growing rectangle (from best to worst,
     // see line marked with ***)
     bool done = false;
     for( int start= 0; start< std::max( parallelismStrategies.size(), loopOrders.size() ) && !done; start++ )
        for( int i = 0; i < parallelismStrategies.size() && !done; i++)
        {
           for( int j = 0; j < loopOrders.size() && !done; j++)
           {
              if( i > start || j > start || (i != start && j != start) ) continue; //these are already done ***
  
              auto numThreadsAtLoop = parallelismStrategies[i];
              auto loopOrder = loopOrders[j];
              auto plan = std::make_shared<Plan>(loopOrder, numThreadsAtLoop );
              const int numTasks = plan->getNumTasks();
  
  #ifdef _OPENMP
  #pragma omp parallel for num_threads(numThreads_) if(numThreads_ > 1)
  #endif
              for( int taskId = 0; taskId < numTasks; taskId++)
              {
                 ComputeNode *currentNode = plan->getRootNode(taskId);
  
                 int numThreadsPerComm = numTasks; //global communicator // e.g., 6
                 int taskIdComm = taskId; // e.g., 0,1,2,3,4,5
                 // divide each loop-level l, corresponding to index loopOrder[l], into numThreadsAtLoop[index] chunks
                 for(int l=0; l < dim_; ++l){
                    const int index = loopOrder[l];
                    currentNode->inc = this->getIncrement( index );
  
                    const int numTasksAtLevel = numThreadsAtLoop[index]; //  e.g., 3
                    const int numParallelismAvailable = (sizeA_[index] + currentNode->inc - 1) / currentNode->inc; // e.g., 5
                    const int workPerThread = (numParallelismAvailable + numTasksAtLevel -1) / numTasksAtLevel; // ceil(5/3) = 2
  
                    numThreadsPerComm /= numTasksAtLevel; //numThreads in next communicator // 6/3 = 2
                    const int commId = (taskIdComm/numThreadsPerComm); //  = 0,0,1,1,2,2
                    taskIdComm = taskIdComm % numThreadsPerComm; // local taskId in next communicator // 0,1,0,1,0,1
  
                    currentNode->start = std::min( sizeA_[index], commId * workPerThread * currentNode->inc );
                    currentNode->end = std::min( sizeA_[index], (commId+1) * workPerThread * currentNode->inc );
  
                    currentNode->lda = lda_[index];
                    currentNode->ldb = ldb_[findPos(index, perm_)];
  
                    if( perm_[0] != 0 || l != dim_-1 ){
                       currentNode->next = new ComputeNode;
                       currentNode = currentNode->next;
                    }
                 }
  
                 //macro-kernel
                 if( perm_[0] != 0 )
                 {
                    currentNode->start = -1;
                    currentNode->end = -1;
                    currentNode->inc = -1;
                    currentNode->lda = lda_[ posStride1B_inA ];
                    currentNode->ldb = ldb_[ posStride1A_inB ];
                    currentNode->next = nullptr;
                 }
              }
              plans.push_back(plan);
              if( selectionMethod_ == ESTIMATE || 
                  selectionMethod_ == MEASURE && plans.size() > 200 || 
                  selectionMethod_ == PATIENT && plans.size() > 400 || 
                  selectionMethod_ == CRAZY && plans.size() > 800 )
                 done = true;
           }
        }
  }
  
  /**
   * Estimates the time in seconds for the given computeTree
   */
  template<typename floatType>
  float Transpose<floatType>::estimateExecutionTime( const std::shared_ptr<Plan> plan)
  {
     auto startTime = std::chrono::high_resolution_clock::now();
     this->executeEstimate(plan.get());
     double elapsedTime = std::chrono::duration_cast<std::chrono::duration<double, std::milli> >
        (std::chrono::high_resolution_clock::now() - startTime).count();
  
     const double minMeasurementTime = 0.1; // in seconds
  
     // do at least 3 repetitions or spent at least 'minMeasurementTime' seconds for each candidate
     int nRepeat = std::min(3, (int) std::ceil(minMeasurementTime / elapsedTime));
  
     //execute just a few iterations and exterpolate the result
     startTime = std::chrono::high_resolution_clock::now();
     for(int i=0;i < nRepeat ; ++i) //ATTENTION: we are not clearing the caches inbetween runs
        this->executeEstimate( plan.get() );
     elapsedTime = std::chrono::duration_cast<std::chrono::duration<double, std::milli> >
        (std::chrono::high_resolution_clock::now() - startTime).count();
     elapsedTime /= nRepeat;
  
  #ifdef DEBUG
     printf("Estimated time: %.3e ms.\n",elapsedTime); 
  #endif
     return elapsedTime; 
  }
  
  template<typename floatType>
  double Transpose<floatType>::getTimeLimit() const
  {
     if( selectionMethod_ == ESTIMATE )
        return 0.0;
     else if( selectionMethod_ == MEASURE)
        return 10.;   // 10s
     else if( selectionMethod_ == PATIENT)
        return 60.;   // 1m
     else if( selectionMethod_ == CRAZY )
        return 3600.; // 1h
     else{
        fprintf(stderr,"[HPTT] ERROR: selectionMethod unknown.\n");
        exit(-1);
     }
     return -1;
  }
  
  template<typename floatType>
  std::shared_ptr<Plan> Transpose<floatType>::selectPlan( const std::vector<std::shared_ptr<Plan> > &plans)
  {
     if( plans.size() <= 0 ){
        fprintf(stderr,"[HPTT] Internal error: not enough plans generated.\n");
        exit(-1);
     }
     if( selectionMethod_ == ESTIMATE ) // fast return
        return plans[0];
  
     double timeLimit = this->getTimeLimit() * 1000; //in ms
     int maxAutotuningCandidates = plans.size();
     if( maxAutotuningCandidates_ != -1 ){
        maxAutotuningCandidates = maxAutotuningCandidates_;
        timeLimit = 1e9;
     }
  
     float minTime = FLT_MAX;
     int bestPlan_id = 0;
  
     if( plans.size() > 1 )
     {
        int plansEvaluated = 0;
        auto startTime = std::chrono::high_resolution_clock::now();
        for( int plan_id = 0; plan_id < maxAutotuningCandidates; plan_id++ )
        {
           auto p = plans[plan_id];
  
           double elapsedTime = std::chrono::duration_cast<std::chrono::duration<double, std::milli> >
              (std::chrono::high_resolution_clock::now() - startTime).count();
           if( elapsedTime >= timeLimit ) // timelimit reached
              break;
  
           float estimatedTime = this->estimateExecutionTime( p );
           plansEvaluated++;
  
           if( estimatedTime < minTime ){
              bestPlan_id = plan_id;
              minTime = estimatedTime;
           }
           if( this->infoLevel_ > 1 ){
              printf("Plan %d will take roughly %f ms.\n", plan_id, estimatedTime * 1000.);
              plans[plan_id]->print();
           }
        }
        if( this->infoLevel_ > 0 )
           printf("We evaluated %d/%zu candidates and selected candidate %d.\n", plansEvaluated, plans.size(), bestPlan_id); 
     }
     return plans[bestPlan_id];
  }
  
  template class Transpose<float>;
  template class Transpose<double>;
  template class Transpose<FloatComplex>;
  template class Transpose<DoubleComplex>;
  
  template void Transpose<float>::execute_expert<true, true, true>();
  template void Transpose<float>::execute_expert<true, false, true>();
  template void Transpose<float>::execute_expert<false, true, true>();
  template void Transpose<float>::execute_expert<false, false, true>();
  template void Transpose<float>::execute_expert<true, true, false>();
  template void Transpose<float>::execute_expert<true, false, false>();
  template void Transpose<float>::execute_expert<false, true, false>();
  template void Transpose<float>::execute_expert<false, false, false>();
  
  template void Transpose<double>::execute_expert<true,true, true>();
  template void Transpose<double>::execute_expert<false,true, true>();
  template void Transpose<double>::execute_expert<true,false, true>();
  template void Transpose<double>::execute_expert<false,false, true>();
  template void Transpose<double>::execute_expert<true,true, false>();
  template void Transpose<double>::execute_expert<false,true, false>();
  template void Transpose<double>::execute_expert<true,false, false>();
  template void Transpose<double>::execute_expert<false,false, false>();
  
  template void Transpose<FloatComplex>::execute_expert<true,true, true>();
  template void Transpose<FloatComplex>::execute_expert<false,true, true>();
  template void Transpose<FloatComplex>::execute_expert<true,false, true>();
  template void Transpose<FloatComplex>::execute_expert<false,false, true>();
  template void Transpose<FloatComplex>::execute_expert<true,true, false>();
  template void Transpose<FloatComplex>::execute_expert<false,true, false>();
  template void Transpose<FloatComplex>::execute_expert<true,false, false>();
  template void Transpose<FloatComplex>::execute_expert<false,false, false>();
  
  template void Transpose<DoubleComplex>::execute_expert<true,true, true>();
  template void Transpose<DoubleComplex>::execute_expert<false,true, true>();
  template void Transpose<DoubleComplex>::execute_expert<true,false, true>();
  template void Transpose<DoubleComplex>::execute_expert<false,false, true>();
  template void Transpose<DoubleComplex>::execute_expert<true,true, false>();
  template void Transpose<DoubleComplex>::execute_expert<false,true, false>();
  template void Transpose<DoubleComplex>::execute_expert<true,false, false>();
  template void Transpose<DoubleComplex>::execute_expert<false,false, false>();
  
  
  }
  