rottest.c

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

#include "atlas_misc.h"
#include <assert.h>
#include "atlas_tst.h"

int FAx=0, MAx=0, FAy=0, MAy=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);
      }
   }
}
#ifndef TEST_ROT
void ATL_ROT(const int, TYPE *, const int, TYPE *, const int,
             const TYPE, const TYPE);
   #define TEST_ROT(N_, alpha_, X_, ix_, beta_, Y_, iy_) \
      ATL_ROT(N_, X_, ix_, Y_, iy_, alpha_, beta_)
#endif

#ifdef TREAL

static void good_rot0(const int N, TYPE *X, const int incX,
                      TYPE *Y, const int incY, const TYPE c, const TYPE s)

{
   int i;
   TYPE tmp;

   for (i=N; i; i--, Y += incY, X += incX)
   {
      tmp = c * *X + s * *Y;
      *Y = c * *Y - s * *X;
      *X = tmp;
   }
}
void good_rot(const int N, const TYPE alpha, TYPE *X, const int incX,
              const TYPE beta, TYPE *Y, const int incY)
{
   good_rot0(N, X, incX, Y, incY, alpha, beta);
}

int CheckPad(int npad, TYPE padval, int N, TYPE *Y, int incY)
{
   int i, iret=0;

   incY = Mabs(incY);
   for (i=0; i < npad; i++)
   {
      if (Y[i] != padval)
      {
         iret = i;
         fprintf(stderr, "OVERWRITE %f IN PREPAD %d before beginning of Y!!\n",
                 Y[i], npad-i);
      }
   }
   Y += npad;
   if (incY != 1)
   {
      for (i=0; i < N*incY; i++)
      {
         if (i%incY)
         {
            if (Y[i] != padval)
            {
               iret = i;
               fprintf(stderr, "INTERNAL OVERWRITE %f AT POSITION %d!!\n",
                       Y[i], i);
            }
         }
      }
   }
   Y += 1 + (N-1)*incY;
   for (i=0; i < npad; i++)
   {
      if (Y[i] != padval)
      {
         iret = i;
         fprintf(stderr, "OVERWRITE %f IN POSTPAD %d past end of Y!!\n",
                 Y[i], i+1);
      }
   }
   return(iret);
}
#else
static void good_rot0
   (const int N, TYPE *X, const int incx, TYPE *Y, const int incy,
    const TYPE c0, const TYPE s0)
{
   const register TYPE c = c0, s = s0;
   register TYPE rx, ix, ry, iy;
   const int incX = incx<<1, incY = incy<<1;
   int i;

   if (N > 0)
   {
      for (i=N; i; i--, X += incX, Y += incY)  /* maybe compiler unrolls */
      {
         rx = *X;  ix = X[1];
         ry = *Y;  iy = Y[1];
         *X   = c * rx + s * ry;
         X[1] = c * ix + s * iy;
         *Y   = c * ry - s * rx;
         Y[1] = c * iy - s * ix;
      }
   }
}

void good_rot(const int N, const TYPE alpha, TYPE *X, const int incX,
              const TYPE beta, TYPE *Y, const int incY)
{
   good_rot0(N, X, incX, Y, incY, alpha, beta);
}
int CheckPad(int npad, TYPE padval, int N, TYPE *Y, int incY)
{
   int i, n, iret=0;

   incY = Mabs(incY);
   npad *= 2;
   for (i=0; i < npad; i++)
   {
      if (Y[i] != padval)
      {
         iret = i;
         fprintf(stderr, "OVERWRITE %f IN PREPAD %d before beginning of Y!!\n",
                 Y[i], npad-i);
      }
   }
   Y += npad;
   if (incY != 1)
   {
      for (i=0; i < N*incY; i++)
      {
         if (i%incY)
         {
            if (Y[2*i] != padval)
            {
               iret = i;
               fprintf(stderr, "INTERNAL REAL OVERWRITE %f AT POSITION %d!!\n",
                       Y[2*i], i);
            }
            if (Y[2*i+1] != padval)
            {
               iret = i+1;
               fprintf(stderr, "INTERNAL IMAG OVERWRITE %f AT POSITION %d!!\n",
                       Y[2*i+1], i);
            }
         }
      }
   }
   Y += 2 + 2*(N-1)*incY;
   for (i=0; i < npad; i++)
   {
      if (Y[i] != padval)
      {
         iret = i;
         fprintf(stderr, "OVERWRITE %f IN POSTPAD %d past end of Y!!\n",
                 Y[i], i+1);
      }
   }
   return(iret);
}
#endif

static TYPE my_infnrm(const int N, const TYPE *X, const int incX)
{
   int i, incx = incX SHIFT;
   TYPE max=ATL_rzero, tmp;
   for (i=N; i; i--, X += incx)
   {
      #ifdef TREAL
         tmp = Mabs(*X);
      #else
         tmp = Mabs(*X) + Mabs(X[1]);
      #endif
      if (tmp > max) max = tmp;
   }
   return(max);
}

int checkXY(int N, TYPE *Xg, int incXg, TYPE *Xc, int incXc,
            TYPE *Yg, int incYg, TYPE *Yc, int incYc,
            TYPE normX, TYPE normY)
{
   TYPE *diff;
   TYPE normDX, normDY, resid;
   TYPE eps, THRESH=50.0;
   TYPE Mjoin(PATL, epsilon)(void);
   int iret=0;

   eps = Mjoin(PATL,epsilon)();
   diff = malloc(sizeof(TYPE)*(N SHIFT));
   assert(diff);

   Mjoin(PATL,vdiff)(N, Xc, incXc, Xg, incXg, diff, 1);
   normDX = my_infnrm(N, diff, 1);
   Mjoin(PATL,vdiff)(N, Yc, incYc, Yg, incYg, diff, 1);
   normDY = my_infnrm(N, diff, 1);

   resid = Mmax(normDX, normDY) /
           (eps * Mmax(normX, ATL_rone) * Mmax(normY, ATL_rone) * N);
   if (resid > THRESH || resid != resid)
   {
      fprintf(stderr,
              "ROT resid=%e (normX=%e, normY=%e, normDX=%e, normDY=%e)!!!\n",
              resid, normX, normY, normDX, normDY);
      iret = -1;
   }
   free(diff);
   return(iret);
}
static int CheckXY(int npad, TYPE padvalX, TYPE padvalY,
                   int N, TYPE *Xg, TYPE *Yg, TYPE normX, TYPE normY,
                   TYPE *Xt, int incXt, TYPE *Yt, int incYt)
/*
 * Xg & Yg are known to be inc=1, with no padding
 */
{
   TYPE *xt = Xt + (npad SHIFT), *yt = Yt + (npad SHIFT);
   int i0, i1, i2;
   if (incXt < 0) xt -= ((N-1)SHIFT)*incXt;
   if (incYt < 0) yt -= ((N-1)SHIFT)*incYt;
   i0 = checkXY(N, Xg, 1, xt, incXt, Yg, 1, yt, incYt, normX, normY);
   i1 = CheckPad(npad, padvalX, N, Xt, incXt);
   i2 = CheckPad(npad, padvalY, N, Yt, incYt);
   if (i0 || i1 || i2) return(1);
   return(0);
}

static void vecset(int N, TYPE alpha, TYPE *X)
{