mmtst.c

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

      Mjoin(NBmm0,_b1)(m_, n_, k_, *(alp_), a_, lda_, (b_)+1, ldb_, \
                       rone, (c_)+1, ldc_); \
   }
#endif
#if 0
   #undef NBmm
   #define NBmm(m_, n_, k_, alp_, a_, lda_, b_, ldb_, bet_, c_, ldc_) \
   { \
      zgemm_("T", "N", &(m_), &(n_), &(k_), alp_, a_, &(lda_), b_, \
             &(ldb_), bet_, c_, &(ldc_)); \
   }
#endif
#endif

void tst_mm(const int M, const int N, const int K, const SCALAR alpha,
            const TYPE *A, const int lda0, const TYPE *B, const int ldb0,
            const SCALAR beta, TYPE *C, const int ldc0)
{
   int i, j, k;
   int lda = lda0 SHIFT, ldb = ldb0 SHIFT, ldc = ldc0 SHIFT;
   register TYPE c0;
   #ifdef TREAL
      for (j=0; j < N; j++)
      {
         for (i=0; i < M; i++)
         {
            c0 = 0.0;
            for (k=0; k < K; k++)
            {
               #if defined(NoTransA) && defined(NoTransB)
                  c0 += A[i+k*lda] * B[j*ldb+k];
               #elif defined(NoTransA) && defined(TransB)
                  c0 += A[i+k*lda] * B[j+k*ldb];
               #elif defined(TransA) && defined(NoTransB)
                  c0 += A[i*lda+k] * B[j*ldb+k];
               #elif defined(TransA) && defined(TransB)
                  c0 += A[i*lda+k] * B[j+k*ldb];
               #endif
            }
            C[i+j*ldc] = beta*C[i+j*ldc] + alpha*c0;
         }
      }
   #else
      register TYPE cr, ci, ar, ai, br, bi;
/*
 *    If matrix is stored split into real & imaginary parts, allocate some
 *    matrices and intermix them for f77-like imaginary matrices
 */
      #if csA == 1 && csB == 1
         int inc, ldaa, ldbb, lda2=lda/2, ldb2=ldb/2;
         TYPE *aa, *bb;
         aa = malloc(M*K*sizeof(TYPE)*2);
         bb = malloc(K*N*sizeof(TYPE)*2);
         assert(aa && bb);
         #ifdef NoTransA
            ldaa = M*2;
            inc = K*lda2;
            for (k=0; k < K; k++)
            {
               j = k*lda;
               for (i=0; i < M; i++)
               {
                  aa[k*ldaa+2*i]   = A[inc+k*lda2+i];
                  aa[k*ldaa+2*i+1] = A[k*lda2+i];
               }
            }
         #else
            ldaa = K*2;
            inc = M*lda2;
            for (i=0; i < M; i++)
            {
               for (k=0; k < K; k++)
               {
                  aa[i*ldaa+2*k]   = A[inc+i*lda2+k];
                  aa[i*ldaa+2*k+1] = A[i*lda2+k];
               }
            }
         #endif
         #ifdef NoTransB
            inc = N*ldb2;
            ldbb = K*2;
            for (j=0; j < N; j++)
            {
               for (k=0; k < K; k++)
               {
                  bb[j*ldbb+2*k]   = B[inc+j*ldb2+k];
                  bb[j*ldbb+2*k+1] = B[j*ldb2+k];
               }
            }
         #else
            inc = K*ldb2;
            ldbb = N*2;
            for (k=0; k < K; k++)
            {
               for (j=0; j < N; j++)
               {
                  bb[k*ldbb+2*j]   = B[inc+k*lda2+j];
                  bb[k*ldbb+2*j+1] = B[k*lda2+j];
               }
            }
         #endif
         A = (const TYPE *) aa;
         B = (const TYPE *) bb;
         lda = ldaa;
         ldb = ldbb;
      #endif

      for (j=0; j < N; j++)
      {
         for (i=0; i < M; i++)
         {
            cr = ci = 0.0;
            for (k=0; k < K; k++)
            {
               #if defined(NoTransA) && defined(NoTransB)
                  ar = A[2*i+k*lda];
                  ai = A[2*i+k*lda+1];
                  br = B[j*ldb+2*k];
                  bi = B[j*ldb+2*k+1];
               #elif defined(NoTransA) && !defined(NoTransB)
                  ar = A[2*i+k*lda] ;
                  ai = A[2*i+k*lda+1];
                  br = B[2*j+k*ldb];
                  bi = B[2*j+k*ldb+1];
               #elif !defined(NoTransA) && defined(NoTransB)
                  ar = A[i*lda+k*2];
                  ai = A[i*lda+k*2+1];
                  br = B[j*ldb+k*2];
                  bi = B[j*ldb+k*2+1];
               #elif !defined(NoTransA) && !defined(NoTransB)
                  ar = A[i*lda+k*2];
                  ai = A[i*lda+k*2+1];
                  br = B[2*j+k*ldb];
                  bi = B[2*j+k*ldb+1];
               #endif
               #ifdef ConjTransA
                  ai = -ai;
               #endif
               #ifdef ConjTransB
                  bi = -bi;
               #endif
               cr += ar * br - ai * bi;
               ci += ar * bi + ai * br;
            }
/*
 *          Scale by alpha
 */
            ar = *alpha;
            ai = alpha[1];
            br = cr;
            bi = ci;
            cr =  br * ar;
            ci =  bi * ar;
            cr -= bi * ai;
            ci += br * ai;
/*
 *          Scale C by beta
 */
            br = *beta;
            bi = beta[1];
            ar = C[2*i+j*ldc];
            ai = C[2*i+j*ldc+1];
            C[2*i+j*ldc]   = ar*br - ai * bi;
            C[2*i+j*ldc+1] = ai*br + ar * bi;
/*
 *          Store answer back to C
 */
            C[2*i+j*ldc]   += cr;
            C[2*i+j*ldc+1] += ci;

         }
      }
      #if csA == 1 && csB == 1
         free(aa);
         free(bb);
      #endif
   #endif
}
int mmtst(void)
{
   char fnam[80];
#if defined(LDA) && LDA != 0
      const int lda=LDA;
#else
      const int lda=2*LDA2;
#endif
#if defined(LDB) && LDB != 0
   const int ldb=LDB;
#else
   const int ldb=2*LDB2;
#endif
#if defined(LDC) && LDC != 0
   const int ldc=LDC;
#else
   const int ldc=2*LDC2;
#endif
   int nA, nB;
   #ifdef TCPLX
      int inca, incb, incc;
      const TYPE one=1.0, none=(-1.0);
      #if (ALPHA == 1)
         TYPE alpha[2] = {1.0, 0.0};
      #elif (ALPHA == -1)
         TYPE alpha[2] = {-1.0, 0.0};
      #else
         TYPE alpha[2] = {2.3, 0.0};
      #endif
      #if (BETA == 1)
         TYPE beta[2] = {1.0, 0.0};
      #elif (BETA == -1)
         TYPE beta[2] = {-1.0, 0.0};
      #elif (BETA == 0)
         TYPE beta[2] = {0.0, 0.0};
      #else
         TYPE beta[2] = {1.3, 0.0};
      #endif
   #else
      #ifdef ALPHA
         TYPE alpha=ALPHA;
      #else
         TYPE alpha=1.0;
      #endif
      #ifdef BETA
         TYPE beta=BETA;
      #else
         TYPE beta=1.0;
      #endif
   #endif
   const TYPE rone=1.0, rnone=(-1.0);
   void *va=NULL, *vb=NULL, *vc=NULL;
   TYPE *C0, *C1, *A, *B;
   TYPE diff, tmp;
   int i, j, k, n, nerr;
   int M=MB, N=NB, K=KB;
   TYPE ErrBound;

   if (!M) M = MB0;
   if (!N) N = NB0;
   if (!K) K = KB0;
   #ifdef TREAL
      ErrBound = 2.0 * (Mabs(alpha) * 2.0*K*EPS + Mabs(beta) * EPS) + EPS;
   #else
      diff = Mabs(*alpha) + Mabs(alpha[1]);
      tmp = Mabs(*beta) + Mabs(beta[1]);
      ErrBound =  2.0 * (diff*8.0*K*EPS + tmp*EPS) + EPS;
   #endif
   #ifdef NoTransA
      nA = K;
   #else
      nA = M;
   #endif
   #ifdef NoTransB
      nB = N;
   #else
      nB = K;
   #endif
   #ifdef TCPLX
      inca = lda*nA;
      incb = ldb*nB;
   #endif
   #ifdef ATL_MinMMAlign
      va = malloc(ATL_MinMMAlign + lda*nA*ATL_sizeof);
      vb = malloc(ATL_MinMMAlign + ldb*nB*ATL_sizeof);
      vc = C0 = malloc(2*ldc*N*ATL_sizeof);
      assert(va && vb && C0);
      A = (TYPE *) ( ( ((size_t) va)/ATL_MinMMAlign ) * ATL_MinMMAlign
                     + ATL_MinMMAlign );
      B = (TYPE *) ( ( ((size_t) vb)/ATL_MinMMAlign ) * ATL_MinMMAlign
                     + ATL_MinMMAlign );
   #else
      C0 = vc = malloc( (2*ldc*N + lda*nA + ldb*nB) * ATL_sizeof);
      assert(vc);
      A = C1 + (ldc * N SHIFT);
      B = A + (lda * nA SHIFT);
   #endif
   C1 = C0 + (ldc * N SHIFT);
   for (n=lda*nA SHIFT, i=0; i < n; i++) A[i] = dumb_rand();
   for (n=ldb*nB SHIFT, i=0; i < n; i++) B[i] = dumb_rand();
   for (n=ldc*N SHIFT, i=0; i < n; i++) C0[i] = C1[i] = dumb_rand();
   tst_mm(M, N, K, alpha, A, lda, B, ldb, beta, C0, ldc);
   NBmm(M, N, K, alpha, A, lda, B, ldb, beta, C1, ldc);
   nerr = 0;
   for (j=0; j < N; j++)
   {
      for (i=0; i < M SHIFT; i++)
      {
         k = i + j*(ldc SHIFT);
         diff = C0[k] - C1[k];
         if (diff < 0.0) diff = -diff;
         if (diff > ErrBound)
         {
            fprintf(stderr, "C(%d,%d) : expected=%f, got=%f\n",
                    i, j, C0[k], C1[k]);
            nerr++;
         }
      }
   }
   free(vc);
   if (va) free(va);
   if (vb) free(vb);
   return(nerr);
}

main()
{
   int ierr;
   ierr = mmtst();
   if (!ierr) fprintf(stdout, "PASSED TEST\n");
   exit(ierr);
}