//此Smith-Waterman 算法分别用mpi与openmp实现是没问题的，但是两个混合编程的时候就会出各种问题，希望懂的能够给指条明路。。。万分感谢

#include <omp.h>
#include <time.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <mpi.h>
/*File pointers*/
FILE *ptr_file_1, *ptr_file_2,*ptr_file_3;	//test file pointers

/*Definitions*/
#define TRUE 1
#define FALSE 0
#define Match 1
#define MissMatch -1
#define GapPenalty 2
#define MAXLEN 20

/*Global Variables*/
char inputC[5];		//user input character
int inputI;
int StrLen1,StrLen2;
int intcheck = TRUE;

char holder, ch;
int filelen1 = 0;
int filelen2 = 0;
int filelen3 = 0;
int i,j,k,l,m,n,lenA,lenB,compval,top,left,left_top,rows,columns,contrl;
char FASTA1[MAXLEN],FASTA2[MAXLEN]; 
char dash = '-';

char strA[MAXLEN];			//holds 1st string to be aligned in character array		
char strB[MAXLEN];			//holds 2nd string to be aligned in character array 
int HiScore;			//holds value of highest scoring alignment(s).
int HiScorePos[2];		//holds the position of the HiScore
int SWArray[MAXLEN][MAXLEN];	//S-W Matrix

char MaxA[MAXLEN];
char MaxB[MAXLEN];
char OptA[MAXLEN];
char OptB[MAXLEN];

int MaxAcounter = 1;	//MaxA counter	
int MaxBcounter = 1;	//MaxB counter
int cont = TRUE;
int check;


int rank,size;
MPI_Status status;
int buff[3];



/*ALIGNMENT FUNCTION*/
int Align(int PosA,int PosB) {

	/*Function Variables*/
	int relmax = -1;		//hold highest value in sub columns and rows
	int relmaxpos[2];		//holds position of relmax

	if(SWArray[PosA-1][PosB-1] == 0) {
		cont = FALSE;
	}

	while(cont == TRUE) {	//until the diagonal of the current cell has a value of zero
		
		/*Find relmax in sub columns and rows*/
		for(i=PosA; i>0; --i) {
	
			if(relmax < SWArray[i-1][PosB-1]) {
				
				relmax = SWArray[i-1][PosB-1];
				relmaxpos[0] = i-1;
				relmaxpos[1] = PosB-1;
			}
		}

		for(j=PosB; j>0; --j) {

			if(relmax < SWArray[PosA-1][j-1]) {

				relmax = SWArray[PosA-1][j-1];
				relmaxpos[0]=PosA-1;
				relmaxpos[1]=j-1;
			}
		}

		/*Align strings to relmax*/
		if((relmaxpos[0] == PosA-1) && (relmaxpos[1] == PosB-1)) {	//if relmax position is diagonal from current position simply align and increment counters

			MaxA[MaxAcounter] = strA[relmaxpos[0]-1];
			++MaxAcounter;
			MaxB[MaxBcounter] = strB[relmaxpos[1]-1];
			++MaxBcounter;

		}

		else {

			if((relmaxpos[1] == PosB-1) && (relmaxpos[0] != PosA-1)) {	//maxB needs at least one '-'
				
				for(i=PosA-1; i>relmaxpos[0]-1; --i) {	//for all elements of strA between PosA and relmaxpos[0]
				
						MaxA[MaxAcounter]= strA[i-1];
						++MaxAcounter;
				}
				for(j=PosA-1; j>relmaxpos[0]; --j) {	//set dashes to MaxB up to relmax

					MaxB[MaxBcounter] = dash;
					++MaxBcounter;
				}

				MaxB[MaxBcounter] = strB[relmaxpos[1]-1];	//at relmax set pertinent strB value to MaxB
				++MaxBcounter;
			}

			if((relmaxpos[0] == PosA-1) && (relmaxpos[1] != PosB-1)) {	//MaxA needs at least one '-'

				for(j=PosB-1; j>relmaxpos[1]-1; --j) {	//for all elements of strB between PosB and relmaxpos[1]

					MaxB[MaxBcounter] = strB[j-1];
					++MaxBcounter;
				}
				for(i=PosB-1; i>relmaxpos[1]; --i) {		//set dashes to strA

					MaxA[MaxAcounter] = dash;
					++MaxAcounter;
				}

				MaxA[MaxAcounter] = strA[relmaxpos[0]-1];
				++MaxAcounter;
			}
		}

		//printf("(%i,%i)",relmaxpos[0],relmaxpos[1]);
		Align(relmaxpos[0],relmaxpos[1]);
	}

	return(cont);
}


int SWMax(int num1,int num2,int num3)
{
	int max=-1;
	if(num1>num2)
	{
		max = num1;
	}
	else if(num2>num3)
	{
		max = num2;
	}else
	{
		max =num3;
	}
	return max;
}

