#include <math.h>
#include "bench8sc.h"
#define MASK 0x7fffffff           /* A zero followed by 31 ones*/
#define MASK63 0x7fffffffffffffff /* A zero followed by 63 ones */
#define max(A,B) ((A) > (B) ? (A) : (B))

static char cvs_info[] = "BMARKGRP $Date: 2005/01/10 21:04:37 $ $Revision: 1.2 $ $RCSfile: s8sc.c,v $ $Name: rel_5 $"; 


int s8s(int n,
    int k,
    int l,
    int t,
    long *d,
    double *a,
    int32 ia[NONZERO][MATRIXSIZE],
    int penum,
    int npes)
{

void irand(int,int[],int[]);
void brand64(int,int[],long int[]);
void get_my_idata(int, int, int, int, int, int, int*, int*, int*);
void get_my_bdata64(int, int, int, int, int, int, long int[], int*, int*);
void initiran(int,int,int[]);

  unsigned int i,j,m,kk,ii;
  double tmp;
  long itemp, jj, mm;
  int state[632];
  int  gseed = 9911;
  int zero = 0, one = 1, two = 2, ten = 10;
  int datafirst, datalast;
  int myshare, nelems, chunksize, nseqs;
  int index = 0;

  int rdsize = MATRIXSIZE * NONZERO * NMATRICES;
  long *la;
  int maxlen;
  int halft;
  long *randomdata;
  int32 *randint;
  long temp;

  maxlen = max(MATRIXSIZE * NONZERO * NMATRICES/npes, MATRIXSIZE*2);
  randomdata = (long *) malloc(maxlen*8);
  if (randomdata == (long *)0) {
    printf("Can't malloc space for randomdata\n");
    return(-1);
    }

  randint = (int32 *)randomdata;
   /*
    * Generate the D array
    * The D array needs to be replicated on each site so we will generate
    * it in parallel identically at each PE.
    */
    initiran(0, gseed, state); /* initialize the rand#gen identically */

#if defined(__alpha)
   irand(2*t,state,(int *)d);
#elif defined(_CRAYMPP)
   brand64(t,state,(long *)d);
#endif
   /*
     Have to worry about endian-ness.  T3E and SPARC
     are big-endian.  DEC alpha is little-endian.
   */
   for (i=0; i<t; i=i+2)
     {
      d[i] = 1 & d[i];
      d[i+1] = 1 & d[i+1];
#if 0  /* only if d was a short */
#ifdef __alpha
      temp = d[i+1];
      d[i+1] = d[i];
      d[i] = temp;
#endif
#endif
     }

  myshare = n/npes;
  /*printf("npes = %d, rdsize = %d, myshare = %d\n", npes, rdsize, myshare);*/

  /*
   * Generate the A arrays with entries between 0 and 1
   */

    chunksize = l * k;
    nelems = n * chunksize;

    datafirst = MATRIXSIZE * NONZERO * NMATRICES/npes;

    /*check_maxlen( npes, 1, rdsize, chunksize, &maxlen); */

#ifdef FRANDOM
    get_my_bdata64(&penum, &npes, &ten, &rdsize, &chunksize, &gseed, randomdata, &datafirst, &datalast);
#else
    get_my_bdata64(penum, npes, ten, rdsize, chunksize, gseed, randomdata, &datafirst, &datalast);
#endif

   for (kk=0;kk<myshare;kk++)
     for (i=0;i<k;i++)
       for (j=0; j<l; j++)
         {
           la = (long *)&(randomdata[index++]);
           A(i,j,kk) = (MASK63 & *la) * pow(0.5, 63.0);
         }
         
       
   for (i=0; i<k; i++)
     {
       for(kk = 0; kk<myshare; kk++)
         {
           tmp = A(i,0,kk);
           for (j=1; j<l; j++)
             {
               tmp = tmp + A(i,j,kk);
             }
           for (j=0; j<l; j++)
             {
               A(i,j,kk)  = A(i,j,kk)/tmp;
             }
         }
     }

   /*
    *  Generate the array of indices of the non-zero entries in A
    */
   rdsize = 5*l*n;
   for (j=0; j<l; j++)
     for (m = 0; m<n; m++) ia[j][m] = m;
   
   /*
    *  Randomly permute the array IA
    */
#ifdef FRANDOM
    initiran(&zero, &gseed, state); /* initialize the rand#gen identically */
#else
    initiran(0, gseed, state); /* initialize the rand#gen identically */
#endif
L100:
   /*printf("At 100\n"); */
   nelems = n * 2;
   for (i=0; i<5; i++)
     for (j=0; j<l; j++) {

#ifdef __alpha
#ifdef FRANDOM
       irand(&nelems,state,(int *)randomdata);
#else
       irand(nelems,state,(int *)randomdata);
#endif
#endif
#ifdef _CRAYMPP
       halft = nelems/2;
       brand64(halft,state, (long *)randomdata);
       /*brand64(nelems/2,state, (int *)randomdata);*/
#endif
       for (m = 0; m<n; m++)
         {
#ifdef __alpha
           /* little endian */
           jj = MASK & randint[2*m];
           mm = MASK & randint[2*m+1];
#else
           /* big endian */
           mm = MASK & randint[2*m];
           jj = MASK & randint[2*m+1];
#endif
           jj = (jj % l);
           mm = (mm % n);
           /*printf("mod jj = %u, mm = %u\n",jj,mm);  */
           itemp = ia[j][m];
           ia[j][m] = ia[jj][mm];
           ia[jj][mm] = itemp;
         }
       }
/*
  for (i=0;i<4;i++)
    for (ii = 0;ii<1000;ii++)
      printf("%d\n",ia[i][ii]); 
*/

   /*
    * Make final adjustments to prevent repeats in the index in a column
    */
   for (m=0; m<n; m++)
     {
L140:
     /*printf("at 140, m = %d\n",m); */
     for (j=1; j<l; j++)
       for (jj = 0; jj<j; jj++)
         {
           if (ia[j][m] == ia[jj][m])
             {
               if (m != n)
                 {
                   itemp = ia[j][m];
                   ia[j][m] = ia[j][m+1];
                   ia[j][m+1] = itemp;
                   /*printf("j = %d, m = %d start over\n",j,m); */
                   goto L140;
                  }
                else
                  { /*printf("start way over\n");*/ goto L100; }
             }
         }
     }
free(randomdata);
return 0;
}