/****   Copyright (c) 1997-2002 Alexandros V. Gerbessiotis
 ****
 ****  Permission to use, copy, modify, and distribute this 
 ****  software, and to incorporate it, in whole or in part, 
 ****  into other software, is hereby granted without fee, 
 ****  provided that 
 ****  (1) the above copyright notice and this permission 
 ****      notice appear in all copies of the source code, 
 ****      and the above copyright notice appear in clearly 
 ****      visible form on all supporting documentation
 ****      and distribution media, and 
 ****  (2) any redistribution of the software, in original 
 ****      or modified form, be without fee and subject to 
 ****      these same conditions. 
 ****   No guarantee is offered that the code works as 
 ****   advertised or is absent of any, even damaging, 
 ****   side-effects; use this code on your own personal risk 
 ****/
#include "matmul.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 genmatmul(double *A,double *B,double *C,int n,int pi,int pj) 
{
register int	nprocs; /* number of processors */
register int	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 */
  nprocs = AINPROCS(); pid = AIPID();
  ni=n/pi; 
  nj=n/pj;
#ifdef SANITY_CHECK
  check_if_valid((ni*pi),n,"gmatmul.c   ","n/pi must be integer");
  check_if_valid((nj*pj),n,"gmatmul.c   ","n/pj must be integer");
  /* nprocs must be equal to pi*pj or error */
  check_if_valid((pi*pj),nprocs,"gmatmul.c   ","nprocs != pi*pj");
#endif

  /* 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);
#ifdef SANITY_CHECK
  check_if_null((void *)a,"in gmatmul.c   ","a");
#endif
  b = (double *)malloc(sj);
#ifdef SANITY_CHECK
  check_if_null((void *)b,"in gmatmul.c   ","b");
#endif
  tA= (double *)malloc(ni*nj*sizeof(double));
#ifdef SANITY_CHECK
  check_if_null((void *)tA,"in gmatmul.c   ","tA");
#endif


  /* 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 */
  AIREGISTER(a,si);
  AISYNC();  
  for(l=0;l<pj;l++) {
    p_b = ppi;      /* base processor */
    p_o = (ppj+l)%pj;/* round-robin reading */ 
    p_b += p_o*pi;/* pid reads pid,pid+1,..,nprocs-1,0,..,pid-1 */
    for(j=0;j<ni;j++) {
     a_b = j*n + ppj*nj;/* base dest.address of read data in a */
     AIHPPUT(p_b,&tA[j*nj],a,a_b*sizeof(double),nj*sizeof(double));
    }
  }
  AISYNC();  

  /* Read pi blocks of matrix B locally,  as needed */
  AIREGISTER(b,sj);
  AISYNC();  
  for(l=0;l<pi;l++) {
    p_b = (ppj*pi);
    p_o = (ppi+l)%pi;
    p_b += p_o;
    for(j=0;j<nj;j++) {
     a_b = j*n + ppi*ni;
     AIHPPUT(p_b,&B[j*ni],b,a_b*sizeof(double),ni*sizeof(double));
    }
  }
  AISYNC();  

  /* Do the multiplication locally */
  for (j=0,t1=0;j<nj;j++,t1+=n) 
   for (i=0,t2=0;i<ni;i++,t2+=n) {
     for (k=0,tmp=0.0;k<n;k++) 
       tmp += a[t2+k] * b[t1+k]; /*C[i,j]+=tA[k,i]*B[k,j]*/
     C[j*ni+i] = tmp;
   }

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

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

