atl_gemvn_4x4_1.c

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

      z4 += x2 * A2[12];
      z5 += x3 * A3[13];
      z6 += x0 * A0[14];
      z7 += x1 * A1[15];

      z0 += x3 * A3[8];
      z1 += x0 * A0[9];
      z2 += x1 * A1[10];
      z3 += x2 * A2[11];
      z4 += x3 * A3[12];
      z5 += x0 * A0[13];
      z6 += x1 * A1[14];
      z7 += x2 * A2[15];
      y[8] = z0;
      y[9] = z1;
      y[10] = z2;
      y[11] = z3;
      y[12] = z4;
      y[13] = z5;
      y[14] = z6;
      y[15] = z7;
      if (M-M16) gemvMlt8(M-M16, N, A0+16, lda, x, beta, y+16);
   }
   else if (N) gemvMlt8(M, N, A, lda, x, beta, y);
}

static void gemv32x4(const int M, const int N, const TYPE *A, const int lda,
                     const TYPE *X, const SCALAR beta, TYPE *Y)
{
   #ifdef BETA1
      int j;
   #endif
   const int incA = lda<<2;

   if (N >= 4)
   {
      if (M >= 32)
      {
         #ifdef BETA1
            for (j=(N>>2); j; j--, A += incA, X += 4)
               gemvN32x4(M, 4, A, lda, X, ATL_rone, Y);
            if ( (j = N-((N>>2)<<2)) ) gemvNle4(M, j, A, lda, X, ATL_rone, Y);
         #else
            gemvN32x4(M, 4, A, lda, X, beta, Y);
            if (N != 4)
               Mjoin(PATL,gemvN_a1_x1_b1_y1)
                  (M, N-4, ATL_rone, A+incA, lda, X+4, 1, ATL_rone, Y, 1);
         #endif
      }
      else gemvMlt8(M, N, A, lda, X, beta, Y);
   }
   else if (M) gemvNle4(M, N, A, lda, X, beta, Y);
}

static void gemvMlt8(const int M, const int N, const TYPE *A, const int lda,
                     const TYPE *X, const SCALAR beta, TYPE *Y)
{
   int i;
   register TYPE y0;
   for (i=M; i; i--)
   {
      #ifdef BETA0
         y0 = Mjoin(PATL,dot)(N, A, lda, X, 1);
      #else
         Yget(y0, *Y, beta);
         y0 += Mjoin(PATL,dot)(N, A, lda, X, 1);
      #endif
      *Y++ = y0;
      A++;
   }
}
static void gemvNle4(const int M, const int N, const TYPE *A, const int lda,
                     const TYPE *X, const SCALAR beta, TYPE *Y)
{
   int i;
   const TYPE *A0 = A, *A1 = A+lda, *A2 = A1+lda, *A3 = A2+lda;
   register TYPE x0, x1, x2, x3;
   #ifdef BETAX
      const register TYPE bet=beta;
   #endif

   switch(N)
   {
   case 1:
      #if defined(BETA0)
         Mjoin(PATL,cpsc)(M, *X, A, 1, Y, 1);
      #elif defined(BETAX)
         Mjoin(PATL,axpby)(M, *X, A, 1, beta, Y, 1);
      #else
         Mjoin(PATL,axpy)(M, *X, A, 1, Y, 1);
      #endif
      break;
   case 2:
      x0 = *X; x1 = X[1];
      for (i=0; i != M; i++)
      #ifdef BETA0
         Y[i] = A0[i] * x0 + A1[i] * x1;
      #elif defined(BETAX)
         Y[i] = Y[i]*bet + A0[i] * x0 + A1[i] * x1;
      #else
         Y[i] += A0[i] * x0 + A1[i] * x1;
      #endif
      break;
   case 3:
      x0 = *X; x1 = X[1]; x2 = X[2];
      for (i=0; i != M; i++)
      #ifdef BETA0
         Y[i] = A0[i] * x0 + A1[i] * x1 + A2[i] * x2;
      #elif defined(BETAX)
         Y[i] = Y[i]*bet + A0[i] * x0 + A1[i] * x1 + A2[i] * x2;
      #else
         Y[i] += A0[i] * x0 + A1[i] * x1 + A2[i] * x2;
      #endif
      break;
   case 4:
      if (M >= 32) gemv32x4(M, 4, A, lda, X, beta, Y);
      else
      {
         x0 = *X; x1 = X[1]; x2 = X[2]; x3 = X[3];
         for (i=0; i != M; i++)
         #ifdef BETA0
            Y[i] = A0[i] * x0 + A1[i] * x1 + A2[i] * x2 + A3[i] * x3;
         #elif defined(BETAX)
            Y[i] = Y[i]*bet + A0[i] * x0 + A1[i] * x1 + A2[i] * x2 + A3[i] * x3;
         #else
            Y[i] += A0[i] * x0 + A1[i] * x1 + A2[i] * x2 + A3[i] * x3;
         #endif
      }
      break;
   default:
      ATL_assert(!N);
   }
}

