#include <stdio.h>
#include "bsp.h"
#include "matrix.h"

/* It implements a memory inefficient matrix multiplication routine,
 * on the BSP model suggested by Valiant (C.ACM 33(8),pp103-111,Aug 1990)
 * Initially, two global $n \times n$ matrices $A_g $ and $B_g$ are
 * distributed among  $p$ processors. The $p$ processors are divided 
 * into $pj$ groups of $pi$ processors each.  Element $(i,j)$ of $A_g$ 
 * or $B_g$ is stored in the $i/pi$-th processor of the $j/pj$ processor 
 * group, that is, processor $(j/pj)*pi + i/pi$.
 * Function {\tt multiply\_par} requires six arguments, the two input 
 * (distributed among the processors) matrices $A$ and $B$ each
 * of dimension  (n/pi X n/pj), the result (already malloced)
 * of dimension (n/pi X n/pj), the dimension n, pi and pj.
 * $A$, $B$ and $C$ are ANSI-C pointers to a {\tt double} data type.
 * We store matrices in the form of an on dimensional array. This way, 
 * element $(i,j)$ of a two-dimensional matrix is stored in position
 * $j*n+i$ of say, $A$. All indices are in the range $0, \ldots, n-1$.
 *
 * Each processor in one superstep reads the all the block it needs
 * to perform the matrix multiplication, and then it performs it.
 *  Size of memory per processor instead of 3n^2/p is 
 *  4n^2 /p + 2n^2/q, where q=min{pi,pj}.
 *  All matrices are one-dimensional. Thus, element (i,j) of an
 *  nXn matrix X is the (n*j+i)-th element of X  X[n*j+i]; a column
 *  major order is used to store the matrix.
 *  Matrix A is transposed prior to communication and computation; this
 *  may increase efficiency due to locality - caching - issues.
 */

void multiply_par(double *A,double *B,double *C,int n,int pi,int pj) 
{
register int	nprocs=bsp_nprocs(); /* number of processors */
register int	pid=bsp_pid();       /* who am i ? */
double 		*a,*b,*tA;	     /* aux space for A, B and transpose(A) */
register double	tmp;		     /* aux. variable in matrix mult.  */
register int 	i,j,k,l;	     /* Indices */
register int   	t,  		     /* Auxiliary variable  */
		ppi,ppj,  	     /* Index within and processor group of pid */
		ni,nj,               /* A and B are n/pi X n/pj matrices */
		si,sj,               /* Total size in bytes of a and b  */
		p_b,  		     /* Base processor id. in communication  */
		p_o,		     /* Offset from base   in communication  */
		a_b,		     /* Base address of communicated data    */
		t1,t2;		     /* Indices/Addresses in multiplication  */

  /* Dimension of A, B and C ; must be integer */
  ni=n/pi; 
  nj=n/pj;
  check_if_valid((ni*pi),n,"multiply_par","n/pi must be integer");
  check_if_valid((nj*pj),n,"multiply_par","n/pj must be integer");

  /* nprocs must be equal to pi*pj or error */
  check_if_valid((pi*pj),nprocs,"multiply_par","nprocs != pi*pj");

  /* size of a and b (in bytes) respectively */
  si=ni*n*sizeof(double);
  sj=n*nj*sizeof(double);
  ppi=pid % pi;  /* index of pid within its (column block) group */
  ppj=pid / pi;  /* (column-block) group pid belongs to */
  
  /* Allocation and checking */
  a = (double *)malloc(si);
  check_if_null((void *)a,"in multiply_par","a");
  b = (double *)malloc(sj);
  check_if_null((void *)b,"in multiply_par","b");
  tA= (double *)malloc(ni*nj*sizeof(double));
  check_if_null((void *)tA,"in multiply_par","tA");

  bsp_pushregister(a,si);
  bsp_pushregister(b,sj);
  bsp_sync();

  /* Initialize result */
  for (j=0;j<nj;j++) 
    for (i=0,t=j*ni;i<ni;i++) 
      C[t + i]= 0.0;

   /* Transpose A first for efficiency */
      for (j=0,t1=0;j<ni;j++,t1+=nj) 
        for (i=0,t2=0;i<nj;i++,t2+=ni)  {
	    tA[t1+i] = A[t2+j];
        }


  /* read pj  blocks of matrix A locally, as needed */
  for(l=0;l<pj;l++) {
/* Implementation Point 1 */
      }
   }
   bsp_sync();

  /* Read pi blocks of matrix B locally,  as needed */
  for(l=0;l<pi;l++) {
/* Implementation Point 2 */
   }
   bsp_sync();

   /* Do the multiplication locally */
/* Implementation Point 3 */

  /* De-registration */
  bsp_popregister(b) ; 
  bsp_popregister(a) ;

  /* De-allocation */
  free((void *)tA) ;
  free((void *)b) ;
  free((void *)a) ;
}

