#include <stdafx.h>
#include <winsock2.h>
/*****************************************************************
**Matrix Multiplication Testing Functions(contains a main() )
**Copyright (c) 1997,2001-2003  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  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
**/
#undef  WIN2K
#define WIN2K
#include <stdio.h>
#include <memory.h>
#include <math.h>
#ifndef WIN2K
#include <sys/utsname.h>
#include <unistd.h>
#endif
#include <stdlib.h>
#include "allmult.h"
#ifndef BRD_CUTOFF
#define BRD_CUTOFF      1200
#endif
#ifndef PPF_CUTOFF
#define PPF_CUTOFF      200
#endif

#define		SQRT(X)	((int)ceil(sqrt((double) (X))))
#define		SQRT_F(X)	((int)floor(sqrt((double) (X))))
#define         CEIL(X,Y)       ((X != 0)?(1+ (X-1)/Y):0)  
#define         FLOOR(X,Y)      (X/Y)  
#define         LEFT(x,y)       (x-(x/y)*y)

void check_if_null(void *x,char *s,char *t);
void __my_estimate_bsp_values();
void bspbroad(int fromp, char *from, char *to, int size);
void foldmax(double *z, double *x, double *y);
void foldmin(double *z, double *x, double *y);
void bspfold(void (*oper)(void *,void*,void*),char *from,char *to,int size);

#define PRINTMAT(X)   bspfold((void (*)(void *,void *,void *)) foldmin,(char *)&t1,(char *)&t3,sizeof(double)); bspfold((void (*)(void *,void *,void *)) foldmax,(char *)&t1,(char *)&t4,sizeof(double)); if (pid == 0) printf("**"#X" Time %6.2f(%6.2f)-%6.2f(%6.2f) (Mflop %5.2f-%5.2f)\n", t3,(t3)/runs,t4,(t4)/runs, runs*a/(nprocs*t4*1000000.0), runs*a/(nprocs*t3*1000000.0));fflush(stdout);

void pmatmul(double *,double *,double *, int ,int);
void pomatmul(double *, double *, double *,int,int);
void gmatmul(double *, double *, double *,int,int);
void gomatmul(double *,double *, double *,int,int);
void genmatmul();
void init_blocked(double **, int , int ,int ,int);
void print_mat(double *,int,int,char *);
void collect_blocked_matrix(double *,double *,int,int,int);

double __my_estimate_l = (double) 500;
double __my_estimate_g = (double) 0.5;
int    __my_estimate_nprocs = (int) 1;
int    __my_first_call=(int) 0;
int    __my_estimate_brd=(int) BRD_CUTOFF;
int    __my_estimate_ppf=(int) PPF_CUTOFF;


#ifdef   MPILIB
#define AIMPIWINMAX   128
MPI_Comm comm;
MPI_Win  win;
MPI_Win  aiwin[AIMPIWINMAX];
void*    bspmpi[AIMPIWINMAX];
int bspmpi_max = AIMPIWINMAX;
int bspmpi_first=1;
int bspmpicur_max = 0;
int mpi_pid = -1;
int mpi_nprocs = -1;
double mpi_bsp_l = 300.0;
double mpi_bsp_g = 0.5;
double mpi_basetime=0.0;
#endif

#ifdef MPILIB
int aisearch(void * address) {
 register int i;

 for (i=bspmpi_max; i >=0; i-- ){
   if (bspmpi[i] == address) {
        return(i);
   }
 }
 return(-1);
}

int aiinsert(void *address) {
 register int i;
 if (bspmpi_first == 1) {
    bspmpi_first =0;
    for(i=0;i<bspmpi_max;i++) {
          bspmpi[i]=NULL;
    }
 }
 for(i=0; i< bspmpi_max ; ++i){
  if (bspmpi[i] == NULL)  {
    bspmpi[i]=address;
    bspmpicur_max++;
    return (i);
  }
 }
 /* runout of space */
 return(-1);
}

void aidelete(void *address) {
 register int i;
 for(i=0; i< bspmpi_max; ++i) {
   if (bspmpi[i]==address){
     bspmpi[i]= NULL;
     bspmpicur_max --;
   }
 }
}
#endif
void check_if_null(void *x,char *s,char *t)
{
  if (NULL == x) {
     printf("%s:: Null pointer error (%s)\n",s,t);
     AIABORT("check_if_null");
  }
}

void __my_estimate_bsp_values()
{
        __my_estimate_l = AIL();
        __my_estimate_g = AIG();
        __my_estimate_nprocs = AINPROCS();
        __my_estimate_brd = BRD_CUTOFF;
        __my_estimate_ppf = PPF_CUTOFF;
}

