/**************************************************************************
*
*   Subroutine S11 builds the "random" matrices A and B.
*
*   Parameters:
*
*    Provided by calling routine:
*       N        =  The size of the data to be multiplied
*       NUMSIZE  =  The number of words used to store an MP integer
*                   (passed in so S11 can declare arrays with exact dimensions)
*       MODBITS  =  The number of bits in MODULUS
*
*    Returned by this routine:
*       A        =  Matrix of multiple precision numbers
*       B        =  Matrix of multiple precision numbers
*
***************************************************************************/

#include "bench11.h"

/* CVS info   */
/* $Date: 2005/01/10 20:57:47 $     */
/* $Revision: 1.2 $ */
/* $RCSfile: s11.c,v $  */
/* $Name: rel_5 $     */

static char cvs_info[] = "BMARKGRP $Date: 2005/01/10 20:57:47 $ $Revision: 1.2 $ $RCSfile: s11.c,v $ $Name: rel_5 $";


/*
 * Fortran can define A(NUMSIZE,M,M) with variable dimensions.
 * this kludge in C lets us use same notation as if we declared a[m][m][numsize]
 * where m and numsize are arguments of this subr call
 * The 3rd dimension is for the MP integers which are vectors of 64-bit ints
 */
 
#define a(I,J,K)        a[ (K) + numsize * ( (J) + n * (I) ) ]
#define b(I,J,K)        b[ (K) + numsize * ( (J) + n * (I) ) ]

