// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2020 Everton Constantino (everton.constantino@ibm.com)
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef EIGEN_MATRIX_PRODUCT_ALTIVEC_H
#define EIGEN_MATRIX_PRODUCT_ALTIVEC_H

#include "MatrixProductCommon.h"

#if __GNUC__ > 10 || \
    (__GNUC__ == 10 && (__GNUC_MINOR__ > 2 || \
                       (__GNUC_MINOR__ == 2 && \
                        __GNUC_PATCHLEVEL__ >= 1)))
  #define ALTIVEC_MMA_SUPPORT
#endif

#if defined(ALTIVEC_MMA_SUPPORT) && !defined(EIGEN_ALTIVEC_DISABLE_MMA)
  #include "MatrixProductMMA.h"
#endif

/**************************************************************************************************
 * TODO                                                                                           *
 * - Check StorageOrder on lhs_pack (the innermost second loop seems unvectorized when it could). *
 * - Check the possibility of transposing as GETREAL and GETIMAG when needed.                     *
 * - Check if change conjugation to xor instead of mul gains any performance.                     *
 * - Remove IsComplex template argument from complex packing.                                     *
 **************************************************************************************************/
namespace Eigen {

namespace internal {

/**************************
 * Constants and typedefs *
 **************************/
const int QuadRegisterCount = 8;

template<typename Scalar>
struct quad_traits
{
  typedef typename packet_traits<Scalar>::type    vectortype;
  typedef PacketBlock<vectortype, 4>                    type;
  typedef vectortype                                 rhstype;
  enum
  {
    vectorsize = packet_traits<Scalar>::size,
    size = 4,
    rows = 4
  };
};

template<>
struct quad_traits<double>
{
  typedef Packet2d                        vectortype;
  typedef PacketBlock<vectortype, 4>            type;
  typedef PacketBlock<Packet2d,2>            rhstype;
  enum
  {
    vectorsize = packet_traits<double>::size,
    size = 2,
    rows = 4
  };
};

// MatrixProduct decomposes real/imaginary vectors into a real vector and an imaginary vector, this turned out
// to be faster than Eigen's usual approach of having real/imaginary pairs on a single vector. This constants then
// are responsible to extract from convert between Eigen's and MatrixProduct approach.
const static Packet4f p4f_CONJUGATE = {-1.0f, -1.0f, -1.0f, -1.0f};

const static Packet2d p2d_CONJUGATE = {-1.0, -1.0};

const static Packet16uc p16uc_GETREAL32 = {  0,  1,  2,  3,
                                             8,  9, 10, 11,
                                            16, 17, 18, 19,
                                            24, 25, 26, 27};

const static Packet16uc p16uc_GETIMAG32 = {  4,  5,  6,  7,
                                            12, 13, 14, 15,
                                            20, 21, 22, 23,
                                            28, 29, 30, 31};
const static Packet16uc p16uc_GETREAL64 = {  0,  1,  2,  3,  4,  5,  6,  7,
                                            16, 17, 18, 19, 20, 21, 22, 23};

//[a,ai],[b,bi] = [ai,bi]
const static Packet16uc p16uc_GETIMAG64 = {  8,  9, 10, 11, 12, 13, 14, 15,
                                            24, 25, 26, 27, 28, 29, 30, 31};

/*********************************************
 * Single precision real and complex packing *
 * *******************************************/

/**
 * Symm packing is related to packing of symmetric adjoint blocks, as expected the packing leaves
 * the diagonal real, whatever is below it is copied from the respective upper diagonal element and 
 * conjugated. There's no PanelMode available for symm packing.
 *
 * Packing in general is supposed to leave the lhs block and the rhs block easy to be read by gemm using 
 * it's respective rank-update instructions. The float32/64 versions are different because at this moment
 * the size of the accumulator is fixed at 512-bits so you can't have a 4x4 accumulator of 64-bit elements.
 * 
 * As mentioned earlier MatrixProduct breaks complex numbers into a real vector and a complex vector so packing has
 * to take that into account, at the moment, we run pack the real part and then the imaginary part, this is the main
 * reason why packing for complex is broken down into several different parts, also the reason why we endup having a
 * float32/64 and complex float32/64 version.
 **/
template<typename Scalar, typename Index, int StorageOrder>
EIGEN_STRONG_INLINE std::complex<Scalar> getAdjointVal(Index i, Index j, const_blas_data_mapper<std::complex<Scalar>, Index, StorageOrder>& dt)
{
  std::complex<Scalar> v;
  if(i < j)
  {
    v.real(dt(j,i).real());
    v.imag(-dt(j,i).imag());
  } else if(i > j)
  {
    v.real(dt(i,j).real());
    v.imag(dt(i,j).imag());
  } else {
    v.real(dt(i,j).real());
    v.imag((Scalar)0.0f);
  }
  return v;
}

template<typename Scalar, typename Index, int StorageOrder, int N>
EIGEN_STRONG_INLINE void symm_pack_complex_rhs_helper(std::complex<Scalar> *blockB, const std::complex<Scalar>* _rhs, Index rhsStride, Index rows, Index cols, Index k2)
{
  const Index depth = k2 + rows;
  const_blas_data_mapper<std::complex<Scalar>, Index, StorageOrder> rhs(_rhs, rhsStride);
  const int vectorSize = N*quad_traits<Scalar>::vectorsize;
  Scalar* blockBf = reinterpret_cast<Scalar *>(blockB);

  Index ri = 0, j = 0;
  for(; j + vectorSize < cols; j+=vectorSize)
  {
      Index i = k2;
      for(; i < depth; i++)
      {
        for(Index k = 0; k < vectorSize; k++)
        {
          std::complex<Scalar> v = getAdjointVal<Scalar, Index, StorageOrder>(i, j + k, rhs);
          blockBf[ri + k] = v.real();
        }
        ri += vectorSize;
      }

      i = k2;

      for(; i < depth; i++)
      {
        for(Index k = 0; k < vectorSize; k++)
        {
          std::complex<Scalar> v = getAdjointVal<Scalar, Index, StorageOrder>(i, j + k, rhs);
          blockBf[ri + k] = v.imag();
        }
        ri += vectorSize;
      }
  }
  for(Index i = k2; i < depth; i++)
  {
      Index k = j;
      for(; k < cols; k++)
      {
        std::complex<Scalar> v = getAdjointVal<Scalar, Index, StorageOrder>(i, k, rhs);
        blockBf[ri] = v.real();
        ri += 1;
      }
  }
  for(Index i = k2; i < depth; i++)
  {
      Index k = j;
      for(; k < cols; k++)
      {
        std::complex<Scalar> v = getAdjointVal<Scalar, Index, StorageOrder>(i, k, rhs);
        blockBf[ri] = v.imag();
        ri += 1;
      }
  }
}

template<typename Scalar, typename Index, int StorageOrder>
EIGEN_STRONG_INLINE void symm_pack_complex_lhs_helper(std::complex<Scalar> *blockA, const std::complex<Scalar>* _lhs, Index lhsStride, Index cols, Index rows)
{
  const Index depth = cols;
  const_blas_data_mapper<std::complex<Scalar>, Index, StorageOrder> lhs(_lhs, lhsStride);
  const int vectorSize = quad_traits<Scalar>::vectorsize;
  Index ri = 0, j = 0;
  Scalar *blockAf = (Scalar *)(blockA);

  for(; j + vectorSize < rows; j+=vectorSize)
  {
      Index i = 0;

      for(; i < depth; i++)
      {
        for(int k = 0; k < vectorSize; k++)
        {
          std::complex<Scalar> v = getAdjointVal<Scalar, Index, StorageOrder>(j+k, i, lhs);
          blockAf[ri + k] = v.real();
        }
        ri += vectorSize;
      }
      i = 0;
      for(; i < depth; i++)
      {
        for(int k = 0; k < vectorSize; k++)
        {
          std::complex<Scalar> v = getAdjointVal<Scalar, Index, StorageOrder>(j+k, i, lhs);
          blockAf[ri + k] = v.imag();
        }
        ri += vectorSize;
      }
  }

  for(Index i = 0; i < depth; i++)
  {
      Index k = j;
      for(; k < rows; k++)
      {
          std::complex<Scalar> v = getAdjointVal<Scalar, Index, StorageOrder>(k, i, lhs);
          blockAf[ri] = v.real();
          ri += 1;
      }
  }
  for(Index i = 0; i < depth; i++)
  {
      Index k = j;
      for(; k < rows; k++)
      {
          std::complex<Scalar> v = getAdjointVal<Scalar, Index, StorageOrder>(k, i, lhs);
          blockAf[ri] = v.imag();
          ri += 1;
      }
  }
}

template<typename Scalar, typename Index, int StorageOrder, int N>
EIGEN_STRONG_INLINE void symm_pack_rhs_helper(Scalar *blockB, const Scalar* _rhs, Index rhsStride, Index rows, Index cols, Index k2)
{
  const Index depth = k2 + rows;
  const_blas_data_mapper<Scalar, Index, StorageOrder> rhs(_rhs, rhsStride);
  const int vectorSize = quad_traits<Scalar>::vectorsize;

  Index ri = 0, j = 0;
  for(; j + N*vectorSize < cols; j+=N*vectorSize)
  {
      Index i = k2;
      for(; i < depth; i++)
      {
        for(int k = 0; k < N*vectorSize; k++)
        {
          if(i <= j+k)
            blockB[ri + k] = rhs(j+k, i);
          else
            blockB[ri + k] = rhs(i, j+k);
        }
        ri += N*vectorSize;
      }
  }
  for(Index i = k2; i < depth; i++)
  {
      Index k = j;
      for(; k < cols; k++)
      {
        if(k <= i)
          blockB[ri] = rhs(i, k);
        else
          blockB[ri] = rhs(k, i);
        ri += 1;
      }
  }
}

template<typename Scalar, typename Index, int StorageOrder>
EIGEN_STRONG_INLINE void symm_pack_lhs_helper(Scalar *blockA, const Scalar* _lhs, Index lhsStride, Index cols, Index rows)
{
  const Index depth = cols;
  const_blas_data_mapper<Scalar, Index, StorageOrder> lhs(_lhs, lhsStride);
  const int vectorSize = quad_traits<Scalar>::vectorsize;
  Index ri = 0, j = 0;

  for(j = 0; j + vectorSize < rows; j+=vectorSize)
  {
      Index i = 0;

      for(; i < depth; i++)
      {
        for(int k = 0; k < vectorSize; k++)
        {
          if(i <= j+k)
            blockA[ri + k] = lhs(j+k, i);
          else
            blockA[ri + k] = lhs(i, j+k);
        }
        ri += vectorSize;
      }
  }

  for(Index i = 0; i < depth; i++)
  {
      Index k = j;
      for(; k < rows; k++)
      {
          if(i <= k)
            blockA[ri] = lhs(k, i);
          else
            blockA[ri] = lhs(i, k);
          ri += 1;
      }
  }
}

template<typename Index, int nr, int StorageOrder>
struct symm_pack_rhs<std::complex<float>, Index, nr, StorageOrder>
{
  void operator()(std::complex<float>* blockB, const std::complex<float>* _rhs, Index rhsStride, Index rows, Index cols, Index k2)
  {
    symm_pack_complex_rhs_helper<float, Index, StorageOrder, 1>(blockB, _rhs, rhsStride, rows, cols, k2);
  }
};

template<typename Index, int Pack1, int Pack2_dummy, int StorageOrder>
struct symm_pack_lhs<std::complex<float>, Index, Pack1, Pack2_dummy, StorageOrder>
{
  void operator()(std::complex<float>* blockA, const std::complex<float>* _lhs, Index lhsStride, Index cols, Index rows)
  {
   symm_pack_complex_lhs_helper<float, Index, StorageOrder>(blockA, _lhs, lhsStride, cols, rows);
  }
};

// *********** symm_pack std::complex<float64> ***********

template<typename Index, int nr, int StorageOrder>
struct symm_pack_rhs<std::complex<double>, Index, nr, StorageOrder>
{
  void operator()(std::complex<double>* blockB, const std::complex<double>* _rhs, Index rhsStride, Index rows, Index cols, Index k2)
  {
    symm_pack_complex_rhs_helper<double, Index, StorageOrder, 2>(blockB, _rhs, rhsStride, rows, cols, k2);
  }
};

template<typename Index, int Pack1, int Pack2_dummy, int StorageOrder>
struct symm_pack_lhs<std::complex<double>, Index, Pack1, Pack2_dummy, StorageOrder>
{
  void operator()(std::complex<double>* blockA, const std::complex<double>* _lhs, Index lhsStride, Index cols, Index rows)
  {
    symm_pack_complex_lhs_helper<double, Index, StorageOrder>(blockA, _lhs, lhsStride, cols, rows);
  }
};

// *********** symm_pack float32 ***********
template<typename Index, int nr, int StorageOrder>
struct symm_pack_rhs<float, Index, nr, StorageOrder>
{
  void operator()(float* blockB, const float* _rhs, Index rhsStride, Index rows, Index cols, Index k2)
  {
   symm_pack_rhs_helper<float, Index, StorageOrder, 1>(blockB, _rhs, rhsStride, rows, cols, k2);
  }
};

template<typename Index, int Pack1, int Pack2_dummy, int StorageOrder>
struct symm_pack_lhs<float, Index, Pack1, Pack2_dummy, StorageOrder>
{
  void operator()(float* blockA, const float* _lhs, Index lhsStride, Index cols, Index rows)
  {
   symm_pack_lhs_helper<float, Index, StorageOrder>(blockA, _lhs, lhsStride, cols, rows);
  }
};

// *********** symm_pack float64 ***********
template<typename Index, int nr, int StorageOrder>
struct symm_pack_rhs<double, Index, nr, StorageOrder>
{
  void operator()(double* blockB, const double* _rhs, Index rhsStride, Index rows, Index cols, Index k2)
  {
    symm_pack_rhs_helper<double, Index, StorageOrder, 2>(blockB, _rhs, rhsStride, rows, cols, k2);
  }
};

template<typename Index, int Pack1, int Pack2_dummy, int StorageOrder>
struct symm_pack_lhs<double, Index, Pack1, Pack2_dummy, StorageOrder>
{
  void operator()(double* blockA, const double* _lhs, Index lhsStride, Index cols, Index rows)
  {
   symm_pack_lhs_helper<double, Index, StorageOrder>(blockA, _lhs, lhsStride, cols, rows);
  }
};

/**
 * PanelMode
 * Packing might be called several times before being multiplied by gebp_kernel, this happens because 
 * on special occasions it fills part of block with other parts of the matrix. Two variables control
 * how PanelMode should behave: offset and stride. The idea is that those variables represent whatever
 * is going to be the real offset and stride in the future and this is what you should obey. The process
 * is to behave as you would with normal packing but leave the start of each part with the correct offset
 * and the end as well respecting the real stride the block will have. Gebp is aware of both blocks stride
 * and offset and behaves accordingly.
 **/

// General template for lhs complex packing.
template<typename Scalar, bool IsComplex, typename Index, typename DataMapper, typename Packet, typename PacketC, int StorageOrder, bool Conjugate, bool PanelMode>
struct lhs_cpack {
  EIGEN_STRONG_INLINE void operator()(std::complex<Scalar>* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride, Index offset)
  {
    const int vectorSize = quad_traits<Scalar>::vectorsize;
    Index ri = 0, j = 0;
    Scalar *blockAt  = reinterpret_cast<Scalar *>(blockA);
    Packet conj = pset1<Packet>((Scalar)-1.0f);

    for(j = 0; j + vectorSize < rows; j+=vectorSize)
    {
      Index i = 0;

      if(PanelMode) ri += vectorSize*offset;

      for(; i + vectorSize < depth; i+=vectorSize)
      {
        PacketBlock<Packet, 4> block;

        PacketBlock<PacketC, 8> cblock;
        if(StorageOrder == ColMajor)
        {
          cblock.packet[0] = lhs.template loadPacket<PacketC>(j, i + 0);
          cblock.packet[1] = lhs.template loadPacket<PacketC>(j, i + 1);
          cblock.packet[2] = lhs.template loadPacket<PacketC>(j, i + 2);
          cblock.packet[3] = lhs.template loadPacket<PacketC>(j, i + 3);

          cblock.packet[4] = lhs.template loadPacket<PacketC>(j + 2, i + 0);
          cblock.packet[5] = lhs.template loadPacket<PacketC>(j + 2, i + 1);
          cblock.packet[6] = lhs.template loadPacket<PacketC>(j + 2, i + 2);
          cblock.packet[7] = lhs.template loadPacket<PacketC>(j + 2, i + 3);
        } else {
          cblock.packet[0] = lhs.template loadPacket<PacketC>(j + 0, i);
          cblock.packet[1] = lhs.template loadPacket<PacketC>(j + 1, i);
          cblock.packet[2] = lhs.template loadPacket<PacketC>(j + 2, i);
          cblock.packet[3] = lhs.template loadPacket<PacketC>(j + 3, i);

          cblock.packet[4] = lhs.template loadPacket<PacketC>(j + 0, i + 2);
          cblock.packet[5] = lhs.template loadPacket<PacketC>(j + 1, i + 2);
          cblock.packet[6] = lhs.template loadPacket<PacketC>(j + 2, i + 2);
          cblock.packet[7] = lhs.template loadPacket<PacketC>(j + 3, i + 2);
        }

        block.packet[0] = vec_perm(cblock.packet[0].v , cblock.packet[4].v, p16uc_GETREAL32);
        block.packet[1] = vec_perm(cblock.packet[1].v , cblock.packet[5].v, p16uc_GETREAL32);
        block.packet[2] = vec_perm(cblock.packet[2].v , cblock.packet[6].v, p16uc_GETREAL32);
        block.packet[3] = vec_perm(cblock.packet[3].v , cblock.packet[7].v, p16uc_GETREAL32);

        if(StorageOrder == RowMajor) ptranspose(block);

        pstore<Scalar>(blockAt + ri     , block.packet[0]);
        pstore<Scalar>(blockAt + ri +  4, block.packet[1]);
        pstore<Scalar>(blockAt + ri +  8, block.packet[2]);
        pstore<Scalar>(blockAt + ri + 12, block.packet[3]);

        ri += 4*vectorSize;
      }
      for(; i < depth; i++)
      {
        blockAt[ri + 0] = lhs(j + 0, i).real();
        blockAt[ri + 1] = lhs(j + 1, i).real();
        blockAt[ri + 2] = lhs(j + 2, i).real();
        blockAt[ri + 3] = lhs(j + 3, i).real();

        ri += vectorSize;
      }
      if(PanelMode) ri += vectorSize*(stride - offset - depth);
      
      i = 0;

      if(PanelMode) ri += vectorSize*offset;

      for(; i + vectorSize < depth; i+=vectorSize)
      {
        PacketBlock<PacketC, 8> cblock;
        if(StorageOrder == ColMajor)
        {
          cblock.packet[0] = lhs.template loadPacket<PacketC>(j, i + 0);
          cblock.packet[1] = lhs.template loadPacket<PacketC>(j, i + 1);
          cblock.packet[2] = lhs.template loadPacket<PacketC>(j, i + 2);
          cblock.packet[3] = lhs.template loadPacket<PacketC>(j, i + 3);

          cblock.packet[4] = lhs.template loadPacket<PacketC>(j + 2, i + 0);
          cblock.packet[5] = lhs.template loadPacket<PacketC>(j + 2, i + 1);
          cblock.packet[6] = lhs.template loadPacket<PacketC>(j + 2, i + 2);
          cblock.packet[7] = lhs.template loadPacket<PacketC>(j + 2, i + 3);
        } else {
          cblock.packet[0] = lhs.template loadPacket<PacketC>(j + 0, i);
          cblock.packet[1] = lhs.template loadPacket<PacketC>(j + 1, i);
          cblock.packet[2] = lhs.template loadPacket<PacketC>(j + 2, i);
          cblock.packet[3] = lhs.template loadPacket<PacketC>(j + 3, i);

          cblock.packet[4] = lhs.template loadPacket<PacketC>(j + 0, i + 2);
          cblock.packet[5] = lhs.template loadPacket<PacketC>(j + 1, i + 2);
          cblock.packet[6] = lhs.template loadPacket<PacketC>(j + 2, i + 2);
          cblock.packet[7] = lhs.template loadPacket<PacketC>(j + 3, i + 2);
        }

        PacketBlock<Packet, 4> block;
        block.packet[0] = vec_perm(cblock.packet[0].v , cblock.packet[4].v, p16uc_GETIMAG32);
        block.packet[1] = vec_perm(cblock.packet[1].v , cblock.packet[5].v, p16uc_GETIMAG32);
        block.packet[2] = vec_perm(cblock.packet[2].v , cblock.packet[6].v, p16uc_GETIMAG32);
        block.packet[3] = vec_perm(cblock.packet[3].v , cblock.packet[7].v, p16uc_GETIMAG32);

        if(Conjugate)
        {
          block.packet[0] *= conj;
          block.packet[1] *= conj;
          block.packet[2] *= conj;
          block.packet[3] *= conj;
        }

        if(StorageOrder == RowMajor) ptranspose(block);

        pstore<Scalar>(blockAt + ri     , block.packet[0]);
        pstore<Scalar>(blockAt + ri +  4, block.packet[1]);
        pstore<Scalar>(blockAt + ri +  8, block.packet[2]);
        pstore<Scalar>(blockAt + ri + 12, block.packet[3]);

        ri += 4*vectorSize;
      }
      for(; i < depth; i++)
      {
        if(Conjugate)
        {
          blockAt[ri + 0] = -lhs(j + 0, i).imag();
          blockAt[ri + 1] = -lhs(j + 1, i).imag();
          blockAt[ri + 2] = -lhs(j + 2, i).imag();
          blockAt[ri + 3] = -lhs(j + 3, i).imag();
        } else {
          blockAt[ri + 0] = lhs(j + 0, i).imag();
          blockAt[ri + 1] = lhs(j + 1, i).imag();
          blockAt[ri + 2] = lhs(j + 2, i).imag();
          blockAt[ri + 3] = lhs(j + 3, i).imag();
        }

        ri += vectorSize;
      }
      if(PanelMode) ri += vectorSize*(stride - offset - depth);
    }

    if(PanelMode) ri += offset*(rows - j);

    for(Index i = 0; i < depth; i++)
    {
      Index k = j;
      for(; k < rows; k++)
      {
        blockAt[ri] = lhs(k, i).real();
        ri += 1;
      }
    }

    if(PanelMode) ri += (rows - j)*(stride - offset - depth);

    if(PanelMode) ri += offset*(rows - j);

    for(Index i = 0; i < depth; i++)
    {
      Index k = j;
      for(; k < rows; k++)
      {
        if(Conjugate)
          blockAt[ri] = -lhs(k, i).imag();
        else
          blockAt[ri] = lhs(k, i).imag();
        ri += 1;
      }
    }

    if(PanelMode) ri += (rows - j)*(stride - offset - depth);
  }
};

// General template for lhs packing.
template<typename Scalar, typename Index, typename DataMapper, typename Packet, int StorageOrder, bool PanelMode>
struct lhs_pack{
  EIGEN_STRONG_INLINE void operator()(Scalar* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride, Index offset)
  {
    const int vectorSize = quad_traits<Scalar>::vectorsize;
    Index ri = 0, j = 0;

    for(j = 0; j + vectorSize < rows; j+=vectorSize)
    {
      Index i = 0;

      if(PanelMode) ri += vectorSize*offset;

      for(; i + vectorSize < depth; i+=vectorSize)
      {
        PacketBlock<Packet, 4> block;

        if(StorageOrder == RowMajor)
        {
          block.packet[0] = lhs.template loadPacket<Packet>(j + 0, i);
          block.packet[1] = lhs.template loadPacket<Packet>(j + 1, i);
          block.packet[2] = lhs.template loadPacket<Packet>(j + 2, i);
          block.packet[3] = lhs.template loadPacket<Packet>(j + 3, i);

          ptranspose(block);
        } else {
          block.packet[0] = lhs.template loadPacket<Packet>(j, i + 0);
          block.packet[1] = lhs.template loadPacket<Packet>(j, i + 1);
          block.packet[2] = lhs.template loadPacket<Packet>(j, i + 2);
          block.packet[3] = lhs.template loadPacket<Packet>(j, i + 3);
        }

        pstore<Scalar>(blockA + ri     , block.packet[0]);
        pstore<Scalar>(blockA + ri +  4, block.packet[1]);
        pstore<Scalar>(blockA + ri +  8, block.packet[2]);
        pstore<Scalar>(blockA + ri + 12, block.packet[3]);

        ri += 4*vectorSize;
      }
      for(; i < depth; i++)
      {
        if(StorageOrder == RowMajor)
        {
          blockA[ri+0] = lhs(j+0, i);
          blockA[ri+1] = lhs(j+1, i);
          blockA[ri+2] = lhs(j+2, i);
          blockA[ri+3] = lhs(j+3, i);
        } else {
          Packet lhsV = lhs.template loadPacket<Packet>(j, i);
          pstore<Scalar>(blockA + ri, lhsV);
        }

        ri += vectorSize;
      }
      if(PanelMode) ri += vectorSize*(stride - offset - depth);
    }

    if(PanelMode) ri += offset*(rows - j);

    for(Index i = 0; i < depth; i++)
    {
      Index k = j;
      for(; k < rows; k++)
      {
        blockA[ri] = lhs(k, i);
        ri += 1;
      }
    }

    if(PanelMode) ri += (rows - j)*(stride - offset - depth);
  }
};

// General template for rhs complex packing.
template<typename Scalar, typename Index, typename DataMapper, typename Packet, typename PacketC, int StorageOrder, bool Conjugate, bool PanelMode>
struct rhs_cpack
{
  EIGEN_STRONG_INLINE void operator()(std::complex<Scalar>* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset)
  {
    const int vectorSize = quad_traits<Scalar>::vectorsize;
    Scalar *blockBt = reinterpret_cast<Scalar *>(blockB);
    Packet conj = pset1<Packet>((Scalar)-1.0f);

    Index ri = 0, j = 0;
    for(; j + vectorSize < cols; j+=vectorSize)
    {
        Index i = 0;

        if(PanelMode) ri += offset*vectorSize;

        for(; i + vectorSize < depth; i+=vectorSize)
        {
            PacketBlock<PacketC, 8> cblock;
            if(StorageOrder == ColMajor)
            {
              cblock.packet[0] = rhs.template loadPacket<PacketC>(i, j + 0);
              cblock.packet[1] = rhs.template loadPacket<PacketC>(i, j + 1);
              cblock.packet[2] = rhs.template loadPacket<PacketC>(i, j + 2);
              cblock.packet[3] = rhs.template loadPacket<PacketC>(i, j + 3);

              cblock.packet[4] = rhs.template loadPacket<PacketC>(i + 2, j + 0);
              cblock.packet[5] = rhs.template loadPacket<PacketC>(i + 2, j + 1);
              cblock.packet[6] = rhs.template loadPacket<PacketC>(i + 2, j + 2);
              cblock.packet[7] = rhs.template loadPacket<PacketC>(i + 2, j + 3);
            } else {
              cblock.packet[0] = rhs.template loadPacket<PacketC>(i + 0, j);
              cblock.packet[1] = rhs.template loadPacket<PacketC>(i + 1, j);
              cblock.packet[2] = rhs.template loadPacket<PacketC>(i + 2, j);
              cblock.packet[3] = rhs.template loadPacket<PacketC>(i + 3, j);

              cblock.packet[4] = rhs.template loadPacket<PacketC>(i + 0, j + 2);
              cblock.packet[5] = rhs.template loadPacket<PacketC>(i + 1, j + 2);
              cblock.packet[6] = rhs.template loadPacket<PacketC>(i + 2, j + 2);
              cblock.packet[7] = rhs.template loadPacket<PacketC>(i + 3, j + 2);
            }

            PacketBlock<Packet, 4> block;
            block.packet[0] = vec_perm(cblock.packet[0].v , cblock.packet[4].v, p16uc_GETREAL32);
            block.packet[1] = vec_perm(cblock.packet[1].v , cblock.packet[5].v, p16uc_GETREAL32);
            block.packet[2] = vec_perm(cblock.packet[2].v , cblock.packet[6].v, p16uc_GETREAL32);
            block.packet[3] = vec_perm(cblock.packet[3].v , cblock.packet[7].v, p16uc_GETREAL32);

            if(StorageOrder == ColMajor) ptranspose(block);

            pstore<Scalar>(blockBt + ri     , block.packet[0]);
            pstore<Scalar>(blockBt + ri +  4, block.packet[1]);
            pstore<Scalar>(blockBt + ri +  8, block.packet[2]);
            pstore<Scalar>(blockBt + ri + 12, block.packet[3]);

            ri += 4*vectorSize;
        }
        for(; i < depth; i++)
        {
            blockBt[ri+0] = rhs(i, j+0).real();
            blockBt[ri+1] = rhs(i, j+1).real();
            blockBt[ri+2] = rhs(i, j+2).real();
            blockBt[ri+3] = rhs(i, j+3).real();
            ri += vectorSize;
        }

        if(PanelMode) ri += vectorSize*(stride - offset - depth);

        i = 0;

        if(PanelMode) ri += offset*vectorSize;

        for(; i + vectorSize < depth; i+=vectorSize)
        {
          PacketBlock<PacketC, 8> cblock;
          if(StorageOrder == ColMajor)
          {

            cblock.packet[0] = rhs.template loadPacket<PacketC>(i, j + 0);
            cblock.packet[1] = rhs.template loadPacket<PacketC>(i, j + 1);
            cblock.packet[2] = rhs.template loadPacket<PacketC>(i, j + 2);
            cblock.packet[3] = rhs.template loadPacket<PacketC>(i, j + 3);

            cblock.packet[4] = rhs.template loadPacket<PacketC>(i + 2, j + 0);
            cblock.packet[5] = rhs.template loadPacket<PacketC>(i + 2, j + 1);
            cblock.packet[6] = rhs.template loadPacket<PacketC>(i + 2, j + 2);
            cblock.packet[7] = rhs.template loadPacket<PacketC>(i + 2, j + 3);
          } else {
            cblock.packet[0] = rhs.template loadPacket<PacketC>(i + 0, j);
            cblock.packet[1] = rhs.template loadPacket<PacketC>(i + 1, j);
            cblock.packet[2] = rhs.template loadPacket<PacketC>(i + 2, j);
            cblock.packet[3] = rhs.template loadPacket<PacketC>(i + 3, j);

            cblock.packet[4] = rhs.template loadPacket<PacketC>(i + 0, j + 2);
            cblock.packet[5] = rhs.template loadPacket<PacketC>(i + 1, j + 2);
            cblock.packet[6] = rhs.template loadPacket<PacketC>(i + 2, j + 2);
            cblock.packet[7] = rhs.template loadPacket<PacketC>(i + 3, j + 2);
          }

          PacketBlock<Packet, 4> block;
          block.packet[0] = vec_perm(cblock.packet[0].v , cblock.packet[4].v, p16uc_GETIMAG32);
          block.packet[1] = vec_perm(cblock.packet[1].v , cblock.packet[5].v, p16uc_GETIMAG32);
          block.packet[2] = vec_perm(cblock.packet[2].v , cblock.packet[6].v, p16uc_GETIMAG32);
          block.packet[3] = vec_perm(cblock.packet[3].v , cblock.packet[7].v, p16uc_GETIMAG32);

          if(Conjugate)
          {
            block.packet[0] *= conj;
            block.packet[1] *= conj;
            block.packet[2] *= conj;
            block.packet[3] *= conj;
          }

          if(StorageOrder == ColMajor) ptranspose(block);

          pstore<Scalar>(blockBt + ri     , block.packet[0]);
          pstore<Scalar>(blockBt + ri +  4, block.packet[1]);
          pstore<Scalar>(blockBt + ri +  8, block.packet[2]);
          pstore<Scalar>(blockBt + ri + 12, block.packet[3]);

          ri += 4*vectorSize;
        }
        for(; i < depth; i++)
        {
            if(Conjugate)
            {
              blockBt[ri+0] = -rhs(i, j+0).imag();
              blockBt[ri+1] = -rhs(i, j+1).imag();
              blockBt[ri+2] = -rhs(i, j+2).imag();
              blockBt[ri+3] = -rhs(i, j+3).imag();
            } else {
              blockBt[ri+0] = rhs(i, j+0).imag();
              blockBt[ri+1] = rhs(i, j+1).imag();
              blockBt[ri+2] = rhs(i, j+2).imag();
              blockBt[ri+3] = rhs(i, j+3).imag();
            }
            ri += vectorSize;
        }

        if(PanelMode) ri += vectorSize*(stride - offset - depth);
    }

    if(PanelMode) ri += offset*(cols - j);

    for(Index i = 0; i < depth; i++)
    {
        Index k = j;
        for(; k < cols; k++)
        {
            blockBt[ri] = rhs(i, k).real();
            ri += 1;
        }
    }
    if(PanelMode) ri += (cols - j)*(stride - offset - depth);

    if(PanelMode) ri += offset*(cols - j);

    for(Index i = 0; i < depth; i++)
    {
        Index k = j;
        for(; k < cols; k++)
        {
            if(Conjugate)
              blockBt[ri] = -rhs(i, k).imag();
            else
              blockBt[ri] = rhs(i, k).imag();
            ri += 1;
        }
    }
    if(PanelMode) ri += (cols - j)*(stride - offset - depth);
  }
};

// General template for rhs packing.
template<typename Scalar, typename Index, typename DataMapper, typename Packet, int StorageOrder, bool PanelMode>
struct rhs_pack {
  EIGEN_STRONG_INLINE void operator()(Scalar* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset)
  {
    const int vectorSize = quad_traits<Scalar>::vectorsize;
    Index ri = 0, j = 0;
    for(; j + vectorSize < cols; j+=vectorSize)
    {
      Index i = 0;

      if(PanelMode) ri += offset*vectorSize;

      for(; i + vectorSize < depth; i+=vectorSize)
      {
        PacketBlock<Packet, 4> block;
        if(StorageOrder == ColMajor)
        {
          block.packet[0] = rhs.template loadPacket<Packet>(i, j + 0);
          block.packet[1] = rhs.template loadPacket<Packet>(i, j + 1);
          block.packet[2] = rhs.template loadPacket<Packet>(i, j + 2);
          block.packet[3] = rhs.template loadPacket<Packet>(i, j + 3);

          ptranspose(block);
        } else {
          block.packet[0] = rhs.template loadPacket<Packet>(i + 0, j);
          block.packet[1] = rhs.template loadPacket<Packet>(i + 1, j);
          block.packet[2] = rhs.template loadPacket<Packet>(i + 2, j);
          block.packet[3] = rhs.template loadPacket<Packet>(i + 3, j);
        }

        pstore<Scalar>(blockB + ri     , block.packet[0]);
        pstore<Scalar>(blockB + ri +  4, block.packet[1]);
        pstore<Scalar>(blockB + ri +  8, block.packet[2]);
        pstore<Scalar>(blockB + ri + 12, block.packet[3]);

        ri += 4*vectorSize;
      }
      for(; i < depth; i++)
      {
        if(StorageOrder == ColMajor)
        {
          blockB[ri+0] = rhs(i, j+0);
          blockB[ri+1] = rhs(i, j+1);
          blockB[ri+2] = rhs(i, j+2);
          blockB[ri+3] = rhs(i, j+3);
        } else {
          Packet rhsV = rhs.template loadPacket<Packet>(i, j);
          pstore<Scalar>(blockB + ri, rhsV);
        }
        ri += vectorSize;
      }

      if(PanelMode) ri += vectorSize*(stride - offset - depth);
    }

    if(PanelMode) ri += offset*(cols - j);

    for(Index i = 0; i < depth; i++)
    {
      Index k = j;
      for(; k < cols; k++)
      {
        blockB[ri] = rhs(i, k);
        ri += 1;
      }
    }
    if(PanelMode) ri += (cols - j)*(stride - offset - depth);
  }
};

// General template for lhs packing, float64 specialization.
template<typename Index, typename DataMapper, int StorageOrder, bool PanelMode>
struct lhs_pack<double,Index, DataMapper, Packet2d, StorageOrder, PanelMode>
{
  EIGEN_STRONG_INLINE void operator()(double* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride, Index offset)
  {
    const int vectorSize = quad_traits<double>::vectorsize;
    Index ri = 0, j = 0;

    for(j = 0; j + vectorSize < rows; j+=vectorSize)
    {
      Index i = 0;

      if(PanelMode) ri += vectorSize*offset;

      for(; i + vectorSize < depth; i+=vectorSize)
      {
        PacketBlock<Packet2d, 2> block;
        if(StorageOrder == RowMajor)
        {
          block.packet[0] = lhs.template loadPacket<Packet2d>(j + 0, i);
          block.packet[1] = lhs.template loadPacket<Packet2d>(j + 1, i);

          ptranspose(block);
        } else {
          block.packet[0] = lhs.template loadPacket<Packet2d>(j, i + 0);
          block.packet[1] = lhs.template loadPacket<Packet2d>(j, i + 1);
        }

        pstore<double>(blockA + ri    , block.packet[0]);
        pstore<double>(blockA + ri + 2, block.packet[1]);

        ri += 2*vectorSize;
      }
      for(; i < depth; i++)
      {
        if(StorageOrder == RowMajor)
        {
          blockA[ri+0] = lhs(j+0, i);
          blockA[ri+1] = lhs(j+1, i);
        } else {
          Packet2d lhsV = lhs.template loadPacket<Packet2d>(j, i);
          pstore<double>(blockA + ri, lhsV);
        }

        ri += vectorSize;
      }
      if(PanelMode) ri += vectorSize*(stride - offset - depth);
    }

    if(PanelMode) ri += offset*(rows - j);

    for(Index i = 0; i < depth; i++)
    {
      Index k = j;
      for(; k < rows; k++)
      {
        blockA[ri] = lhs(k, i);
        ri += 1;
      }
    }

    if(PanelMode) ri += (rows - j)*(stride - offset - depth);
  }
};

// General template for rhs packing, float64 specialization.
template<typename Index, typename DataMapper, int StorageOrder, bool PanelMode>
struct rhs_pack<double, Index, DataMapper, Packet2d, StorageOrder, PanelMode>
{
  EIGEN_STRONG_INLINE void operator()(double* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset)
  {
    const int vectorSize = quad_traits<double>::vectorsize;
    Index ri = 0, j = 0;
    for(; j + 2*vectorSize < cols; j+=2*vectorSize)
    {
      Index i = 0;

      if(PanelMode) ri += offset*(2*vectorSize);
      for(; i + vectorSize < depth; i+=vectorSize)
      {
        PacketBlock<Packet2d, 4> block;
        if(StorageOrder == ColMajor)
        {
          PacketBlock<Packet2d, 2> block1, block2;
          block1.packet[0] = rhs.template loadPacket<Packet2d>(i, j + 0);
          block1.packet[1] = rhs.template loadPacket<Packet2d>(i, j + 1);
          block2.packet[0] = rhs.template loadPacket<Packet2d>(i, j + 2);
          block2.packet[1] = rhs.template loadPacket<Packet2d>(i, j + 3);

          ptranspose(block1);
          ptranspose(block2);

          pstore<double>(blockB + ri    , block1.packet[0]);
          pstore<double>(blockB + ri + 2, block2.packet[0]);
          pstore<double>(blockB + ri + 4, block1.packet[1]);
          pstore<double>(blockB + ri + 6, block2.packet[1]);
        } else {
          block.packet[0] = rhs.template loadPacket<Packet2d>(i + 0, j + 0); //[a1 a2]
          block.packet[1] = rhs.template loadPacket<Packet2d>(i + 0, j + 2); //[a3 a4]
          block.packet[2] = rhs.template loadPacket<Packet2d>(i + 1, j + 0); //[b1 b2]
          block.packet[3] = rhs.template loadPacket<Packet2d>(i + 1, j + 2); //[b3 b4]

          pstore<double>(blockB + ri    , block.packet[0]);
          pstore<double>(blockB + ri + 2, block.packet[1]);
          pstore<double>(blockB + ri + 4, block.packet[2]);
          pstore<double>(blockB + ri + 6, block.packet[3]);
        }

        ri += 4*vectorSize;
      }
      for(; i < depth; i++)
      {
        if(StorageOrder == ColMajor)
        {
          blockB[ri+0] = rhs(i, j+0);
          blockB[ri+1] = rhs(i, j+1);

          ri += vectorSize;

          blockB[ri+0] = rhs(i, j+2);
          blockB[ri+1] = rhs(i, j+3);
        } else {
          Packet2d rhsV = rhs.template loadPacket<Packet2d>(i, j);
          pstore<double>(blockB + ri, rhsV);

          ri += vectorSize;

          rhsV = rhs.template loadPacket<Packet2d>(i, j + 2);
          pstore<double>(blockB + ri, rhsV);
        }
        ri += vectorSize;
      }

      if(PanelMode) ri += (2*vectorSize)*(stride - offset - depth);
    }

    if(PanelMode) ri += offset*(cols - j);

    for(Index i = 0; i < depth; i++)
    {
      Index k = j;
      for(; k < cols; k++)
      {
        blockB[ri] = rhs(i, k);
        ri += 1;
      }
    }
    if(PanelMode) ri += (cols - j)*(stride - offset - depth);
  }
};

// General template for lhs complex packing, float64 specialization.
template<bool IsComplex, typename Index, typename DataMapper, typename Packet, typename PacketC, int StorageOrder, bool Conjugate, bool PanelMode>
struct lhs_cpack<double, IsComplex, Index, DataMapper, Packet, PacketC, StorageOrder, Conjugate, PanelMode>
{
  EIGEN_STRONG_INLINE void operator()(std::complex<double>* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride, Index offset)
  {
    const int vectorSize = quad_traits<double>::vectorsize;
    Index ri = 0, j = 0;
    double *blockAt  = reinterpret_cast<double *>(blockA);
    Packet conj = pset1<Packet>(-1.0);

    for(j = 0; j + vectorSize < rows; j+=vectorSize)
    {
      Index i = 0;

      if(PanelMode) ri += vectorSize*offset;

      for(; i + vectorSize < depth; i+=vectorSize)
      {
        PacketBlock<Packet, 2> block;

        PacketBlock<PacketC, 4> cblock;
        if(StorageOrder == ColMajor)
        {
          cblock.packet[0] = lhs.template loadPacket<PacketC>(j, i + 0); //[a1 a1i]
          cblock.packet[1] = lhs.template loadPacket<PacketC>(j, i + 1); //[b1 b1i]

          cblock.packet[2] = lhs.template loadPacket<PacketC>(j + 1, i + 0); //[a2 a2i]
          cblock.packet[3] = lhs.template loadPacket<PacketC>(j + 1, i + 1); //[b2 b2i]

          block.packet[0] = vec_perm(cblock.packet[0].v , cblock.packet[2].v, p16uc_GETREAL64); //[a1 a2]
          block.packet[1] = vec_perm(cblock.packet[1].v , cblock.packet[3].v, p16uc_GETREAL64); //[b1 b2]
        } else {
          cblock.packet[0] = lhs.template loadPacket<PacketC>(j + 0, i); //[a1 a1i]
          cblock.packet[1] = lhs.template loadPacket<PacketC>(j + 1, i); //[a2 a2i]

          cblock.packet[2] = lhs.template loadPacket<PacketC>(j + 0, i + 1); //[b1 b1i]
          cblock.packet[3] = lhs.template loadPacket<PacketC>(j + 1, i + 1); //[b2 b2i]

          block.packet[0] = vec_perm(cblock.packet[0].v , cblock.packet[1].v, p16uc_GETREAL64); //[a1 a2]
          block.packet[1] = vec_perm(cblock.packet[2].v , cblock.packet[3].v, p16uc_GETREAL64); //[b1 b2]
        }

        pstore<double>(blockAt + ri     , block.packet[0]);
        pstore<double>(blockAt + ri +  2, block.packet[1]);

        ri += 2*vectorSize;
      }
      for(; i < depth; i++)
      {
        blockAt[ri + 0] = lhs(j + 0, i).real();
        blockAt[ri + 1] = lhs(j + 1, i).real();
        ri += vectorSize;
      }
      if(PanelMode) ri += vectorSize*(stride - offset - depth);
      
      i = 0;

      if(PanelMode) ri += vectorSize*offset;

      for(; i + vectorSize < depth; i+=vectorSize)
      {
        PacketBlock<Packet, 2> block;

        PacketBlock<PacketC, 4> cblock;
        if(StorageOrder == ColMajor)
        {
          cblock.packet[0] = lhs.template loadPacket<PacketC>(j, i + 0);
          cblock.packet[1] = lhs.template loadPacket<PacketC>(j, i + 1);

          cblock.packet[2] = lhs.template loadPacket<PacketC>(j + 1, i + 0);
          cblock.packet[3] = lhs.template loadPacket<PacketC>(j + 1, i + 1);

          block.packet[0] = vec_perm(cblock.packet[0].v , cblock.packet[2].v, p16uc_GETIMAG64);
          block.packet[1] = vec_perm(cblock.packet[1].v , cblock.packet[3].v, p16uc_GETIMAG64);
        } else {
          cblock.packet[0] = lhs.template loadPacket<PacketC>(j + 0, i);
          cblock.packet[1] = lhs.template loadPacket<PacketC>(j + 1, i);

          cblock.packet[2] = lhs.template loadPacket<PacketC>(j + 0, i + 1);
          cblock.packet[3] = lhs.template loadPacket<PacketC>(j + 1, i + 1);

          block.packet[0] = vec_perm(cblock.packet[0].v , cblock.packet[1].v, p16uc_GETIMAG64);
          block.packet[1] = vec_perm(cblock.packet[2].v , cblock.packet[3].v, p16uc_GETIMAG64);
        }

        if(Conjugate)
        {
          block.packet[0] *= conj;
          block.packet[1] *= conj;
        }

        pstore<double>(blockAt + ri     , block.packet[0]);
        pstore<double>(blockAt + ri +  2, block.packet[1]);

        ri += 2*vectorSize;
      }
      for(; i < depth; i++)
      {
        if(Conjugate)
        {
          blockAt[ri + 0] = -lhs(j + 0, i).imag();
          blockAt[ri + 1] = -lhs(j + 1, i).imag();
        } else {
          blockAt[ri + 0] = lhs(j + 0, i).imag();
          blockAt[ri + 1] = lhs(j + 1, i).imag();
        }

        ri += vectorSize;
      }
      if(PanelMode) ri += vectorSize*(stride - offset - depth);
    }

    if(PanelMode) ri += offset*(rows - j);

    for(Index i = 0; i < depth; i++)
    {
      Index k = j;
      for(; k < rows; k++)
      {
        blockAt[ri] = lhs(k, i).real();
        ri += 1;
      }
    }

    if(PanelMode) ri += (rows - j)*(stride - offset - depth);

    if(PanelMode) ri += offset*(rows - j);

    for(Index i = 0; i < depth; i++)
    {
      Index k = j;
      for(; k < rows; k++)
      {
        if(Conjugate)
          blockAt[ri] = -lhs(k, i).imag();
        else
          blockAt[ri] = lhs(k, i).imag();
        ri += 1;
      }
    }

    if(PanelMode) ri += (rows - j)*(stride - offset - depth);
  }
};

// General template for rhs complex packing, float64 specialization.
template<typename Index, typename DataMapper, typename Packet, typename PacketC, int StorageOrder, bool Conjugate, bool PanelMode>
struct rhs_cpack<double, Index, DataMapper, Packet, PacketC, StorageOrder, Conjugate, PanelMode>
{
  EIGEN_STRONG_INLINE void operator()(std::complex<double>* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset)
  {
    const int vectorSize = quad_traits<double>::vectorsize;
    double *blockBt = reinterpret_cast<double *>(blockB);
    Packet conj = pset1<Packet>(-1.0);

    Index ri = 0, j = 0;
    for(; j + 2*vectorSize < cols; j+=2*vectorSize)
    {
      Index i = 0;

      if(PanelMode) ri += offset*(2*vectorSize);

      for(; i < depth; i++)
      {
        PacketBlock<PacketC, 4> cblock;
        PacketBlock<Packet, 2> block;

        cblock.packet[0] = rhs.template loadPacket<PacketC>(i, j + 0);
        cblock.packet[1] = rhs.template loadPacket<PacketC>(i, j + 1);
        cblock.packet[2] = rhs.template loadPacket<PacketC>(i, j + 2);
        cblock.packet[3] = rhs.template loadPacket<PacketC>(i, j + 3);

        block.packet[0] = vec_perm(cblock.packet[0].v , cblock.packet[1].v, p16uc_GETREAL64);
        block.packet[1] = vec_perm(cblock.packet[2].v , cblock.packet[3].v, p16uc_GETREAL64);

        pstore<double>(blockBt + ri    , block.packet[0]);
        pstore<double>(blockBt + ri + 2, block.packet[1]);

        ri += 2*vectorSize;
      }

      if(PanelMode) ri += (2*vectorSize)*(stride - offset - depth);

      i = 0;

      if(PanelMode) ri += offset*(2*vectorSize);

      for(; i < depth; i++)
      {
        PacketBlock<PacketC, 4> cblock;
        PacketBlock<Packet, 2> block;

        cblock.packet[0] = rhs.template loadPacket<PacketC>(i, j + 0); //[a1 a1i]
        cblock.packet[1] = rhs.template loadPacket<PacketC>(i, j + 1); //[b1 b1i]
        cblock.packet[2] = rhs.template loadPacket<PacketC>(i, j + 2); //[c1 c1i]
        cblock.packet[3] = rhs.template loadPacket<PacketC>(i, j + 3); //[d1 d1i]

        block.packet[0] = vec_perm(cblock.packet[0].v , cblock.packet[1].v, p16uc_GETIMAG64);
        block.packet[1] = vec_perm(cblock.packet[2].v , cblock.packet[3].v, p16uc_GETIMAG64);

        if(Conjugate)
        {
          block.packet[0] *= conj;
          block.packet[1] *= conj;
        }

        pstore<double>(blockBt + ri     , block.packet[0]);
        pstore<double>(blockBt + ri +  2, block.packet[1]);

        ri += 2*vectorSize;
      }
      if(PanelMode) ri += (2*vectorSize)*(stride - offset - depth);
    }

    if(PanelMode) ri += offset*(cols - j);

    for(Index i = 0; i < depth; i++)
    {
      Index k = j;
      for(; k < cols; k++)
      {
        blockBt[ri] = rhs(i, k).real();
        ri += 1;
      }
    }
    if(PanelMode) ri += (cols - j)*(stride - offset - depth);

    if(PanelMode) ri += offset*(cols - j);

    for(Index i = 0; i < depth; i++)
    {
      Index k = j;
      for(; k < cols; k++)
      {
        if(Conjugate)
          blockBt[ri] = -rhs(i, k).imag();
        else
          blockBt[ri] = rhs(i, k).imag();
        ri += 1;
      }
    }
    if(PanelMode) ri += (cols - j)*(stride - offset - depth);
  }
};

/**************
 * GEMM utils *
 **************/

// 512-bits rank1-update of acc. It can either positive or negative accumulate (useful for complex gemm).
template<typename Scalar, typename Packet, bool NegativeAccumulate>
EIGEN_STRONG_INLINE void pger(PacketBlock<Packet, 4> *acc, const Scalar* lhs, const Scalar* rhs)
{
  Packet lhsV = *((Packet *) lhs);
  Packet rhsV1 = pset1<Packet>(rhs[0]);
  Packet rhsV2 = pset1<Packet>(rhs[1]);
  Packet rhsV3 = pset1<Packet>(rhs[2]);
  Packet rhsV4 = pset1<Packet>(rhs[3]);

  if(NegativeAccumulate)
  {
    acc->packet[0] -= lhsV*rhsV1;
    acc->packet[1] -= lhsV*rhsV2;
    acc->packet[2] -= lhsV*rhsV3;
    acc->packet[3] -= lhsV*rhsV4;
  } else {
    acc->packet[0] += lhsV*rhsV1;
    acc->packet[1] += lhsV*rhsV2;
    acc->packet[2] += lhsV*rhsV3;
    acc->packet[3] += lhsV*rhsV4;
  }
}

// 512-bits rank1-update of complex acc. It takes decoupled accumulators as entries. It also takes cares of mixed types real * complex and complex * real.
template<typename Scalar, typename Packet, bool ConjugateLhs, bool ConjugateRhs, bool LhsIsReal, bool RhsIsReal>
EIGEN_STRONG_INLINE void pgerc(PacketBlock<Packet, 4>& accReal, PacketBlock<Packet,4>& accImag, const Scalar *rhs_ptr, const Scalar *rhs_ptr_imag, const Scalar *lhs_ptr, const Scalar* lhs_ptr_imag, Packet& conj)
{
  Packet lhsV  = *((Packet *) lhs_ptr);
  Packet rhsV1 = pset1<Packet>(rhs_ptr[0]);
  Packet rhsV2 = pset1<Packet>(rhs_ptr[1]);
  Packet rhsV3 = pset1<Packet>(rhs_ptr[2]);
  Packet rhsV4 = pset1<Packet>(rhs_ptr[3]);

  Packet lhsVi;
  if(!LhsIsReal) lhsVi = *((Packet *) lhs_ptr_imag);
  Packet rhsV1i, rhsV2i, rhsV3i, rhsV4i;
  if(!RhsIsReal)
  {
    rhsV1i = pset1<Packet>(rhs_ptr_imag[0]);
    rhsV2i = pset1<Packet>(rhs_ptr_imag[1]);
    rhsV3i = pset1<Packet>(rhs_ptr_imag[2]);
    rhsV4i = pset1<Packet>(rhs_ptr_imag[3]);
  }

  if(ConjugateLhs && !LhsIsReal) lhsVi = pmul<Packet>(lhsVi,conj);
  if(ConjugateRhs && !RhsIsReal)
  {
    rhsV1i = pmul<Packet>(rhsV1i,conj);
    rhsV2i = pmul<Packet>(rhsV2i,conj);
    rhsV3i = pmul<Packet>(rhsV3i,conj);
    rhsV4i = pmul<Packet>(rhsV4i,conj);
  }

  if(LhsIsReal)
  {
    accReal.packet[0] = pmadd<Packet>(rhsV1, lhsV, accReal.packet[0]);
    accReal.packet[1] = pmadd<Packet>(rhsV2, lhsV, accReal.packet[1]);
    accReal.packet[2] = pmadd<Packet>(rhsV3, lhsV, accReal.packet[2]);
    accReal.packet[3] = pmadd<Packet>(rhsV4, lhsV, accReal.packet[3]);

    accImag.packet[0] = pmadd<Packet>(rhsV1i, lhsV, accImag.packet[0]);
    accImag.packet[1] = pmadd<Packet>(rhsV2i, lhsV, accImag.packet[1]);
    accImag.packet[2] = pmadd<Packet>(rhsV3i, lhsV, accImag.packet[2]);
    accImag.packet[3] = pmadd<Packet>(rhsV4i, lhsV, accImag.packet[3]);
  } else if(RhsIsReal) {
    accReal.packet[0] = pmadd<Packet>(rhsV1, lhsV, accReal.packet[0]);
    accReal.packet[1] = pmadd<Packet>(rhsV2, lhsV, accReal.packet[1]);
    accReal.packet[2] = pmadd<Packet>(rhsV3, lhsV, accReal.packet[2]);
    accReal.packet[3] = pmadd<Packet>(rhsV4, lhsV, accReal.packet[3]);

    accImag.packet[0] = pmadd<Packet>(rhsV1, lhsVi, accImag.packet[0]);
    accImag.packet[1] = pmadd<Packet>(rhsV2, lhsVi, accImag.packet[1]);
    accImag.packet[2] = pmadd<Packet>(rhsV3, lhsVi, accImag.packet[2]);
    accImag.packet[3] = pmadd<Packet>(rhsV4, lhsVi, accImag.packet[3]);
  } else {
    accReal.packet[0] = pmadd<Packet>(rhsV1, lhsV, accReal.packet[0]);
    accReal.packet[1] = pmadd<Packet>(rhsV2, lhsV, accReal.packet[1]);
    accReal.packet[2] = pmadd<Packet>(rhsV3, lhsV, accReal.packet[2]);
    accReal.packet[3] = pmadd<Packet>(rhsV4, lhsV, accReal.packet[3]);

    accImag.packet[0] = pmadd<Packet>(rhsV1i, lhsV, accImag.packet[0]);
    accImag.packet[1] = pmadd<Packet>(rhsV2i, lhsV, accImag.packet[1]);
    accImag.packet[2] = pmadd<Packet>(rhsV3i, lhsV, accImag.packet[2]);
    accImag.packet[3] = pmadd<Packet>(rhsV4i, lhsV, accImag.packet[3]);

    accReal.packet[0] = psub<Packet>(accReal.packet[0], pmul<Packet>(rhsV1i, lhsVi));
    accReal.packet[1] = psub<Packet>(accReal.packet[1], pmul<Packet>(rhsV2i, lhsVi));
    accReal.packet[2] = psub<Packet>(accReal.packet[2], pmul<Packet>(rhsV3i, lhsVi));
    accReal.packet[3] = psub<Packet>(accReal.packet[3], pmul<Packet>(rhsV4i, lhsVi));

    accImag.packet[0] = pmadd<Packet>(rhsV1, lhsVi, accImag.packet[0]);
    accImag.packet[1] = pmadd<Packet>(rhsV2, lhsVi, accImag.packet[1]);
    accImag.packet[2] = pmadd<Packet>(rhsV3, lhsVi, accImag.packet[2]);
    accImag.packet[3] = pmadd<Packet>(rhsV4, lhsVi, accImag.packet[3]);
  }
}

template<typename Scalar, typename Packet>
EIGEN_STRONG_INLINE Packet ploadLhs(const Scalar *lhs)
{
    return *((Packet *)lhs);
}

// Zero the accumulator on PacketBlock.
template<typename Scalar, typename Packet>
EIGEN_STRONG_INLINE void bsetzero(PacketBlock<Packet,4>& acc)
{
  acc.packet[0] = pset1<Packet>((Scalar)0);
  acc.packet[1] = pset1<Packet>((Scalar)0);
  acc.packet[2] = pset1<Packet>((Scalar)0);
  acc.packet[3] = pset1<Packet>((Scalar)0);
}

// Scale the PacketBlock vectors by alpha.
template<typename Packet>
EIGEN_STRONG_INLINE void bscale(PacketBlock<Packet,4>& acc, PacketBlock<Packet,4>& accZ, const Packet& pAlpha)
{
  acc.packet[0] = pmadd(pAlpha,accZ.packet[0], acc.packet[0]);
  acc.packet[1] = pmadd(pAlpha,accZ.packet[1], acc.packet[1]);
  acc.packet[2] = pmadd(pAlpha,accZ.packet[2], acc.packet[2]);
  acc.packet[3] = pmadd(pAlpha,accZ.packet[3], acc.packet[3]);
}

// Complex version of PacketBlock scaling.
template<typename Packet>
EIGEN_STRONG_INLINE void bscalec(PacketBlock<Packet,4>& aReal, PacketBlock<Packet,4>& aImag, const Packet& bReal, const Packet& bImag, PacketBlock<Packet,4>& cReal, PacketBlock<Packet,4>& cImag)
{
  cReal.packet[0] = pmul<Packet>(aReal.packet[0], bReal);
  cReal.packet[1] = pmul<Packet>(aReal.packet[1], bReal);
  cReal.packet[2] = pmul<Packet>(aReal.packet[2], bReal);
  cReal.packet[3] = pmul<Packet>(aReal.packet[3], bReal);

  cImag.packet[0] = pmul<Packet>(aImag.packet[0], bReal);
  cImag.packet[1] = pmul<Packet>(aImag.packet[1], bReal);
  cImag.packet[2] = pmul<Packet>(aImag.packet[2], bReal);
  cImag.packet[3] = pmul<Packet>(aImag.packet[3], bReal);

  cReal.packet[0] = psub<Packet>(cReal.packet[0], pmul<Packet>(aImag.packet[0], bImag));
  cReal.packet[1] = psub<Packet>(cReal.packet[1], pmul<Packet>(aImag.packet[1], bImag));
  cReal.packet[2] = psub<Packet>(cReal.packet[2], pmul<Packet>(aImag.packet[2], bImag));
  cReal.packet[3] = psub<Packet>(cReal.packet[3], pmul<Packet>(aImag.packet[3], bImag));

  cImag.packet[0] = pmadd<Packet>(aReal.packet[0], bImag, cImag.packet[0]);
  cImag.packet[1] = pmadd<Packet>(aReal.packet[1], bImag, cImag.packet[1]);
  cImag.packet[2] = pmadd<Packet>(aReal.packet[2], bImag, cImag.packet[2]);
  cImag.packet[3] = pmadd<Packet>(aReal.packet[3], bImag, cImag.packet[3]);
}

// Load a PacketBlock, the N parameters make tunning gemm easier so we can add more accumulators as needed.
template<typename DataMapper, typename Packet, typename Index, int N>
EIGEN_STRONG_INLINE void bload(PacketBlock<Packet,4>& acc, const DataMapper& res, Index row, Index col, Index accCols)
{
  acc.packet[0] = res.template loadPacket<Packet>(row + N*accCols, col + 0);
  acc.packet[1] = res.template loadPacket<Packet>(row + N*accCols, col + 1);
  acc.packet[2] = res.template loadPacket<Packet>(row + N*accCols, col + 2);
  acc.packet[3] = res.template loadPacket<Packet>(row + N*accCols, col + 3);
}

// An overload of bload when you have a PacketBLock with 8 vectors.
template<typename DataMapper, typename Packet, typename Index, int N>
EIGEN_STRONG_INLINE void bload(PacketBlock<Packet,8>& acc, const DataMapper& res, Index row, Index col, Index accCols)
{
  acc.packet[0] = res.template loadPacket<Packet>(row + N*accCols, col + 0);
  acc.packet[1] = res.template loadPacket<Packet>(row + N*accCols, col + 1);
  acc.packet[2] = res.template loadPacket<Packet>(row + N*accCols, col + 2);
  acc.packet[3] = res.template loadPacket<Packet>(row + N*accCols, col + 3);
  acc.packet[4] = res.template loadPacket<Packet>(row + (N+1)*accCols, col + 0);
  acc.packet[5] = res.template loadPacket<Packet>(row + (N+1)*accCols, col + 1);
  acc.packet[6] = res.template loadPacket<Packet>(row + (N+1)*accCols, col + 2);
  acc.packet[7] = res.template loadPacket<Packet>(row + (N+1)*accCols, col + 3);
}


// PEEL loop factor.
#define PEEL 10

/****************
 * GEMM kernels *
 * **************/
template<typename Scalar, typename Index, typename Packet, typename RhsPacket, typename DataMapper>
EIGEN_STRONG_INLINE void gemm(const DataMapper& res, const Scalar* blockA, const Scalar* blockB,
          Index rows, Index depth, Index cols, Scalar alpha, Index strideA, Index strideB, Index offsetA, Index offsetB, const int accRows, const int accCols)
{
      const Index remaining_rows = rows % accCols;
      const Index remaining_cols = cols % accRows;

      if( strideA == -1 ) strideA = depth;
      if( strideB == -1 ) strideB = depth;

      const Packet pAlpha = pset1<Packet>(alpha);
      Index col = 0;
      for(; col + accRows <= cols; col += accRows)
      {
        const Scalar *rhs_base = blockB + ( col/accRows     )*strideB*accRows;
        const Scalar *lhs_base = blockA;

        Index row = 0;
        for(; row + 6*accCols <= rows; row += 6*accCols)
        {
#define MICRO() \
            pger<Scalar, Packet, false>(&accZero1, lhs_ptr1, rhs_ptr); \
            lhs_ptr1 += accCols; \
            pger<Scalar, Packet, false>(&accZero2, lhs_ptr2, rhs_ptr); \
            lhs_ptr2 += accCols; \
            pger<Scalar, Packet, false>(&accZero3, lhs_ptr3, rhs_ptr); \
            lhs_ptr3 += accCols; \
            pger<Scalar, Packet, false>(&accZero4, lhs_ptr4, rhs_ptr); \
            lhs_ptr4 += accCols; \
            pger<Scalar, Packet, false>(&accZero5, lhs_ptr5, rhs_ptr); \
            lhs_ptr5 += accCols; \
            pger<Scalar, Packet, false>(&accZero6, lhs_ptr6, rhs_ptr); \
            lhs_ptr6 += accCols; \
            rhs_ptr += accRows;

          const Scalar *rhs_ptr  = rhs_base;
          const Scalar *lhs_ptr1 = lhs_base + (row/accCols)*strideA*accCols;
          const Scalar *lhs_ptr2 = lhs_base + ((row/accCols) + 1)*strideA*accCols;
          const Scalar *lhs_ptr3 = lhs_base + ((row/accCols) + 2)*strideA*accCols;
          const Scalar *lhs_ptr4 = lhs_base + ((row/accCols) + 3)*strideA*accCols;
          const Scalar *lhs_ptr5 = lhs_base + ((row/accCols) + 4)*strideA*accCols;
          const Scalar *lhs_ptr6 = lhs_base + ((row/accCols) + 5)*strideA*accCols;

          PacketBlock<Packet,4> acc1, accZero1;
          PacketBlock<Packet,4> acc2, accZero2;
          PacketBlock<Packet,4> acc3, accZero3;
          PacketBlock<Packet,4> acc4, accZero4;
          PacketBlock<Packet,4> acc5, accZero5;
          PacketBlock<Packet,4> acc6, accZero6;

          bload<DataMapper, Packet, Index, 0>(acc1, res, row, col, accCols);
          bsetzero<Scalar, Packet>(accZero1);
          bload<DataMapper, Packet, Index, 1>(acc2, res, row, col, accCols);
          bsetzero<Scalar, Packet>(accZero2);
          bload<DataMapper, Packet, Index, 2>(acc3, res, row, col, accCols);
          bsetzero<Scalar, Packet>(accZero3);
          bload<DataMapper, Packet, Index, 3>(acc4, res, row, col, accCols);
          bsetzero<Scalar, Packet>(accZero4);
          bload<DataMapper, Packet, Index, 4>(acc5, res, row, col, accCols);
          bsetzero<Scalar, Packet>(accZero5);
          bload<DataMapper, Packet, Index, 5>(acc6, res, row, col, accCols);
          bsetzero<Scalar, Packet>(accZero6);

          lhs_ptr1 += accCols*offsetA;
          lhs_ptr2 += accCols*offsetA;
          lhs_ptr3 += accCols*offsetA;
          lhs_ptr4 += accCols*offsetA;
          lhs_ptr5 += accCols*offsetA;
          lhs_ptr6 += accCols*offsetA;
          rhs_ptr += accRows*offsetB;

          Index k = 0;
          for(; k + PEEL < depth; k+= PEEL)
          {
            prefetch(rhs_ptr);
            prefetch(lhs_ptr1);
            prefetch(lhs_ptr2);
            prefetch(lhs_ptr3);
            prefetch(lhs_ptr4);
            prefetch(lhs_ptr5);
            prefetch(lhs_ptr6);
            MICRO();
            MICRO();
            MICRO();
            MICRO();
            MICRO();
            MICRO();
            MICRO();
            MICRO();
#if PEEL > 8
            MICRO();
            MICRO();
#endif
          }
          for(; k < depth; k++)
          {
            MICRO();
          }

          bscale<Packet>(acc1,accZero1, pAlpha);
          bscale<Packet>(acc2,accZero2, pAlpha);
          bscale<Packet>(acc3,accZero3, pAlpha);
          bscale<Packet>(acc4,accZero4, pAlpha);
          bscale<Packet>(acc5,accZero5, pAlpha);
          bscale<Packet>(acc6,accZero6, pAlpha);

          res.template storePacketBlock<Packet, 4>(row + 0*accCols, col, acc1);
          res.template storePacketBlock<Packet, 4>(row + 1*accCols, col, acc2);
          res.template storePacketBlock<Packet, 4>(row + 2*accCols, col, acc3);
          res.template storePacketBlock<Packet, 4>(row + 3*accCols, col, acc4);
          res.template storePacketBlock<Packet, 4>(row + 4*accCols, col, acc5);
          res.template storePacketBlock<Packet, 4>(row + 5*accCols, col, acc6);
#undef MICRO
        }
        for(; row + 5*accCols <= rows; row += 5*accCols)
        {
#define MICRO() \
            pger<Scalar, Packet, false>(&accZero1, lhs_ptr1, rhs_ptr); \
            lhs_ptr1 += accCols; \
            pger<Scalar, Packet, false>(&accZero2, lhs_ptr2, rhs_ptr); \
            lhs_ptr2 += accCols; \
            pger<Scalar, Packet, false>(&accZero3, lhs_ptr3, rhs_ptr); \
            lhs_ptr3 += accCols; \
            pger<Scalar, Packet, false>(&accZero4, lhs_ptr4, rhs_ptr); \
            lhs_ptr4 += accCols; \
            pger<Scalar, Packet, false>(&accZero5, lhs_ptr5, rhs_ptr); \
            lhs_ptr5 += accCols; \
            rhs_ptr += accRows;

          const Scalar *rhs_ptr  = rhs_base;
          const Scalar *lhs_ptr1 = lhs_base + (row/accCols      )*strideA*accCols;
          const Scalar *lhs_ptr2 = lhs_base + ((row/accCols) + 1)*strideA*accCols;
          const Scalar *lhs_ptr3 = lhs_base + ((row/accCols) + 2)*strideA*accCols;
          const Scalar *lhs_ptr4 = lhs_base + ((row/accCols) + 3)*strideA*accCols;
          const Scalar *lhs_ptr5 = lhs_base + ((row/accCols) + 4)*strideA*accCols;

          PacketBlock<Packet,4> acc1, accZero1;
          PacketBlock<Packet,4> acc2, accZero2;
          PacketBlock<Packet,4> acc3, accZero3;
          PacketBlock<Packet,4> acc4, accZero4;
          PacketBlock<Packet,4> acc5, accZero5;

          bload<DataMapper, Packet, Index, 0>(acc1, res, row, col, accCols);
          bsetzero<Scalar, Packet>(accZero1);
          bload<DataMapper, Packet, Index, 1>(acc2, res, row, col, accCols);
          bsetzero<Scalar, Packet>(accZero2);
          bload<DataMapper, Packet, Index, 2>(acc3, res, row, col, accCols);
          bsetzero<Scalar, Packet>(accZero3);
          bload<DataMapper, Packet, Index, 3>(acc4, res, row, col, accCols);
          bsetzero<Scalar, Packet>(accZero4);
          bload<DataMapper, Packet, Index, 4>(acc5, res, row, col, accCols);
          bsetzero<Scalar, Packet>(accZero5);

          lhs_ptr1 += accCols*offsetA;
          lhs_ptr2 += accCols*offsetA;
          lhs_ptr3 += accCols*offsetA;
          lhs_ptr4 += accCols*offsetA;
          lhs_ptr5 += accCols*offsetA;
          rhs_ptr += accRows*offsetB;
          Index k = 0;

          for(; k + PEEL < depth; k+= PEEL)
          {
            prefetch(rhs_ptr);
            prefetch(lhs_ptr1);
            prefetch(lhs_ptr2);
            prefetch(lhs_ptr3);
            prefetch(lhs_ptr4);
            prefetch(lhs_ptr5);

            MICRO();
            MICRO();
            MICRO();
            MICRO();
            MICRO();
            MICRO();
            MICRO();
            MICRO();
#if PEEL > 8
            MICRO();
            MICRO();
#endif
          }
          for(; k < depth; k++)
          {
            MICRO();
          }

          bscale<Packet>(acc1,accZero1, pAlpha);
          bscale<Packet>(acc2,accZero2, pAlpha);
          bscale<Packet>(acc3,accZero3, pAlpha);
          bscale<Packet>(acc4,accZero4, pAlpha);
          bscale<Packet>(acc5,accZero5, pAlpha);

          res.template storePacketBlock<Packet, 4>(row + 0*accCols, col, acc1);
          res.template storePacketBlock<Packet, 4>(row + 1*accCols, col, acc2);
          res.template storePacketBlock<Packet, 4>(row + 2*accCols, col, acc3);
          res.template storePacketBlock<Packet, 4>(row + 3*accCols, col, acc4);
          res.template storePacketBlock<Packet, 4>(row + 4*accCols, col, acc5);
#undef MICRO
        }
        for(; row + 4*accCols <= rows; row += 4*accCols)
        {
#define MICRO() \
            pger<Scalar, Packet, false>(&accZero1, lhs_ptr1, rhs_ptr); \
            lhs_ptr1 += accCols; \
            pger<Scalar, Packet, false>(&accZero2, lhs_ptr2, rhs_ptr); \
            lhs_ptr2 += accCols; \
            pger<Scalar, Packet, false>(&accZero3, lhs_ptr3, rhs_ptr); \
            lhs_ptr3 += accCols; \
            pger<Scalar, Packet, false>(&accZero4, lhs_ptr4, rhs_ptr); \
            lhs_ptr4 += accCols; \
            rhs_ptr += accRows;

          const Scalar *rhs_ptr  = rhs_base;
          const Scalar *lhs_ptr1 = lhs_base + (row/accCols      )*strideA*accCols;
          const Scalar *lhs_ptr2 = lhs_base + ((row/accCols) + 1)*strideA*accCols;
          const Scalar *lhs_ptr3 = lhs_base + ((row/accCols) + 2)*strideA*accCols;
          const Scalar *lhs_ptr4 = lhs_base + ((row/accCols) + 3)*strideA*accCols;

          PacketBlock<Packet,4> acc1, accZero1;
          PacketBlock<Packet,4> acc2, accZero2;
          PacketBlock<Packet,4> acc3, accZero3;
          PacketBlock<Packet,4> acc4, accZero4;

          bload<DataMapper, Packet, Index, 0>(acc1, res, row, col, accCols);
          bsetzero<Scalar, Packet>(accZero1);
          bload<DataMapper, Packet, Index, 1>(acc2, res, row, col, accCols);
          bsetzero<Scalar, Packet>(accZero2);
          bload<DataMapper, Packet, Index, 2>(acc3, res, row, col, accCols);
          bsetzero<Scalar, Packet>(accZero3);
          bload<DataMapper, Packet, Index, 3>(acc4, res, row, col, accCols);
          bsetzero<Scalar, Packet>(accZero4);

          lhs_ptr1 += accCols*offsetA;
          lhs_ptr2 += accCols*offsetA;
          lhs_ptr3 += accCols*offsetA;
          lhs_ptr4 += accCols*offsetA;
          rhs_ptr += accRows*offsetB;
          Index k = 0;

          for(; k + PEEL < depth; k+= PEEL)
          {
            prefetch(rhs_ptr);
            prefetch(lhs_ptr1);
            prefetch(lhs_ptr2);
            prefetch(lhs_ptr3);
            prefetch(lhs_ptr4);
            MICRO();
            MICRO();
            MICRO();
            MICRO();
            MICRO();
            MICRO();
            MICRO();
            MICRO();
#if PEEL > 8
            MICRO();
            MICRO();
#endif
          }
          for(; k < depth; k++)
          {
            MICRO();
          }

          bscale<Packet>(acc1,accZero1, pAlpha);
          bscale<Packet>(acc2,accZero2, pAlpha);
          bscale<Packet>(acc3,accZero3, pAlpha);
          bscale<Packet>(acc4,accZero4, pAlpha);

          res.template storePacketBlock<Packet, 4>(row + 0*accCols, col, acc1);
          res.template storePacketBlock<Packet, 4>(row + 1*accCols, col, acc2);
          res.template storePacketBlock<Packet, 4>(row + 2*accCols, col, acc3);
          res.template storePacketBlock<Packet, 4>(row + 3*accCols, col, acc4);
#undef MICRO
        }
        for(; row + 3*accCols <= rows; row += 3*accCols)
        {
#define MICRO() \
            pger<Scalar, Packet, false>(&accZero1, lhs_ptr1, rhs_ptr); \
            lhs_ptr1 += accCols; \
            pger<Scalar, Packet, false>(&accZero2, lhs_ptr2, rhs_ptr); \
            lhs_ptr2 += accCols; \
            pger<Scalar, Packet, false>(&accZero3, lhs_ptr3, rhs_ptr); \
            lhs_ptr3 += accCols; \
            rhs_ptr += accRows;

          const Scalar *rhs_ptr  = rhs_base;
          const Scalar *lhs_ptr1 = lhs_base + (row/accCols      )*strideA*accCols;
          const Scalar *lhs_ptr2 = lhs_base + ((row/accCols) + 1)*strideA*accCols;
          const Scalar *lhs_ptr3 = lhs_base + ((row/accCols) + 2)*strideA*accCols;

          PacketBlock<Packet,4> acc1, accZero1;
          PacketBlock<Packet,4> acc2, accZero2;
          PacketBlock<Packet,4> acc3, accZero3;

          bload<DataMapper, Packet, Index, 0>(acc1, res, row, col, accCols);
          bsetzero<Scalar, Packet>(accZero1);
          bload<DataMapper, Packet, Index, 1>(acc2, res, row, col, accCols);
          bsetzero<Scalar, Packet>(accZero2);
          bload<DataMapper, Packet, Index, 2>(acc3, res, row, col, accCols);
          bsetzero<Scalar, Packet>(accZero3);

          lhs_ptr1 += accCols*offsetA;
          lhs_ptr2 += accCols*offsetA;
          lhs_ptr3 += accCols*offsetA;
          rhs_ptr += accRows*offsetB;
          Index k = 0;
          for(; k + PEEL < depth; k+= PEEL)
          {
            prefetch(rhs_ptr);
            prefetch(lhs_ptr1);
            prefetch(lhs_ptr2);
            prefetch(lhs_ptr3);

            MICRO();
            MICRO();
            MICRO();
            MICRO();
            MICRO();
            MICRO();
            MICRO();
            MICRO();
#if PEEL > 8
            MICRO();
            MICRO();
#endif
          }
          for(; k < depth; k++)
          {
            MICRO();
          }

          bscale<Packet>(acc1,accZero1, pAlpha);
          bscale<Packet>(acc2,accZero2, pAlpha);
          bscale<Packet>(acc3,accZero3, pAlpha);

          res.template storePacketBlock<Packet, 4>(row + 0*accCols, col, acc1);
          res.template storePacketBlock<Packet, 4>(row + 1*accCols, col, acc2);
          res.template storePacketBlock<Packet, 4>(row + 2*accCols, col, acc3);
#undef MICRO
        }
        for(; row + 2*accCols <= rows; row += 2*accCols)
        {
#define MICRO() \
            pger<Scalar, Packet, false>(&accZero1, lhs_ptr1, rhs_ptr); \
            lhs_ptr1 += accCols; \
            pger<Scalar, Packet, false>(&accZero2, lhs_ptr2, rhs_ptr); \
            lhs_ptr2 += accCols; \
            rhs_ptr += accRows;

          const Scalar *rhs_ptr  = rhs_base;
          const Scalar *lhs_ptr1 = lhs_base + (row/accCols      )*strideA*accCols;
          const Scalar *lhs_ptr2 = lhs_base + ((row/accCols) + 1)*strideA*accCols;
          
          PacketBlock<Packet,4> acc1, accZero1;
          PacketBlock<Packet,4> acc2, accZero2;

          bload<DataMapper, Packet, Index, 0>(acc1, res, row, col, accCols);
          bsetzero<Scalar, Packet>(accZero1);
          bload<DataMapper, Packet, Index, 1>(acc2, res, row, col, accCols);
          bsetzero<Scalar, Packet>(accZero2);

          lhs_ptr1 += accCols*offsetA;
          lhs_ptr2 += accCols*offsetA;
          rhs_ptr += accRows*offsetB;
          Index k = 0;
          for(; k + PEEL < depth; k+= PEEL)
          {
            prefetch(rhs_ptr);
            prefetch(lhs_ptr1);
            prefetch(lhs_ptr2);

            MICRO();
            MICRO();
            MICRO();
            MICRO();
            MICRO();
            MICRO();
            MICRO();
            MICRO();
#if PEEL > 8
            MICRO();
            MICRO();
#endif
          }
          for(; k < depth; k++)
          {
            MICRO();
          }

          bscale<Packet>(acc1,accZero1, pAlpha);
          bscale<Packet>(acc2,accZero2, pAlpha);

          res.template storePacketBlock<Packet, 4>(row + 0*accCols, col, acc1);
          res.template storePacketBlock<Packet, 4>(row + 1*accCols, col, acc2);
#undef MICRO
        }

        for(; row + accCols <= rows; row += accCols)
        {
#define MICRO() \
            pger<Scalar, Packet, false>(&accZero1, lhs_ptr1, rhs_ptr); \
            lhs_ptr1 += accCols; \
            rhs_ptr += accRows;

          const Scalar *rhs_ptr  = rhs_base;
          const Scalar *lhs_ptr1 = lhs_base + (row/accCols)*strideA*accCols;

          PacketBlock<Packet,4> acc1, accZero1;

          bload<DataMapper, Packet, Index, 0>(acc1, res, row, col, accCols);
          bsetzero<Scalar, Packet>(accZero1);

          lhs_ptr1 += accCols*offsetA;
          rhs_ptr += accRows*offsetB;
          Index k = 0;
          for(; k + PEEL < depth; k+= PEEL)
          {
            prefetch(rhs_ptr);
            prefetch(lhs_ptr1);

            MICRO();
            MICRO();
            MICRO();
            MICRO();
            MICRO();
            MICRO();
            MICRO();
            MICRO();
#if PEEL > 8
            MICRO();
            MICRO();
#endif
          }
          for(; k < depth; k++)
          {
            MICRO();
          }

          bscale<Packet>(acc1,accZero1, pAlpha);

          res.template storePacketBlock<Packet, 4>(row, col, acc1);
#undef MICRO
        }
        if(remaining_rows > 0)
        {
          const Scalar *rhs_ptr  = rhs_base;
          const Scalar *lhs_ptr = lhs_base + (row/accCols)*strideA*accCols;

          lhs_ptr += remaining_rows*offsetA;
          rhs_ptr += accRows*offsetB;
          for(Index k = 0; k < depth; k++)
          {
              for(Index arow = 0; arow < remaining_rows; arow++)
              {
                  for(Index acol = 0; acol < accRows; acol++ )
                  {
                    res(row + arow, col + acol) += alpha*lhs_ptr[arow]*rhs_ptr[acol];
                  }
              }
              rhs_ptr += accRows;
              lhs_ptr += remaining_rows;
          }
        }
    }

    if(remaining_cols > 0)
    {
      const Scalar *rhs_base = blockB + (col/accRows)*strideB*accRows;
      const Scalar *lhs_base = blockA;

      Index row = 0;
      for(; row + accCols <= rows; row += accCols)
      {
        const Scalar *rhs_ptr = rhs_base;
        const Scalar *lhs_ptr = lhs_base + (row/accCols)*strideA*accCols;

        lhs_ptr += accCols*offsetA;
        rhs_ptr += remaining_cols*offsetB;
        for(Index k = 0; k < depth; k++)
        {
          for(Index arow = 0; arow < accCols; arow++)
          {
            for(Index acol = 0; acol < remaining_cols; acol++ )
            {
              res(row + arow, col + acol) += alpha*lhs_ptr[arow]*rhs_ptr[acol];
            }
          }
          rhs_ptr += remaining_cols;
          lhs_ptr += accCols;
        }
      }
      
      if(remaining_rows > 0 )
      {
        const Scalar *rhs_ptr  = rhs_base;
        const Scalar *lhs_ptr = lhs_base + (row/accCols)*strideA*accCols;

        lhs_ptr += remaining_rows*offsetA;
        rhs_ptr += remaining_cols*offsetB;
        for(Index k = 0; k < depth; k++)
        {
            for(Index arow = 0; arow < remaining_rows; arow++)
            {
                for(Index acol = 0; acol < remaining_cols; acol++ )
                {
                  res(row + arow, col + acol) += alpha*lhs_ptr[arow]*rhs_ptr[acol];
                }
            }
            rhs_ptr += remaining_cols;
            lhs_ptr += remaining_rows;
        }
      }
    }
}

template<typename LhsScalar, typename RhsScalar, typename Scalarc, typename Scalar, typename Index, typename Packet, typename Packetc, typename RhsPacket, typename DataMapper, bool ConjugateLhs, bool ConjugateRhs, bool LhsIsReal, bool RhsIsReal>
EIGEN_STRONG_INLINE void gemm_complex(const DataMapper& res, const LhsScalar* blockAc, const RhsScalar* blockBc,
          Index rows, Index depth, Index cols, Scalarc alpha, Index strideA, Index strideB, Index offsetA, Index offsetB, const int accRows, const int accCols)
{
      const int remaining_rows = rows % accCols;
      const int remaining_cols = cols % accRows;
      const int accColsC = accCols / 2;
      int advanceCols = 2;
      int advanceRows = 2;

      if(LhsIsReal) advanceRows = 1;
      if(RhsIsReal) advanceCols = 1;

      if( strideA == -1 ) strideA = depth;
      if( strideB == -1 ) strideB = depth;

      const Packet pAlphaReal = pset1<Packet>(alpha.real());
      const Packet pAlphaImag = pset1<Packet>(alpha.imag());

      const Scalar *blockA = (Scalar *) blockAc;
      const Scalar *blockB = (Scalar *) blockBc;

      Packet conj = pset1<Packet>((Scalar)-1.0f);

      Index col = 0;
      for(; col + accRows <= cols; col += accRows)
      {
        const Scalar *rhs_base = blockB + ( (advanceCols*col)/accRows     )*strideB*accRows;
        const Scalar *lhs_base = blockA;

        Index row = 0;
        for(; row + accCols <= rows; row += accCols)
        {
#define MICRO() \
            pgerc<Scalar, Packet, ConjugateLhs, ConjugateRhs, LhsIsReal, RhsIsReal>(accReal1, accImag1, rhs_ptr, rhs_ptr_imag, lhs_ptr1, lhs_ptr_imag1, conj); \
            lhs_ptr1 += accCols; \
            rhs_ptr += accRows; \
            if(!LhsIsReal) \
              lhs_ptr_imag1 += accCols; \
            if(!RhsIsReal) \
              rhs_ptr_imag += accRows;

          const Scalar *rhs_ptr  = rhs_base;
          const Scalar *rhs_ptr_imag = rhs_ptr + accRows*strideB;
          const Scalar *lhs_ptr1 = lhs_base + ((advanceRows*row)/accCols)*strideA*accCols;
          const Scalar *lhs_ptr_imag1 = lhs_ptr1 + accCols*strideA;

          PacketBlock<Packet,4> accReal1, accImag1;
          bsetzero<Scalar, Packet>(accReal1);
          bsetzero<Scalar, Packet>(accImag1);

          lhs_ptr1 += accCols*offsetA;
          if(!LhsIsReal)
            lhs_ptr_imag1 += accCols*offsetA;
          rhs_ptr += accRows*offsetB;
          if(!RhsIsReal)
            rhs_ptr_imag += accRows*offsetB;
          Index k = 0;
          for(; k + PEEL < depth; k+=PEEL)
          {
            prefetch(rhs_ptr);
            prefetch(rhs_ptr_imag);
            prefetch(lhs_ptr1);
            prefetch(lhs_ptr_imag1);
            MICRO();
            MICRO();
            MICRO();
            MICRO();
            MICRO();
            MICRO();
            MICRO();
            MICRO();
#if PEEL > 8
            MICRO();
            MICRO();
#endif
          }
          for(; k < depth; k++)
          {
            MICRO();
          }
          PacketBlock<Packet,4> taccReal, taccImag;
          bscalec<Packet>(accReal1, accImag1, pAlphaReal, pAlphaImag, taccReal, taccImag);

          PacketBlock<Packetc, 8> tRes;
          bload<DataMapper, Packetc, Index, 0>(tRes, res, row, col, accColsC);

          PacketBlock<Packetc, 4> acc1, acc2;
          bcouple<Packet, Packetc>(taccReal, taccImag, tRes, acc1, acc2);

          res.template storePacketBlock<Packetc, 4>(row + 0, col, acc1);
          res.template storePacketBlock<Packetc, 4>(row + accColsC, col, acc2);
#undef MICRO
      }
          if(remaining_rows > 0)
          {
            const Scalar *rhs_ptr  = rhs_base;
            const Scalar *rhs_ptr_imag = rhs_ptr + accRows*strideB;
            const Scalar *lhs_ptr = lhs_base + ((advanceRows*row)/accCols)*strideA*accCols;
            const Scalar *lhs_ptr_imag = lhs_ptr + remaining_rows*strideA;

            lhs_ptr += remaining_rows*offsetA;
            if(!LhsIsReal)
              lhs_ptr_imag += remaining_rows*offsetA;
            rhs_ptr += accRows*offsetB;
            if(!RhsIsReal)
              rhs_ptr_imag += accRows*offsetB;
            for(Index k = 0; k < depth; k++)
            {
              for(Index arow = 0; arow < remaining_rows; arow++)
              {
                Scalar lhs_real = lhs_ptr[arow];
                Scalar lhs_imag;
                if(!LhsIsReal) lhs_imag = lhs_ptr_imag[arow];

                Scalarc lhsc;

                lhsc.real(lhs_real);
                if(!LhsIsReal)
                {
                  if(ConjugateLhs)
                    lhsc.imag(-lhs_imag);
                  else
                    lhsc.imag(lhs_imag);
                } else {
                  //Lazy approach for now
                  lhsc.imag((Scalar)0);
                }

                for(int acol = 0; acol < accRows; acol++ )
                {
                  Scalar rhs_real = rhs_ptr[acol];
                  Scalar rhs_imag;
                  if(!RhsIsReal) rhs_imag = rhs_ptr_imag[acol];
                  Scalarc rhsc;

                  rhsc.real(rhs_real);
                  if(!RhsIsReal)
                  {
                    if(ConjugateRhs)
                      rhsc.imag(-rhs_imag);
                    else
                      rhsc.imag(rhs_imag);
                  } else {
                    //Lazy approach for now
                    rhsc.imag((Scalar)0);
                  }
                  res(row + arow, col + acol) += alpha*lhsc*rhsc;
                }
              }
              rhs_ptr += accRows;
              lhs_ptr += remaining_rows;
              if(!LhsIsReal)
                lhs_ptr_imag += remaining_rows;
              if(!RhsIsReal)
                rhs_ptr_imag += accRows;
            }
          }
      }

      if(remaining_cols > 0)
      {
        const Scalar *rhs_base = blockB + ( (advanceCols*col)/accRows     )*strideB*accRows;
        const Scalar *lhs_base = blockA;
        Index row = 0;

        for(; row + accCols <= rows; row += accCols)
        {
          const Scalar *rhs_ptr  = rhs_base;
          const Scalar *rhs_ptr_imag = rhs_ptr + remaining_cols*strideB;
          const Scalar *lhs_ptr = lhs_base + ((advanceRows*row)/accCols)*strideA*accCols;
          const Scalar *lhs_ptr_imag = lhs_ptr + accCols*strideA;

          lhs_ptr += accCols*offsetA;
          if(!LhsIsReal)
            lhs_ptr_imag += accCols*offsetA;
          rhs_ptr += remaining_cols*offsetB;
          if(!RhsIsReal)
            rhs_ptr_imag += remaining_cols*offsetB;
          Scalarc scalarAcc[4][4];
          for(Index arow = 0; arow < 4; arow++ )
          {
            for(Index acol = 0; acol < 4; acol++ )
            {
              scalarAcc[arow][acol].real((Scalar)0.0f);
              scalarAcc[arow][acol].imag((Scalar)0.0f);
            }
          }
          for(Index k = 0; k < depth; k++)
          {
            for(Index arow = 0; arow < accCols; arow++)
            {
              Scalar lhs_real = lhs_ptr[arow];
              Scalar lhs_imag;
              if(!LhsIsReal)
              {
                lhs_imag = lhs_ptr_imag[arow];

                if(ConjugateLhs)
                  lhs_imag *= -1;
              } else {
                lhs_imag = (Scalar)0;
              }

              for(int acol = 0; acol < remaining_cols; acol++ )
              {
                Scalar rhs_real = rhs_ptr[acol];
                Scalar rhs_imag;
                if(!RhsIsReal)
                {
                  rhs_imag = rhs_ptr_imag[acol];

                  if(ConjugateRhs)
                    rhs_imag *= -1;
                } else {
                  rhs_imag = (Scalar)0;
                }

                scalarAcc[arow][acol].real(scalarAcc[arow][acol].real() + lhs_real*rhs_real - lhs_imag*rhs_imag);
                scalarAcc[arow][acol].imag(scalarAcc[arow][acol].imag() + lhs_imag*rhs_real + lhs_real*rhs_imag);
              }
            }
            rhs_ptr += remaining_cols;
            lhs_ptr += accCols;
            if(!RhsIsReal)
              rhs_ptr_imag += remaining_cols;
            if(!LhsIsReal)
              lhs_ptr_imag += accCols;
          }
          for(int arow = 0; arow < accCols; arow++ )
          {
            for(int acol = 0; acol < remaining_cols; acol++ )
            {
              Scalar accR = scalarAcc[arow][acol].real();
              Scalar accI = scalarAcc[arow][acol].imag();
              Scalar aR = alpha.real();
              Scalar aI = alpha.imag();
              Scalar resR = res(row + arow, col + acol).real();
              Scalar resI = res(row + arow, col + acol).imag();

              res(row + arow, col + acol).real(resR + accR*aR - accI*aI);
              res(row + arow, col + acol).imag(resI + accR*aI + accI*aR);
            }
          }
        }

        if(remaining_rows > 0)
        {
          const Scalar *rhs_ptr  = rhs_base;
          const Scalar *rhs_ptr_imag = rhs_ptr + remaining_cols*strideB;
          const Scalar *lhs_ptr = lhs_base + ((advanceRows*row)/accCols)*strideA*accCols;
          const Scalar *lhs_ptr_imag = lhs_ptr + remaining_rows*strideA;

          lhs_ptr += remaining_rows*offsetA;
          if(!LhsIsReal)
            lhs_ptr_imag += remaining_rows*offsetA;
          rhs_ptr += remaining_cols*offsetB;
          if(!RhsIsReal)
            rhs_ptr_imag += remaining_cols*offsetB;
          for(Index k = 0; k < depth; k++)
          {
            for(Index arow = 0; arow < remaining_rows; arow++)
            {
              Scalar lhs_real = lhs_ptr[arow];
              Scalar lhs_imag;
              if(!LhsIsReal) lhs_imag = lhs_ptr_imag[arow];
              Scalarc lhsc;

              lhsc.real(lhs_real);
              if(!LhsIsReal)
              {
                if(ConjugateLhs)
                  lhsc.imag(-lhs_imag);
                else
                  lhsc.imag(lhs_imag);
              } else {
                lhsc.imag((Scalar)0);
              }

              for(Index acol = 0; acol < remaining_cols; acol++ )
              {
                Scalar rhs_real = rhs_ptr[acol];
                Scalar rhs_imag;
                if(!RhsIsReal) rhs_imag = rhs_ptr_imag[acol];
                Scalarc rhsc;

                rhsc.real(rhs_real);
                if(!RhsIsReal)
                {
                  if(ConjugateRhs)
                    rhsc.imag(-rhs_imag);
                  else
                    rhsc.imag(rhs_imag);
                } else {
                  rhsc.imag((Scalar)0);
                }
                res(row + arow, col + acol) += alpha*lhsc*rhsc;
              }
            }
            rhs_ptr += remaining_cols;
            lhs_ptr += remaining_rows;
            if(!LhsIsReal)
              lhs_ptr_imag += remaining_rows;
            if(!RhsIsReal)
              rhs_ptr_imag += remaining_cols;
          }
        }
      }
}

/************************************
 * ppc64le template specializations *
 * **********************************/
template<typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
struct gemm_pack_lhs<double, Index, DataMapper, Pack1, Pack2, Packet, ColMajor, Conjugate, PanelMode>
{
  void operator()(double* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride=0, Index offset=0);
};

template<typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
void gemm_pack_lhs<double, Index, DataMapper, Pack1, Pack2, Packet, ColMajor, Conjugate, PanelMode>
  ::operator()(double* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride, Index offset)
{
    lhs_pack<double, Index, DataMapper, Packet2d, ColMajor, PanelMode> pack;
    pack(blockA, lhs, depth, rows, stride, offset);
}

template<typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
struct gemm_pack_lhs<double, Index, DataMapper, Pack1, Pack2, Packet, RowMajor, Conjugate, PanelMode>
{
  void operator()(double* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride=0, Index offset=0);
};

template<typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
void gemm_pack_lhs<double, Index, DataMapper, Pack1, Pack2, Packet, RowMajor, Conjugate, PanelMode>
  ::operator()(double* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride, Index offset)
{
    lhs_pack<double, Index, DataMapper, Packet2d, RowMajor, PanelMode> pack;
    pack(blockA, lhs, depth, rows, stride, offset);
}

template<typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
struct gemm_pack_rhs<double, Index, DataMapper, nr, ColMajor, Conjugate, PanelMode>
{
  void operator()(double* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride=0, Index offset=0);
};

template<typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
void gemm_pack_rhs<double, Index, DataMapper, nr, ColMajor, Conjugate, PanelMode>
  ::operator()(double* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset)
{
  rhs_pack<double, Index, DataMapper, Packet2d, ColMajor, PanelMode> pack;
  pack(blockB, rhs, depth, cols, stride, offset);
}

template<typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
struct gemm_pack_rhs<double, Index, DataMapper, nr, RowMajor, Conjugate, PanelMode>
{
  void operator()(double* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride=0, Index offset=0);
};

template<typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
void gemm_pack_rhs<double, Index, DataMapper, nr, RowMajor, Conjugate, PanelMode>
  ::operator()(double* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset)
{
  rhs_pack<double, Index, DataMapper, Packet2d, RowMajor, PanelMode> pack;
  pack(blockB, rhs, depth, cols, stride, offset);
}

template<typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
struct gemm_pack_lhs<float, Index, DataMapper, Pack1, Pack2, Packet, RowMajor, Conjugate, PanelMode>
{
  void operator()(float* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride=0, Index offset=0);
};

template<typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
void gemm_pack_lhs<float, Index, DataMapper, Pack1, Pack2, Packet, RowMajor, Conjugate, PanelMode>
  ::operator()(float* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride, Index offset)
{
  lhs_pack<float, Index, DataMapper, Packet4f, RowMajor, PanelMode> pack;
  pack(blockA, lhs, depth, rows, stride, offset);
}

template<typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
struct gemm_pack_lhs<float, Index, DataMapper, Pack1, Pack2, Packet, ColMajor, Conjugate, PanelMode>
{
  void operator()(float* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride=0, Index offset=0);
};

template<typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
void gemm_pack_lhs<float, Index, DataMapper, Pack1, Pack2, Packet, ColMajor, Conjugate, PanelMode>
  ::operator()(float* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride, Index offset)
{
  lhs_pack<float, Index, DataMapper, Packet4f, ColMajor, PanelMode> pack;
  pack(blockA, lhs, depth, rows, stride, offset);
}
template<typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
struct gemm_pack_lhs<std::complex<float>, Index, DataMapper, Pack1, Pack2, Packet, RowMajor, Conjugate, PanelMode>
{
  void operator()(std::complex<float>* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride=0, Index offset=0);
};

template<typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
void gemm_pack_lhs<std::complex<float>, Index, DataMapper, Pack1, Pack2, Packet, RowMajor, Conjugate, PanelMode>
  ::operator()(std::complex<float>* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride, Index offset)
{
  lhs_cpack<float, true, Index, DataMapper, Packet4f, Packet2cf, RowMajor, Conjugate, PanelMode> pack;
  pack(blockA, lhs, depth, rows, stride, offset);
}

template<typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
struct gemm_pack_lhs<std::complex<float>, Index, DataMapper, Pack1, Pack2, Packet, ColMajor, Conjugate, PanelMode>
{
  void operator()(std::complex<float>* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride=0, Index offset=0);
};

template<typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
void gemm_pack_lhs<std::complex<float>, Index, DataMapper, Pack1, Pack2, Packet, ColMajor, Conjugate, PanelMode>
  ::operator()(std::complex<float>* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride, Index offset)
{
  lhs_cpack<float, true, Index, DataMapper, Packet4f, Packet2cf, ColMajor, Conjugate, PanelMode> pack;
  pack(blockA, lhs, depth, rows, stride, offset);
}

template<typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
struct gemm_pack_rhs<float, Index, DataMapper, nr, ColMajor, Conjugate, PanelMode>
{
  void operator()(float* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride=0, Index offset=0);
};

template<typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
void gemm_pack_rhs<float, Index, DataMapper, nr, ColMajor, Conjugate, PanelMode>
  ::operator()(float* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset)
{
  rhs_pack<float, Index, DataMapper, Packet4f, ColMajor, PanelMode> pack;
  pack(blockB, rhs, depth, cols, stride, offset);
}

template<typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
struct gemm_pack_rhs<float, Index, DataMapper, nr, RowMajor, Conjugate, PanelMode>
{
  void operator()(float* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride=0, Index offset=0);
};

template<typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
void gemm_pack_rhs<float, Index, DataMapper, nr, RowMajor, Conjugate, PanelMode>
  ::operator()(float* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset)
{
  rhs_pack<float, Index, DataMapper, Packet4f, RowMajor, PanelMode> pack;
  pack(blockB, rhs, depth, cols, stride, offset);
}

template<typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
struct gemm_pack_rhs<std::complex<float>, Index, DataMapper, nr, ColMajor, Conjugate, PanelMode>
{
  void operator()(std::complex<float>* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride=0, Index offset=0);
};

template<typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
void gemm_pack_rhs<std::complex<float>, Index, DataMapper, nr, ColMajor, Conjugate, PanelMode>
  ::operator()(std::complex<float>* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset)
{
  rhs_cpack<float, Index, DataMapper, Packet4f, Packet2cf, ColMajor, Conjugate, PanelMode> pack;
  pack(blockB, rhs, depth, cols, stride, offset);
}

template<typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
struct gemm_pack_rhs<std::complex<float>, Index, DataMapper, nr, RowMajor, Conjugate, PanelMode>
{
  void operator()(std::complex<float>* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride=0, Index offset=0);
};

template<typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
void gemm_pack_rhs<std::complex<float>, Index, DataMapper, nr, RowMajor, Conjugate, PanelMode>
  ::operator()(std::complex<float>* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset)
{
  rhs_cpack<float, Index, DataMapper, Packet4f, Packet2cf, RowMajor, Conjugate, PanelMode> pack;
  pack(blockB, rhs, depth, cols, stride, offset);
}

template<typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
struct gemm_pack_lhs<std::complex<double>, Index, DataMapper, Pack1, Pack2, Packet, RowMajor, Conjugate, PanelMode>
{
  void operator()(std::complex<double>* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride=0, Index offset=0);
};

template<typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
void gemm_pack_lhs<std::complex<double>, Index, DataMapper, Pack1, Pack2, Packet, RowMajor, Conjugate, PanelMode>
  ::operator()(std::complex<double>* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride, Index offset)
{
  lhs_cpack<double, true, Index, DataMapper, Packet2d, Packet1cd, RowMajor, Conjugate, PanelMode> pack;
  pack(blockA, lhs, depth, rows, stride, offset);
}

template<typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
struct gemm_pack_lhs<std::complex<double>, Index, DataMapper, Pack1, Pack2, Packet, ColMajor, Conjugate, PanelMode>
{
  void operator()(std::complex<double>* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride=0, Index offset=0);
};

template<typename Index, typename DataMapper, int Pack1, int Pack2, typename Packet, bool Conjugate, bool PanelMode>
void gemm_pack_lhs<std::complex<double>, Index, DataMapper, Pack1, Pack2, Packet, ColMajor, Conjugate, PanelMode>
  ::operator()(std::complex<double>* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride, Index offset)
{
  lhs_cpack<double, true, Index, DataMapper, Packet2d, Packet1cd, ColMajor, Conjugate, PanelMode> pack;
  pack(blockA, lhs, depth, rows, stride, offset);
}

template<typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
struct gemm_pack_rhs<std::complex<double>, Index, DataMapper, nr, ColMajor, Conjugate, PanelMode>
{
  void operator()(std::complex<double>* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride=0, Index offset=0);
};

template<typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
void gemm_pack_rhs<std::complex<double>, Index, DataMapper, nr, ColMajor, Conjugate, PanelMode>
  ::operator()(std::complex<double>* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset)
{
  rhs_cpack<double, Index, DataMapper, Packet2d, Packet1cd, ColMajor, Conjugate, PanelMode> pack;
  pack(blockB, rhs, depth, cols, stride, offset);
}

template<typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
struct gemm_pack_rhs<std::complex<double>, Index, DataMapper, nr, RowMajor, Conjugate, PanelMode>
{
  void operator()(std::complex<double>* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride=0, Index offset=0);
};

template<typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
void gemm_pack_rhs<std::complex<double>, Index, DataMapper, nr, RowMajor, Conjugate, PanelMode>
  ::operator()(std::complex<double>* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset)
{
  rhs_cpack<double, Index, DataMapper, Packet2d, Packet1cd, RowMajor, Conjugate, PanelMode> pack;
  pack(blockB, rhs, depth, cols, stride, offset);
}

// ********* gebp specializations *********
template<typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs>
struct gebp_kernel<float, float, Index, DataMapper, mr, nr, ConjugateLhs, ConjugateRhs>
{
  typedef typename quad_traits<float>::vectortype   Packet;
  typedef typename quad_traits<float>::rhstype      RhsPacket;

