#include "mpi.h"
#include "math.h"
#include "stdio.h"

#define  M 32
#define  N 128

int main( int argc, char *argv[])
{
    double starttime, endtime;
	double arr[M][N],p[8][8],q[8][8],c[8][8],a[8][8];
	double  ans1[8][8],DCT[M][N],a1[M][N],ans[M][N];
	int myid ,master,numprocs,len;
	int i,j,m,n,k,e,f,i1,j1,x,y;
	int numsent;
	char processor_name[MPI_MAX_PROCESSOR_NAME];
	MPI_Datatype anewtype;

	MPI_Init(&argc,&argv);
	MPI_Comm_rank(MPI_COMM_WORLD,&myid);
    MPI_Comm_size(MPI_COMM_WORLD,&numprocs);
	MPI_Get_processor_name(processor_name,&len);


	printf("processor %d  on%s\n",myid,processor_name);

	//由進(jìn)程0讀取image1.dat,xishu1.dat,xishu2.dat中的數(shù)據(jù)，并將其分別存入arr，q，p數(shù)組中。
	if (myid==0)
	{
		starttime=MPI_Wtime();
		FILE *fp; 
		FILE *fp1;
        FILE *fp2;
		if((fp=fopen("image1.dat","rb"))==NULL) //文件名是image1.dat 
			printf("Read error!!"); 
		
		for(i=0;i<M;i++) 
		{
			for(j=0;j<N;j++)
			{
				arr[i][j] = 0 ; //初始化數(shù)組		
			}
		}
		
		for(i=0;i<M;i++) 
		{
			for(j=0;j<N;j++)
			{
				fread(&arr[i][j],8,1,fp); //讀取數(shù)據(jù)		
			}
		}
		fclose(fp);
		for(i=0;i<M;i++) 
		{ 
			for(j=0;j<N;j++) 
				printf("arr %f \t",arr[i][j]); 
		}
		
	
		
		if((fp1=fopen("xishu2.dat","rb"))==NULL) //文件名是xishu2.dat 
			printf("Read error!!"); 
		for(i=0;i<8;i++) 
		{
			for(j=0;j<8;j++)
			{
				p[i][j]=0; //初始化數(shù)組		
			}
		}
		for(i=0;i<8;i++) 
		{
			for(j=0;j<8;j++)
			{
				fread(&p[i][j],8,1,fp1); //讀取數(shù)據(jù)		
			}
		}
		fclose(fp1);
		
		puts("The arrange is");    //打印該數(shù)組
		for(i=0;i<8;i++) 
		{ 
			for(j=0;j<8;j++) 
				printf("p %f \t",p[i][j]); 
		}
		
		if((fp2=fopen("xishu1.dat","rb"))==NULL) //文件名是xishu1.dat 
			printf("Read error!!"); 
		for(i=0;i<8;i++) 
		{
			for(j=0;j<8;j++)
			{
				q[i][j]=0; //初始化數(shù)組		
			}
		}
		for(i=0;i<8;i++) 
		{
			for(j=0;j<8;j++)
			{
				fread(&q[i][j],8,1,fp2); //讀取數(shù)據(jù)		
			}
		}
		fclose(fp2);

		puts("The arrange is");   //打印該數(shù)組
		for(i=0;i<8;i++) 
		{ 
			for(j=0;j<8;j++) 
				printf("q %f \t",q[i][j]); 
		}
		
	}





	
	

master=0;
numsent=M/numprocs;

if (numsent*numprocs<M) numsent=numsent+1;
// 定義a的行向量類型anewtype并提交
MPI_Type_vector(N,1,1,MPI_DOUBLE,&anewtype);
MPI_Type_commit(&anewtype);

//廣播矩陣p，q
MPI_Bcast(p,64,MPI_DOUBLE,master,MPI_COMM_WORLD);
MPI_Bcast(q,64,MPI_DOUBLE,master,MPI_COMM_WORLD);
//分發(fā)矩陣ayy
MPI_Scatter(arr,numsent,anewtype,a1,numsent,anewtype,master,MPI_COMM_WORLD);
//把矩陣的每個(gè)行塊劃分成8*8；對每塊進(jìn)行矩陣乘法
for(i=0;i<M/numprocs;i++) 
{ 
	for(j=0;j<N;j++) 
		printf("a1 %f \t",a1[i][j]); 
}
for(i=0;i<8;i++) 
{ 
	for(j=0;j<8;j++) 
		a[i][j]=0;
}

	


 for(i=0;i<numsent/8;i++)
	{
	 for(j=0;j<N/8;j++){
                    e=-1;
			for(x=(i*8);x<i*8+8;x++){
		        e=e+1; f=0;
				for(y=j*8;y<j*8+8;y++)
			
				{	a[e][f]=a1[x][y];
					f=f+1;
				}

			}
				for(m=0;m<8;m++) 
				{ 
					for(n=0;n<8;n++) 
						printf("a %f \t",a[m][n]); 
				}
				
	
				Multiply(8,p,a,c);
				for(m=0;m<8;m++) 
				{ 
					for(n=0;n<8;n++) 
						printf("c %f \t",c[m][n]); 
				}
				
	
				Multiply(8,c,q,ans1);

				for(m=0;m<8;m++) 
				{ 
					for(n=0;n<8;n++) 
						printf("ans1 %f \t",ans1[m][n]); 
				}
				

				for (i1=0;i1<8;i1++)
					for(j1=0;j1<8;j1++)

					{ 
						ans[i*8+i1][j*8+j1]=ans1[i1][j1];
					}
					
	 }		
					
	
}
 for(m=0;m<numsent;m++) 
 { 
	 for(j=0;n<N;n++) 
		 printf("ans %f \t",ans[m][n]); 
 }
//收集計(jì)算結(jié)果
  MPI_Gather(ans,numsent,anewtype,DCT,numsent,anewtype,master,MPI_COMM_WORLD);
  if(myid==0){
	  endtime=MPI_Wtime();
	   printf("time = %f\n",endtime-starttime);
	  printf("  The results are an fllows: \n");
	  for(i=0;i<M;i++){
		  for(j=0;j<N;j++)
			  printf("DCT %f \t",DCT[i][j]);
		  printf("\n");
	  }
  }
  MPI_Finalize();
}

Multiply(int n,double *a,double *b ,double *c)
{   
	int i,j,k;
	for(i=0;i<n;i++)
		for (j=0;j<n;j++)
			for(k=0;k<n;k++)
				c[i*n+j]+=a[i*n+k]*b[k*n+j];
}