/*****************************************************************
**Matrix Multiplication Testing Functions(contains a main() )
**Copyright (c) 1997,2001  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
**/

#include <stdio.h>
#include <sys/utsname.h>
#include <unistd.h>
#include "matmul.h"
#include "aimpi.h"

#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));

void pmatmul();
void pomatmul();
void gmatmul();
void gomatmul();
void genmatmul();
void init_blocked();
void print_mat();
void collect_blocked_matrix();

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 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");
  }
}



/* 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();
 for (i = 0; i < nprocs;i++) {
   if (pid==fromp)
        AIHPPUT(i,from,to,0,size);
 }
 AISYNC();
}

/* 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 (*operator)(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();
 }
 /* each proc sends its from to every other proc; results stored in list */
 for (i = 0; i < nprocs;i++) {  
   AIHPPUT(i,from,list,pid*size,size);
 }
 AISYNC();
 /* apply operator to local list */
 for (i = 0; i < nprocs-1;i++) {
  operator(to,(list+i*size),(list+(i+1)*size));
  memcpy((list+(i+1)*size),to,size);
 }
    /* free memory and deregister */
 
    AIDEREGISTER(list);
    free((void *)list);
}


/* MAIN function */

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

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

 gethostname(hostname,sizeof(hostname));
 /* 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);
   }
 }
 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));
   AISYNC();
      bspbroad(0,(void *)&n ,(void *)&n ,sizeof(int));
   AISYNC();
   AIDEREGISTER(&n );

   AIREGISTER(&pi,sizeof(int));
   AISYNC();
      bspbroad(0,(void *)&pi ,(void *)&pi ,sizeof(int));
   AISYNC();
   AIDEREGISTER(&pi );

   AIREGISTER(&pj,sizeof(int));
   AISYNC();
      bspbroad(0,(void *)&pj ,(void *)&pj ,sizeof(int));
   AISYNC();
   AIDEREGISTER(&pj );

   AIREGISTER(&debug,sizeof(int));
   AISYNC();
      bspbroad(0,(void *)&debug,(void *)&debug  ,sizeof(int));
   AISYNC();
   AIDEREGISTER(&debug);

   AIREGISTER(&runs,sizeof(int));
   AISYNC();
      bspbroad(0,(void *)&runs,(void *)&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");
 }


  /*************/
  t1=0.0;
  AIBARRIER();
  t3=AITIME();  
   for(k=0;k<runs;k++) {
    pmatmul(mat1,mat2,mat,n,pi);
   }
  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);
   }
  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)");
  /*************/
  /*************/
  t1=0.0;
  AIBARRIER();
  t3=AITIME();  
   for(k=0;k<runs;k++) {
    gmatmul(mat1,mat2,mat,n,pi);
   }
  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);
   }
  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);
*/
   }
  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();
}

 




/* 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%10));
   else        
/*
     mat[k] = (double) (i);
*/
     mat[k] = (double) (i+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));
  AISYNC();

  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);
}


