axpbysrch.c

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

{
   int i;
   for (i=0; i < n; i++)
      if (ix[i] == 0 && iy[i] == 0 && ia[i] == AlphaX && ib[i] == AlphaX)
         return(bp[i]);
   return(NULL);
}

int KillSpecCases(int n, int *ix, int *iy, int *ia, int *ib, FILENODE **bp)
/*
 * Kills all special cases that are actually handled by general case
 */
{
   FILENODE *gp;
   int i, j;
   gp = FindGen(n, ix, iy, ia, ib, bp);
   if (gp == NULL)
   {
      fprintf(stderr, "NO GENERAL CASE SURVIVED!!  ABORTING!!\n");
      exit(-1);
   }
   for (i=0; i < n; i++)
   {
      if (bp[i] == gp)
      {
         if (ix[i] != 0 || iy[i] != 0 || ia[i] != AlphaX || ib[i] != AlphaX)
            bp[i] = NULL;
      }
   }
/*
 * eliminate entries that are NULL
 */
   for (j=i=0; i < n; i++)
   {
      if (bp[i] != NULL)
      {
         ix[j] = ix[i];
         iy[j] = iy[i];
         ia[j] = ia[i];
         ib[j] = ib[i];
         bp[j] = bp[i];
         j++;
      }
   }
   return(j);
}

void SwapCases(int i, int j, int *ix, int *iy, int *ia, int *ib, FILENODE **bp)
{
   int itmp;
   FILENODE *ftmp;
   if (i != j)
   {
      ftmp = bp[i];
      bp[i] = bp[j];
      bp[j] = ftmp;
      itmp = ix[i];
      ix[i] = ix[j];
      ix[j] = itmp;
      itmp = iy[i];
      iy[i] = iy[j];
      iy[j] = itmp;
      itmp = ia[i];
      ia[i] = ia[j];
      ia[j] = itmp;
      itmp = ib[i];
      ib[i] = ib[j];
      ib[j] = itmp;
   }
}

void SortCases(int n, int *ix, int *iy, int *ia, int *ib, FILENODE **bp)
/*
 * sorts (BFI sort) cases so they are in this order:
 * incX=1, incY=1, alpha=X, beta=X case first
 * other specific incX & incY cases next
 * specific incY, general incX cases next
 * specific incX, general incY cases next
 * general case last
 */
{
/*
 * put specific incX and incY cases first
 */
   int i, j, itmp;
   FILENODE *ftmp;
   for (i=0; i < n; i++)
   {
      if (ix[i] != 0 && iy[i] != 0) continue;
      for (j=i+1; j < n; j++) /* search for other spec cases */
         if (ix[j] != 0 && iy[j] != 0) break;  /* found another case */
      if (j >= n) break; /* no more cases */
      SwapCases(i, j, ix, iy, ia, ib, bp);
   }
/*
 * Put incX=1, incY=1, alpha=X, beta=X case first, if it exists
 */
   for (j=0; j < i; j++)
   {
      if (ix[j] == 1 && iy[j] == 1 && ia[i] == AlphaX && ib[i] == AlphaX)
      {
         SwapCases(0, j, ix, iy, ia, ib, bp);
         break;
      }
   }
/*
 * Put specific incY next
 */
   for (i=0; i < n; i++)  /* find first non-spec case */
      if (ix[i] == 0 || iy[i] == 0) break;
   for (; i < n; i++)
   {
      if (iy[i] != 0) continue;
      for (j=i+1; j < n; j++) /* find other spec incY cases */
         if (iy[j] != 0) break;
      if (j >= n) break; /* no more spec incY cases */
      SwapCases(i, j, ix, iy, ia, ib, bp);
   }
/*
 * Put specific incX next
 */
   for (i=0; i < n; i++)  /* find first non-spec case */
      if (ix[i] == 0 || iy[i] == 0) break;
   for (; i < n; i++) if (iy[i] == 0) break; /* find first gen incY */
   for (; i < n; i++)
   {
      if (ix[i] != 0) continue;
      for (j=i+1; j < n; j++) /* find other spec incX cases */
         if (ix[j] != 0) break;
      if (j >= n) break; /* no more spec incX cases */
      SwapCases(i, j, ix, iy, ia, ib, bp);
   }
/*
 * General case will now have been forced to last by previous sorts
 */
}