  void operator()(const DataMapper& res, const float* blockA, const float* blockB,
                  Index rows, Index depth, Index cols, float alpha,
                  Index strideA=-1, Index strideB=-1, Index offsetA=0, Index offsetB=0);
};

template<typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs>
void gebp_kernel<float, float, Index, DataMapper, mr, nr, ConjugateLhs, ConjugateRhs>
  ::operator()(const DataMapper& res, const float* blockA, const float* blockB,
               Index rows, Index depth, Index cols, float alpha,
               Index strideA, Index strideB, Index offsetA, Index offsetB)
  {
    const int accRows = quad_traits<float>::rows;
    const int accCols = quad_traits<float>::size;
    void (*gemm_function)(const DataMapper&, const float*, const float*, Index, Index, Index, float, Index, Index, Index, Index, const int, const int);

    #ifdef EIGEN_ALTIVEC_MMA_ONLY
      //generate with MMA only
      gemm_function = &Eigen::internal::gemmMMA<float, Index, Packet, RhsPacket, DataMapper>;
    #elif defined(ALTIVEC_MMA_SUPPORT) && !defined(EIGEN_ALTIVEC_DISABLE_MMA)
      if (__builtin_cpu_supports ("arch_3_1") && __builtin_cpu_supports ("mma")){
        gemm_function = &Eigen::internal::gemmMMA<float, Index, Packet, RhsPacket, DataMapper>;
      }
      else{
        gemm_function = &Eigen::internal::gemm<float, Index, Packet, RhsPacket, DataMapper>;
      }
    #else
      gemm_function = &Eigen::internal::gemm<float, Index, Packet, RhsPacket, DataMapper>;
    #endif
      gemm_function(res, blockA, blockB, rows, depth, cols, alpha, strideA, strideB, offsetA, offsetB, accRows, accCols);
  }

template<typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs>
struct gebp_kernel<std::complex<float>, std::complex<float>, Index, DataMapper, mr, nr, ConjugateLhs, ConjugateRhs>
{
  typedef Packet4f   Packet;
  typedef Packet2cf  Packetc;
  typedef Packet4f   RhsPacket;

