atl_gemvn_8x32_2.c

来自「基于Blas CLapck的.用过的人知道是干啥的」· C语言 代码 · 共 448 行 · 第 1/2 页

      A23++;
      A24++;
      A25++;
      A26++;
      A27++;
      A28++;
      A29++;
      A30++;
      A31++;
   }
   while (Y != stY);
}
static void gemv8x16(const int M, const int N, const TYPE *A, const int lda,
                     const TYPE *X, const SCALAR beta, TYPE *Y)
/*
 * Try to let compiler do all the work
 */
{
   const TYPE *A0=A, *A1=A+lda, *A2=A1+lda, *A3=A2+lda;
   const TYPE *A4=A3+lda, *A5=A4+lda, *A6=A5+lda, *A7=A6+lda;
   const TYPE *A8=A7+lda, *A9=A8+lda, *A10=A9+lda, *A11=A10+lda;
   const TYPE *A12=A11+lda, *A13=A12+lda, *A14=A13+lda, *A15=A14+lda;
   const TYPE x0=*X, x1=X[1], x2=X[2], x3=X[3];
   const TYPE x4=X[4], x5=X[5], x6=X[6], x7=X[7];
   const TYPE x8=X[8], x9=X[9], x10=X[10], x11=X[11];
   const TYPE x12=X[12], x13=X[13], x14=X[14], x15=X[15];
   TYPE *stY = Y + ((M>>3)<<3);

   ATL_assert(M >= 8 && (N == 16));
   do
   {
      Yass(*Y) *A0 * x0 + *A1 * x1 + *A2 * x2 + *A3 * x3 +
               *A4 * x4 + *A5 * x5 + *A6 * x6 + *A7 * x7 +
               *A8 * x8 + *A9 * x9 + *A10 * x10 + *A11 * x11 +
               *A12 * x12 + *A13 * x13 + *A14 * x14 + *A15 * x15;
      Yass(Y[1]) A0[1]*x0 + A1[1]*x1 + A2[1]*x2 + A3[1]*x3 +
                 A4[1]*x4 + A5[1]*x5 + A6[1]*x6 + A7[1]*x7 +
                 A8[1]*x8 + A9[1]*x9 + A10[1]*x10 + A11[1]*x11 +
                 A12[1]*x12 + A13[1]*x13 + A14[1]*x14 + A15[1]*x15;
      Yass(Y[2]) A0[2]*x0 + A1[2]*x1 + A2[2]*x2 + A3[2]*x3 +
                 A4[2]*x4 + A5[2]*x5 + A6[2]*x6 + A7[2]*x7 +
                 A8[2]*x8 + A9[2]*x9 + A10[2]*x10 + A11[2]*x11 +
                 A12[2]*x12 + A13[2]*x13 + A14[2]*x14 + A15[2]*x15;
      Yass(Y[3]) A0[3]*x0 + A1[3]*x1 + A2[3]*x2 + A3[3]*x3 +
                 A4[3]*x4 + A5[3]*x5 + A6[3]*x6 + A7[3]*x7 +
                 A8[3]*x8 + A9[3]*x9 + A10[3]*x10 + A11[3]*x11 +
                 A12[3]*x12 + A13[3]*x13 + A14[3]*x14 + A15[3]*x15;
      Yass(Y[4]) A0[4]*x0 + A1[4]*x1 + A2[4]*x2 + A3[4]*x3 +
                 A4[4]*x4 + A5[4]*x5 + A6[4]*x6 + A7[4]*x7 +
                 A8[4]*x8 + A9[4]*x9 + A10[4]*x10 + A11[4]*x11 +
                 A12[4]*x12 + A13[4]*x13 + A14[4]*x14 + A15[4]*x15;
      Yass(Y[5]) A0[5]*x0 + A1[5]*x1 + A2[5]*x2 + A3[5]*x3 +
                 A4[5]*x4 + A5[5]*x5 + A6[5]*x6 + A7[5]*x7 +
                 A8[5]*x8 + A9[5]*x9 + A10[5]*x10 + A11[5]*x11 +
                 A12[5]*x12 + A13[5]*x13 + A14[5]*x14 + A15[5]*x15;
      Yass(Y[6]) A0[6]*x0 + A1[6]*x1 + A2[6]*x2 + A3[6]*x3 +
                 A4[6]*x4 + A5[6]*x5 + A6[6]*x6 + A7[6]*x7 +
                 A8[6]*x8 + A9[6]*x9 + A10[6]*x10 + A11[6]*x11 +
                 A12[6]*x12 + A13[6]*x13 + A14[6]*x14 + A15[6]*x15;
      Yass(Y[7]) A0[7]*x0 + A1[7]*x1 + A2[7]*x2 + A3[7]*x3 +
                 A4[7]*x4 + A5[7]*x5 + A6[7]*x6 + A7[7]*x7 +
                 A8[7]*x8 + A9[7]*x9 + A10[7]*x10 + A11[7]*x11 +
                 A12[7]*x12 + A13[7]*x13 + A14[7]*x14 + A15[7]*x15;
      A0 += 8;
      A1 += 8;
      A2 += 8;
      A3 += 8;
      A4 += 8;
      A5 += 8;
      A6 += 8;
      A7 += 8;
      A8 += 8;
      A9 += 8;
      A10 += 8;
      A11 += 8;
      A12 += 8;
      A13 += 8;
      A14 += 8;
      A15 += 8;
      Y += 8;
   }
   while (Y != stY);
}
static void gemv8x8(const int M, const int N, const TYPE *A, const int lda,
                    const TYPE *X, const SCALAR beta, TYPE *Y)
