#include<mpi.h>
#include<stdlib.h>
#define MAXNPES 16
#define MATRIXSIZE 8
#define NUMENTRY 3
typedef struct{
    int nsnbrs,*spes;
    int nrnbrs,*rpes;
    int *sendptr,*recvptr,*sendind,*recvind;
    double *sendbuf,*recvbuf;
  }CommInfoType;
main(int argc, char *argv[])
{
  int myrank,npes,i,j,r;
  int nlocal,*rowptr,*colind;
  double *values,*b,*y,*gb,*y2;
  double Dotproduct(int nlocal, double *x1, double *x2);
  CommInfoType cinfo;
  MPI_Init(&argc, &argv);
  MPI_Comm_size(MPI_COMM_WORLD, &npes);
  MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
  nlocal=MATRIXSIZE/npes;
  y=(double*)malloc((nlocal)*sizeof(double));
  y2=(double*)malloc((nlocal)*sizeof(double));
  rowptr=(int*)malloc((nlocal+1)*sizeof(int));
  b=(double*)malloc(nlocal*sizeof(double));
  srand(myrank);
  srand48((myrank)*1.0);
  rowptr[0]=0;
  for(i=0;i<nlocal;i++){
    r=rand()%(NUMENTRY)+1;
    rowptr[i+1]=rowptr[i]+r;
  }
  colind=(int*)malloc(rowptr[nlocal]*sizeof(int));
  values=(double*)malloc(rowptr[nlocal]*sizeof(double));
  
  
  SpareMatrixGenerator(MATRIXSIZE,NUMENTRY,rowptr,colind,values,b,
		       MPI_COMM_WORLD);

  MPI_Barrier(MPI_COMM_WORLD);

  
  SparseMatVec(nlocal,rowptr,colind, values,b,y, MPI_COMM_WORLD);

  for(i=0;i<nlocal;i++)printf("myrank %d,%f\n",myrank,y[i]);
  MPI_Barrier(MPI_COMM_WORLD);

  gb=(double*)malloc(nlocal*npes*sizeof(double));

  MPI_Allgather(b,nlocal,MPI_DOUBLE,gb,nlocal,MPI_DOUBLE,MPI_COMM_WORLD);

  SerialSpareMatVec(nlocal,rowptr,colind,values,gb,y2);

  for(i=0;i<nlocal;i++)printf("myrank %d,%f\n",myrank,y2[i]);

  printf("dot product of y and y2 %f\n",Dotproduct(nlocal, y, y2));

  free(rowptr);
  free(colind);
  free(values);
  free(y);
  free(b);
  free(gb);
  free(y2);
  MPI_Finalize();
}

double Dotproduct(int nlocal, double *x1, double *x2)
{
  int i;
  double localdot,globaldot;
  localdot=0.0;
  for(i=0;i<nlocal;i++)
    localdot+=x1[i]*x2[i];
  MPI_Barrier(MPI_COMM_WORLD);
  MPI_Allreduce((void *)&localdot,(void *)&globaldot,1,MPI_DOUBLE,
		MPI_SUM,MPI_COMM_WORLD);
  printf("iiiii %f\n",globaldot);
  return(globaldot);
}

SpareMatrixGenerator(int matrixsize,int numentry,
		     int *rowptr,int *colind,double *values,double *b,
		     MPI_Comm comm)
{
 
  int i,j,k,nlocal,r,d,numproc,mype;
  int IncOrder(const void *e1,const void *e2);
  MPI_Comm_rank(comm,&mype);
  MPI_Comm_size(comm,&numproc);
  nlocal=matrixsize/numproc;
  /*rowptr=(int*)malloc((nlocal+1)*sizeof(int));
  b=(double*)malloc(nlocal*sizeof(double));
  srand(mype);
  srand48((mype)*1.0);
  rowptr[0]=0;
  for(i=0;i<nlocal;i++){
    r=rand()%(numentry)+1;
    rowptr[i+1]=rowptr[i]+r;
  }
  colind=(int*)malloc(rowptr[nlocal]*sizeof(int));
  values=(double*)malloc(rowptr[nlocal]*sizeof(double));*/

  /* generate entries */

  for(i=0;i<rowptr[nlocal];i++ )
	values[i]=drand48();

  for(j=0;j<nlocal;j++){
    for(i=rowptr[j];i<rowptr[j+1];i++){
    regenerate:
    d=rand()%matrixsize;
    for(k=rowptr[j];k<i;++k)
      if(d==colind[k]) goto regenerate;
    colind[i]=d;
    }
    qsort(colind+rowptr[j],rowptr[j+1]-rowptr[j],sizeof(int),IncOrder);
  }
  /* generate b */
  
  for(i=0;i<nlocal;i++) b[i]=drand48();

  /*printf("rowptr %d ",mype);
  for(i=0;i<nlocal+1;i++)printf("%d ",rowptr[i]);
  for(j=0;j<nlocal;j++){
    printf("\n colind %d  ",mype);
    for(i=rowptr[j];i<rowptr[j+1];i++)printf("%d ",colind[i]);
  }
  printf("\n values %d ",mype);
  for(i=0;i<rowptr[nlocal];i++)printf("%f ",values[i]);
  printf("\n b %d ",mype);
  for(i=0;i<nlocal;i++)printf("%f ",b[i]);
  printf("\n");*/
}