  void operator()(const DataMapper& res, const std::complex<float>* blockA, const std::complex<float>* blockB,
                  Index rows, Index depth, Index cols, std::complex<float> alpha,
                  Index strideA=-1, Index strideB=-1, Index offsetA=0, Index offsetB=0);
};

template<typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs>
void gebp_kernel<std::complex<float>, std::complex<float>, Index, DataMapper, mr, nr, ConjugateLhs, ConjugateRhs>
  ::operator()(const DataMapper& res, const std::complex<float>* blockA, const std::complex<float>* blockB,
               Index rows, Index depth, Index cols, std::complex<float> alpha,
               Index strideA, Index strideB, Index offsetA, Index offsetB)
  {
    const int accRows = quad_traits<float>::rows;
    const int accCols = quad_traits<float>::size;
    void (*gemm_function)(const DataMapper&, const std::complex<float>*, const std::complex<float>*,
          Index, Index, Index, std::complex<float>, Index, Index , Index, Index, const int, const int);

    #ifdef EIGEN_ALTIVEC_MMA_ONLY
       //generate with MMA only
       gemm_function = &Eigen::internal::gemm_complexMMA<std::complex<float>, std::complex<float>, std::complex<float>, float, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, false>;
     #elif defined(ALTIVEC_MMA_SUPPORT) && !defined(EIGEN_ALTIVEC_DISABLE_MMA)
       if (__builtin_cpu_supports ("arch_3_1") && __builtin_cpu_supports ("mma")){
         gemm_function = &Eigen::internal::gemm_complexMMA<std::complex<float>, std::complex<float>, std::complex<float>, float, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, false>;
       }
       else{
         gemm_function = &Eigen::internal::gemm_complex<std::complex<float>, std::complex<float>, std::complex<float>, float, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, false>;
       }
     #else
       gemm_function = &Eigen::internal::gemm_complex<std::complex<float>, std::complex<float>, std::complex<float>, float, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, false>;
     #endif
       gemm_function(res, blockA, blockB, rows, depth, cols, alpha, strideA, strideB, offsetA, offsetB, accRows, accCols);
  }

template<typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs>
struct gebp_kernel<float, std::complex<float>, Index, DataMapper, mr, nr, ConjugateLhs, ConjugateRhs>
{
  typedef Packet4f   Packet;
  typedef Packet2cf  Packetc;
  typedef Packet4f   RhsPacket;