/*MAIN FUNCTIONS*/
int main(int argc,char ** argv)
{

	MPI_Init(&argc,&argv);//MPI初始化
	MPI_Comm_rank(MPI_COMM_WORLD,&rank);//得到当前进程标识
	MPI_Comm_size(MPI_COMM_WORLD,&size);//总的进程个数

	printf("size=%d\n",size);


	/*open first file*/
	ptr_file_1 = fopen("d:/mpi/str1.fa", "r");
	/*check to see if it opened okay*/
	if(ptr_file_1 == NULL) {
		printf("Error opening 'str1.fa'\n");
		system("PAUSE");
		exit(8);
	}

	/*open second file*/
	ptr_file_2 = fopen("d:/mpi/str2.fa", "r");
	/*check to see if it opened okay*/
	if(ptr_file_2 == NULL) {
		printf("Error opening 'str2.fa'\n");
		system("PAUSE");	
		exit(8);
	} 
	/*measure file1 length*/
	while(holder != EOF) {
		holder = fgetc(ptr_file_1);
		if(filelen1 < MAXLEN && holder !=EOF) {
			strA[filelen1] = holder;
			++filelen1;
		}		
	}
	strA[filelen1] = -1;
	holder = '0';
	/*measure file2 length*/
	while(holder != EOF) {
		holder = fgetc(ptr_file_2);
		if(filelen2 < MAXLEN && holder !=EOF) {
			strB[filelen2] = holder;
			++filelen2;
		}
	}
	strB[filelen2] = -1;
	fclose(ptr_file_1);
	fclose(ptr_file_2);								
							
	lenA = strlen(strA)-1;
	lenB = strlen(strB)-1;



	
	
	time_t  t1 = time(NULL);
	printf("t1=%d\n",t1);

	/*----------------------------问题代码处*---------------------------*/
	for(contrl=0;contrl<lenA+lenB+1;contrl++)
	{
		if(contrl<=lenA)
		{
			rows = contrl;
			j = 0;
		}
		else 
		{
			rows = lenA;
			j = contrl-lenA;
		}

		if(rank==0)
		{

			//#pragma omp parallel for
			for(columns = j;columns <= lenB;columns++)
			{
				//printf("%d\n",columns);
				if(rows == 0||columns == 0)
				{
					SWArray[rows][columns]=0;
				}
				else
				{
					left = SWArray[rows][columns-1] - GapPenalty;
					top = SWArray[rows-1][columns] - GapPenalty;
					if(strA[rows-1] == strB[columns-1]) 
					{				
						left_top = (SWArray[rows-1][columns-1] + Match);	
					}else
					{
						left_top = SWArray[rows-1][columns-1] + MissMatch;
					}
					compval = SWMax(left,top,left_top);
					if(compval<0) compval = 0;
					SWArray[rows][columns] = compval;
					buff[0]=compval;
					buff[1]=rows;
					buff[2]=columns;
					MPI_Send(buff,4,MPI_INT,rank+1,99,MPI_COMM_WORLD);
					MPI_Recv(buff,4,MPI_INT,rank+1,99,MPI_COMM_WORLD,&status);
					SWArray[buff[1]][buff[2]]=buff[0];
					compval = 0;
				}	
			}
		}
		else
		{
			//#pragma omp parallel for
			for(columns = j;columns <= lenB;columns++)
			{
				//printf("%d\n",columns);
				if(rows == 0||columns == 0)
				{
					SWArray[rows][columns]=0;
				}
				else
				{
					left = SWArray[rows][columns-1] - GapPenalty;
					top = SWArray[rows-1][columns] - GapPenalty;
					if(strA[rows-1] == strB[columns-1]) 
					{				
						left_top = (SWArray[rows-1][columns-1] + Match);	
					}else
					{
						left_top = SWArray[rows-1][columns-1] + MissMatch;
					}
					compval = SWMax(left,top,left_top);
					if(compval<0) compval = 0;
					SWArray[rows][columns] = compval;
					buff[0]=compval;
					buff[1]=rows;
					buff[2]=columns;
					MPI_Send(buff,4,MPI_INT,rank-1,99,MPI_COMM_WORLD);
					MPI_Recv(buff,4,MPI_INT,rank-1,99,MPI_COMM_WORLD,&status);
					SWArray[buff[1]][buff[2]]=buff[0];
					compval = 0;
				}	
			}
			
		}
		MPI_Barrier(MPI_COMM_WORLD);
	}


	/*PRINT S-W Table*/
	ptr_file_3 = fopen("str3.fa", "w+");
	if(ptr_file_3 == NULL) {
		printf("Error opening 'str3.fa'\n");
		system("PAUSE");	
		exit(8);
	}
	fprintf(ptr_file_3,"   0");
	for(i = 0; i <= lenB; ++i) {
		fprintf(ptr_file_3,"  %c",strB[i]);
	}
	fprintf(ptr_file_3,"\n");

	for(i = 0; i <= lenA; ++i) {
		if(i < 1) {
			fprintf(ptr_file_3,"0");
		}

		if(i > 0) {
			fprintf(ptr_file_3,"%c",strA[i-1]);
		}

		for(j = 0; j <= lenB; ++j) {
				fprintf(ptr_file_3,"%3i",SWArray[i][j]);
		}
		fprintf(ptr_file_3,"\n");
	}
	fclose(ptr_file_3);
	/*MAKE ALIGNMENTT*/
	for(i=0; i<=lenA; ++i) {	//find highest score in matrix: this is the starting point of an optimal local alignment

		for(j=0; j<=lenB; ++j) {
	
			if(SWArray[i][j] > HiScore) {

				HiScore = SWArray[i][j];
				HiScorePos[0] = i;
				HiScorePos[1] = j;
			}
		}
	}

	//send Position to alignment function
	MaxA[0] = strA[HiScorePos[0]-1];
	MaxB[0] = strB[HiScorePos[1]-1];

	check = Align(HiScorePos[0],HiScorePos[1]);

	//in the end reverse Max A and B
	k=0;
	for(i = strlen(MaxA)-1; i > -1; --i) {
		OptA[k] = MaxA[i];
		++k;
	}

	k=0;
	for(j=strlen(MaxB)-1; j > -1; --j) {
		OptB[k] = MaxB[j];
		++k;
	}
 
	printf("%s\n%s	\n",OptA,OptB);


	time_t  t2 = time(NULL);
	printf("t2=%d\n",t2);
	printf("Time-consuming is:%ds\n",(t2-t1));
	
	

	MPI_Finalize();

	return(0);

}