SparseMatVec(int nlocal, int *rowptr, int *colind, double *values,
	    double *b,double *y, MPI_Comm comm)
{
  int npes;
  CommInfoType cinfo;
  double *gb;
  MPI_Comm_size(comm,&npes);
  gb=(double *)malloc(nlocal*npes*sizeof(double));

  SetUpCommInfo(nlocal,rowptr,colind,&cinfo, comm);

  CommInterfaceValues(nlocal,b,&cinfo,gb,comm);
  
  SerialSpareMatVec(nlocal,rowptr,colind,values,gb,y);

  
  free(gb);
}


SetUpCommInfo(int nlocal,int *rowptr, int *colind, CommInfoType *cinfo,
	      MPI_Comm comm)
{ 
  int i,j,k;
  int npes,mype,*pes,*ptr;
  int nsend,nrecv,*rentries;
  int *receives,*sends;
  int *gvec;
  int IncOrder(const void *e1,const void *e2);
  MPI_Request sreq[MAXNPES],rreq[MAXNPES];
  MPI_Status statuses[MAXNPES];
  MPI_Comm_size(comm,&npes);
  MPI_Comm_rank(comm,&mype);

  pes=(int *)malloc(npes*sizeof(int));
  ptr=(int *)malloc((npes+1)*sizeof(int));
  gvec=(int *)malloc(nlocal*(npes)*sizeof(int));

  rentries=(int *)malloc((rowptr[nlocal]+1)*sizeof(int));
  for(i=0;i<npes*nlocal;i++)
    gvec[i]=0;
  for(i=mype*nlocal;i<(mype+1)*nlocal;i++)
    gvec[i]=1;
  for(nrecv=i=0;i<nlocal;i++){
    for(j=rowptr[i];j<rowptr[i+1];j++){
      k=colind[j];
      if(gvec[k]==0){
	gvec[k]=1;
	rentries[nrecv++]=k;
      }
    }
  }
  rentries[nrecv]=npes*nlocal;

  qsort(rentries,nrecv,sizeof(int),IncOrder);

  cinfo->nrnbrs=ptr[0]=0;
  for(j=i=0;i<npes;i++){
    for(k=j;rentries[j]<(i+1)*nlocal;j++);
    if(j>k){
      pes[cinfo->nrnbrs++]=i;
      ptr[cinfo->nrnbrs]=j;
    }
  }

  cinfo->rpes    = (int *)malloc(cinfo->nrnbrs*sizeof(int));
  cinfo->recvptr = (int *)malloc((cinfo->nrnbrs+1)*sizeof(int));
  cinfo->recvind = (int *)malloc(nrecv*sizeof(int));
  cinfo->recvbuf = (double  *)malloc(nrecv*sizeof(double));
  
  for(i=0;i<cinfo->nrnbrs;i++) cinfo->rpes[i]=pes[i];
  for(i=0;i<cinfo->nrnbrs+1;i++) cinfo->recvptr[i]=ptr[i];
  for(i=0;i<nrecv;i++) cinfo->recvind[i]=rentries[i];

  receives=(int *)malloc(npes*sizeof(int));
  sends=(int *)malloc(npes*sizeof(int));

  for(i=0;i<npes;i++) receives[i]=0;
  for(i=0;i<cinfo->nrnbrs;i++)
    receives[cinfo->rpes[i]]=cinfo->recvptr[i+1]-cinfo->recvptr[i];

  MPI_Alltoall((void *)receives,1,MPI_INT,(void *)sends,1,MPI_INT,comm);

  cinfo->nsnbrs=ptr[0]=0;
  for(i=0;i<npes;i++){
    if(sends[i]>0){
      pes[cinfo->nsnbrs++]=i;
      ptr[cinfo->nsnbrs]=ptr[cinfo->nsnbrs-1]+sends[i];
    }
  }
  nsend=ptr[cinfo->nsnbrs];

  cinfo->spes=(int *)malloc(cinfo->nsnbrs*sizeof(int));
  cinfo->sendptr=(int *)malloc((cinfo->nsnbrs+1)*sizeof(int));
  cinfo->sendind=(int *)malloc(nsend*sizeof(int));
  cinfo->sendbuf=(double  *)malloc(nsend*sizeof(double));

  for(i=0;i<cinfo->nsnbrs;i++) cinfo->spes[i]=pes[i];
  for(i=0;i<cinfo->nsnbrs+1;i++) cinfo->sendptr[i]=ptr[i];

  for(i=0;i<cinfo->nsnbrs;i++)
    MPI_Irecv((void *)(cinfo->sendind+cinfo->sendptr[i]),
	      cinfo->sendptr[i+1]-cinfo->sendptr[i],MPI_INT,
	      cinfo->spes[i],1,comm,rreq+i);

  for(i=0;i<cinfo->nrnbrs;i++)
    MPI_Isend((void *)(cinfo->recvind+cinfo->recvptr[i]),
	      cinfo->recvptr[i+1]-cinfo->recvptr[i],MPI_INT,
	      cinfo->rpes[i],1,comm,sreq+i);

  MPI_Waitall(cinfo->nsnbrs,rreq,statuses);
  MPI_Waitall(cinfo->nrnbrs,sreq,statuses);

  free(rentries);
  free(pes);
  free(ptr);
  free(receives);
  free(sends);


  /* if(mype==2){
    printf("my rank %d\n",mype);
    printf("nsnbrs %d",cinfo->nsnbrs);
    printf("\n spes ");
    for(i=0;i<cinfo->nsnbrs;i++)printf("%d",cinfo->spes[i]);
    printf("\n sendptr ");
    for(i=0;i<cinfo->nsnbrs+1;i++)printf("%d",cinfo->sendptr[i]);
    printf("\n sendind ");
    for(i=0;i<nsend;i++)printf("%d",cinfo->sendind[i]);
  }*/
  /*for(i=0;i<npes;i++) printf("%d,  %d ",sends[i], mype);*/
  
}