  void operator()(const DataMapper& res, const float* blockA, const std::complex<float>* blockB,
                  Index rows, Index depth, Index cols, std::complex<float> alpha,
                  Index strideA=-1, Index strideB=-1, Index offsetA=0, Index offsetB=0);
};

template<typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs>
void gebp_kernel<float, std::complex<float>, Index, DataMapper, mr, nr, ConjugateLhs, ConjugateRhs>
  ::operator()(const DataMapper& res, const float* blockA, const std::complex<float>* blockB,
               Index rows, Index depth, Index cols, std::complex<float> alpha,
               Index strideA, Index strideB, Index offsetA, Index offsetB)
  {
    const int accRows = quad_traits<float>::rows;
    const int accCols = quad_traits<float>::size;
    void (*gemm_function)(const DataMapper&, const float*, const std::complex<float>*,
          Index, Index, Index, std::complex<float>, Index, Index , Index, Index, const int, const int);
    #ifdef EIGEN_ALTIVEC_MMA_ONLY
       //generate with MMA only
       gemm_function = &Eigen::internal::gemm_complexMMA<float, std::complex<float>, std::complex<float>, float, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, true, false>;
     #elif defined(ALTIVEC_MMA_SUPPORT) && !defined(EIGEN_ALTIVEC_DISABLE_MMA)
       if (__builtin_cpu_supports ("arch_3_1") && __builtin_cpu_supports ("mma")){
         gemm_function = &Eigen::internal::gemm_complexMMA<float, std::complex<float>, std::complex<float>, float, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, true, false>;
       }
       else{
         gemm_function = &Eigen::internal::gemm_complex<float, std::complex<float>, std::complex<float>, float, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, true, false>;
       }
     #else
       gemm_function = &Eigen::internal::gemm_complex<float, std::complex<float>, std::complex<float>, float, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, true, false>;
     #endif
       gemm_function(res, blockA, blockB, rows, depth, cols, alpha, strideA, strideB, offsetA, offsetB, accRows, accCols);
  }

template<typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs>
struct gebp_kernel<std::complex<float>, float, Index, DataMapper, mr, nr, ConjugateLhs, ConjugateRhs>
{
  typedef Packet4f   Packet;
  typedef Packet2cf  Packetc;
  typedef Packet4f   RhsPacket;

