mvtest.c

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

/*
 *             Automatically Tuned Linear Algebra Software v3.8.0
 *                    (C) Copyright 2000 R. Clint Whaley
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *   1. Redistributions of source code must retain the above copyright
 *      notice, this list of conditions and the following disclaimer.
 *   2. Redistributions in binary form must reproduce the above copyright
 *      notice, this list of conditions, and the following disclaimer in the
 *      documentation and/or other materials provided with the distribution.
 *   3. The name of the ATLAS group or the names of its contributers may
 *      not be used to endorse or promote products derived from this
 *      software without specific written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE ATLAS GROUP OR ITS CONTRIBUTORS
 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 *
 */
#include "atlas_misc.h"
#include "atlas_tst.h"
#include "atlas_lvl2.h"
int FAx=0, MAx=0, FAy=0, MAy=0, FAa=0, MAa=0;
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>

struct FA_allocs
{
   void *mem, *memA;
   struct FA_allocs *next;
} *allocQ=NULL;

struct FA_allocs *NewAlloc(size_t size, struct FA_allocs *next,
                           int align, int misalign)
/*
 * Allocates size allocation that is aligned to [align], but not aligned
 * to [misalign].  Therefore, misalign > align.  Align must minimally be sizeof
 * while misalign may be 0 if we don't need to avoid a particular alignment.
 */
{
   void *vp;
   char *cp;
   struct FA_allocs *ap;
   int n, i;
   const int malign = align >= misalign ? align : misalign;

   n = size + align + malign;
   i = (n >> 3)<<3;
   if (n != i)
      n += n - i;
   cp = malloc(n + sizeof(struct FA_allocs));
   assert(cp);
   ap = (struct FA_allocs *) (cp + n);
   ap->mem = cp;
/*
 * Align to min alignment
 */
   ap->memA = align ? (void*) ((((size_t) cp)/align)*align + align) : cp;
/*
 * Misalign to misalign
 */
   if (misalign)
   {
      if (((size_t)ap->memA)%misalign == 0)
         ap->memA = ((char*)ap->memA) + align;
   }
   ap->next = next;
   return(ap);
}

/*
 * no-align malloc free retaining system default behavior
 */
void *NA_malloc(size_t size)
{
   return(malloc(size));
}
void *NA_calloc(size_t n, size_t size)
{
   return(calloc(n, size));
}
void NA_free(void *ptr)
{
   free(ptr);
}


/*
 * malloc/free pair that aligns data to align, but not to misalign
 */
void *FA_malloc(size_t size, int align, int misalign)
{
   if ((!misalign && align <= 8) || !size)
      return(malloc(size));
   else
   {
      allocQ = NewAlloc(size, allocQ, align, misalign);
      return(allocQ->memA);
   }
}
void *FA_calloc(size_t n, size_t size, int align, int misalign)
{
   char *cp;
   int *ip;
   double *dp;
   size_t i;
   size_t tsize;
   tsize = n * size;
   cp = FA_malloc(tsize, align, misalign);
   if (size == sizeof(int))
      for (ip=(int*)cp,i=0; i < n; i++)
        ip[i] = 0;
   else if (size == sizeof(double))
      for (dp=(double*)cp,i=0; i < n; i++)
        dp[i] = 0.0;
   else
      for (i=0; i < tsize; i++)
        cp[i] = 0;
   return(cp);
}

void FA_free(void *ptr, int align, int misalign)
/*
 * Part of malloc/free pair that aligns data to FALIGN
 */
{
   struct FA_allocs *ap, *prev;
   if (ptr)
   {
      if ((!misalign && align <= 8))
         free(ptr);
      else
      {
         for (ap=allocQ; ap && ap->memA != ptr; ap = ap->next) prev = ap;
         if (!ap)
         {
            fprintf(stderr, "Couldn't find mem=%ld\nmemQ=\n", ptr);
            for (ap=allocQ; ap; ap = ap->next)
               fprintf(stderr, "   %ld, %ld\n", ap->memA, ap->mem);
         }
         assert(ap);
         if (ap == allocQ)
            allocQ = allocQ->next;
         else
            prev->next = ap->next;
         free(ap->mem);
      }
   }
}