static void gemv4x4(const int M, const int N, const TYPE *A, const int lda,
                    const TYPE *X, const SCALAR beta, TYPE *Y)
{
   const int M4 = (M>>2)<<2, N4 = (N>>2)<<2, nr = N-N4;
   const int incA = lda << 2, incAm = 4 - N4*lda;
   const TYPE *stX = X + N4, *x;
   const TYPE *A0 = A, *A1 = A + lda, *A2 = A1 + lda, *A3 = A2 + lda;
   TYPE *stY = Y + M4;
   register TYPE a00, a10, a20, a30, a01, a11, a21, a31;
   register TYPE a02, a12, a22, a32, a03, a13, a23, a33;
   register TYPE y0, y1, y2, y3, yy0, yy1, yy2, yy3;
   register TYPE x0, x1, x2, x3;

   if (N > 4)
   {
      if (M4)
      {
         do
         {
            a00 = *A0; a01 = *A1;  a02 = *A2; a03 = *A3;
            a10 = A0[1]; a11 = A1[1]; a12 = A2[1]; a13 = A3[1];
            a20 = A0[2]; a21 = A1[2]; a22 = A2[2]; a23 = A3[2];
            a30 = A0[3]; a31 = A1[3]; a32 = A2[3]; a33 = A3[3];
            A0 += incA; A1 += incA; A2 += incA; A3 += incA;
            #ifdef BETA0
               yy0 = yy1 = yy2 = yy3 = y0 = y1 = y2 = y3 = ATL_rzero;
            #elif defined BETAX
               y0 = beta;
               yy0 = *Y; yy1 = Y[1]; yy2 = Y[2];  yy3 = Y[3];
               yy0 *= y0; yy1 *= y0; yy2 *= y0; yy3 *= y0;
               y0 = y1 = y2 = y3 = ATL_rzero;
            #else
               yy0 = *Y; yy1 = Y[1]; yy2 = Y[2];  yy3 = Y[3];
               y0 = y1 = y2 = y3 = ATL_rzero;
            #endif
            x0 = *X; x1 = X[1]; x2 = X[2]; x3 = X[3];
            x = X + 4;
            if (N4 != 4)
            {
               do
               {
                  y0  += x0 * a00; a00 = *A0;
                  yy0 += x1 * a01; a01 = *A1;
                  y1  += x0 * a10; a10 = A0[1];
                  yy1 += x1 * a11; a11 = A1[1];
                  y2  += x0 * a20; a20 = A0[2];
                  yy2 += x1 * a21; a21 = A1[2];
                  y3  += x0 * a30; a30 = A0[3]; A0 += incA; x0 = *x;
                  yy3 += x1 * a31; a31 = A1[3]; A1 += incA; x1 = x[1];

                  y0  += x2 * a02; a02 = *A2;
                  yy0 += x3 * a03; a03 = *A3;
                  y1  += x2 * a12; a12 = A2[1];
                  yy1 += x3 * a13; a13 = A3[1];
                  y2  += x2 * a22; a22 = A2[2];
                  yy2 += x3 * a23; a23 = A3[2];
                  y3  += x2 * a32; a32 = A2[3]; A2 += incA; x2 = x[2];
                  yy3 += x3 * a33; a33 = A3[3]; A3 += incA; x3 = x[3]; x += 4;
               }
               while(x != stX);
            }
/*
 *          Drain pipe
 */
            y0  += x0 * a00;
            yy0 += x1 * a01;
            y1  += x0 * a10;
            yy1 += x1 * a11;
            y2  += x0 * a20;
            yy2 += x1 * a21;
            y3  += x0 * a30;
            yy3 += x1 * a31;

            y0  += x2 * a02;
            yy0 += x3 * a03;
            y1  += x2 * a12;
            yy1 += x3 * a13;
            y2  += x2 * a22;
            yy2 += x3 * a23;
            y3  += x2 * a32;
            yy3 += x3 * a33;
/*
 *          loop cleanup
 */
            switch(nr)
            {
            case 1:
               x0 = *x;
               y0 += x0 * *A0;
               y1 += x0 * A0[1];
               y2 += x0 * A0[2];
               y3 += x0 * A0[3];
               break;
            case 2:
               x0 = *x;
               x1 = x[1];
               y0  += x0 * *A0;
               yy0 += x1 * *A1;
               y1  += x0 * A0[1];
               yy1 += x1 * A1[1];
               y2  += x0 * A0[2];
               yy2 += x1 * A1[2];
               y3  += x0 * A0[3];
               yy3 += x1 * A1[3];
               break;
            case 3:
               x0 = *x; x1 = x[1]; x2 = x[2];
               y0  += x0 * *A0;
               yy0 += x1 * *A1;
               y1  += x0 * A0[1];
               yy1 += x1 * A1[1];
               y2  += x0 * A0[2];
               yy2 += x1 * A1[2];
               y3  += x0 * A0[3];
               yy3 += x1 * A1[3];
               y0 += x2 * *A2;
               y1 += x2 * A2[1];
               y2 += x2 * A2[2];
               y3 += x2 * A2[3];
               break;
            default:;
            }
            A0 += incAm;
            y0 += yy0;
            A1 += incAm;
            y1 += yy1;
            A2 += incAm;
            y2 += yy2;
            A3 += incAm;
            y3 += yy3;
            *Y = y0;
            Y[1] = y1;
            Y[2] = y2;
            Y[3] = y3;
            Y += 4;
         }
         while (Y != stY);
      }
      if (M-M4) gemvMlt8(M-M4, N, A0, lda, X, beta, Y);
   }
   else if (M) gemvNle4(M, N, A, lda, X, beta, Y);
}

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)
{
   gemv4x4(M, N, A, lda, X, beta, Y);
}