  void operator()(const DataMapper& res, const std::complex<float>* blockA, const float* blockB,
                  Index rows, Index depth, Index cols, std::complex<float> alpha,
                  Index strideA=-1, Index strideB=-1, Index offsetA=0, Index offsetB=0);
};

template<typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs>
void gebp_kernel<std::complex<float>, float, Index, DataMapper, mr, nr, ConjugateLhs, ConjugateRhs>
  ::operator()(const DataMapper& res, const std::complex<float>* blockA, const float* blockB,
               Index rows, Index depth, Index cols, std::complex<float> alpha,
               Index strideA, Index strideB, Index offsetA, Index offsetB)
  {
    const int accRows = quad_traits<float>::rows;
    const int accCols = quad_traits<float>::size;
    void (*gemm_function)(const DataMapper&, const std::complex<float>*, const float*,
          Index, Index, Index, std::complex<float>, Index, Index , Index, Index, const int, const int);
    #ifdef EIGEN_ALTIVEC_MMA_ONLY
       //generate with MMA only
       gemm_function = &Eigen::internal::gemm_complexMMA<std::complex<float>, float, std::complex<float>, float, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, true>;
     #elif defined(ALTIVEC_MMA_SUPPORT) && !defined(EIGEN_ALTIVEC_DISABLE_MMA)
       if (__builtin_cpu_supports ("arch_3_1") && __builtin_cpu_supports ("mma")){
         gemm_function = &Eigen::internal::gemm_complexMMA<std::complex<float>, float, std::complex<float>, float, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, true>;
       }
       else{
         gemm_function = &Eigen::internal::gemm_complex<std::complex<float>, float, std::complex<float>, float, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, true>;
       }
     #else
       gemm_function = &Eigen::internal::gemm_complex<std::complex<float>, float, std::complex<float>, float, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, true>;
     #endif
       gemm_function(res, blockA, blockB, rows, depth, cols, alpha, strideA, strideB, offsetA, offsetB, accRows, accCols);
  }

template<typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs>
struct gebp_kernel<double, double, Index, DataMapper, mr, nr, ConjugateLhs, ConjugateRhs>
{
  typedef typename quad_traits<double>::vectortype  Packet;
  typedef typename quad_traits<double>::rhstype     RhsPacket;