FILENODE *KillFN(FILENODE *fbase, FILENODE *fk)
/*
 * Finds fk in Q, and removes it
 */
{
   FILENODE *fp, *pr=NULL;
   if (fk == fbase)
   {
      fp = fbase->next;
      fbase->next = NULL;
      KillAllFN(fbase);
      fbase = fp;
   }
   else
   {
      for (fp=fbase; fp != fk && fp; fp = fp->next) pr = fp;
      assert(fp && pr);
      pr->next = fp->next;
      fp->next = NULL;
      KillAllFN(fp);
   }
   return(fbase);
}

FILENODE *TstFile(char pre, FILENODE *fbase, int *N)
{
   FILENODE *fp, *fn;
   char ln[2048];
   int i, k=0;

   for (fp=fbase; fp; fp = fn)
   {
      fn = fp->next;
      i = sprintf(ln, "make %caxpbytest urout=%s opt=\"", pre, fp->rout);
      if (fp->incX == 0) i += sprintf(ln+i, " -X 4 1 -1 2 -3");
      if (fp->incY == 0) i += sprintf(ln+i, " -Y 4 1 -1 3 -2");
      if (pre == 'c' && pre == 'z') /* complex alphas */
      {
         if (fp->alpha == 0) /* imag = 0 */
            i += sprintf(ln+i, " -a 3 1.0 0.0 -1.0 0.0 0.9 0.0");
         else if (fp->alpha ==  1) i += sprintf(ln+i, " -a 1 1.0 0.0");
         else if (fp->alpha == -1) i += sprintf(ln+i, " -a 1 -1.0 0.0");
      }
      else if (fp->alpha ==  1) i += sprintf(ln+i, " -a 1 1.0");
      else if (fp->alpha == -1) i += sprintf(ln+i, " -a 1 -1.0");
      if (pre == 'c' && pre == 'z') /* complex betas */
      {
         if (fp->beta == 0) /* imag = 0 */
            i += sprintf(ln+i, " -b 3 1.0 0.0 -1.0 0.0 0.8 0.0");
         else if (fp->beta ==  1) i += sprintf(ln+i, " -b 1 1.0 0.0");
         else if (fp->beta == -1) i += sprintf(ln+i, " -b 1 -1.0 0.0");
      }
      else if (fp->beta ==  1) i += sprintf(ln+i, " -b 1 1.0");
      else if (fp->beta == -1) i += sprintf(ln+i, " -b 1 -1.0");
      i += sprintf(ln+i, "\" ");
      if (fp->cc) i += sprintf(ln+i, "%cUCC=\"%s\" ", pre, fp->cc);
      if (fp->ccflags) i += sprintf(ln+i, "%cUCCFLAGS=\"%s\" ",pre,fp->ccflags);
      sprintf(ln+i, "\n");
      fprintf(stderr, "TST: %s", ln);
      #if 1
         if (system(ln)) fbase = KillFN(fbase, fp);
         else k++;
      #else
         k++;
      #endif
   }
   *N = k;
   return(fbase);
}

int NumOfAlpha(int alpha, int n, int *ia)
/*
 * returns number of occurences of alpha in ia
 */
{
   int i, j;
   for (i=j=0; i < n; i++) if (ia[i] == alpha) j++;
   return(j);
}

char *GetNam(char pre, int alpha, int beta, int incX, int incY)
{
   static char nam[128];
   char sx, sy, ca, cb;
   char *rout="axpby";

   if (alpha == 1) ca = '1';
   else if (alpha == -1) ca = 'n';
   else if (alpha == 0 && (pre=='c' || pre=='z')) ca = 'r';
   else ca = 'X';
   if (beta == 1) cb = '1';
   else if (beta == -1) cb = 'n';
   else if (beta == 0 && (pre=='c' || pre=='z')) cb = 'r';
   else cb = 'X';
   if (incX < 0) { sx = 'n'; incX = -incX; }
   else sx = 'p';
   if (incY < 0) { sy = 'n'; incY = -incY; }
   else sy = 'p';
   sprintf(nam, "ATL_%c%s_x%c%dy%c%da%cb%c",
           pre, rout, sx, incX, sy, incY, ca, cb);
   return(nam);
}

int HaveScalCase(int scal, int n, int *is)
{
   int i;
   for (i=0; i < n; i++)
      if (is[i] == scal) return(i+1);
   return(0);
}
int HaveABcase(int alpha, int beta, int n, int *ia, int *ib)
{
   int i;

   for (i=0; i < n; i++)
      if (alpha == ia[i] && beta == ib[i]) return(i+1);
   return(0);
}

