biology.cpp

生物领域的一个并行计算代码

#include "mpi.h"
#include <fstream>
#include<stdio.h>
#include <cmath>
#include<fstream>
#include<strstream>
#include <sstream>
#include <time.h>
using namespace std;

int deg(int **a,int i, int num);
double max_value(double x,double y);


void parseString(char* pString, int *a, int num)
{
	std::istrstream inp(pString,0);
//	int aNumber;
//	int balance;
	for(int i = 0; i < num ; i++)
	{
		inp>>a[i];
	}
//	char* pBuffer = new char[128]; 此处的代码已经没有使用,可以屏蔽
//	ostrstream outp(pBuffer,128);
//	outp<<"aNumber="<<aNumber
//		<<"balance="<<balance<<endl;
//	return pBuffer;
}

void parseString1(char* pString1, double *b, int num1)
{
	std::istrstream inp(pString1,0);
	for(int i = 0; i < num1 ; i++)
	{
		inp>>b[i];
	}
//	char* pBuffer = new char[128]; 此处的代码已经没有使用,可以屏蔽
//	ostrstream outp(pBuffer,128);
//	return pBuffer;
}

int deg(int **a,int i, int num)
{
	int deg = 0;
	for(int j = 0; j < num; j++)
	{
		deg = deg + a[i][j];
	}
	return deg;
}

double max_value(double x,double y)
{
	if(x>y) return x;
	else return y;
}

const int BUFFERSIZE = 1024;
	