  void operator()(const DataMapper& res, const double* blockA, const double* blockB,
                  Index rows, Index depth, Index cols, double alpha,
                  Index strideA=-1, Index strideB=-1, Index offsetA=0, Index offsetB=0);
};

template<typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs>
void gebp_kernel<double, double, Index, DataMapper, mr, nr, ConjugateLhs, ConjugateRhs>
  ::operator()(const DataMapper& res, const double* blockA, const double* blockB,
               Index rows, Index depth, Index cols, double alpha,
               Index strideA, Index strideB, Index offsetA, Index offsetB)
  {
    const int accRows = quad_traits<double>::rows;
    const int accCols = quad_traits<double>::size;
    void (*gemm_function)(const DataMapper&, const double*, const double*, Index, Index, Index, double, Index, Index, Index, Index, const int, const int);

    #ifdef EIGEN_ALTIVEC_MMA_ONLY
      //generate with MMA only
      gemm_function = &Eigen::internal::gemmMMA<double, Index, Packet, RhsPacket, DataMapper>;
    #elif defined(ALTIVEC_MMA_SUPPORT) && !defined(EIGEN_ALTIVEC_DISABLE_MMA)
      if (__builtin_cpu_supports ("arch_3_1") && __builtin_cpu_supports ("mma")){
        gemm_function = &Eigen::internal::gemmMMA<double, Index, Packet, RhsPacket, DataMapper>;
      }
      else{
        gemm_function = &Eigen::internal::gemm<double, Index, Packet, RhsPacket, DataMapper>;
      }
    #else
      gemm_function = &Eigen::internal::gemm<double, Index, Packet, RhsPacket, DataMapper>;
    #endif
      gemm_function(res, blockA, blockB, rows, depth, cols, alpha, strideA, strideB, offsetA, offsetB, accRows, accCols);
  }

template<typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs>
struct gebp_kernel<std::complex<double>, std::complex<double>, Index, DataMapper, mr, nr, ConjugateLhs, ConjugateRhs>
{
  typedef quad_traits<double>::vectortype   Packet;
  typedef Packet1cd  Packetc;
  typedef quad_traits<double>::rhstype   RhsPacket;

