test.cpp

本程序采用VC2005结合openmp实现并行矩阵相乘的Fox算法

#define p 16
#include "omp.h"
#include "conio.h"
#include <iostream>
#include <iomanip>
#include<time.h>
#include "stdlib.h"

class CMatrix
{
public:
	double m_data[4][4];
	void SetZero()
	{
		for(int i=0;i<4;i++)
			for(int j=0;j<4;j++)
				m_data[i][j]=0;
	}
};
CMatrix operator *(CMatrix &m1,CMatrix &m2)
{
	CMatrix m;
	int i, j, k;

	for (i=0; i<4; i++)
	{
		for (j=0; j<4; j++)
		{
			m.m_data[i][j] = 0;
			for (k=0; k<4; k++)
				m.m_data[i][j] += m1.m_data[i][k]*m2.m_data[k][j];
		}
	}
	return m;
}
CMatrix operator +(CMatrix &m1,CMatrix &m2)
{
	CMatrix m;
	int i, j;

	for (i=0; i<4; i++)
	{
		for (j=0; j<4; j++)
		{
			m.m_data[i][j] = m1.m_data[i][j]+m2.m_data[i][j];
		}
	}
	return m;
}

void UpMove(CMatrix A[4][4])
{
	int i,j;
	CMatrix m;

	#pragma omp parallel sections shared(A) private(i,j)
	{
		#pragma omp section
		{
			for(j=0;j<2;j++)
			for(i=0;i<3;i++)
			{
				m=A[i][j];
				A[i][j]=A[(i+1)%4][j];
				A[(i+1)%4][j]=m;
			}
		}

		#pragma omp section
		{
			for(j=2;j<4;j++)
			for(i=0;i<3;i++)
			{
				m=A[i][j];
				A[i][j]=A[(i+1)%4][j];
				A[(i+1)%4][j]=m;
			}
		}
	}
}
void main()
{
	clock_t begin,end,begin_normal,end_normal;
	begin=clock();
	int Aold[16][16],Bold[16][16],Result[16][16],text[16][16];
	int i,j,row;
	FILE *fpA = fopen( "A.txt", "w+" );
	FILE *fpB = fopen( "B.txt", "w+" );
	for(i=0;i<16;i++)   //对Aold,Bold进行初始化
	{
		for(j=0;j<16;j++)
		{
			if(i==j)
			{
			Aold[i][j]=2.0;
			}
			else Aold[i][j]=1.0;
			Bold[i][j]=(double)(i+j+1);
			fprintf(fpA,"\nA[%d][%d}=%d",i,j,Aold[i][j]);
			fprintf(fpB,"\nB[%d][%d]=%d",i,j,Bold[i][j]);

		}
	}
	fclose(fpA);
	fclose(fpB);

	FILE *fpD = fopen( "D.txt", "w+" );	

	begin_normal=clock();

	for(i = 0;i<16;i++)
		for(j = 0;j<16;j++)
		{
			text[i][j] = 0;
			for(int m = 0;m<16;m++)
				text[i][j]+=Aold[i][m]*Bold[m][j];
		}

    for( i=0; i<16; i++ )
		for( j=0; j<16; j++ )
		{
			//printf( "\n C[%d][%d] = %d ", i, j, C[i][j] );
			fprintf( fpD, "\n D[%d][%d] = %d ", i, j, text[i][j] );
			printf("\nC[%d][%d] = %d",i,j,text[i][j]);
		}
		fclose(fpD);

	end_normal=clock();
	printf("using time %f.",(double)end_normal-(double)begin_normal);


	CMatrix  A[4][4];
	CMatrix  B[4][4];
	CMatrix  C[4][4];
	CMatrix  W[4];
	for(i=0;i<16;i++)
		for(j=0;j<16;j++)
		{
			A[i/4][j/4].m_data[i%4][j%4]=Aold[i][j];
		}

	for(i=0;i<16;i++)
		for(j=0;j<16;j++)
		{
			B[i/4][j/4].m_data[i%4][j%4]=Bold[i][j];
		}
	for(i=0;i<4;i++)
		for(j=0;j<4;j++)
		{
			C[i][j].SetZero();
		}

	for(int k=0;k<4;k++)
	{
		for(row=0;row<4;row++)
			W[row]=A[row][(row+k)%4];  //取行中广播的元素
	#pragma omp parallel sections shared(A,C,W,B) private(i,j)
	{
		#pragma omp section
		{
			for(i=0;i<2;i++)
			for(j=0;j<4;j++)
			{
				C[i][j]=C[i][j]+W[i]*B[i][j];
			}
		}

		#pragma omp section
		{
			for(i=2;i<4;i++)
			for(j=0;j<4;j++)
			{
				C[i][j]=C[i][j]+W[i]*B[i][j];
			}
		}
	}
		UpMove(B);  //B向上循环移位*
	}
//将块阵重新组合为Result矩阵	
	for(i=0;i<4;i++)
		for(j=0;j<4;j++)
			for(int k=0;k<4;k++)
				for(int q=0;q<4;q++)
		{
			Result[i*4+k][j*4+q]=C[i][j].m_data[k][q];
		}
	
	end=clock();
	printf("\n Time = %f ", ((double)end-(double)begin));

			// 显示本次迭代结果
		FILE *fpC = fopen( "C.txt", "w+" );
		for( i=0; i<16; i++ )
		for( j=0; j<16; j++ )
		{
			//printf( "\n C[%d][%d] = %d ", i, j, C[i][j] );
			fprintf( fpC, "\n C[%d][%d] = %d ", i, j, Result[i][j] );
			printf("\nC[%d][%d] = %d",i,j,Result[i][j]);
		}
		fclose(fpC);

		
		getch();
}