int IncOrder(const void *e1,const void *e2)
{
  return(*((int *)e1)-*((int *)e2));
}

CommInterfaceValues(int nlocal,double *b,CommInfoType *cinfo,
		    double *gb, MPI_Comm comm)
{
  int i,j,mype,firstrow;
  MPI_Request sreq[MAXNPES],rreq[MAXNPES];
  MPI_Status statuses[MAXNPES];
  

  MPI_Comm_rank(comm,&mype);
  for(i=0;i<cinfo->nrnbrs;i++)
    MPI_Irecv((void *)(cinfo->recvbuf+cinfo->recvptr[i]),
	      cinfo->recvptr[i+1]-cinfo->recvptr[i],MPI_DOUBLE,
	      cinfo->rpes[i],1,comm,rreq+i);

  firstrow=mype*nlocal;
  for(i=0;i<cinfo->sendptr[cinfo->nsnbrs];i++)
	cinfo->sendbuf[i]=b[cinfo->sendind[i]-firstrow];

  for(i=0;i<cinfo->nsnbrs;i++)
    MPI_Isend((void *)(cinfo->sendbuf+cinfo->sendptr[i]),
	      cinfo->sendptr[i+1]-cinfo->sendptr[i],MPI_DOUBLE,
	      cinfo->spes[i],1,comm,sreq+i);

  MPI_Waitall(cinfo->nrnbrs,rreq,statuses);
  MPI_Waitall(cinfo->nsnbrs,sreq,statuses); 

  for(i=0;i<cinfo->recvptr[cinfo->nrnbrs];i++)
    gb[cinfo->recvind[i]]=cinfo->recvbuf[i];
  for(i=0;i<nlocal;i++)
    gb[firstrow+i]=b[i];
}
SerialSpareMatVec(int n, int *rowptr, int *colind,
		       double *values, double *b, double *y)
{
  int i,j;
  for(i=0;i<n;i++){
    y[i]=0.0;
    for(j=rowptr[i];j<rowptr[i+1];j++)
      y[i]+=values[j]*b[colind[j]];
  }
}