/*
 * Try to let compiler do all the work
 */
{
   const TYPE *A0=A, *A1=A+lda, *A2=A1+lda, *A3=A2+lda;
   const TYPE *A4=A3+lda, *A5=A4+lda, *A6=A5+lda, *A7=A6+lda;
   const TYPE x0=*X, x1=X[1], x2=X[2], x3=X[3];
   const TYPE x4=X[4], x5=X[5], x6=X[6], x7=X[7];
   TYPE *stY = Y + ((M>>3)<<3);

   ATL_assert(M >= 8 && (N == 8));
   do
   {
      Yass(*Y) *A0 * x0 + *A1 * x1 + *A2 * x2 + *A3 * x3 +
               *A4 * x4 + *A5 * x5 + *A6 * x6 + *A7 * x7;
      Yass(Y[1]) A0[1]*x0 + A1[1]*x1 + A2[1]*x2 + A3[1]*x3 +
                 A4[1]*x4 + A5[1]*x5 + A6[1]*x6 + A7[1]*x7;
      Yass(Y[2]) A0[2]*x0 + A1[2]*x1 + A2[2]*x2 + A3[2]*x3 +
                 A4[2]*x4 + A5[2]*x5 + A6[2]*x6 + A7[2]*x7;
      Yass(Y[3]) A0[3]*x0 + A1[3]*x1 + A2[3]*x2 + A3[3]*x3 +
                 A4[3]*x4 + A5[3]*x5 + A6[3]*x6 + A7[3]*x7;
      Yass(Y[4]) A0[4]*x0 + A1[4]*x1 + A2[4]*x2 + A3[4]*x3 +
                 A4[4]*x4 + A5[4]*x5 + A6[4]*x6 + A7[4]*x7;
      Yass(Y[5]) A0[5]*x0 + A1[5]*x1 + A2[5]*x2 + A3[5]*x3 +
                 A4[5]*x4 + A5[5]*x5 + A6[5]*x6 + A7[5]*x7;
      Yass(Y[6]) A0[6]*x0 + A1[6]*x1 + A2[6]*x2 + A3[6]*x3 +
                 A4[6]*x4 + A5[6]*x5 + A6[6]*x6 + A7[6]*x7;
      Yass(Y[7]) A0[7]*x0 + A1[7]*x1 + A2[7]*x2 + A3[7]*x3 +
                 A4[7]*x4 + A5[7]*x5 + A6[7]*x6 + A7[7]*x7;
      A0 += 8;
      A1 += 8;
      A2 += 8;
      A3 += 8;
      A4 += 8;
      A5 += 8;
      A6 += 8;
      A7 += 8;
      Y += 8;
   }
   while (Y != stY);
}

static void gemv8x4(const int M, const int N, const TYPE *A, const int lda,
                    const TYPE *X, const SCALAR beta, TYPE *Y)
