Smith-waterman算法 openmp+mpi实现

//此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 characterint 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 arraychar 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 HiScoreint SWArray[MAXLEN][MAXLEN];//S-W Matrixchar MaxA[MAXLEN];char MaxB[MAXLEN];char OptA[MAXLEN];char OptB[MAXLEN];int MaxAcounter = 1;//MaxA counterint MaxBcounter = 1;//MaxB counterint 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 rowsint relmaxpos[2];//holds position of relmaxif(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 countersMaxA[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 relmaxMaxB[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 strAMaxA[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 forfor(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 forfor(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 alignmentfor(j=0; j<=lenB; ++j) {if(SWArray[i][j] > HiScore) {HiScore = SWArray[i][j];HiScorePos[0] = i;HiScorePos[1] = j;}}}//send Position to alignment functionMaxA[0] = strA[HiScorePos[0]-1];MaxB[0] = strB[HiScorePos[1]-1];check = Align(HiScorePos[0],HiScorePos[1]);//in the end reverse Max A and Bk=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);}