int main(int argc,char *argv[])
{
	int i,j,m,n,u,v;
	double max,localMax;
	double starttime, endtime;
	int myid,numprocs;
	std::ifstream inFile;
	std::ofstream outFile;
	char *str;
	str = new char[BUFFERSIZE];
	inFile.open("testdata.txt",ios::in);
	if(!inFile)
	{
		printf("File1.txt can't open.\n");
		abort();
	}
	inFile.getline(str, BUFFERSIZE);
	inFile.getline(str, BUFFERSIZE);
	parseString(str, &m, 1);
	inFile.getline(str, BUFFERSIZE);
	inFile.getline(str, BUFFERSIZE);
	parseString(str, &n, 1);
	printf("%d %d\n",m,n);
	inFile.getline(str, BUFFERSIZE);

	int **A;     //分配动态数组A；
	A = new int * [m];
	A[0] = new int [m*m];
	for (i=0;i<m;i++)
	{
		A[i] = A[0] + m*i;
	}
	
	for(j = 0; j < m; j ++)			//读取A数组信息；
	{
		inFile.getline(str,BUFFERSIZE);
		parseString(str, A[j], m);
			
	}
	
	//cout<<A[0][6];
	int **AA;     //分配动态数组AA；
	AA = new int * [n];
	AA[0] = new int [n*n];
	for (i=0;i<n;i++)
	{
		AA[i] = AA[0] + n*i;
	}
	inFile.getline(str, BUFFERSIZE); //此处应过滤一行
	for(i = 0; i < n; i ++)			//读取AA数组信息；
	{
		inFile.getline(str,BUFFERSIZE);
		parseString(str, AA[i], n);
			
	}
	//cout<<AA[4][1];


	double **s;     //分配动态数组s；
	s = new double * [m];
	s[0] = new double [m*n];
	for (i=0;i<m;i++) //此处应为m 因为是m行
	{
		s[i] = s[0] + n*i;
	}
	inFile.getline(str, BUFFERSIZE); //此处应过滤一行
	for(i = 0; i < m; i ++)	//此处应为m 因为是m行		//读取s数组信息；
	{
		inFile.getline(str,BUFFERSIZE);
		parseString1(str, s[i], n);
			
	}
	inFile.close();
	
	
	double **sk0;     
	sk0 = new double * [m];
	sk0[0] = new double [m*n];
	for (i=0;i<m;i++) //此处应为m 因为是m行
	{
		sk0[i] = sk0[0] + n*i;
	}
	double **sk1;     
	sk1 = new double * [m];
	sk1[0] = new double [m*n]; //此处应为sk1 不是sk0
	for (i=0;i<m;i++) //此处应为m 因为是m行
	{
		sk1[i] = sk1[0] + n*i;
	}
	
	for(i=0;i<m;i++)
	{
		for(j=0;j<n;j++)
		{
			sk0[i][j] = s[i][j];
			sk1[i][j] = 0;
		}
		printf("\n");
	}
double deltaMax;

printf("%d ",AA[7][0]);

	MPI_Init(&argc,&argv);				//启动MPI环境
	MPI_Comm_size(MPI_COMM_WORLD,&numprocs);	//确定进程数
	MPI_Comm_rank(MPI_COMM_WORLD,&myid);
	MPI_Status status;

	int aveNum = m/numprocs;
	int modNum = m%numprocs;
	int myNum, myN, beginRow,bRow;
	if(myid < modNum) 
	{
		myNum = aveNum + 1;
		beginRow = myNum * myid;
	}
	else 
	{
		myNum = aveNum;
		beginRow = modNum * (aveNum + 1) + myNum * (myid - modNum);
	}
	int tag = 0;
	
	starttime = MPI_Wtime();
do
{
	tag++;
	localMax = 0;
	max = 0;
	
	MPI_Bcast(sk0[beginRow], m*n, MPI_DOUBLE, 0, MPI_COMM_WORLD);  //广播出去
	printf("one");
	
	for(i=beginRow;i < myNum;i++)
		for(j=0;j<n;j++)
		{
			double N1=0;
			double N2=0;
			for(u=0;u<m;u++)
			if(u!=i)
			{
				for(v=0;v<n;v++)
				{
					double n1=0;
					double n2=0;
					
					if(v!=j)
					{
						if(A[i][u]==1)
						{	
							if(AA[j][v]==1)
							{
								n1=sk0[u][v];
							}
						}
						else
						{
							if(AA[j][v]==0)
							{
								n2=sk0[u][v];
							}
						}
						if((deg(A,i,m)==0)&&(deg(AA,j,n)==0))
						{
							N1=N1 + n1/(m*n);
						}
						else if((deg(A,i,m)!=0)&&(deg(AA,j,n)!=0))
						{
							N1=N1 + n1/(deg(A,i,m)*deg(AA,j,n));
						}
						if((deg(A,i,m)!=m)&&(deg(AA,j,n)!=n))
						{
							N2=N2 + n2/((m-deg(A,i,m))*(n-deg(AA,j,n)));
						}
						else if(((m-deg(A,i,m))==m)&&(n-deg(AA,j,n)==n))
						{
							N2=N2 + n2/(m*n);
						}
					}//if()
				}//for v
			}//for u
			
			sk1[i][j]=(N1+N2)*s[i][j]/2;
			localMax = max_value(sk1[i][j], localMax); //求最大值放到此处，可以少一个二重循环
		}// for i,j
	
	MPI_Reduce(&localMax, &max, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD); 
	printf("two");
	if(myid == 0)
	{
		for(int p = 1; p < numprocs; p ++)
		{
			if(p < modNum) 
			{
				myN = aveNum + 1;  //myNum 不要和上面的重复,此处为局部使用
				bRow = myNum * p;
			}
			else 
			{
				myN = aveNum;
				bRow = modNum * (aveNum + 1) + myNum * (p - modNum);
			}

			MPI_Recv(sk1[bRow], n*myN, MPI_DOUBLE, p, tag, MPI_COMM_WORLD, &status);
			printf("three");
		}
	}
	else
	{
		MPI_Send(sk1[beginRow], n*myNum, MPI_DOUBLE,0, tag, MPI_COMM_WORLD);
		printf("four");
	}

	

	deltaMax = 0;
	if(myid == 0)
	{
	for(int k=0;k<m;k++)
		for(int z=0;z<n;z++)
		{
			sk1[k][z]=sk1[k][z]/max;
			if(deltaMax < fabs(fabs(sk1[k][z]) - fabs(sk0[k][z]))) deltaMax = fabs(fabs(sk1[k][z]) - fabs(sk0[k][z]));
			//sk0[k][z] = sk1[k][z];
		}
	//这里是指针复制，不是数组的数据复制
	double **temp;
	temp = sk0;
	sk0 = sk1;
	sk1 = temp;
	}

	MPI_Bcast(&deltaMax, 1,MPI_DOUBLE,0, MPI_COMM_WORLD);
	printf("five");

}while(deltaMax>0.01);

endtime = MPI_Wtime();
	printf("%f\n", endtime-starttime);

if(myid==0)
{
	outFile.open("file1.txt",ios::out);
	if(!outFile)
		{
		printf("File2.txt can't open.\n");
		abort();
		}
	for(i=0;i<m;i++)
	{
		for(j=0;j<n;j++)
			{		
			
				outFile<<sk0[i][j]<<" ";  //输出数据不是指针
			}
		outFile<<endl;
	}
	outFile.close();

	
	delete [] A[0];
	delete [] A;

	delete [] AA[0];
	delete [] AA;


	delete [] sk0[0];
	delete [] sk0;

	delete [] sk1[0];
	delete [] sk1;

	delete []s[0];
	delete []s;
	//这里对A AA sk0 sk1 s进行空间的释放。
}
	MPI_Finalize();		//结束MPI

}