#include <stdio.h>
#include "matmul.h"

/* It implements a memory efficient matrix multiplication routine.
 * Matrices A.B and C are stored in terms of block of size
 * dimension (n/sqrt(p) X n/sqrt(p)), p # of processors.
 * Algorithm requires n/sqrt(p) to be an integer.
 * Each processor in sqrt(p) supersteps reads one block from A and one
 * from B, and computes a partial result of the multiplication  C=AxB
 *  Size of memory per processor instead is  3n^2 /p
 *  All amtrices 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 in +omatmul is transposed prior to computation; this
 *  may increase efficiency due to locality - caching - effects.
 */

/*  Copyright (c) 2001-2002  Alexandros V. Gerbessiotis  */
void pmatmul(double *A,double *B,double *C,int n,int p_sqrt) 
{
  register int	pid;    /* who am i? */
  double        *a,*b;  /* local storage  */
  register double tmp;  /* temporary variable */
  register int 	i,j,k,l; /* indices */
  register int  t,	 /* temporary index */
  		ppi,ppj, /* index and processor group of pid */
		ni,	 /* n/sqrt(p) */
		size;	 /* size of a and b */

  /* Dimension of A and B, a and b */
  pid=AIPID(); 
  ni=n/p_sqrt; 
  size=ni*ni*sizeof(double); 
  ppi=pid % p_sqrt; /*index of pid within its (column group */
  ppj=pid / p_sqrt; /* (column-block) group pid belongs to */
  
  /* Allocation and checking */
  a = (double *)malloc(size);
#ifdef SANITY_CHECK
  check_if_null((void *)a,"in pmatmul      ","a");
#endif
  b = (double *)malloc(size);
#ifdef SANITY_CHECK
  check_if_null((void *)b,"in pmatmul      ","b");
#endif

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

  AIOREGISTER(a,size); 
  AIOREGISTER(b,size); 
  AIOSYNC();
for(l=0;l<p_sqrt;++l) {
  AIOHPPUT((((2*p_sqrt-ppi-ppj+l)%p_sqrt)*p_sqrt+ppi),A,a,0,size);
  AIOSYNC();

  AIOHPPUT((((2*p_sqrt-ppj-ppi+l)%p_sqrt)+ppj*p_sqrt),B,b,0,size);
  AIOSYNC();

  /* Computation */
  for(j=0;j<ni;++j) 
    for(i=0,t=j*ni;i<ni;++i) {
     for(k=0,tmp=0.0;k<ni;++k) 
      tmp += (b[t+k]*a[k*ni+i]); 
      C[t+i] += tmp;
    }
}

 AIODEREGISTER(b) ;
 AIODEREGISTER(a) ;

 free((void *)b) ;
 free((void *)a) ;
}


/*  Copyright (c) 2001  Alexandros V. Gerbessiotis  */
/* Similar to the previous routine except for a transposition of A.
 * Faster by 50% on the average.
 * Modification by Alexandros V. Gerbessiotis, April 1997, April 2001
 */

void pomatmul(double *A,double *B,double *C,int n,int p_sqrt ) 
{
  register int    pid;    /* who am i? */
  double          *a,*b;	/* local storage  */
  register double tmp;	/* temporary variable */
  register int    i,j,k,l;/* indices */
  register int    t1,t2,	/* temporary indices */
		ppi,ppj,/* index and processor group of pid */
		ni,	/* n/sqrt(p) */
                size;   

  /* Dimension of A and B, a and b */
  pid=AIPID();
  ni=n/p_sqrt; 
  ppi=pid % p_sqrt; /*index of pid within its group */
  ppj=pid / p_sqrt; /* proc group pid belongs to */
  
  size= ni*ni*sizeof(double);
  /* Allocation and checking */
  a = (double *)malloc(size);
#ifdef SANITY_CHECK
  check_if_null((void *)a,"in pomatmul     ","a");
#endif
  b = (double *)malloc(size);
#ifdef SANITY_CHECK
  check_if_null((void *)b,"in pomatmul     ","b");
#endif

  /* Initialization */
  for (j=0;j<ni;++j) 
    for (i=0,t1=j*ni;i<ni;++i) 
      C[t1 + i]= 0.0;

    AIOREGISTER(a,  size);
    AIOREGISTER(b,  size);
    AIOSYNC();
  for (l=0;l<p_sqrt;++l) {
    /* Efficient non-conflicting communication */
    /* A 2n^2/p relation is communicated */
    AIOHPPUT((((2*p_sqrt-ppj-ppi+l)%p_sqrt)*p_sqrt+ppi),A,a,0,size);
    AIOSYNC();

    AIOHPPUT((((2*p_sqrt-ppj-ppi+l)%p_sqrt)+ppj*p_sqrt),B,b,0,size);
    AIOSYNC();
    /* Transpose a */
    for(i=0,t2=0;i<ni;++i,t2+=ni) 
     for(j=i,t1=i*ni;j<ni;++j,t1+=ni)  {
      tmp=a[t2+j];
      a[t2+j]=a[t1+i];
      a[t1+i]=tmp;
     }

    /* Multiplication */
    for (j=0;j<ni;++j) { 
     t1=j*ni;
     for (i=0;i<ni;++i) {
       t2=i*ni;
       tmp=0.0;
       for (k=0;k<ni;++k) 
        tmp += (b[t1+k]*a[t2+k]); 
       C[t1+i] += tmp;
     }
    }
  }

  AIODEREGISTER(b) ;
  AIODEREGISTER(a) ;

  free((void *)b) ;
  free((void *)a) ;
}