/* One superstep broadcasting to communicate command line args */
void bspbroad(int fromp, char *from, char *to, int size)
{
 int pid;
 int nprocs;
 int i;

 pid=AIPID(); /* processor id */
 nprocs=AINPROCS(); /* number of processors */
 AIREGISTER(to,size);
 AISYNC();
 AICOMMIT();
 for (i = 0; i < nprocs;i++) {
   if (pid==fromp)
        AIPUT(i,from,to,0,size);
 }
 AISYNC();
 AIDEREGISTER((void*)to);
}

/* operators for bspfold */

void foldmax(double *z, double *x, double *y)
{
  *z=(((*x)>=(*y))?(*x):(*y));
}

void foldmin(double *z, double *x, double *y)
{
  *z=(((*x)<(*y))?(*x):(*y));
}



/* One superstep folding operation to compute global min/max of timing data */
void 
  bspfold(void (*oper)(void *,void*,void*),char *from,char *to,int size)
{
 int  pid; 
 int  nprocs; 
 int  i;
 char *list;
 int  lsize;

 pid=AIPID(); 
 nprocs=AINPROCS(); 
 if (1==nprocs)  /* if local, copy */
    memcpy(to,from,size*sizeof(char));
 else {    /* if global, allocate memory O(p) per processor */
    lsize = nprocs*size; 
    list=(char *) malloc(lsize*sizeof(char));
    check_if_null(list,"bspfold","list");
    AIREGISTER(list,lsize);
    AISYNC();
    AICOMMIT();
 }
 /* each proc sends its from to every other proc; results stored in list */
 for (i = 0; i < nprocs;i++) {  
   AIPUT(i,from,list,pid*size,size);
 }
 AISYNC();
 /* apply operator to local list */
 for (i = 0; i < nprocs-1;i++) {
  oper(to,(list+i*size),(list+(i+1)*size));
  memcpy((list+(i+1)*size),to,size);
 }
    /* free memory and deregister */
 
    AIDEREGISTER(list);
    free((void *)list);
}
void check_if_valid(int x,int y, char *s,char *t)
{
  if (x != y) {
    printf("%s::Integer parameter %d should be equal to %d (%s)\n",s,x,y,t);
    AIABORT("check_if_valid");
  }
}

/* MAIN function */