  void operator()(const DataMapper& res, const std::complex<double>* blockA, const std::complex<double>* blockB,
                  Index rows, Index depth, Index cols, std::complex<double> alpha,
                  Index strideA=-1, Index strideB=-1, Index offsetA=0, Index offsetB=0);
};

template<typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs>
void gebp_kernel<std::complex<double>, std::complex<double>, Index, DataMapper, mr, nr, ConjugateLhs, ConjugateRhs>
  ::operator()(const DataMapper& res, const std::complex<double>* blockA, const std::complex<double>* blockB,
               Index rows, Index depth, Index cols, std::complex<double> alpha,
               Index strideA, Index strideB, Index offsetA, Index offsetB)
  {
    const int accRows = quad_traits<double>::rows;
    const int accCols = quad_traits<double>::size;
    void (*gemm_function)(const DataMapper&, const std::complex<double>*, const std::complex<double>*,
          Index, Index, Index, std::complex<double>, Index, Index , Index, Index, const int, const int);
    #ifdef EIGEN_ALTIVEC_MMA_ONLY
       //generate with MMA only
       gemm_function = &Eigen::internal::gemm_complexMMA<std::complex<double>, std::complex<double>, std::complex<double>, double, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, false>;
     #elif defined(ALTIVEC_MMA_SUPPORT) && !defined(EIGEN_ALTIVEC_DISABLE_MMA)
       if (__builtin_cpu_supports ("arch_3_1") && __builtin_cpu_supports ("mma")){
         gemm_function = &Eigen::internal::gemm_complexMMA<std::complex<double>, std::complex<double>, std::complex<double>, double, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, false>;
       }
       else{
         gemm_function = &Eigen::internal::gemm_complex<std::complex<double>, std::complex<double>, std::complex<double>, double, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, false>;
       }
     #else
       gemm_function = &Eigen::internal::gemm_complex<std::complex<double>, std::complex<double>, std::complex<double>, double, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, false>;
     #endif
       gemm_function(res, blockA, blockB, rows, depth, cols, alpha, strideA, strideB, offsetA, offsetB, accRows, accCols);
  }

template<typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs>
struct gebp_kernel<std::complex<double>, double, Index, DataMapper, mr, nr, ConjugateLhs, ConjugateRhs>
{
  typedef quad_traits<double>::vectortype   Packet;
  typedef Packet1cd  Packetc;
  typedef quad_traits<double>::rhstype   RhsPacket;

