fully adapt the gebp kernel and optimize it for CPU with only 8 registers

This commit is contained in:
Gael Guennebaud 2010-02-22 16:35:05 +01:00
parent dd569c7d0f
commit fc4a85ecd5

View File

@ -27,6 +27,12 @@
#ifndef EIGEN_EXTERN_INSTANTIATIONS
#ifdef EIGEN_HAS_FUSE_CJMADD
#define CJMADD(A,B,C,T) C = cj.pmadd(A,B,C);
#else
#define CJMADD(A,B,C,T) T = A; T = cj.pmul(T,B); C = ei_padd(C,T);
#endif
// optimized GEneral packed Block * packed Panel product kernel
template<typename Scalar, int mr, int nr, typename Conj>
struct ei_gebp_kernel
@ -74,133 +80,111 @@ struct ei_gebp_kernel
const Scalar* blB = &blockB[j2*strideB*PacketSize+offsetB*nr];
for(int k=0; k<peeled_kc; k+=4)
{
PacketType B0, B1, B2, B3, A0, A1;
if(nr==2)
{
PacketType B0, T0, A0, A1;
#if 0
A0 = ei_pload(&blA[0*PacketSize]);
A1 = ei_pload(&blA[1*PacketSize]);
B0 = ei_pload(&blB[0*PacketSize]);
B1 = ei_pload(&blB[1*PacketSize]);
C0 = cj.pmadd(A0, B0, C0);
if(nr==4) B2 = ei_pload(&blB[2*PacketSize]);
C4 = cj.pmadd(A1, B0, C4);
if(nr==4) B3 = ei_pload(&blB[3*PacketSize]);
B0 = ei_pload(&blB[(nr==4 ? 4 : 2)*PacketSize]);
C1 = cj.pmadd(A0, B1, C1);
C5 = cj.pmadd(A1, B1, C5);
B1 = ei_pload(&blB[(nr==4 ? 5 : 3)*PacketSize]);
if(nr==4) C2 = cj.pmadd(A0, B2, C2);
if(nr==4) C6 = cj.pmadd(A1, B2, C6);
if(nr==4) B2 = ei_pload(&blB[6*PacketSize]);
if(nr==4) C3 = cj.pmadd(A0, B3, C3);
CJMADD(A0,B0,C0,T0);
CJMADD(A1,B0,C4,T0);
B0 = ei_pload(&blB[1*PacketSize]);
CJMADD(A0,B0,C1,T0);
CJMADD(A1,B0,C5,T0);
A0 = ei_pload(&blA[2*PacketSize]);
if(nr==4) C7 = cj.pmadd(A1, B3, C7);
A1 = ei_pload(&blA[3*PacketSize]);
if(nr==4) B3 = ei_pload(&blB[7*PacketSize]);
C0 = cj.pmadd(A0, B0, C0);
C4 = cj.pmadd(A1, B0, C4);
B0 = ei_pload(&blB[(nr==4 ? 8 : 4)*PacketSize]);
C1 = cj.pmadd(A0, B1, C1);
C5 = cj.pmadd(A1, B1, C5);
B1 = ei_pload(&blB[(nr==4 ? 9 : 5)*PacketSize]);
if(nr==4) C2 = cj.pmadd(A0, B2, C2);
if(nr==4) C6 = cj.pmadd(A1, B2, C6);
if(nr==4) B2 = ei_pload(&blB[10*PacketSize]);
if(nr==4) C3 = cj.pmadd(A0, B3, C3);
B0 = ei_pload(&blB[2*PacketSize]);
CJMADD(A0,B0,C0,T0);
CJMADD(A1,B0,C4,T0);
B0 = ei_pload(&blB[3*PacketSize]);
CJMADD(A0,B0,C1,T0);
CJMADD(A1,B0,C5,T0);
A0 = ei_pload(&blA[4*PacketSize]);
if(nr==4) C7 = cj.pmadd(A1, B3, C7);
A1 = ei_pload(&blA[5*PacketSize]);
if(nr==4) B3 = ei_pload(&blB[11*PacketSize]);
B0 = ei_pload(&blB[4*PacketSize]);
CJMADD(A0,B0,C0,T0);
CJMADD(A1,B0,C4,T0);
B0 = ei_pload(&blB[5*PacketSize]);
CJMADD(A0,B0,C1,T0);
CJMADD(A1,B0,C5,T0);
C0 = cj.pmadd(A0, B0, C0);
C4 = cj.pmadd(A1, B0, C4);
B0 = ei_pload(&blB[(nr==4 ? 12 : 6)*PacketSize]);
C1 = cj.pmadd(A0, B1, C1);
C5 = cj.pmadd(A1, B1, C5);
B1 = ei_pload(&blB[(nr==4 ? 13 : 7)*PacketSize]);
if(nr==4) C2 = cj.pmadd(A0, B2, C2);
if(nr==4) C6 = cj.pmadd(A1, B2, C6);
if(nr==4) B2 = ei_pload(&blB[14*PacketSize]);
if(nr==4) C3 = cj.pmadd(A0, B3, C3);
A0 = ei_pload(&blA[6*PacketSize]);
if(nr==4) C7 = cj.pmadd(A1, B3, C7);
A1 = ei_pload(&blA[7*PacketSize]);
if(nr==4) B3 = ei_pload(&blB[15*PacketSize]);
C0 = cj.pmadd(A0, B0, C0);
C4 = cj.pmadd(A1, B0, C4);
C1 = cj.pmadd(A0, B1, C1);
C5 = cj.pmadd(A1, B1, C5);
if(nr==4) C2 = cj.pmadd(A0, B2, C2);
if(nr==4) C6 = cj.pmadd(A1, B2, C6);
if(nr==4) C3 = cj.pmadd(A0, B3, C3);
if(nr==4) C7 = cj.pmadd(A1, B3, C7);
#else
B0 = ei_pload(&blB[6*PacketSize]);
CJMADD(A0,B0,C0,T0);
CJMADD(A1,B0,C4,T0);
B0 = ei_pload(&blB[7*PacketSize]);
CJMADD(A0,B0,C1,T0);
CJMADD(A1,B0,C5,T0);
}
else
{
PacketType B0, B1, B2, B3, A0, A1;
PacketType T0, T1;
#define MADD(A,B,C,T) { T = A; T = ei_pmul(T,B); C = ei_padd(C,T); }
A0 = ei_pload(&blA[0*PacketSize]);
A1 = ei_pload(&blA[1*PacketSize]);
B0 = ei_pload(&blB[0*PacketSize]);
B1 = ei_pload(&blB[1*PacketSize]);
MADD(A0,B0,C0,T0); // C0 = cj.pmadd(A0, B0, C0);
CJMADD(A0,B0,C0,T0);
if(nr==4) B2 = ei_pload(&blB[2*PacketSize]);
MADD(A1,B0,C4,T1); // C4 = cj.pmadd(A1, B0, C4);
CJMADD(A1,B0,C4,T1);
if(nr==4) B3 = ei_pload(&blB[3*PacketSize]);
B0 = ei_pload(&blB[(nr==4 ? 4 : 2)*PacketSize]);
MADD(A0,B1,C1,T0); // C1 = cj.pmadd(A0, B1, C1);
MADD(A1,B1,C5,T1); // C5 = cj.pmadd(A1, B1, C5);
CJMADD(A0,B1,C1,T0);
CJMADD(A1,B1,C5,T1);
B1 = ei_pload(&blB[(nr==4 ? 5 : 3)*PacketSize]);
if(nr==4) { MADD(A0,B2,C2,T0); }// C2 = cj.pmadd(A0, B2, C2);
if(nr==4) { MADD(A1,B2,C6,T1); }// C6 = cj.pmadd(A1, B2, C6);
if(nr==4) { CJMADD(A0,B2,C2,T0); }
if(nr==4) { CJMADD(A1,B2,C6,T1); }
if(nr==4) B2 = ei_pload(&blB[6*PacketSize]);
if(nr==4) { MADD(A0,B3,C3,T0); }// C3 = cj.pmadd(A0, B3, C3);
if(nr==4) { CJMADD(A0,B3,C3,T0); }
A0 = ei_pload(&blA[2*PacketSize]);
if(nr==4) { MADD(A1,B3,C7,T1); }// C7 = cj.pmadd(A1, B3, C7);
if(nr==4) { CJMADD(A1,B3,C7,T1); }
A1 = ei_pload(&blA[3*PacketSize]);
if(nr==4) B3 = ei_pload(&blB[7*PacketSize]);
MADD(A0,B0,C0,T0); // C0 = cj.pmadd(A0, B0, C0);
MADD(A1,B0,C4,T1); // C4 = cj.pmadd(A1, B0, C4);
CJMADD(A0,B0,C0,T0);
CJMADD(A1,B0,C4,T1);
B0 = ei_pload(&blB[(nr==4 ? 8 : 4)*PacketSize]);
MADD(A0,B1,C1,T0); // C1 = cj.pmadd(A0, B1, C1);
MADD(A1,B1,C5,T1); // C5 = cj.pmadd(A1, B1, C5);
CJMADD(A0,B1,C1,T0);
CJMADD(A1,B1,C5,T1);
B1 = ei_pload(&blB[(nr==4 ? 9 : 5)*PacketSize]);
if(nr==4) { MADD(A0,B2,C2,T0); }// C2 = cj.pmadd(A0, B2, C2);
if(nr==4) { MADD(A1,B2,C6,T1); }// C6 = cj.pmadd(A1, B2, C6);
if(nr==4) { CJMADD(A0,B2,C2,T0); }
if(nr==4) { CJMADD(A1,B2,C6,T1); }
if(nr==4) B2 = ei_pload(&blB[10*PacketSize]);
if(nr==4) { MADD(A0,B3,C3,T0); } // C3 = cj.pmadd(A0, B3, C3);
if(nr==4) { CJMADD(A0,B3,C3,T0); }
A0 = ei_pload(&blA[4*PacketSize]);
if(nr==4) { MADD(A1,B3,C7,T1); }// C7 = cj.pmadd(A1, B3, C7);
if(nr==4) { CJMADD(A1,B3,C7,T1); }
A1 = ei_pload(&blA[5*PacketSize]);
if(nr==4) B3 = ei_pload(&blB[11*PacketSize]);
MADD(A0,B0,C0,T0); // C0 = cj.pmadd(A0, B0, C0);
MADD(A1,B0,C4,T1); // C4 = cj.pmadd(A1, B0, C4);
CJMADD(A0,B0,C0,T0);
CJMADD(A1,B0,C4,T1);
B0 = ei_pload(&blB[(nr==4 ? 12 : 6)*PacketSize]);
MADD(A0,B1,C1,T0); // C1 = cj.pmadd(A0, B1, C1);
MADD(A1,B1,C5,T1); // C5 = cj.pmadd(A1, B1, C5);
CJMADD(A0,B1,C1,T0);
CJMADD(A1,B1,C5,T1);
B1 = ei_pload(&blB[(nr==4 ? 13 : 7)*PacketSize]);
if(nr==4) { MADD(A0,B2,C2,T0); }// C2 = cj.pmadd(A0, B2, C2);
if(nr==4) { MADD(A1,B2,C6,T1); }// C6 = cj.pmadd(A1, B2, C6);
if(nr==4) { CJMADD(A0,B2,C2,T0); }
if(nr==4) { CJMADD(A1,B2,C6,T1); }
if(nr==4) B2 = ei_pload(&blB[14*PacketSize]);
if(nr==4) { MADD(A0,B3,C3,T0); } // C3 = cj.pmadd(A0, B3, C3);
if(nr==4) { CJMADD(A0,B3,C3,T0); }
A0 = ei_pload(&blA[6*PacketSize]);
if(nr==4) { MADD(A1,B3,C7,T1); } // C7 = cj.pmadd(A1, B3, C7);
if(nr==4) { CJMADD(A1,B3,C7,T1); }
A1 = ei_pload(&blA[7*PacketSize]);
if(nr==4) B3 = ei_pload(&blB[15*PacketSize]);
MADD(A0,B0,C0,T0); // C0 = cj.pmadd(A0, B0, C0);
MADD(A1,B0,C4,T1); // C4 = cj.pmadd(A1, B0, C4);
MADD(A0,B1,C1,T0); // C1 = cj.pmadd(A0, B1, C1);
MADD(A1,B1,C5,T1); // C5 = cj.pmadd(A1, B1, C5);
if(nr==4) { MADD(A0,B2,C2,T0); }// C2 = cj.pmadd(A0, B2, C2);
if(nr==4) { MADD(A1,B2,C6,T1); }//C6 = cj.pmadd(A1, B2, C6);
if(nr==4) { MADD(A0,B3,C3,T0); }//C3 = cj.pmadd(A0, B3, C3);
if(nr==4) { MADD(A1,B3,C7,T1); }//C7 = cj.pmadd(A1, B3, C7);
#endif
CJMADD(A0,B0,C0,T0);
CJMADD(A1,B0,C4,T1);
CJMADD(A0,B1,C1,T0);
CJMADD(A1,B1,C5,T1);
if(nr==4) { CJMADD(A0,B2,C2,T0); }
if(nr==4) { CJMADD(A1,B2,C6,T1); }
if(nr==4) { CJMADD(A0,B3,C3,T0); }
if(nr==4) { CJMADD(A1,B3,C7,T1); }
}
blB += 4*nr*PacketSize;
blA += 4*mr;
@ -208,22 +192,40 @@ struct ei_gebp_kernel
// process remaining peeled loop
for(int k=peeled_kc; k<depth; k++)
{
PacketType B0, B1, B2, B3, A0, A1;
if(nr==2)
{
PacketType B0, T0, A0, A1;
A0 = ei_pload(&blA[0*PacketSize]);
A1 = ei_pload(&blA[1*PacketSize]);
B0 = ei_pload(&blB[0*PacketSize]);
CJMADD(A0,B0,C0,T0);
CJMADD(A1,B0,C4,T0);
B0 = ei_pload(&blB[1*PacketSize]);
CJMADD(A0,B0,C1,T0);
CJMADD(A1,B0,C5,T0);
}
else
{
PacketType B0, B1, B2, B3, A0, A1, T0, T1;
A0 = ei_pload(&blA[0*PacketSize]);
A1 = ei_pload(&blA[1*PacketSize]);
B0 = ei_pload(&blB[0*PacketSize]);
B1 = ei_pload(&blB[1*PacketSize]);
C0 = cj.pmadd(A0, B0, C0);
CJMADD(A0,B0,C0,T0);
if(nr==4) B2 = ei_pload(&blB[2*PacketSize]);
C4 = cj.pmadd(A1, B0, C4);
CJMADD(A1,B0,C4,T1);
if(nr==4) B3 = ei_pload(&blB[3*PacketSize]);
C1 = cj.pmadd(A0, B1, C1);
C5 = cj.pmadd(A1, B1, C5);
if(nr==4) C2 = cj.pmadd(A0, B2, C2);
if(nr==4) C6 = cj.pmadd(A1, B2, C6);
if(nr==4) C3 = cj.pmadd(A0, B3, C3);
if(nr==4) C7 = cj.pmadd(A1, B3, C7);
B0 = ei_pload(&blB[(nr==4 ? 4 : 2)*PacketSize]);
CJMADD(A0,B1,C1,T0);
CJMADD(A1,B1,C5,T1);
if(nr==4) { CJMADD(A0,B2,C2,T0); }
if(nr==4) { CJMADD(A1,B2,C6,T1); }
if(nr==4) { CJMADD(A0,B3,C3,T0); }
if(nr==4) { CJMADD(A1,B3,C7,T1); }
}
blB += nr*PacketSize;
blA += mr;
@ -257,63 +259,112 @@ struct ei_gebp_kernel
const Scalar* blB = &blockB[j2*strideB*PacketSize+offsetB*nr];
for(int k=0; k<peeled_kc; k+=4)
{
if(nr==2)
{
PacketType B0, T0, A0;
A0 = ei_pload(&blA[0*PacketSize]);
B0 = ei_pload(&blB[0*PacketSize]);
CJMADD(A0,B0,C0,T0);
B0 = ei_pload(&blB[1*PacketSize]);
CJMADD(A0,B0,C1,T0);
A0 = ei_pload(&blA[1*PacketSize]);
B0 = ei_pload(&blB[2*PacketSize]);
CJMADD(A0,B0,C0,T0);
B0 = ei_pload(&blB[3*PacketSize]);
CJMADD(A0,B0,C1,T0);
A0 = ei_pload(&blA[2*PacketSize]);
B0 = ei_pload(&blB[4*PacketSize]);
CJMADD(A0,B0,C0,T0);
B0 = ei_pload(&blB[5*PacketSize]);
CJMADD(A0,B0,C1,T0);
A0 = ei_pload(&blA[3*PacketSize]);
B0 = ei_pload(&blB[6*PacketSize]);
CJMADD(A0,B0,C0,T0);
B0 = ei_pload(&blB[7*PacketSize]);
CJMADD(A0,B0,C1,T0);
}
else
{
PacketType B0, B1, B2, B3, A0;
PacketType T0, T1;
A0 = ei_pload(&blA[0*PacketSize]);
B0 = ei_pload(&blB[0*PacketSize]);
B1 = ei_pload(&blB[1*PacketSize]);
C0 = cj.pmadd(A0, B0, C0);
CJMADD(A0,B0,C0,T0);
if(nr==4) B2 = ei_pload(&blB[2*PacketSize]);
if(nr==4) B3 = ei_pload(&blB[3*PacketSize]);
B0 = ei_pload(&blB[(nr==4 ? 4 : 2)*PacketSize]);
C1 = cj.pmadd(A0, B1, C1);
CJMADD(A0,B1,C1,T1);
B1 = ei_pload(&blB[(nr==4 ? 5 : 3)*PacketSize]);
if(nr==4) C2 = cj.pmadd(A0, B2, C2);
if(nr==4) { CJMADD(A0,B2,C2,T0); }
if(nr==4) B2 = ei_pload(&blB[6*PacketSize]);
if(nr==4) C3 = cj.pmadd(A0, B3, C3);
if(nr==4) { CJMADD(A0,B3,C3,T1); }
A0 = ei_pload(&blA[1*PacketSize]);
if(nr==4) B3 = ei_pload(&blB[7*PacketSize]);
C0 = cj.pmadd(A0, B0, C0);
CJMADD(A0,B0,C0,T0);
B0 = ei_pload(&blB[(nr==4 ? 8 : 4)*PacketSize]);
C1 = cj.pmadd(A0, B1, C1);
CJMADD(A0,B1,C1,T1);
B1 = ei_pload(&blB[(nr==4 ? 9 : 5)*PacketSize]);
if(nr==4) C2 = cj.pmadd(A0, B2, C2);
if(nr==4) { CJMADD(A0,B2,C2,T0); }
if(nr==4) B2 = ei_pload(&blB[10*PacketSize]);
if(nr==4) C3 = cj.pmadd(A0, B3, C3);
if(nr==4) { CJMADD(A0,B3,C3,T1); }
A0 = ei_pload(&blA[2*PacketSize]);
if(nr==4) B3 = ei_pload(&blB[11*PacketSize]);
C0 = cj.pmadd(A0, B0, C0);
CJMADD(A0,B0,C0,T0);
B0 = ei_pload(&blB[(nr==4 ? 12 : 6)*PacketSize]);
C1 = cj.pmadd(A0, B1, C1);
CJMADD(A0,B1,C1,T1);
B1 = ei_pload(&blB[(nr==4 ? 13 : 7)*PacketSize]);
if(nr==4) C2 = cj.pmadd(A0, B2, C2);
if(nr==4) { CJMADD(A0,B2,C2,T0); }
if(nr==4) B2 = ei_pload(&blB[14*PacketSize]);
if(nr==4) C3 = cj.pmadd(A0, B3, C3);
if(nr==4) { CJMADD(A0,B3,C3,T1); }
A0 = ei_pload(&blA[3*PacketSize]);
if(nr==4) B3 = ei_pload(&blB[15*PacketSize]);
C0 = cj.pmadd(A0, B0, C0);
C1 = cj.pmadd(A0, B1, C1);
if(nr==4) C2 = cj.pmadd(A0, B2, C2);
if(nr==4) C3 = cj.pmadd(A0, B3, C3);
CJMADD(A0,B0,C0,T0);
CJMADD(A0,B1,C1,T1);
if(nr==4) { CJMADD(A0,B2,C2,T0); }
if(nr==4) { CJMADD(A0,B3,C3,T1); }
}
blB += 4*nr*PacketSize;
blA += 4*PacketSize;
}
// process remaining peeled loop
for(int k=peeled_kc; k<depth; k++)
{
if(nr==2)
{
PacketType B0, T0, A0;
A0 = ei_pload(&blA[0*PacketSize]);
B0 = ei_pload(&blB[0*PacketSize]);
CJMADD(A0,B0,C0,T0);
B0 = ei_pload(&blB[1*PacketSize]);
CJMADD(A0,B0,C1,T0);
}
else
{
PacketType B0, B1, B2, B3, A0;
PacketType T0, T1;
A0 = ei_pload(&blA[0*PacketSize]);
B0 = ei_pload(&blB[0*PacketSize]);
B1 = ei_pload(&blB[1*PacketSize]);
C0 = cj.pmadd(A0, B0, C0);
if(nr==4) B2 = ei_pload(&blB[2*PacketSize]);
if(nr==4) B3 = ei_pload(&blB[3*PacketSize]);
C1 = cj.pmadd(A0, B1, C1);
if(nr==4) C2 = cj.pmadd(A0, B2, C2);
if(nr==4) C3 = cj.pmadd(A0, B3, C3);
CJMADD(A0,B0,C0,T0);
CJMADD(A0,B1,C1,T1);
if(nr==4) { CJMADD(A0,B2,C2,T0); }
if(nr==4) { CJMADD(A0,B3,C3,T1); }
}
blB += nr*PacketSize;
blA += PacketSize;
@ -335,18 +386,33 @@ struct ei_gebp_kernel
Scalar C0(0), C1(0), C2(0), C3(0);
const Scalar* blB = &blockB[j2*strideB*PacketSize+offsetB*nr];
for(int k=0; k<depth; k++)
{
if(nr==2)
{
Scalar B0, T0, A0;
A0 = blA[0*PacketSize];
B0 = blB[0*PacketSize];
CJMADD(A0,B0,C0,T0);
B0 = blB[1*PacketSize];
CJMADD(A0,B0,C1,T0);
}
else
{
Scalar B0, B1, B2, B3, A0;
Scalar T0, T1;
A0 = blA[k];
B0 = blB[0*PacketSize];
B1 = blB[1*PacketSize];
C0 = cj.pmadd(A0, B0, C0);
if(nr==4) B2 = blB[2*PacketSize];
if(nr==4) B3 = blB[3*PacketSize];
C1 = cj.pmadd(A0, B1, C1);
if(nr==4) C2 = cj.pmadd(A0, B2, C2);
if(nr==4) C3 = cj.pmadd(A0, B3, C3);
CJMADD(A0,B0,C0,T0);
CJMADD(A0,B1,C1,T1);
if(nr==4) { CJMADD(A0,B2,C2,T0); }
if(nr==4) { CJMADD(A0,B3,C3,T1); }
}
blB += nr*PacketSize;
}
@ -378,13 +444,13 @@ struct ei_gebp_kernel
const Scalar* blB = &blockB[j2*strideB*PacketSize+offsetB];
for(int k=0; k<depth; k++)
{
PacketType B0, A0, A1;
PacketType B0, A0, A1, T0, T1;
A0 = ei_pload(&blA[0*PacketSize]);
A1 = ei_pload(&blA[1*PacketSize]);
B0 = ei_pload(&blB[0*PacketSize]);
C0 = cj.pmadd(A0, B0, C0);
C4 = cj.pmadd(A1, B0, C4);
CJMADD(A0,B0,C0,T0);
CJMADD(A1,B0,C4,T1);
blB += PacketSize;
blA += mr;
@ -431,6 +497,8 @@ struct ei_gebp_kernel
}
};
#undef CJMADD
// pack a block of the lhs
// The travesal is as follow (mr==4):
// 0 4 8 12 ...