void GenIncCases(char pre, char *spc, char *ifp, FILE *fpout, int alpha,
                 int beta, int n, int *ix, int *iy, int *ia, int *ib)
{
   int i;
   for (i=0; i < n; i++)
   {
      if (ia[i] == alpha && ib[i] == beta)
      {
         if (ix[i] != 0 && iy[i] != 0)
            fprintf(fpout, "%s%sif (incx == %d && incy == %d)\n%s{\n",
                    spc, ifp, ix[i], iy[i], spc);
         else if (ix[i] == 0 && iy[i])
            fprintf(fpout, "%s%sif (incy == %d)\n%s{\n", spc, ifp, iy[i], spc);
         else if (iy[i] == 0 && ix[i])
            fprintf(fpout, "%s%sif (incx == %d)\n%s{\n", spc, ifp, ix[i], spc);
         else if (!ifp[0]) /* only general case */
         {
            fprintf(fpout, "%s%s(N, alpha, X, incx, beta, Y, incy);\n", spc,
                    GetNam(pre, alpha, beta, 0, 0));
            fprintf(fpout, "%sreturn;\n", spc);
            return;
         }
         else fprintf(fpout, "%s%selse\n%s{\n", spc, spc);
         fprintf(fpout, "%s   %s(N, alpha, X, incx, beta, Y, incy);\n", spc,
                 GetNam(pre, alpha, beta, ix[i], iy[i]));
         fprintf(fpout, "%s   return;\n%s}\n", spc, spc);
      }
      ifp = "else ";
   }
}