int main(int argc, char *argv[])
{  
 int     k,n,nprocs,runs,debug,pid,pi,pj,size;
 double  *mat,*mat1,*mat2,*mat3;
 double t1,t2,t3,t4,a;
 char hostname[512];

 AIBEGIN(AINPROCS());
 pid = AIPID(); nprocs = AINPROCS();

 gethostname(hostname,sizeof(hostname));
#ifdef SANITY_CHECK
 /* Show hostnames of processors; */
 for(i=0;i<nprocs;i++) {
   if (i==pid) {
      printf("Processor %d of %d [%s]\n",pid,nprocs,hostname);
      fflush(stdout);
   }
 }
#endif
 AIBARRIER();
  /* Command line processing */
 n=8;debug=0;runs=1; pi=pj=1;
 if (0 == pid ){
 if (6 != argc){
   printf("Usage: %s dimension pi pj runs debug\n",argv[0]);
   fflush(stdout);
   AIABORT("exit");
 } 
 else {
  n=atoi(argv[1]);
  pi=atoi(argv[2]);
  pj=atoi(argv[3]);
  runs=atoi(argv[4]);
  debug=atoi(argv[5]);
 }
    if (pi != pj) {
      printf(" pi must be equalt to pj (pi=%d pj=%d)\n",pi,pj);
      AIABORT("exit pi!=pj");
    }
   }
	
   /* Parallel Command line processing */
   AIREGISTER(&n,sizeof(int));
   AIINIT();AICOMMIT();
      bspbroad(0,(char *)&n ,(char *)&n ,sizeof(int));
   AISYNC();
   AIDEREGISTER(&n );

   AIREGISTER(&pi,sizeof(int));
   AIINIT();AICOMMIT();
      bspbroad(0,(char *)&pi ,(char *)&pi ,sizeof(int));
   AISYNC();
   AIDEREGISTER(&pi );

   AIREGISTER(&pj,sizeof(int));
   AIINIT();AICOMMIT();
      bspbroad(0,(char *)&pj ,(char *)&pj ,sizeof(int));
   AISYNC();
   AIDEREGISTER(&pj );

   AIREGISTER(&debug,sizeof(int));
   AIINIT();AICOMMIT();
      bspbroad(0,(char *)&debug,(char *)&debug  ,sizeof(int));
   AISYNC();
   AIDEREGISTER(&debug);

   AIREGISTER(&runs,sizeof(int));
   AIINIT();AICOMMIT();
      bspbroad(0,(char *)&runs,(char *)&runs  ,sizeof(int));
   AISYNC();
   AIDEREGISTER(&runs);

   if ((nprocs-1)==pid) {
     printf("%s:Proc %2d of %d: n=%d (pi,pj)=(%d,%d) debug=%d runs=%d\n",
          argv[0],pid, nprocs,n,pi,pj,debug,runs);
     fflush(stdout);
   }

 /* Initialization */
  t1=t2=0.0;
  size = CEIL(n,pi)*CEIL(n,pj)*sizeof(double);

 /* Memory allocation */
  mat = (double *) malloc(size);
#ifdef SANITY_CHECK
  check_if_null((void *) mat,"in main of mmmain.c","mat");
#endif
 if ( 1 == debug) {
   mat3 = (double *) malloc(n*n*sizeof(double));
#ifdef SANITY_CHECK
   check_if_null((void *)mat3,"in main of mmmain.c","mat3");
#endif
 }
 /* Read/Initialize Test Matrix A */
 init_blocked(&mat1,n,pi,pj,1);
 if ((1)== debug) {
      collect_blocked_matrix(mat3,mat1,n,pi,pj);
      print_mat(mat3,n,1,"mat1");
 }
 /* Read/Initialize Test Matrix B */
 init_blocked(&mat2,n,pi,pj,0);
 if ((1)== debug) {
   collect_blocked_matrix(mat3,mat2,n,pi,pj);
   print_mat(mat3,n,1,"mat2");
 }


  /*************/
#ifdef ALEXG
#endif
  t1=0.0;
  AIBARRIER();
  t3=AITIME();  
   for(k=0;k<runs;k++) {
    pmatmul(mat1,mat2,mat,n,pi);
    AIBARRIER();
   }
  t4=AITIME(); 
  t1 +=(t4-t3);

  if ((1)== debug) {
   collect_blocked_matrix(mat3,mat,n,pi,pj);
   print_mat(mat3,n,1,"Put");
  }
  a= (2.0*n-1.0)*n*n;
  PRINTMAT("Matrix Mult.1(puts)");
  /*************/

  /*************/
  t1=0.0;
  AIBARRIER();
  t3=AITIME();  
   for(k=0;k<runs;k++) {
    pomatmul(mat1,mat2,mat,n,pi);
    AIBARRIER();
   }
  t4=AITIME(); 
  t1 +=(t4-t3);

  if ((1)== debug) {
   collect_blocked_matrix(mat3,mat,n,pi,pj);
   print_mat(mat3,n,1,"Puto");
  }
  a= (2.0*n-1.0)*n*n;
  PRINTMAT("Matrix Mult.2(puts)");
  /*************/
  /*************/
#ifdef ALEXG
#endif
  t1=0.0;
  AIBARRIER();
  t3=AITIME();  
   for(k=0;k<runs;k++) {
    gmatmul(mat1,mat2,mat,n,pi);
    AIBARRIER();
   }
  t4=AITIME(); 
  t1 +=(t4-t3);

  if ((1)== debug) {
   collect_blocked_matrix(mat3,mat,n,pi,pj);
   print_mat(mat3,n,1,"Get");
  }
  a= (2.0*n-1.0)*n*n;
  PRINTMAT("Matrix Mult.1(gets)");
  /*************/
  /*************/
  t1=0.0;
  AIBARRIER();
  t3=AITIME();  
   for(k=0;k<runs;k++) {
    gomatmul(mat1,mat2,mat,n,pi);
    AIBARRIER();
   }
  t4=AITIME(); 
  t1 +=(t4-t3);

  if ((1)== debug) {
   collect_blocked_matrix(mat3,mat,n,pi,pj);
   print_mat(mat3,n,1,"Geto");
  }
  a= (2.0*n-1.0)*n*n;
  PRINTMAT("Matrix Mult.2(gets)");
  /*************/

/*
  t1=0.0;
  AIBARRIER();
  t3=AITIME();  
   for(k=0;k<runs;k++) {
    genmatmul(mat1,mat2,mat,n,pi,pj);
    AIBARRIER();
   }
  t4=AITIME(); 
  t1 +=(t4-t3);

  if ((1)== debug) {
   collect_blocked_matrix(mat3,mat,n,pi,pj);
   print_mat(mat3,n,1,"mat=mat1Xmat2");
  }
  a= (2.0*n-1.0)*n*n;
  PRINTMAT("Matrix Mult.3(gen )");

*/

  free ((void *) mat2);
  free ((void *) mat1);
  if (1 == debug ) free ((void *) mat3);
  free ((void *) mat);

  AIEND();
  return(0);
}

 




