cpi.c

mpi 结合vc编程用模拟退火法算一个最小路径的值

#include<C:\Program Files\MPICH\SDK\Include\mpi.h>
#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#include<math.h>
#define  maxn 200
#define  maxnp 500
#define  maxlp  500

//产生[0,m)随机数的函数
int myrandom(int m)
{
	time_t s;
	return abs(((int)time(&s)*rand())%m);
}
//产生[0,1)随机数的函数
float myrandom0_1()
{
	return ((float)myrandom(32767)/32767.f);
}
float dist(int x1,int y1,int x2,int y2)
{
	return (float)sqrt((float)(x1-x2)*(float)(x1-x2)+(float)(y1-y2)*(float)(y1-y2));
}
void init_d(int *x,int *y,float d[maxn][maxn],int *path,int n)
{
	int i,j;
	for(i=1;i<=n;i++)
	{
		path[i]=i;
		for(j=1;j<=i;j++)d[i][j]=dist(x[i],y[i],x[j],y[j]);
	}
	for(i=1;i<=n;i++)for(j=i;j<=n;j++)d[i][j]=d[j][i];
}
void generate(int n,int *r,int *m,int *c)
{
	int i;
	do
	{
		c[1]=1+myrandom(n);
		c[2]=1+myrandom(n-1);
		if(c[2]==c[1])c[2]++;
		i=1+(c[1]-c[2]+n-1)%n;
	}while(i<3);
	c[3]=1+(c[1]+n-2)%n;
	c[4]=1+c[2]%n;
	c[5]=0;
	c[6]=0;
	*r=2+myrandom(2);
	*m=4;
	if(*r==3)
	{
		c[5]=1+(c[2]+myrandom(abs(i-2)))%n;
		c[6]=1+c[5]%n;
		*m=6;
	}
}
float delta_length(int r,int m,int *c0,int *p,float d[maxn][maxn])
{
	float df;
	int j,c[7];
	for(j=1;j<=m;j++)c[j]=p[c0[j]];
	if(r==2)
		df=d[c[1]][c[4]]+d[c[2]][c[3]]-d[c[1]][c[3]]-d[c[2]][c[4]];
	else 
		df=d[c[1]][c[5]]+d[c[2]][c[6]]+d[c[3]][c[4]]-d[c[1]][c[3]]-d[c[2]][c[4]]-d[c[5]][c[6]];
	return(df);
}
int accept(float t,float df)
{
	return(df<0.0)||((df/t<88.)&&(exp(-df/t)>myrandom0_1()));

}
void twochain(int n,int *c,int *p)
{
	int i,j,u,v,w;
	i=(1+(c[2]-c[1]+n)%n)/2;
	for(j=1;j<=i;j++)
	{
		u=1+(c[1]+j-2)%n;
		v=1+(c[2]-j+n)%n;
		w=p[u];
		p[u]=p[v];
		p[v]=w;

	}
}
void threechain(int n,int *c,int *p)
{
	int p0[maxn],i,j,m1,m2,m3,w;
	m1=1+(c[2]-c[1]+n)%n;m2=1+(c[3]-c[6]+n)%n;m3=1+(c[5]-c[4]+n)%n;
	i=1;
	for(j=1;j<=m1;j++){w=1+(j+c[1]-2)%n; p0[i]=p[w]; i++;}
	for(j=1;j<=m2;j++){w=1+(j+c[6]-2)%n; p0[i]=p[w]; i++;}
	for(j=1;j<=m3;j++){w=1+(j+c[4]-2)%n; p0[i]=p[w]; i++;}
	for(j=1;j<=n;j++)p[j]=p0[j];
}
//模拟退火法求解TSP的主程序
float annealing(int n,int lp,int s1,float t,float dt,int *p,float d[maxn][maxn])
{
	float df,f=0;
	int c[7],k,r,m,a,s0=0;
	for(k=1;k<n;k++)f=f+d[p[k]][p[k+1]];
	f=f+d[p[n]][p[1]];
	do
	{
		a=0;
		for(k=1;k<=lp;k++)
		{ 
			generate(n,&r,&m,c); 
			df=delta_length(r,m,c,p,d);
			if(accept(t,df)==1)
			{
				if(r==2)twochain(n,c,p);
				else
					threechain(n,c,p);
				f=f+df;
				a=1;
			}
		}
		t=t*dt;
		if(a==0)s0++;
		else
			s0=0;
	}while(s0<s1);
	return(f);
}
//模拟退火法求解TSP144主程序