void GenRealIf(char pre, char *spc, FILE *fpout,
               int n, int *ix, int *iy, int *ia, int *ib)
{
   char *ifp="";
   fprintf(fpout, "%sif (beta == ATL_rone)\n%s{\n", spc, spc);
   fprintf(fpout, "%s   ATL_%caxpy(N, alpha, X, incx, Y, incy);\n", spc, pre);
   fprintf(fpout, "%s   return;\n%s}\n", spc, spc);
   fprintf(fpout, "%selse if (beta == ATL_rzero)\n%s{\n", spc, spc);
   fprintf(fpout, "%s   ATL_%ccpsc(N, alpha, X, incx, Y, incy);\n", spc, pre);
   fprintf(fpout, "%s   return;\n%s}\n", spc, spc);
   fprintf(fpout, "%selse if (alpha == ATL_rzero)\n%s{\n", spc, spc);
   fprintf(fpout, "%s   ATL_%cscal(N, beta, Y, incy);\n", spc, pre);
   fprintf(fpout, "%s   return;\n%s}\n", spc, spc);
   if ( HaveScalCase(1, n, ia) )
   {
      fprintf(fpout, "%selse if (alpha == ATL_rone)\n%s{\n", spc, spc);
      spc -= 3;
      if (HaveABcase(1, -1, n, ia, ib))
      {
         fprintf(fpout, "%sif (beta == ATL_rnone)\n%s{\n", spc, spc);
         GenIncCases(pre, spc-3, "", fpout, 1, -1, n, ix, iy, ia, ib);
         fprintf(fpout, "%s}\n", spc);
      }
      if (HaveABcase(1, AlphaX, n, ia, ib))
      {
         GenIncCases(pre, spc, "", fpout, 1, AlphaX, n, ix, iy, ia, ib);
      }
      spc += 3;
      fprintf(fpout, "%s}\n", spc);
   }
   if ( HaveScalCase(-1, n, ia) )
   {
      fprintf(fpout, "%selse if (alpha == ATL_rnone)\n%s{\n", spc, spc);
      spc -= 3;
      if (HaveABcase(-1, -1, n, ia, ib))
      {
         fprintf(fpout, "%sif (beta == ATL_rnone)\n%s{\n", spc, spc);
         GenIncCases(pre, spc-3, "", fpout, -1, -1, n, ix, iy, ia, ib);
         fprintf(fpout, "%s}\n", spc);
      }
      if (HaveABcase(-1, AlphaX, n, ia, ib))
      {
         GenIncCases(pre, spc, "", fpout, -1, AlphaX, n, ix, iy, ia, ib);
      }
      spc += 3;
      fprintf(fpout, "%s}\n", spc);
   }
   if ( HaveABcase(AlphaX, -1, n, ia, ib) )
   {
      fprintf(fpout, "%sif (beta == ATL_rnone)\n%s{\n", spc, spc);
      GenIncCases(pre, spc-3, "", fpout, AlphaX, -1, n, ix, iy, ia, ib);
      fprintf(fpout, "%s}\n", spc);
   }
}
void GenCplxIf(char pre, char *spc, FILE *fpout,
               int n, int *ix, int *iy, int *ia, int *ib)
{ /* HERE HERE HERE */
   fprintf(fpout, "%sif (ibeta == ATL_rzero) /* real beta */\n%s{\n",
           spc, spc);
   fprintf(fpout, "%s   if (rbeta == ATL_rone)\n%s   {\n", spc, spc);
   fprintf(fpout, "%s      ATL_%caxpy(N, alpha, X, incx, Y, incy);\n",
           spc, pre);
   fprintf(fpout, "%s      return;\n%s   }\n", spc, spc);
   fprintf(fpout, "%s   else if (rbeta == ATL_rzero)\n%s   {\n", spc, spc);
   fprintf(fpout, "%s      ATL_%ccpsc(N, alpha, X, incx, Y, incy);\n",
           spc, pre);
   fprintf(fpout, "%s      return;\n%s   }\n", spc, spc);
   fprintf(fpout,
      "%s   else if (ialpha == ATL_rzero) /* both scalars real */\n%s   {\n",
           spc, spc);
   fprintf(fpout, "%s      if (incX == 1 && incY == 1)\n%s      {\n", spc, spc);
   fprintf(fpout,
           "%s         Mjoin(PATLU,axpby)(N+N, ralpha, X, 1, rbeta, Y, 1);\n",
           spc);
   fprintf(fpout, "%s         return;\n%s      }\n", spc, spc);
   if ( HaveABcase(0, 0, n, ia, ib) )
   {
      fprintf(fpout, "%s      else\n%s      {\n", spc, spc);
      GenIncCases(pre, spc-9, "", fpout, 0, 0, n, ix, iy, ia, ib);
      fprintf(fpout, "%s      }\n", spc);
   }

   if ( HaveABcase(AlphaX, 0, n, ia, ib) )
   {
      fprintf(fpout, "%s   }\n%s   else /* cplx alpha, real beta cases\n",
              spc, spc);
      GenIncCases(pre, spc-6, "", fpout, AlphaX, 0, n, ix, iy, ia, ib);
   }
   fprintf(fpout, "%s   }\n", spc);
   fprintf(fpout,
   "%s}\n%selse if (ialpha == ATL_rzero) /* real alpha, cplx beta */\n%s{\n",
           spc, spc, spc);
   fprintf(fpout, "%s   if (ralpha == ATL_rzero)\n%s   {\n", spc, spc);
   fprintf(fpout, "%s      ATL_%cscal(N, beta, Y, incy);\n", spc, pre);
   fprintf(fpout, "%s      return;\n%s   }\n", spc, spc);

   if ( HaveABcase(1, AlphaX, n, ia, ib) )
   {
      fprintf(fpout, "%s   else if (ralpha == ATL_rone)\n%s   {\n", spc, spc);
      GenIncCases(pre, spc-6, "", fpout, 1, AlphaX, n, ix, iy, ia, ib);
      fprintf(fpout, "%s   }\n", spc);
   }
   if ( HaveABcase(-1, AlphaX, n, ia, ib) )
   {
      fprintf(fpout, "%s   else if (ralpha == ATL_rnone)\n%s   {\n", spc, spc);
      GenIncCases(pre, spc-6, "", fpout, -1, AlphaX, n, ix, iy, ia, ib);
      fprintf(fpout, "%s   }\n", spc);
   }
   if ( HaveABcase(0, AlphaX, n, ia, ib) )
   {
      fprintf(fpout, "%s   else\n%s   {\n", spc, spc);
      GenIncCases(pre, spc-6, "", fpout, 0, AlphaX, n, ix, iy, ia, ib);
      fprintf(fpout, "%s   }\n", spc);
   }
   fprintf(fpout, "%s}\n", spc);
}
void GenAlpBetCase(char pre, char *spc, FILE *fpout,
                   int n, int *ix, int *iy, int *ia, int *ib)
{
   if (pre == 'c' || pre == 'z') GenCplxIf(pre, spc, fpout, n, ix, iy, ia, ib);
   else GenRealIf(pre, spc, fpout, n, ix, iy, ia, ib);
}
void MangleIncs(FILE *fpout, char *spc, char pre)
/*
 * makes incs positive if possible, if one must be negative, make sure
 * it is not stride 1; if both are stride 1, or neither,
 * make incY positive over incX
 * NOTE: incX and incY will therefore never both be neg,
 *       incY will be negative only if incX = -1