final_ok_comment.cpp

使用“Biconjugate Gradient Method”并行程序实现复杂的线性方程求解

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <mpi.h>
#include <conio.h>

struct MatElement{
	int col, row;
	double value;
};

struct MatComplex{
	int col, row;
	double real, imag;
};

void main(int argc, char** argv)
{
	int nMyRank;
	int nCommSize;
	const double CONV_MIN = 1.0e-3;	//The convergence value
	const double MIN = 1.0e-9;
	bool failed=false;

	MPI_Init(&argc, &argv);
	MPI_Comm_rank(MPI_COMM_WORLD, &nMyRank);
	MPI_Comm_size(MPI_COMM_WORLD, &nCommSize);

	MatComplex *mtxA;
	MatElement *mtxM;
	double *b;
	double *r, *rbar;
	double *p, *pbar;
	double *oldp, *oldpbar;
	double *q, *qbar;
	double *oldq, *oldqbar;
	double *z, *zbar;
	double *x;
	int rowCount, colCount;
	int oldrowCount, oldcolCount;
	int nnz;
	int row_index;
	double alpha;
	FILE *fa, *fb, *finfo;
	FILE *fout;
	int i;
	int iteration=0;
	int xcount, xstart;
	
	////////////FOR DEBUG ON PC//////////////
	char fn_a[20], fn_info[20];
	sprintf(fn_a, "A%d.dat", nMyRank);
	sprintf(fn_info, "INFO%d.dat", nMyRank);
	/////////////////////////////////////////

	//Read A,b from file and construct r
	if ((fa=fopen(fn_a, "r+"))==NULL
		|| (fb=fopen("B.dat", "r+"))==NULL
		|| (finfo=fopen(fn_info, "r+"))==NULL)
	{
		printf("File read error\n");
		exit(1);
	}
	//Read general information from
	fscanf(finfo, "%d %d %d %d", &oldrowCount, &oldcolCount, &row_index, &nnz);
	rowCount = oldrowCount << 1;	//expand from complex matrix to real matrix
	colCount = oldcolCount << 1;
	b = new double[colCount];	//b is the full length, is this so???
	r = new double[colCount+1];
	rbar = new double[colCount+1];

	p = new double[colCount+1];
	pbar = new double[colCount+1];
	oldp = new double[colCount+1];
	oldpbar = new double[colCount+1];

	q = new double[colCount+1];
	qbar = new double[colCount+1];
	oldq = new double[colCount+1];
	oldqbar = new double[colCount+1];

	z = new double[colCount];
	zbar = new double[colCount];
	
	xcount = colCount/nCommSize;
	xstart = xcount * nMyRank;
	x = new double[colCount];
	mtxA = new MatComplex[nnz];		//just-enough space to store

	double& rou = p[colCount];	//A reference locates in last space of p
	double& rouresult = oldp[colCount];

	double& pq = q[colCount];	//A reference locates in last space of q
	double& pqresult = oldq[colCount];

	double sum;
	double sumresult;
	//Read matrix A data
	for(i=0; i<nnz; i++)
	{		
		fscanf(fa, "%d %d %lf %lf", 
			&mtxA[i].row, &mtxA[i].col, &mtxA[i].real, &mtxA[i].imag);
//		printf("[%d %d] %lg + j %lg\n", 
//			mtxA[i].row, mtxA[i].col, mtxA[i].real, mtxA[i].imag);
	}
	//Read matrix B data
	for(i=0; i<oldcolCount; i++)
	{
		fscanf(fb, "%lf %lf", &b[i], &b[i+oldcolCount]);
//		if (nMyRank==0)
//			printf("[%d] b = %lg + j %lg\n", i, b[i], b[i+oldcolCount]);
	} 

	//Generate pseudo precondition matrix M^-1
	int mnnz = colCount / nCommSize; 
	int mstart = mnnz * nMyRank;
//	printf("[%d] mnnz = %d          mstart = %d\n", nMyRank, mnnz, mstart);
	mtxM = new MatElement[mnnz+1];
	for (i=0; i<mnnz; i++)
	{
		mtxM[i].col = mtxM[i].row = mstart+i;
		mtxM[i].value=1.0;		//Just generate an M
//		printf("[%d]...M[%d %d] = %lg\n", nMyRank, mstart+i, mstart+i, 1.0);
	}

	//Initialization
	memset(x, 0, colCount*sizeof(double));
	memmove(r, b, colCount*sizeof(double));
	memmove(rbar, b, colCount*sizeof(double));

//	for (i=0; i<colCount; i++)
//	{
//		if (nMyRank==1)
//			printf("[%d] r[%d] = %lg\n", nMyRank, i, r[i]);
//	}
		
	//Iteration cycles
	while(true)
	{
		iteration++;
//		if (nMyRank==0)
//		{
//			printf("-----iteration %d--------\n", iteration);
//			getch();
//		}

		//Solve precondition function M*z=r
		//Here I suppose M^-1 is already in mtxM[]
		memset(z, 0, colCount*sizeof(double));
		memset(zbar, 0, colCount*sizeof(double));
		for(i=0; i<mnnz; i++)
		{
			z[mtxM[i].row] += mtxM[i].value * r[mtxM[i].col];
//			if (iteration==1) 
//				printf("!!![%d]      r[%d] = %lg       z[%d]=%lg!!!\n",
//					nMyRank, mtxM[i].col, r[mtxM[i].col],
//					mtxM[i].row, z[mtxM[i].row]);
			zbar[mtxM[i].col] += mtxM[i].value * rbar[mtxM[i].row];
		}
/*
		for (i=0; i<colCount; i++)
		{
			printf("[%d](it%d) z[%d] = %lg\n", nMyRank, iteration, i, z[i]);
			printf("[%d](it%d) zbar[%d] = %lg\n", nMyRank, iteration, i, z[i]);
		}
*/
		rou = 0;
//		printf("\n==================\n");
		for(i=0; i<colCount; i++)
		{
			rou += z[i] * rbar[i];
//			printf ("rou=%lg   add    %lg * %lg\n",
//					rou, z[i], rbar[i]);
			p[i] = z[i];
			pbar[i] = zbar[i];
		}
//		printf("\n==================\n");
		if (iteration!=1)
		{
			for(i=0; i<colCount; i++)
			{
				p[i] += oldp[i]*rou/rouresult;
				pbar[i] += oldpbar[i]*rou/rouresult;
			}
		}
//		printf("[%d] {it%d} rou = %lg\n", nMyRank, iteration, rou);
		//Calculate the sum of p in all processes and save value in oldp
		//Notice rou is also calculated
		MPI_Allreduce(p, oldp, colCount+1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
		MPI_Allreduce(pbar, oldpbar, colCount+1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
/*
		if (nMyRank==0)
		{
			for(i=0;i<colCount;i++)	printf("p[%d]=%lg\n", i, oldp[i]);
			printf("\n");
			for(i=0;i<colCount;i++)	printf("p~[%d]=%lg\n", i, oldpbar[i]);
		}
		printf("[%d] rouresult= %lg\n", nMyRank, rouresult);
*/
		if (rouresult<MIN && rouresult>-MIN)
		{
			failed=true;
			break;	//BiCG method failed
		}
/*
		if (nMyRank==1)
			for (i=0; i<nnz; i++)
				printf("mtxA[%d, %d]= %lg + %lg\n", mtxA[i].row, mtxA[i].col, mtxA[i].real, mtxA[i].imag);
*/
		//Calculate q value
		memset(q, 0, colCount*sizeof(double));
		memset(qbar, 0, colCount*sizeof(double));
		for(i=0; i<nnz; i++)
		{
			q[mtxA[i].row] += mtxA[i].real * oldp[mtxA[i].col];
//			if (nMyRank==1)
//				printf("adding.......q[%d](%lg) = real(%lg) + p[%d](%lg)\n",
//					q[mtxA[i].row], mtxA[i].row, mtxA[i].real, mtxA[i].col, p[mtxA[i].col]);
			q[mtxA[i].row+oldrowCount] += mtxA[i].real * oldp[mtxA[i].col+oldcolCount];
			q[mtxA[i].row+oldrowCount] += mtxA[i].imag * oldp[mtxA[i].col];
			q[mtxA[i].row] -= mtxA[i].imag * oldp[mtxA[i].col+oldcolCount];

			qbar[mtxA[i].col] += mtxA[i].real * oldpbar[mtxA[i].row];
			qbar[mtxA[i].col+oldcolCount] += mtxA[i].real * oldpbar[mtxA[i].row+oldrowCount];
			qbar[mtxA[i].col] += mtxA[i].imag * oldpbar[mtxA[i].row+oldrowCount]; 
			qbar[mtxA[i].col+oldcolCount] -= mtxA[i].imag * oldpbar[mtxA[i].row];
		}

		pq = 0;
		for(i=0; i<colCount; i++)
			pq += oldpbar[i] * q[i];	//pq locates in the last element of q[]
		MPI_Allreduce(q, oldq, colCount+1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
		MPI_Allreduce(qbar, oldqbar, colCount+1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
		alpha = rouresult / pqresult;	//These two figures are worked out by MPI_Allreduce()
/*
		if (nMyRank==0)
		{
			for(i=0;i<colCount;i++)	printf("q[%d]=%lg\n", i, oldq[i]);
			printf("\n");
			for(i=0;i<colCount;i++)	printf("q~[%d]=%lg\n", i, oldqbar[i]);
			printf("\n");
			printf("alpha = %lg\nrou = %lg\npq = %lg\n", alpha, rouresult, pqresult);
		}
*/
		//Update x
		for(i=0; i<xcount; i++)
			x[xstart+i] += alpha * oldp[xstart+i];

//		if (nMyRank==0)
//			for(i=0;i<xcount;i++) printf("x[%d]=%lg\n", i+xstart, x[i+xstart]);

		//Go on to the next iteration
		//Update r value, each process calculates a part of r
		for(i=0; i<xcount; i++)
		{
			r[xstart+i] -= alpha * oldq[xstart+i];
			rbar[xstart+i] -= alpha * oldqbar[xstart+i];
			sum = r[xstart+i] *  r[xstart+i];
		}
		MPI_Allreduce(&sum, &sumresult, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);

		//Test convergence of iteration
//		if (nMyRank==0)
//			printf("converge value = %lg\n", iteration, sumresult);
		if (sumresult < CONV_MIN) //if small enough
			break;

		//Update r in every process
		MPI_Allgather(r+xstart, xcount, MPI_DOUBLE, r, xcount, MPI_DOUBLE, MPI_COMM_WORLD);

//		if (nMyRank==0){
//			printf("=======================================\n");
//			getch();
//		}		
	}

	//Collect x value and output answer
	MPI_Gather(x+xstart, xcount, MPI_DOUBLE, x, xcount, MPI_DOUBLE, 0, MPI_COMM_WORLD);

	if (!failed)
	{
		if (nMyRank==0)
		{
			fout = fopen("result.dat", "w");
			for(i=0; i<oldcolCount; i++)
			{
				fprintf(fout, "x[%d] = %lg + j %lg\n", i, x[i], x[i+oldcolCount]);
				printf("x[%d] = %lg + j %lg\n", i, x[i], x[i+oldcolCount]);
			}
			printf("Converged after %d iterations\n",iteration);
			printf("Finished.\n");
		}
	}
	else
		printf("Method failed.\n");

	fclose(fa);
	fclose(fb);
	fclose(finfo);
	fclose(fout);
	MPI_Finalize();
}