void s11 ( mplong a[], mplong b[], int n, int numsize, int modbits )
{
   int abcnumsize;
   int fullwords, extrabits;
   int skipbits, bitctr, wordctr;
   mplong modsize, asize, bsize, extrabitmask;

   /* Array to hold enough bits for one row (N long) of A and B */
   unsigned mplong * tempbits;

   unsigned mplong temp;

   /* Other variables */

   int irow, myfirst, mylast, ranwds;
   int i, j, k;
   int one = 1, seed = -99911;

   /* this is set in M11 */
   extern int verbose;

   abcnumsize = numsize - MPHL;

   if (verbose)
   {
      printf( "CALL S11(%d, modulus, -, -) \n", n );
      printf( "           (Set up the matrices) \n" );
   }

   modsize = numsize - MPHL;

   /* Set up counters and mask needed to fill A and B */

   fullwords = modbits / BITSPERWORD;
   extrabits = modbits % BITSPERWORD;
   if ( extrabits > 0 )
#if (BITSPERWORD == 64)
      extrabitmask = ((unsigned mplong)(-1L)) >> (BITSPERWORD - extrabits);
#else  /* BPW == 32 */
      extrabitmask = BitMask32 >> (BITSPERWORD - extrabits);
#endif

   if (verbose)
   {
      printf( "          Modulus has %d bits \n", modbits );
      printf( "           Fill the matrices \n" );
   }

/**************************************************************************
*
*  Use the random number generator to fill the matrices A and B.
*
*  Place the random bits into the multiple precision numbers which
*  make up A and B.
*
**************************************************************************/

/* Initialize the structures for the multi-precision integers
   Just the local submatrices */

   for ( i = 0; i < n; i++ )
   {
      for ( j = 0; j < n; j++ )
      {
#if (GMP32 == 0)
         a(i,j,MPALLOC) = abcnumsize;
         a(i,j,MPSIZE)  = modsize;
#else
         a(i,j,MPALLOC) = abcnumsize + (modsize << 32);
#endif
#ifdef GMP
         a(i,j,MPPTR)   = (mplong) & a(i,j,MPDATA);
#endif
#if (GMP32 == 0)
         b(i,j,MPALLOC) = abcnumsize;
         b(i,j,MPSIZE)  = modsize;
#else
         b(i,j,MPALLOC) = abcnumsize + (modsize << 32);
#endif
#ifdef GMP
         b(i,j,MPPTR)   = (mplong) & b(i,j,MPDATA);
#endif
      }
   }

/* RANWDS 64-bit random binary words will fill a row of A and B */

   ranwds = ( (2*n*modbits) + 63 ) / 64;

   tempbits = (unsigned mplong *) malloc( 8 * ranwds + 16 );

   for ( i = 0; i < n; i++ )
   {
      /* Generate random bits for this row # of A and B in the array tempbits */

      irow = 1 + i;  /* add 1 to agree with Fortran subr */

      get_my_bdata64( irow, 1, 1, ranwds, 1, seed, tempbits, &myfirst, &mylast);

      /* These counters keep track of position in tempbits */
      wordctr = 0;
      bitctr =  0;

      for ( j = 0; j < n; j++ )
      {
         /* Fill the low order words of A(I,J) */

         for ( k = MPDATA; k < MPDATA + fullwords; k++ )
         {
#if (BITSPERWORD==64)
            /* need this test in case << 64 is a no-op (ex: Alpha) */
            if ( bitctr == 0 )
               temp = tempbits[wordctr];
            else
	       temp =
               ( ( tempbits[wordctr]   >> bitctr ) |
                 ( tempbits[wordctr+1] << 64-bitctr ) );
            a(i,j,k) = temp;
            wordctr++;
#else
            if ( bitctr == 0 )
               temp = tempbits[wordctr];
            else
	       temp =
               ( ( tempbits[wordctr]   >> bitctr ) |
                 ( tempbits[wordctr+1] << 64-bitctr ) );
            a(i,j,k) = temp & BitMask;
            bitctr += BITSPERWORD;
            if ( bitctr >= 64 )
            {
               wordctr++;
               bitctr -= 64;
            }
#endif
         }

         /* Partially fill the last word of A(I,J) with EXTRABITS bits */

         if (extrabits)
         {
            k = MPDATA + fullwords;

            if ( bitctr == 0 )
               temp = tempbits[wordctr];
            else
	       temp =
               ( ( tempbits[wordctr]   >> bitctr ) |
                 ( tempbits[wordctr+1] << 64-bitctr ) );
            a(i,j,k) = temp & extrabitmask;
            bitctr += extrabits;
            if ( bitctr >= 64 )
            {
               wordctr++;
               bitctr -= 64;
            }
         }

         /* Fill the low order words of B(I,J) */

         for ( k = MPDATA; k < MPDATA + fullwords; k++ )
         {
#if (BITSPERWORD==64)
            if ( bitctr == 0 )
               temp = tempbits[wordctr];
            else
	       temp =
               ( ( tempbits[wordctr]   >> bitctr ) |
                 ( tempbits[wordctr+1] << 64-bitctr ) );
            b(i,j,k) = temp;
            wordctr++;
#else
            if ( bitctr == 0 )
               temp = tempbits[wordctr];
            else
	       temp =
               ( ( tempbits[wordctr]   >> bitctr ) |
                 ( tempbits[wordctr+1] << 64-bitctr ) );
            b(i,j,k) = temp & BitMask;
            bitctr += BITSPERWORD;
            if ( bitctr >= 64 )
            {
               wordctr++;
               bitctr -= 64;
            }
#endif
         }

         /* Partially fill the last word of B(I,J) with EXTRABITS bits */

         if (extrabits)
         {
            k = MPDATA + fullwords;

            if ( bitctr == 0 )
               temp = tempbits[wordctr];
            else
	       temp =
               ( ( tempbits[wordctr]   >> bitctr ) |
                 ( tempbits[wordctr+1] << 64-bitctr ) );
            b(i,j,k) = temp & extrabitmask;
            bitctr += extrabits;
            if ( bitctr >= 64 )
            {
               wordctr++;
               bitctr -= 64;
            }
         }
      } /* loop on j */

   } /* loop on i */

   free ( tempbits );

/*
 * Go through the entries of A and B.
 * If the top order word(s) = zero, reset the usage counter.
 * This normalizes the form of the multiple precision number.
 */
   for ( i = 0; i < n; i++ )
   for ( j = 0; j < n; j++ )
   {
#if (GMP32 == 0)
      asize = a(i,j,MPSIZE);
#else
      asize = a(i,j,MPALLOC) >> 32;
#endif
      if (asize)
      {
         /* work backwards from highest-order limb */
         for ( k = numsize - 1; k >= MPDATA; k-- )
         {
            if ( a(i,j,k) == 0 )
               asize--;
            else
               break;
         }
#if (GMP32 == 0)
         a(i,j,MPSIZE) = asize;
#else
         a(i,j,MPALLOC) = (a(i,j,MPALLOC) & BitMask32) + (asize << 32);
#endif
      }

#if (GMP32 == 0)
      bsize = b(i,j,MPSIZE);
#else
      bsize = b(i,j,MPALLOC) >> 32;
#endif
      if (bsize)
      {
         for ( k = numsize - 1; k >= MPDATA; k-- )
         {
            if ( b(i,j,k) == 0 )
               bsize--;
            else
               break;
         }
#if (GMP32 == 0)
         b(i,j,MPSIZE) = bsize;
#else
         b(i,j,MPALLOC) = (b(i,j,MPALLOC) & BitMask32) + (bsize << 32);
#endif
      }

   } /* end loop on i, j */

   /* if very verbose, test for consistency by printing the filled matrices */
   if ( verbose > 1 )
   {
      printf( "S11 Consistency Check - Matrix A: \n" );
      for ( i = 0; i < n; i++ )
      {
         printf( "Row %d \n", i );
         for ( j = 0; j < n; j++ )
         {
            mpwrite( &a(i,j,0) );
            printf ( "\n" );
         }
      }

      printf( "S11 Consistency Check - Matrix B: \n" );
      for ( i = 0; i < n; i++ )
      {
         printf( "Row %d \n", i );
         for ( j = 0; j < n; j++ )
         {
            mpwrite( &b(i,j,0) );
            printf ( "\n" );
         }
      }

      /* for debugging, may just want to see the matrices and exit */
      /* exit(0); */

   } /* verbose>1 */

}