#include	<stdio.h>
#include	<math.h>
#include	"bsp.h"
#include	"matrix.h"

void multiply_par();

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

	n=8;debug=0;nprocs=1;runs=1;pi=pj=1;
	if (1 == argc ) {
	  printf("Usage: %s dimension nprocs pi pj debug runs \n",argv[0]);
	  printf("Requires, currently, that pi=pj \n");
	  return(1);
	} else  
          if (2 == argc){
	    n=atoi(argv[1]);
	  } else  
            if (3 == argc){
	      n=atoi(argv[1]); nprocs=atoi(argv[2]);
	    } else  
              if (4 == argc){
	        n=atoi(argv[1]); nprocs=atoi(argv[2]); 
                pi=atoi(argv[3]);
	      } if (5 == argc){
		n=atoi(argv[1]); nprocs=atoi(argv[2]);
                pi=atoi(argv[3]); pj=atoi(argv[4]);
	        } 
		else {
		      n=atoi(argv[1]); nprocs=atoi(argv[2]);
		      pi=atoi(argv[3]); pj=atoi(argv[4]);
		      debug=atoi(argv[5]); runs=atoi(argv[6]);
		     }

/* Parallel execution begins here... 
   It makes the implicit assumption that all processors get the
   argument values. You may need to broadcast, otherwise, these values.
   For checking purposes nprocs-1 processor prints command line 
 */
   bsp_begin(nprocs);
	pid=bsp_pid(); 	nprocs=bsp_nprocs();

	if ((nprocs-1) == pid) {
	     for(i=0;i<argc;i++) 
               printf("%s ",argv[i]); 
           printf("\n");
         }

	t1=0.0; t2=0.0;
	for(k=0;k<runs;k++) {

	  size=CEIL(n,pi)*CEIL(n,pj)*sizeof(double);
	  mat = (double *) malloc(size);
	  check_if_null((void *) mat,"in main of mult.c","mat");
	  
	  if ((1)==debug) {
	     mat3 = (double *) malloc(n*n*sizeof(double));
	     check_if_null((void *)mat3,"in main of mult.c","mat3");
	  }

	  init_blocked(&mat1,n,pi,pj,1);
	  if ((1)== debug) {
	     collect_blocked_matrix(mat3,mat1,n,pi,pj);
             print_mat(mat3,n,1,"mat1");
	  }
	  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;
	  bsp_sync(); t3=bsp_time();
	  	multiply_par(mat1,mat2,mat,n,pi,pj);
    	  bsp_sync(); t4=bsp_time(); t1 += (t4-t3);

	  if ((1)== debug) {
	     collect_blocked_matrix(mat3,mat,n,pi,pj);
             print_mat(mat3,n,1,"mat=mat1Xmat2");
	  }

	  if (0 == pid) {
            printf("Matrix Mult.-1 (general): %8.5f (elapsed %8.5f)\n",t1,t1); 
            a= (2.0*n-1.0)*n*n;
            printf("Mflop is      : %8.5f\n",(a)/(nprocs*(t1)*1000000.0)); 
          }


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

   bsp_end();
}
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);
     bsp_abort("check_if_valid");
  }
}


void check_if_null(void *x,char *s,char *t)
{
  if (NULL == x) {
     printf("%s::Malloc error (%s)\n",s,t);
     bsp_abort("check_if_null");
  }
}

void    init_blocked(double **mat_r,int n,int pi,int pj,int option)
{
register int    pid=bsp_pid();
register int    nprocs=bsp_nprocs();
register int    i,j,k,ni,nj,ppi,ppj,nki,nkj,size;
double          *mat;

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

  mat= (double *) malloc(size*sizeof(double));
  check_if_null((void *)mat,"init_blocked","mat");

  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%10)+1);
           else        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 =bsp_pid();

        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=bsp_pid();
register int    nprocs=bsp_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
   */
  ni= CEIL(n,pi); nj= CEIL(n,pj); size=ni*nj;
  ppi= pid % pi; ppj= pid / pi;

  bsp_pushregister(mat,n*n*sizeof(double));
  bsp_sync();

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