main(int argc,char **argv)

{
	
	int r,np,myid,root=0,no[2],loop,Lp,i,j,k,n,lp,s1=0;
	time_t s;
	int x[maxn],y[maxn];
	int path[maxn],best_path[maxn],npath[maxnp][maxn],keep[maxlp][maxn];
	float d[maxn][maxn],ngood[maxnp],nlit,nbig,best,maybe_min,t,dt;
	long   start_time,end_time;
	MPI_Comm comm;
//	char *my_pointer;


	for(r=0;r<maxn;r++)
	{
		x[r]=myrandom(1000);
		y[r]=myrandom(1000);
	
	}
//	FILE *f;

	MPI_Init(&argc,&argv);
	MPI_Comm_size(MPI_COMM_WORLD,&np);
	MPI_Comm_rank(MPI_COMM_WORLD,&myid);
	comm=MPI_COMM_WORLD;
	
	/*if(np<2)
	{ 
		printf("\n np must be set>1");
		goto abandon;
	}*/
	if(myid==root)
	{                             //根进程实施TSP实例的初始化
		n=199;
		init_d(x,y,d,path,n);//根进程实施TSP实例距离矩阵的生成
		printf("\ninput number of loop =(<=500)");
		//scanf("%d",&Lp);
		Lp=115;
		start_time=time(&s);//TSP实例初始化完成
	}
	MPI_Bcast(&n,1,MPI_INT,root,comm);//TSP实例数据发布
	MPI_Bcast(&Lp,1,MPI_INT,root,comm);
	MPI_Bcast(d,maxn*maxn,MPI_FLOAT,root,comm);
	MPI_Bcast(path,maxn,MPI_INT,root,comm);
	for(loop=0;loop<Lp;loop++)
	{                        //并行实施模拟退火法求解TSP的循环

		lp=1000;
		t=100.0f;
		dt=0.95f;
		s1=1;
		maybe_min=annealing(n,lp,s1,t,dt,path,d); //各进程结果因初始路径与随机数值而异
		MPI_Gather(&maybe_min,1,MPI_FLOAT,ngood,1,MPI_FLOAT,root,comm);
		MPI_Gather(path,maxn,MPI_INT,npath,maxn,MPI_INT,root,comm);
		MPI_Bcast(npath,np*maxn,MPI_INT,root,comm);
		if(myid==root)//根进程分析本轮迭代各进程的优化结果，找出其中最短/最长的路径
		{
			if(loop==0)best=maybe_min;
			nlit=maybe_min;
			nbig=nlit;
			i=0;
			j=0;
			for(k=1;k<np;k++)
			{
				if(ngood[k]<nlit){nlit=ngood[k];i=k;}
				if(ngood[k]>nbig){nbig=ngood[k];j=k;}
			}
			no[0]=i;
			no[1]=j;
			keep[loop][0]=(int)nlit;//一般不能用整型量空闲单元存储nlit(双精度)
			for(k=1;k<=n;k++)keep[loop][k]=path[k];
			if(best>nlit)
			{
				best=nlit;
				for(k=0;k<=n;k++)best_path[k]=keep[loop][k];
			}
			printf("\n lp=%4d;best:%9.2f  i=%2d:%9.2f(min)j=%2d:%9.2f(max)",loop,best,i,ngood[i],j,ngood[j]);
		}
		MPI_Bcast(no,2,MPI_INT,root,comm);//公布根进程的分析结果
		if(myid<=np/3)k=no[0];
		else
			k=no[1];//结点采用不同策略选择下一轮的初始路径
		for(i=0;i<maxn;i++)path[i]=npath[k][i];
		MPI_Barrier(comm);
	}//并行实施模拟退火法求TSP循环结束
	//跟进程输出结果
	if(myid==root)
	{
		end_time=time(&s);
		printf("\nTSP-%d,time=%ld  best=%10.2f\n",n,end_time-start_time,best);
	//	f=fopen("TSP_144.dat","w");
		//结果存盘
	//	fclose(f);
	}
    MPI_Finalize();/*abandon:*/



}