/* Init_blocked */
void    init_blocked(double **mat_r,int n,int pi,int pj,int option) {

 register int    pid;
 register int    nprocs;
 register int    i,j,k,ni,nj,ppi,ppj,nki,nkj,size;
 double          *mat;

 pid=AIPID(); nprocs=AINPROCS();
 ni= CEIL(n,pi); nj= CEIL(n,pj); size=ni*nj;
 ppi= pid % pi; ppj= pid / pi;

 mat= (double *) malloc(size*sizeof(double));
#ifdef SANITY_CHECK
 check_if_null((void *)mat,"init_blocked","mat");
#endif

 for(k=0; k <size; k++) {
   nki = k % ni;
   nkj = k / ni;
   j  = nj * ppj + nkj;
   i  = ni * ppi + nki;
   if (1==option) 
/*
     mat[k] = (double) (j);
*/
     mat[k] = (double) ((j+1));
   else        
/*
     mat[k] = (double) (i);
*/
     mat[k] = (double) ((1+i+2*j));
 }
 *mat_r = mat;
}

void    print_mat(double *mat,int n, int debug,char *s) {
 register int i,j,N;
 int pid;

 pid=AIPID();
 if (n > 16) 
   N=16;
 else     
   N=n;

 if (0 == pid ) {
  if (1==debug){
   printf("%s \n",s);
   for(i=0;i<N;i++)
     for(j=0;j<N;j++)
/*
       printf("%10.3f%c",mat[i+j*n],((j<(N-1))?' ':'\n'));
*/
       printf("%5.2f%c",mat[i+j*n],((j<(N-1))?' ':'\n'));
  }
 }
}


void    collect_blocked_matrix(double *mat,double *mat1,int n,int pi,int pj) {

register int    pid;
register int    nprocs;
register int    i,j,k,ni,nj,ppi,ppj,nki,nkj,size;

  /* _blocked assumes a block distribution in which an nXn array is split into
   * p=pi*pj blocks each one of dimension (n/pi)X(n/pj)
   * It assumes pi/n and pj/n
   */
  pid=AIPID(); nprocs =AINPROCS();
  ni= CEIL(n,pi); nj= CEIL(n,pj); size=ni*nj;
  ppi= pid % pi; ppj= pid / pi;

  AIREGISTER(mat,n*n*sizeof(double));
  AIINIT();AICOMMIT();

  for(k=0; k <size; k+= CEIL(n,pi)) {
     nki = k % ni;
     nkj = k / ni;
     j  = nj * ppj + nkj;
     i  = ni * ppi;
     AIPUT(0,&mat1[k],&mat[0],(j*n+i)*sizeof(double),
             CEIL(n,pi)*sizeof(double));
  }
  AISYNC();  
  AIDEREGISTER(mat);
}


/* 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,2003  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,s,	 /* 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); 
  AIOINIT(a);
  AIOREGISTER(b,size); 
  AIOINIT(b);
  AIOCOMMIT();
for(l=0;l<p_sqrt;++l) {
  AIOPUT((((2*p_sqrt-ppi-ppj+l)%p_sqrt)*p_sqrt+ppi),A,a,0,size); /*1*/
  AIOSYNC(a);

  AIOPUT((((2*p_sqrt-ppj-ppi+l)%p_sqrt)+ppj*p_sqrt),B,b,0,size);
  AIOSYNC(b);
  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;
    }
*/
}

 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);
    AIOINIT(a);
    AIOREGISTER(b,  size);
    AIOINIT(b);
    AIOCOMMIT();
  for (l=0;l<p_sqrt;++l) {
    /* Efficient non-conflicting communication */
    /* A 2n^2/p relation is communicated */
    AIOPUT((((2*p_sqrt-ppj-ppi+l)%p_sqrt)*p_sqrt+ppi),A,a,0,size);/*2*/
    AIOSYNC(a);

    AIOPUT((((2*p_sqrt-ppj-ppi+l)%p_sqrt)+ppj*p_sqrt),B,b,0,size);
    AIOSYNC(b);
    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;
     }
    }
/*
    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,s,	 /* 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);
    AIOINIT(A);
    AIOREGISTER(B,size);
    AIOINIT(B);
    AIOCOMMIT();
  for (l=0;l<p_sqrt;++l) {
    AIOGET((((ppi+ppj+l)%p_sqrt)*p_sqrt+ppi),A,0,a,size);/*3*/
    AIOSYNC(A);
    AIOGET((((ppi+ppj+l)%p_sqrt)+ppj*p_sqrt),B,0,b,size);
    AIOSYNC(B);
    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;
    }
    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);
    AIOINIT(A);
    AIOREGISTER(B,size);
    AIOINIT(B);
    AIOCOMMIT();
  for (l=0;l<p_sqrt;++l) {
    /* Efficient non-conflicting communication */
    /* A 2n^2/p relation is communicated */
    AIOGET((((ppi+ppj+l)%p_sqrt)*p_sqrt+ppi),A,0,a,size);/*4*/
    AIOSYNC(A);
    AIOGET((((ppi+ppj+l)%p_sqrt)+ppj*p_sqrt),B,0,b,size);
    AIOSYNC(B);
    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;
     }
    }
/*
    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;
     }
    }
    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) ;
}
