#include <stdio.h>
#include "matmul.h"
/****   Copyright (c) 2001-2003  
 ****   Seung-Yeop Lee for the MPI code and 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
 ****/

/* 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) 2003    */
void Newpmatmul(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,s,	 /* temporary index */
  		ppi,ppj, /* index and processor group of pid */
		ni,	 /* n/sqrt(p) */
		size;	 /* size of a and b */
  MPI_Request   request;
  MPI_Status    status;
  int           temp,temp1;
  register int  Num;

  /* 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;

for(l=0;l<p_sqrt;++l) {

  Num = l+p_sqrt;
  temp = ((2*p_sqrt-ppi-ppj+l)%p_sqrt)*p_sqrt+ppi;
  temp1 = ((2*p_sqrt-ppj-ppi+l)%p_sqrt)+ppj*p_sqrt;

  MPI_Irecv(a,ni*ni,MPI_DOUBLE,MPI_ANY_SOURCE,l,MPI_COMM_WORLD,&request); 
  MPI_Send(A,ni*ni,MPI_DOUBLE,temp,l,MPI_COMM_WORLD); 
  MPI_Wait(&request, &status);
  AIOBARRIER();

  MPI_Irecv(b,ni*ni,MPI_DOUBLE,MPI_ANY_SOURCE,Num,MPI_COMM_WORLD,&request);
  MPI_Send(B,ni*ni,MPI_DOUBLE,temp1,Num,MPI_COMM_WORLD); 
  MPI_Wait(&request, &status);
  AIOBARRIER();

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

/*
  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;
    }
*/
}

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

void Newpomatmul(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;   
  MPI_Request   request;
  MPI_Status    status;
  int   	temp,temp1;	
  register int  Num;

  /* 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;
  for (l=0;l<p_sqrt;++l) {
    /* Efficient non-conflicting communication */
    /* A 2n^2/p relation is communicated */
    Num = l+p_sqrt;    
    temp = ((2*p_sqrt-ppj-ppi+l)%p_sqrt)*p_sqrt+ppi;
    temp1 = ((2*p_sqrt-ppj-ppi+l)%p_sqrt)+ppj*p_sqrt;

    MPI_Irecv(a,ni*ni,MPI_DOUBLE,MPI_ANY_SOURCE,l,MPI_COMM_WORLD,&request); 
    MPI_Send(A,ni*ni,MPI_DOUBLE,temp,l,MPI_COMM_WORLD); 
    MPI_Wait(&request, &status);
    AIOBARRIER();

    MPI_Irecv(b,ni*ni,MPI_DOUBLE,MPI_ANY_SOURCE,Num,MPI_COMM_WORLD,&request);
    MPI_Send(B,ni*ni,MPI_DOUBLE,temp1,Num,MPI_COMM_WORLD); 
    MPI_Wait(&request, &status);
    AIOBARRIER();

    /* 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,t1=0;j<ni;++j,t1+=ni) { 
     for (i=0,t2=0;i<ni;++i,t2==ni) {
       tmp=0.0;
       for (k=0;k<ni;++k) 
        tmp += (b[t1+k]*a[t2+k]); 
       C[t1+i] += tmp;
     }
    }
  }

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