/*  Copyright (c) 1997  Alexandros V. Gerbessiotis and
 *                      Constantinos J. Siniolakis
 */
void gmatmul(double *A,double *B,double *C,int n,int p_sqrt) 
{
  register int	pid;    /* who am i? */
  double        *a,*b; /* local storage  */
  register double tmp;  /* temporary variable */
  register int 	i,j,k,l; /* indices */
  register int  t,	 /* temporary index */
		ppi,ppj, /* index and processor group of pid */
		ni,	 /* n/sqrt(p) */
		size;	 /* size of a and b */

  /* Dimension of A and B, a and b */
  pid=AIPID(); 
  ni=n/p_sqrt; 
  size=ni*ni*sizeof(double); 
  ppi=pid % p_sqrt; /*index of pid within its proc. group */
  ppj=pid / p_sqrt; /* (column-block) group pid belongs to */
  
  /* Allocation and checking */
  a = (double *)malloc(size);
#ifdef SANITY_CHECK
  check_if_null((void *)a,"in gmatmul      ","a");
#endif
  b = (double *)malloc(size);
#ifdef SANITY_CHECK
  check_if_null((void *)b,"in gmatmul      ","b");
#endif

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

    AIOREGISTER(A,size);
    AIOREGISTER(B,size);
    AIOSYNC();
  for (l=0;l<p_sqrt;++l) {
    AIOHPGET((((ppi+ppj+l)%p_sqrt)*p_sqrt+ppi),A,0,a,size);
    AIOSYNC();
    AIOHPGET((((ppi+ppj+l)%p_sqrt)+ppj*p_sqrt),B,0,b,size);
    AIOSYNC();

    /* Computation */
    for (j=0;j<ni;++j) 
      for (i=0,t=j*ni;i<ni;++i) {
        for (k=0,tmp=0.0;k<ni;++k) 
          tmp += (b[t+k]*a[k*ni+i]); 
        C[t+i] += tmp;
      }
  }

  AIODEREGISTER(B);
  AIODEREGISTER(A);

  free((void *)b);
  free((void *)a);
}

/*  Copyright (c) 1997  Alexandros V. Gerbessiotis and
 *                      Constantinos J. Siniolakis
 */
/* Similar to the previous routine except for a transposition of A.
 * Faster by 50% on the average.
 * Modification by Alexandros V. Gerbessiotis, April 1997
 */

void gomatmul(double *A,double *B,double *C,int n,int p_sqrt ) 
{
  register int  pid;    /* who am i? */
  double        *a,*b;	/* local storage  */
  register double tmp;	/* temporary variable */
  register int  i,j,k,l;/* indices */
  register int  t1,t2,	/* temporary indices */
		ppi,ppj,/* index and processor group of pid */
		ni,	/* n/sqrt(p) */
                size;   

  /* Dimension of A and B, a and b */
  pid=AIPID();
  ni=n/p_sqrt; 
  ppi=pid % p_sqrt; /*index of pid within its group */
  ppj=pid / p_sqrt; /* (column-block) group pid belongs to */
  
  size= ni*ni*sizeof(double);
  /* Allocation and checking */
  a = (double *)malloc(size);
#ifdef SANITY_CHECK
  check_if_null((void *)a,"in gomatmul     ","a");
#endif
  b = (double *)malloc(size);
#ifdef SANITY_CHECK
  check_if_null((void *)b,"in gomatmul     ","b");
#endif


  /* Initialization */
  for (j=0;j<ni;++j) 
    for (i=0,t1=j*ni;i<ni;++i) 
      C[t1 + i]= 0.0;

    AIOREGISTER(A,  size);
    AIOREGISTER(B,size);
    AIOSYNC();
  for (l=0;l<p_sqrt;++l) {
    /* Efficient non-conflicting communication */
    /* A 2n^2/p relation is communicated */
    AIOHPGET((((ppi+ppj+l)%p_sqrt)*p_sqrt+ppi),A,0,a,size);
    AIOSYNC();
    AIOHPGET((((ppi+ppj+l)%p_sqrt)+ppj*p_sqrt),B,0,b,size);
    AIOSYNC();

    /* Transpose a */
    for (i=0,t2=0;i<ni;++i,t2+=ni) 
      for (j=i,t1=i*ni;j<ni;++j,t1+=ni)  {
        tmp=a[t2+j];
        a[t2+j]=a[t1+i];
        a[t1+i]=tmp;
      }

    /* Multiplication */
    for (j=0;j<ni;++j) { 
      t1=j*ni;
      for (i=0;i<ni;++i) {
        t2=i*ni;
        tmp=0.0;
        for (k=0;k<ni;++k) 
          tmp += (b[t1+k]*a[t2+k]); 
        C[t1+i] += tmp;
      }
   }
  }

  AIODEREGISTER(B) ;
  AIODEREGISTER(A) ;

  free((void *)b) ;
  free((void *)a) ;
}