/*
 * Try to let compiler do all the work
 */
{
   const TYPE *A0=A, *A1=A+lda, *A2=A1+lda, *A3=A2+lda;
   const TYPE x0=*X, x1=X[1], x2=X[2], x3=X[3];
   TYPE *stY = Y + ((M>>3)<<3);

   ATL_assert(M >= 8 && (N == 4));
   do
   {
      Yass(*Y) *A0 * x0 + *A1 * x1 + *A2 * x2 + *A3 * x3;
      Yass(Y[1]) A0[1]*x0 + A1[1]*x1 + A2[1]*x2 + A3[1]*x3;
      Yass(Y[2]) A0[2]*x0 + A1[2]*x1 + A2[2]*x2 + A3[2]*x3;
      Yass(Y[3]) A0[3]*x0 + A1[3]*x1 + A2[3]*x2 + A3[3]*x3;
      Yass(Y[4]) A0[4]*x0 + A1[4]*x1 + A2[4]*x2 + A3[4]*x3;
      Yass(Y[5]) A0[5]*x0 + A1[5]*x1 + A2[5]*x2 + A3[5]*x3;
      Yass(Y[6]) A0[6]*x0 + A1[6]*x1 + A2[6]*x2 + A3[6]*x3;
      Yass(Y[7]) A0[7]*x0 + A1[7]*x1 + A2[7]*x2 + A3[7]*x3;
      A0 += 8;
      A1 += 8;
      A2 += 8;
      A3 += 8;
      Y += 8;
   }
   while (Y != stY);
}

static void gemv1x1(const int M, const int N, const TYPE *A, const int lda,
                    const TYPE *X, const SCALAR beta, TYPE *Y)
{
   int i;
   for (i=0; i < M; i++) Y[i] += Mjoin(PATL,dot)(N, X, 1, A+i, lda);
}
#endif

static void gemv(const int M, const int N, const TYPE *A, const int lda,
                 const TYPE *X, const SCALAR beta, TYPE *Y)
{
#ifdef BETA0
   Mjoin(PATL,zero)(M, Y, 1);
   Mjoin(PATL,gemvN_a1_x1_b1_y1)(M, N, A, lda, X, beta, Y);
#elif defined(BETAX)
   Mjoin(PATL,scal)(M, beta, Y, 1);
   Mjoin(PATL,gemvN_a1_x1_b1_y1)(M, N, A, lda, X, beta, Y);
#else
   int i, n;
   const int n32 = N>>5, M8 = (M>>3)<<3, mr = M-M8, incA = lda<<5;

   #ifndef BETA1
      ATL_assert(beta == ATL_rone);
   #endif
   for (i=n32; i; i--)
   {
      if (M8) gemv8x32(M8, 32, A, lda, X, beta, Y);
      if (mr) gemv1x32(mr, 32, A+M8, lda, X, beta, Y+M8);
      A += incA;
      X += 32;
   }
   n = N - (n32<<5);
   if (n >= 16)
   {
      n -= 16;
      if (M8) gemv8x16(M8, 16, A, lda, X, beta, Y);
      if (mr) gemv1x1(mr, 16, A+M8, lda, X, beta, Y+M8);
      A += incA>>1;
      X += 16;
   }
   if (n >= 8)
   {
      n -= 8;
      if (M8) gemv8x8(M8, 8, A, lda, X, beta, Y);
      if (mr) gemv1x1(mr, 8, A+M8, lda, X, beta, Y+M8);
      A += incA>>2;
      X += 8;
   }
   if (n >= 4)
   {
      n -= 4;
      if (M8) gemv8x4(M8, 4, A, lda, X, beta, Y);
      if (mr) gemv1x1(mr, 4, A+M8, lda, X, beta, Y+M8);
      A += incA>>3;
      X += 4;
   }
   if (n) gemv1x1(M, n, A, lda, X, beta, Y);
#endif
}

void Mjoin(Mjoin(Mjoin(Mjoin(Mjoin(PATL,gemvN),NM),_x1),BNM),_y1)
   (const int M, const int N, const SCALAR alpha, const TYPE *A, const int lda,
    const TYPE *X, const int incX, const SCALAR beta, TYPE *Y, const int incY)
{
   gemv(M, N, A, lda, X, beta, Y);
}