int CheckAns(int M, int N, SCALAR beta, TYPE *Y, TYPE *Y0)
{
   TYPE eps, diff;
   #ifdef TCPLX
      const int M2 = M<<1, mul=8;
   #else
      const int M2 = M, mul=4;
   #endif
   int i, ierr=0;

   eps = Mjoin(PATL,epsilon)();
   for (i=0; i < M2; i++)
   {
      diff = Y0[i] - Y[i];
      if (diff < 0.0) diff = -diff;
      if (diff > N*mul*eps + 2.0*eps)
      {
         ierr = i+1;
         fprintf(stderr, "Y[%d] = %f, expected %f\n", i, Y[i], Y0[i]);
      }
   }
   return(ierr);
}

#ifdef TREAL
static void dumb_gemv(enum ATLAS_TRANS TA, int M, int N, SCALAR alpha,
                      TYPE *A, int lda, TYPE *X, SCALAR beta, TYPE *Y)
{
   register TYPE y0;
   int i, j;

   ATL_assert(alpha == ATL_rone);
   if (beta == ATL_rzero) Mjoin(PATL,zero)(M, Y, 1);
   else if (beta != ATL_rone) Mjoin(PATL,scal)(M, beta, Y, 1);

   if (TA == AtlasNoTrans)
   {
      for (j=0; j < N; j++) Mjoin(PATL,axpy)(M, X[j], A+j*lda, 1, Y, 1);
   }
   else
   {
      for (i=0; i < M; i++)
      {
         y0 = Y[i];
         for (j=0; j < N; j++) y0 += A[i*lda+j] * X[j];
         Y[i] = y0;
      }
   }
}
#else
static void dumb_gemv(enum ATLAS_TRANS TA, int M, int N, SCALAR alpha,
                      TYPE *A, int lda, TYPE *X, SCALAR beta, TYPE *Y)
{
   register TYPE ry, iy, rx, ix, ra, ia;
   int i, j;
   const int M2 = M<<1, N2 = N<<1;

   ATL_assert(SCALAR_IS_ONE(alpha));
   if ( SCALAR_IS_ZERO(beta) ) Mjoin(PATL,zero)(M, Y, 1);
   else if ( !SCALAR_IS_ONE(beta) ) Mjoin(PATL,scal)(M, beta, Y, 1);

   if (TA == AtlasNoTrans)
      for (j=0; j < N; j++) Mjoin(PATL,axpy)(M, X+2*j, A+j*2*lda, 1, Y, 1);
   else if (TA == AtlasConj)
      for (j=0; j < N; j++) Mjoin(PATL,axpyConj)(M, X+2*j, A+j*2*lda, 1, Y, 1);
   else if (TA == AtlasTrans)
   {
      for (i=0; i < M2; i += 2)
      {
         ry = Y[i];
         iy = Y[i+1];
         for (j=0; j < N2; j += 2)
         {
            ra = A[i*lda+j]; ia = A[i*lda+j+1];
            rx = X[j];       ix = X[j+1];
            ry += rx * ra - ix * ia;
            iy += rx * ia + ix * ra;
         }
         Y[i] = ry;
         Y[i+1] = iy;
      }
   }
   else if (TA == AtlasConjTrans)
   {
      for (i=0; i < M2; i += 2)
      {
         ry = Y[i];
         iy = Y[i+1];
         for (j=0; j < N2; j += 2)
         {
            ra = A[i*lda+j]; ia = -A[i*lda+j+1];
            rx = X[j];       ix = X[j+1];
            ry += rx * ra - ix * ia;
            iy += rx * ia + ix * ra;
         }
         Y[i] = ry;
         Y[i+1] = iy;
      }
   }
}
#endif

int RunTest(enum ATLAS_TRANS TA, int M, int N, int lda, TYPE *beta0)
{
   #ifdef TCPLX
      TYPE alpha[2] = {ATL_rone, ATL_rzero};
      #define beta beta0
   #else
      TYPE alpha= ATL_rone;
      TYPE beta = *beta0;
   #endif
   TYPE *A, *X, *Y, *Y0;
   int mA, nA, ierr;
   char cta, ln[32];

   if (TA == AtlasNoTrans) cta = 'N';
   else if (TA == AtlasTrans) cta = 'T';
   else if (TA == AtlasConjTrans) cta = 'C';
   else cta = '-';
   if (TA == AtlasNoTrans || TA == AtlasConj) { mA = M; nA = N; }
   else { mA = N; nA = M; }

   #ifdef TREAL
      fprintf(stdout, "   TEST TA=%c, M=%d, N=%d, lda=%d, beta=%f STARTED\n",
              cta, M, N, lda, beta);
   #else
      fprintf(stdout,
              "   TEST TA=%c, M=%d, N=%d, lda=%d, beta=(%f,%f) STARTED\n",