/**************************************************************************
*
*   Subroutine P11 performs the matrix multiplication C = A x B
*
*   Parameters:
*
*    Provided by calling routine:
*       A       =  Matrix of multiple precision numbers
*       B       =  Matrix of multiple precision numbers
*	MODULUS =  Large integer modulus, expressed as an array of integers
*       Q       =  The size of the grid of PEs, Q x Q
*       M       =  The size of the sub-matrices on each PE
*       NUMSIZE =  The number of words used to store an MP integer
*
*    Returned by this routine:
*       C       =  Matrix of multiple precision numbers
*
***************************************************************************
*
*    This version for non-shared-memory systems uses Fox`s algorithm for
*    matrix multiplication.
*
**************************************************************************/
#include "bench11.h"

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

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

#define TIMEIT 0

/* use macros to address 3-dim arrays with variable dimension */

#define c(I,J,K)	c[ (K) + numsize * ( (J) + n * (I) ) ]
#define a(I,J,K)	a[ (K) + numsize * ( (J) + n * (I) ) ]
#define b(I,J,K)	b[ (K) + numsize * ( (J) + n * (I) ) ]
#define c(I,J,K)	c[ (K) + numsize * ( (J) + n * (I) ) ]


void p11 ( mplong a[], mplong b[], mplong c[], mplong modulus[],
	   int n, int numsize )
{

   /* Multiply into tempmp, accumulate sums in tempc */
   mplong tempmp [ 2*maxnumsize ];
   mplong tempc  [ 2*maxnumsize ];

   int i, j, k, abcsize, tempsize;
   double wall(), T0, multtime = 0;  /* used if TIMEIT is 1 */

   void mult();

   abcsize    = numsize - MPHL;
   tempsize   = 2*maxnumsize - MPHL;

   /* initialize c to hold product a*b */

   for ( i = 0; i < n; i++ )
   {
      for ( j = 0; j < n; j++ )
      {
         c(i,j,MPALLOC) = abcsize;
#if (GMP32 == 0)
         c(i,j,MPSIZE)  = 0;
#endif
#ifdef GMP
         c(i,j,MPPTR)   = (mplong) & c(i,j,MPDATA);
#endif
      }
   }

   /* initialize temporaries */

   tempmp[MPALLOC] = tempsize;
#if (GMP32 == 0)
   tempmp[MPSIZE]  = 0;
#endif
   tempc[MPALLOC]  = tempsize;
#ifdef GMP
   tempmp[MPPTR] = (mplong) & tempmp[MPDATA];
   tempc[MPPTR]  = (mplong) & tempc[MPDATA];
#endif

   /* do the matrix multiplication */

/* extern int VERBOSE; */

   for ( i = 0; i < n; i++ )
      for ( j = 0; j < n; j++ )
      {
         /* zero tempc - accumulate sum during loop */
#if (GMP32 == 0)
         tempc[MPSIZE] = 0;
#else
         /* erase upper 32 bits of dual-use word 0 */
         tempc[MPALLOC] = tempc[MPALLOC] & BitMask32;
#endif

/* VERBOSE = 2; */

         for ( k = 0; k < n; k++ )
         {
            mpmult ( tempmp, &a(i,k,0), &b(k,j,0) );
            mpadd  ( tempc, tempc, tempmp );
         }

         mpmod ( tempc, tempc, modulus );

#ifdef GMP
         mpset ( &c(i,j,0), tempc );
#else
         c(i,j,MPSIZE) = tempc[MPSIZE];
         for ( k = MPDATA; k < MPDATA + tempc[MPSIZE]; k++ )
            c(i,j,k) = tempc[k];
#endif
      }
}