  void operator()(const DataMapper& res, const std::complex<double>* blockA, const double* blockB,
                  Index rows, Index depth, Index cols, std::complex<double> alpha,
                  Index strideA=-1, Index strideB=-1, Index offsetA=0, Index offsetB=0);
};

template<typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs>
void gebp_kernel<std::complex<double>, double, Index, DataMapper, mr, nr, ConjugateLhs, ConjugateRhs>
  ::operator()(const DataMapper& res, const std::complex<double>* blockA, const double* blockB,
               Index rows, Index depth, Index cols, std::complex<double> alpha,
               Index strideA, Index strideB, Index offsetA, Index offsetB)
  {
    const int accRows = quad_traits<double>::rows;
    const int accCols = quad_traits<double>::size;
    void (*gemm_function)(const DataMapper&, const std::complex<double>*, const double*,
          Index, Index, Index, std::complex<double>, Index, Index , Index, Index, const int, const int);
    #ifdef EIGEN_ALTIVEC_MMA_ONLY
       //generate with MMA only
       gemm_function = &Eigen::internal::gemm_complexMMA<std::complex<double>, double, std::complex<double>, double, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, true>;
     #elif defined(ALTIVEC_MMA_SUPPORT) && !defined(EIGEN_ALTIVEC_DISABLE_MMA)
       if (__builtin_cpu_supports ("arch_3_1") && __builtin_cpu_supports ("mma")){
         gemm_function = &Eigen::internal::gemm_complexMMA<std::complex<double>, double, std::complex<double>, double, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, true>;
       }
       else{
         gemm_function = &Eigen::internal::gemm_complex<std::complex<double>, double, std::complex<double>, double, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, true>;
       }
     #else
       gemm_function = &Eigen::internal::gemm_complex<std::complex<double>, double, std::complex<double>, double, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, true>;
     #endif
       gemm_function(res, blockA, blockB, rows, depth, cols, alpha, strideA, strideB, offsetA, offsetB, accRows, accCols);
  }

template<typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs>
struct gebp_kernel<double, std::complex<double>, Index, DataMapper, mr, nr, ConjugateLhs, ConjugateRhs>
{
  typedef quad_traits<double>::vectortype   Packet;
  typedef Packet1cd  Packetc;
  typedef quad_traits<double>::rhstype   RhsPacket;

  void operator()(const DataMapper& res, const double* blockA, const std::complex<double>* blockB,
                  Index rows, Index depth, Index cols, std::complex<double> alpha,
                  Index strideA=-1, Index strideB=-1, Index offsetA=0, Index offsetB=0);
};

template<typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs>
void gebp_kernel<double, std::complex<double>, Index, DataMapper, mr, nr, ConjugateLhs, ConjugateRhs>
  ::operator()(const DataMapper& res, const double* blockA, const std::complex<double>* blockB,
               Index rows, Index depth, Index cols, std::complex<double> alpha,
               Index strideA, Index strideB, Index offsetA, Index offsetB)
  {
    const int accRows = quad_traits<double>::rows;
    const int accCols = quad_traits<double>::size;
    void (*gemm_function)(const DataMapper&, const double*, const std::complex<double>*,
          Index, Index, Index, std::complex<double>, Index, Index , Index, Index, const int, const int);
    #ifdef EIGEN_ALTIVEC_MMA_ONLY
       //generate with MMA only
       gemm_function = &Eigen::internal::gemm_complexMMA<double, std::complex<double>, std::complex<double>, double, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, true, false>;
     #elif defined(ALTIVEC_MMA_SUPPORT) && !defined(EIGEN_ALTIVEC_DISABLE_MMA)
       if (__builtin_cpu_supports ("arch_3_1") && __builtin_cpu_supports ("mma")){
         gemm_function = &Eigen::internal::gemm_complexMMA<double, std::complex<double>, std::complex<double>, double, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, true, false>;
       }
       else{
         gemm_function = &Eigen::internal::gemm_complex<double, std::complex<double>, std::complex<double>, double, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, true, false>;
       }
     #else
       gemm_function = &Eigen::internal::gemm_complex<double, std::complex<double>, std::complex<double>, double, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, true, false>;
     #endif
       gemm_function(res, blockA, blockB, rows, depth, cols, alpha, strideA, strideB, offsetA, offsetB, accRows, accCols);
  }
} // end namespace internal

} // end namespace Eigen

#endif // EIGEN_MATRIX_PRODUCT_ALTIVEC_H