Hitmap 1.3
 All Data Structures Namespaces Files Functions Variables Typedefs Friends Macros Groups Pages
SWpar_ref.c
Go to the documentation of this file.
1 
28 /*
29  * <license>
30  *
31  * Hitmap v1.2
32  *
33  * This software is provided to enhance knowledge and encourage progress in the scientific
34  * community. It should be used only for research and educational purposes. Any reproduction
35  * or use for commercial purpose, public redistribution, in source or binary forms, with or
36  * without modifications, is NOT ALLOWED without the previous authorization of the copyright
37  * holder. The origin of this software must not be misrepresented; you must not claim that you
38  * wrote the original software. If you use this software for any purpose (e.g. publication),
39  * a reference to the software package and the authors must be included.
40  *
41  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER AND CONTRIBUTORS "AS IS" AND ANY
42  * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
43  * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
44  * THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
45  * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
46  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
47  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
48  * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
49  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
50  *
51  * Copyright (c) 2007-2015, Trasgo Group, Universidad de Valladolid.
52  * All rights reserved.
53  *
54  * More information on http://trasgo.infor.uva.es/
55  *
56  * </license>
57 */
58 
59 // Include the common variables and utilities for all the versions
60 #include "SWcommon.h"
61 
62 // Include the common variables and utilities for the C reference versions
63 #include "SWcommon_ref.h"
64 
65 // MPI include
66 #include <mpi.h>
67 
68 
72 #define BEGIN(coord,dims,size) (coord * size / dims)
73 
77 #define END(coord,dims,size) ((( coord + 1 ) * size / dims)-1)
78 
82 #define SIZE(coord,dims,size) (END(coord,dims,size)-BEGIN(coord,dims,size)+1)
83 
87 void distributeProtein(aa_t * protein, aa_t *lprotein, int dim);
88 
95 void recv_Back(int dir, int * pos);
96 
101 void send_Back(int dir, int pos);
102 
103 
107 void composeOutput(aa_t * out, int * outpos, int * outsize, aa_t * lout, int loutsize, int * gBegin, int * gEnd);
108 
109 
110 
111 
112 
113 // **********************************************************************
114 // Declarations of the different phases
115 // **********************************************************************
116 
120 void phase_read_sequences();
121 
126 
130 void phase_recv_hmatrix();
131 
135 void phase_comp_hmatrix();
136 
140 void phase_send_hmatrix();
141 
145 void phase_recv_back();
146 
150 void phase_comp_back();
151 
155 void phase_send_back();
156 
161 
162 
166 void debug_show_matrices();
167 
172 
173 
174 
175 
176 // **********************************************************************
177 // Parallel related global variables
178 // **********************************************************************
179 
181 int rank;
183 int nProcs;
185 int lsize[2];
187 int dims[2] = {0,0};
189 int lcoords[2];
191 int lbegin[2];
193 MPI_Comm comm;
194 
196 MPI_Op opHMaxLoc;
197 
198 
199 
200 // **********************************************************************
201 // Global root variables
202 // **********************************************************************
203 
204 // Protein sequences
207 
208 // Output proteins with the alignment
211 
212 // X and Y coordinates of the maximum value
213 int xMax;
214 int yMax;
215 
216 // Begin coordinates of the match
217 int xBegin;
218 int yBegin;
219 
220 // Match lengths (total, sequence 1, sequence 2)
224 
225 // Output match begin
227 
228 
229 // **********************************************************************
230 // Global (C) local (parallel) variables
231 // **********************************************************************
232 
233 // Protein sequences
236 
237 // Local output proteins with the alignment
240 
241 // Similarity matrix
242 h_t * H;
243 #define h(i,j) (H[((i)*(lsize[1]+1))+(j)])
244 
245 // Traceback matrices
248 #define xback(i,j) (xTraceback[((i)*(lsize[1]+1))+(j)])
249 #define yback(i,j) (yTraceback[((i)*(lsize[1]+1))+(j)])
250 
251 // Local Max value
253 
254 // X and Y coordinates of the maximum value
255 int lxMax;
256 int lyMax;
257 
258 // Begin coordinates of the match
261 
262 // Match lengths (total, sequence 1, sequence 2)
266 
267 // Output match begin
269 
270 
271 
275 int main(int argc, char * argv[]){
276 
277  // **********************************************************************
278  // Init MPI
279  // **********************************************************************
280  MPI_Init(&argc,&argv);
281  MPI_Comm_rank(MPI_COMM_WORLD,&rank);
282  MPI_Comm_size(MPI_COMM_WORLD,&nProcs);
283 
284 
285  // MPI Datatype for the double_3int struct.
286  double_3int struct_d3i;
287 
288  int blocks[4] = {1,1,1,1};
289  MPI_Aint displ[5];
290  MPI_Datatype types[4] = {MPI_DOUBLE, MPI_INT, MPI_INT, MPI_INT};
291 
292  MPI_Address( &struct_d3i.val, &displ[0] );
293  MPI_Address( &struct_d3i.rank, &displ[1] );
294  MPI_Address( &struct_d3i.i, &displ[2] );
295  MPI_Address( &struct_d3i.j, &displ[3] );
296  MPI_Address( &struct_d3i, &displ[4] );
297 
298  displ[0] -= displ[4];
299  displ[1] -= displ[4];
300  displ[2] -= displ[4];
301  displ[3] -= displ[4];
302 
303  MPI_Type_struct(4, blocks, displ, types, &MPI_DOUBLE_3INT);
304  MPI_Type_commit( &MPI_DOUBLE_3INT );
305 
306  // Create the operation to reduce the double_3int values.
307  MPI_Op_create( HMaxLoc, 1, &opHMaxLoc );
308 
309 
310 
311  // **********************************************************************
312  // Check input parameters
313  // **********************************************************************
314  input_parameters(argc, argv);
315 
316  // **********************************************************************
317  // Layout
318  // **********************************************************************
319 
320  int periods[2] = {0,0};
321  MPI_Dims_create(nProcs, 2, dims);
322  MPI_Cart_create(MPI_COMM_WORLD,2,dims,periods,0,&comm);
323  MPI_Cart_coords(comm,rank,2,lcoords);
324 
325  // Print processor info
326  printfroot("=== Parallel ===\n"
327  " Processors: %d\n"
328  " Processors: %d x %d\n\n",nProcs,dims[0], dims[1]);
329 
330  // Begin, end
331  int k;
332 
333  for(k=0; k<2; k++){
334  lbegin[k] = BEGIN(lcoords[k],dims[k],size[k]);
335  lsize[k] = SIZE(lcoords[k],dims[k],size[k]);
336 
337  if(lsize[k] < 1){
338  fprintf(stderr,"There are inactive processors\n");
339  fflush(stderr);
340  MPI_Abort(comm,EXIT_FAILURE);
341 
342  }
343  }
344 
345  // Print info of each processor
346 #ifdef DEBUG
347  ALL_SEQUENTIAL_BEGIN(seq1);
348  printf("Processor: %d/%d, ",rank,nProcs);
349  printf(" Coords %d x %d\n",lcoords[0],lcoords[1]);
350  printf(" Dim1: Block begin %2d- (size %2d)\n",lbegin[0], lsize[0]);
351  printf(" Dim2: Block begin %2d- (size %2d)\n",lbegin[1], lsize[1]);
352  printf("\n");
353  ALL_SEQUENTIAL_END(seq1);
354 #endif
355 
356  // **********************************************************************
357  // Timers
358  // **********************************************************************
359  double time_initS;
360  double time_initE;
361  double time_compS;
362  double time_compE;
364  double time_phaseS;
365  double time_phaseE;
367  double time_comp_read = 0;
368  double time_comp_hmatrix = 0;
369  double time_comp_back = 0;
371  double time_comp_distr = 0;
372  double time_comp_syncH = 0;
373  double time_comp_syncB = 0;
374  double time_comp_compose = 0;
377  // **********************************************************************
378  // INIT PHASE: Structure allocation and PAM matrix load
379  // **********************************************************************
380 
381  // Init time
382  MPI_Barrier(comm);
383  time_initS = currentTime();
384 
385  // Every one do this, the pam matrix is small and all procs
386  // need it, we don't bother to read it and distribute it.
387  initPAM(pam_name);
388 
389  // Open the proteins files
392 
393  // Allocate memory for the proteins (Only root has them)
394  if(is_root()){
395  p1 = malloc(sizeof(aa_t)*((size_t)(size[0]+1)));
396  p2 = malloc(sizeof(aa_t)*((size_t)(size[1]+1)));
397  CHECK_NULL(p1);
398  CHECK_NULL(p2);
399  }
400 
401  // Local proteins
402  lp1 = malloc((size_t) ( lsize[0]+1 ) * sizeof(aa_t));
403  lp2 = malloc((size_t) ( lsize[1]+1 ) * sizeof(aa_t));
404  CHECK_NULL(lp1);
405  CHECK_NULL(lp2);
406 
407  // Local Similarity matrix
408  H = malloc(sizeof(h_t)*(size_t)((lsize[0]+1)*(lsize[1]+1)));
409  CHECK_NULL(H);
410 
411  // Traceback matrices
412  xTraceback = malloc(sizeof(trace_t)*(size_t)((lsize[0]+1)*(lsize[1]+1)));
413  yTraceback = malloc(sizeof(trace_t)*(size_t)((lsize[0]+1)*(lsize[1]+1)));
416 
417  // Allocate size of the output (local)
418  lout1 = malloc(sizeof(aa_t)*((size_t)(lsize[0]+lsize[1]+1)));
419  lout2 = malloc(sizeof(aa_t)*((size_t)(lsize[0]+lsize[1]+1)));
420  CHECK_NULL(lout1);
421  CHECK_NULL(lout2);
422 
423  // Allocate size of the output (global, only root)
424  if(is_root()){
425  out1 = malloc(sizeof(aa_t)*((size_t)(size[0]+size[1]+1)));
426  out2 = malloc(sizeof(aa_t)*((size_t)(size[0]+size[1]+1)));
427  CHECK_NULL(out1);
428  CHECK_NULL(out2);
429  }
430 
431  // Init time
432  time_initE = currentTime();
433 
434 
435  // **********************************************************************
436  // WHOLE COMPUTATION PHASE: iteration of: Read, H matrix, Backtracking
437  // **********************************************************************
438 
439  // Computation time
440  time_compS = currentTime();
441 
442  int iter;
443  for(iter=0;iter<iterations;iter++){
444 
445 #ifdef DEBUG
446  printfroot("=== Iteration: %d/%d ===\n\n",iter,iterations);
447 #endif
448 
449  // **********************************************************************
450  // READ PHASE
451  // **********************************************************************
452  time_phaseS = currentTime();
454  time_phaseE = currentTime();
455  time_comp_read += (time_phaseE-time_phaseS);
456 
457 
458  // **********************************************************************
459  // DISTRIBUTION PHASE
460  // **********************************************************************
461  time_phaseS = currentTime();
463  time_phaseE = currentTime();
464  time_comp_distr += (time_phaseE-time_phaseS);
465 
466 
467  // **********************************************************************
468  // RECV H PHASE
469  // **********************************************************************
470  time_phaseS = currentTime();
471  // Recv the needed H values from the neighbors
473  time_phaseE = currentTime();
474  time_comp_syncH += (time_phaseE-time_phaseS);
475 
476 
477  // **********************************************************************
478  // H MATRIX PHASE
479  // **********************************************************************
480  time_phaseS = currentTime();
482  time_phaseE = currentTime();
483  time_comp_hmatrix += (time_phaseE-time_phaseS);
484 
485 
486  // **********************************************************************
487  // SEND H PHASE
488  // **********************************************************************
489  time_phaseS = currentTime();
490  // Send the needed H values to the neighbors
492  time_phaseE = currentTime();
493  time_comp_syncH += (time_phaseE-time_phaseS);
494 
495 #ifdef DEBUG
497 #endif
498 
499 
500  // **********************************************************************
501  // RECV BACKTRAING PHASE
502  // **********************************************************************
503  time_phaseS = currentTime();
504  phase_recv_back();
505  time_phaseE = currentTime();
506  time_comp_syncB += (time_phaseE-time_phaseS);
507 
508 
509  // **********************************************************************
510  // BACKTRACKING PHASE
511  // **********************************************************************
512  time_phaseS = currentTime();
513  phase_comp_back();
514  time_phaseE = currentTime();
515  time_comp_back += (time_phaseE-time_phaseS);
516 
517 
518  // **********************************************************************
519  // SEND BACKTRAING PHASE
520  // **********************************************************************
521  time_phaseS = currentTime();
522  phase_send_back();
523  time_phaseE = currentTime();
524  time_comp_syncB += (time_phaseE-time_phaseS);
525 
526 
527  // **********************************************************************
528  // COMPOSE FINAL PROTEIN PHASE
529  // **********************************************************************
530  time_phaseS = currentTime();
532  time_phaseE = currentTime();
533  time_comp_compose += (time_phaseE-time_phaseS);
534 
535 #ifdef DEBUG
537 #endif
538 
539 
540  } // Iterations Loop
541 
542  // Computation time
543  MPI_Barrier(comm);
544  time_compE = currentTime();
545 
546 
547  // **********************************************************************
548  // Show results
549  // **********************************************************************
550 
551 
552  // Reduce the timers
553  // [local,MIN,AVG,MAX][READ,HMAT,BACK,DIST,SYNH,SYNB,COMP]
554 
555 #define T_LOCAL 0
556 #define T_MIN 1
557 #define T_AVG 2
558 #define T_MAX 3
559 
560 #define T_READ 0
561 #define T_HMAT 1
562 #define T_BACK 2
563 #define T_DIST 3
564 #define T_SYNH 4
565 #define T_SYNB 5
566 #define T_COMP 6
567 
568 
569  double timers[4][7];
570  timers[T_LOCAL][0] = time_comp_read;
571  timers[T_LOCAL][1] = time_comp_hmatrix;
572  timers[T_LOCAL][2] = time_comp_back;
573  timers[T_LOCAL][3] = time_comp_distr;
574  timers[T_LOCAL][4] = time_comp_syncH;
575  timers[T_LOCAL][5] = time_comp_syncB;
576  timers[T_LOCAL][6] = time_comp_compose;
577 
578 #define allTimers(timer) timers[T_MIN][timer],timers[T_AVG][timer],timers[T_MAX][timer]
579 
580  MPI_Reduce(timers[0], timers[1], 7, MPI_DOUBLE, MPI_MIN, 0, comm);
581  MPI_Reduce(timers[0], timers[2], 7, MPI_DOUBLE, MPI_SUM, 0, comm);
582  MPI_Reduce(timers[0], timers[3], 7, MPI_DOUBLE, MPI_MAX, 0, comm);
583 
584  int i;
585  for(i=0;i<7;i++) timers[T_AVG][i] /= nProcs; // AVG
586 
587  // Final result.
588  if(is_root()){
589 
592 
593  printf("=== Result ===\n");
594  printf("Init time: %f\n",time_initE-time_initS);
595  printf("Comp time: %f\n",time_compE-time_compS);
596  printf(" Read protein: min: %f, avg: %f, max: %f \n",allTimers(T_READ));
597  printf(" Distribute: min: %f, avg: %f, max: %f \n",allTimers(T_DIST));
598  printf(" H matrix: min: %f, avg: %f, max: %f \n",allTimers(T_HMAT));
599  printf(" sync H mat: min: %f, avg: %f, max: %f \n",allTimers(T_SYNH));
600  printf(" Backtracking: min: %f, avg: %f, max: %f \n",allTimers(T_BACK));
601  printf(" Sync Back: min: %f, avg: %f, max: %f \n",allTimers(T_SYNB));
602  printf(" Compose: min: %f, avg: %f, max: %f \n",allTimers(T_COMP));
603 
604  printf("Last Match length of A: %d\n", match_length1);
605  printf("Last Match length of B: %d\n", match_length2);
606  printf("Last Alignment length: %d\n", match_length);
607 
608 
609  // Print the result for small inputsets
610  if(match_length < 50){
611  printf("\n");
612  printf("MCH1: "); printProteinMatch(p1,size[0],xBegin,xMax);
613  printf("MCH2: "); printProteinMatch(p2,size[1],yBegin,yMax);
614 
615  printf("\n");
616  printf("OUT1: "); printProtein(&out1[match_pos],match_length);
617  printf("OUT2: "); printProtein(&out2[match_pos],match_length);
618  } else {
619  printf("\nInputset too big to print to stdout\n");
620  }
621  }
622 
623 
624  // Free all the resources
625  free(lp1);
626  free(lp2);
627  free(xTraceback);
628  free(yTraceback);
629  free(H);
630  free(lout1);
631  free(lout2);
632  MPI_Op_free(&opHMaxLoc);
633  MPI_Comm_free(&comm);
634  MPI_Type_free(&MPI_DOUBLE_3INT);
635 
636  if(is_root()){
637  free(p1);
638  free(p2);
639  free(out1);
640  free(out2);
641  }
642 
643  MPI_Finalize();
644  return EXIT_SUCCESS;
645 
646 }
647 
648 
649 
650 void phase_read_sequences(){
651 
652  if(! is_root() ) return;
653 
654  // Read the proteins
655  readProtein(&pfile1, p1, size[0]);
656  readProtein(&pfile2, p2, size[1]);
657 
658 #ifdef DEBUG
659  printf("PRT1: "); printProtein(p1,size[0]);
660  printf("PRT2: "); printProtein(p2,size[1]);
661  printf("\n");
662 #endif
663 
664 
665 }
666 
667 
669 
670  // Send and Recv the proteins from the root proc.
673 
674  // Print info of each processor
675 #ifdef DEBUG
676  ALL_SEQUENTIAL_BEGIN(seq2)
677  printf("Processor: %d/%d, ",rank,nProcs);
678  printf("Coords %d x %d\n",lcoords[0],lcoords[1]);
679  printf("Local protein 1: ");
680  printProtein(lp1,lsize[0]);
681  printf("Local protein 2: ");
682  printProtein(lp2,lsize[1]);
683  printf("\n");
684  ALL_SEQUENTIAL_END(seq2)
685 #endif
686 
687 }
688 
689 
690 
691 void phase_recv_hmatrix(){
692 
693  int pcoords[2];
694  int prank;
695 
696  // Arrays to receive
697  h_t array1[lsize[1]];
698  h_t array2[lsize[0]+1];
699 
700  // X-1
701  pcoords[0] = lcoords[0]-1;
702  pcoords[1] = lcoords[1];
703 
704  if(pcoords[0] >= 0){
705 
706  // Get the rank of the neighbor
707  MPI_Cart_rank(comm,pcoords,&prank);
708 
709  // Receive the array
710  MPI_Recv(array1, lsize[1], MPI_H, prank, 0, comm, MPI_STATUS_IGNORE);
711 
712  int j;
713  for(j=0; j<lsize[1]; j++){
714  h(0,j+1) = array1[j];
715  }
716 
717  }
718 
719 
720  // Y-1
721  pcoords[0] = lcoords[0];
722  pcoords[1] = lcoords[1]-1;
723 
724 
725  if(pcoords[1] >= 0){
726 
727  // Get the rank of the neighbor
728  MPI_Cart_rank(comm,pcoords,&prank);
729 
730  // Receive the array
731  MPI_Recv(array2, lsize[0]+1, MPI_H, prank, 0,comm, MPI_STATUS_IGNORE);
732 
733  int i;
734  for(i=0; i<lsize[0]+1; i++){
735  h(i,0) = array2[i];
736  }
737  }
738 
739 
740 }
741 
742 
743 
744 void phase_send_hmatrix(){
745 
746  int pcoords[2];
747  int prank;
748 
749  // Arrays to send
750  h_t array1[lsize[1]];
751  h_t array2[lsize[0]+1];
752 
753  // X-1
754  pcoords[0] = lcoords[0]+1;
755  pcoords[1] = lcoords[1];
756 
757  if(pcoords[0] < dims[0]){
758 
759  // Get the rank of the neighbor
760  MPI_Cart_rank(comm,pcoords,&prank);
761 
762  int j;
763  for(j=0; j<lsize[1]; j++){
764  array1[j] = h(lsize[0],j+1);
765  }
766 
767  MPI_Send(array1, lsize[1], MPI_H, prank, 0, comm);
768  }
769 
770 
771  // Y-1
772  pcoords[0] = lcoords[0];
773  pcoords[1] = lcoords[1]+1;
774 
775  if(pcoords[1] < dims[1]){
776 
777  // Get the rank of the neighbor
778  MPI_Cart_rank(comm,pcoords,&prank);
779 
780  int i;
781  for(i=0; i<lsize[0]+1; i++){
782  array2[i] = h(i,lsize[1]);
783  }
784  MPI_Send(array2, lsize[0]+1, MPI_H, prank, 0,comm);
785  }
786 
787 }
788 
789 
790 
791 
792 
793 
794 void phase_comp_hmatrix(){
795 
796  // Compute the Similarity Matrix (H)
797 
798  // Initialize values
799  int i,j;
800 
801  if(lcoords[1] == 0){
802  for (i=0;i<=lsize[0];i++){
803  xback(i,0) = -1;
804  yback(i,0) = -1;
805  h(i,0)=0;
806  }
807  }
808 
809  if(lcoords[0] == 0){
810  for (j=0;j<=lsize[1];j++){
811  xback(0,j) = -1;
812  yback(0,j) = -1;
813  h(0,j)=0;
814  }
815  }
816 
817  // Maximum values to retrieve the alignment
818  lMax = -1;
819  lxMax = -1;
820  lyMax = -1;
821 
822  // Perform wavefront algorithm
823  for(i=1; i<lsize[0]+1; i++){
824  for(j=1; j<lsize[1]+1; j++){
825 
826  // Get the scores
827  h_t diag = h(i-1,j-1) + pam(lp1[i],lp2[j]);
828  h_t down = h(i-1,j ) + gapPenalty;
829  h_t right = h(i ,j-1) + gapPenalty;
830 
831  // Calculate the maximum
832  int idx;
833  h_t max = MAX4(diag,down,right,0,&idx);
834  h(i,j) = max;
835 
836  // Set the back trace variables.
837  if ( idx == 0 ){ // max == diag
838  xback(i,j) = (trace_t) (i-1);
839  yback(i,j) = (trace_t) (j-1);
840  } else if ( idx == 1 ) { // max == down
841  xback(i,j) = (trace_t) (i-1);
842  yback(i,j) = (trace_t) j;
843  } else if ( idx == 2 ) { // max == right
844  xback(i,j) = (trace_t) i;
845  yback(i,j) = (trace_t) (j-1);
846  } else { // max == 0
847  xback(i,j) = -1;
848  yback(i,j) = -1;
849  }
850 
851  // Maximum
852  if(max >= lMax){
853  lxMax = i;
854  lyMax = j;
855  lMax = max;
856  }
857 
858  } // End loop j;
859  } // End loop i;
860 
861 }
862 
863 
864 
865 void debug_show_matrices(){
866 
867  // Print info of each processor
868  #ifdef DEBUG
869 
870  int i,j;
871 
872  ALL_SEQUENTIAL_BEGIN(seq3);
873 
874  printf("Processor: %d/%d, ",rank,nProcs);
875  printf("Coords %d x %d\n",lcoords[0],lcoords[1]);
876 
877  // Print
878  printf(" ");
879  for (j=0;j<=lsize[1];j++){
880  if(j==0)
881  printf(" -");
882  else
883  printf("%6c",AA2char(lp2[j]));
884  }
885  printf("\n");
886  for (i=0;i<=lsize[0];i++){
887 
888  if(i==0)
889  printf("- |");
890  else
891  printf("%c |",AA2char(lp1[i]));
892 
893  for (j=0;j<=lsize[1];j++){
894  printf("%6.1f",(double)h(i,j));
895  }
896  printf("\n");
897  }
898  printf("\n");
899 
900  ALL_SEQUENTIAL_END(seq3);
901  #endif
902 
903 
904 
905  #ifdef DEBUG
906 
907  ALL_SEQUENTIAL_BEGIN(seq4)
908 
909  printf("Processor: %d/%d, ",rank,nProcs);
910  printf("coords %d x %d\n",lcoords[0],lcoords[1]);
911 
912  // Print
913  printf(" ");
914  for (j=0;j<=lsize[1];j++){
915  if(j==0)
916  printf(" -");
917  else
918  printf("%8c",AA2char(lp2[j]));
919  }
920  printf("\n");
921  for (i=0;i<=lsize[0];i++){
922 
923  if(i==0)
924  printf("- |");
925  else
926  printf("%c |",AA2char(lp1[i]));
927 
928  for (j=0;j<=lsize[1];j++){
929  printf("(%2d,%2d) ",xback(i,j),yback(i,j));
930  }
931  printf("\n");
932  }
933  printf("\n");
934 
935  ALL_SEQUENTIAL_END(seq4)
936 
937  #endif
938 
939 
940 }
941 
942 
943 
944 void phase_recv_back(){
945 
946  // Determine who has the maximum value
947  double_3int Max, outMax;
948 
949  // Fill Max with our values.
950  Max.val = lMax;
951  Max.rank = rank;
952  Max.i = lbegin[0] + lxMax;
953  Max.j = lbegin[1] + lyMax;
954 
955  // Reduce the values.
956  MPI_Allreduce(&Max, &outMax, 1, MPI_DOUBLE_3INT, opHMaxLoc, comm);
957 
958 #ifdef DEBUG
959  // Print the maximum
960  printfroot("Max %f in proc: %d\n\n",outMax.val, outMax.rank);
961 #endif
962 
963  // Receive the point to continue the backtracking algorithm.
964  if(rank == outMax.rank){
965 
966  // Avoid possible deadlock, if we are the processor with
967  // the maximum value, we receive the messages but we ignore them.
968  int pos;
969  recv_Back(0, &pos);
970  recv_Back(1, &pos);
971 
972  // Get the position from a neighbor
973  } else {
974 
975  // By default we do not have the path.
976  lxMax=-1;
977  lyMax=-1;
978 
979  // Receive the path from the neigh from dim 0
980  int pos = -1;
981  recv_Back(0, &pos);
982 
983  if(pos != -1){
984  lxMax = lsize[0];
985  lyMax = pos;
986  }
987 
988  // Receive the path from the neigh from dim 1
989  pos = -1;
990  recv_Back(1, &pos);
991 
992  if(pos != -1){
993  lxMax = pos;
994  lyMax = lsize[1];
995  }
996  }
997 
998 
999 }
1000 
1001 
1002 
1003 
1004 
1005 void phase_send_back(){
1006 
1007  int i = lxBegin;
1008  int j = lyBegin;
1009 
1010  // Send a message to the neighbors to continue or stop
1011  // backtracking algorithm
1012  if( i == -1 || j == -1 ){
1013  // The current processors do not have any point of the path of the alignment
1014  send_Back(0, -1);
1015  send_Back(1, -1);
1016  } else if(h(i,j) > 0 && i == 0 ){
1017  // Neighbor in dim 0 should continue the algorithm.
1018  send_Back(0, j);
1019  send_Back(1, -1);
1020  } else if (h(i,j) > 0 && j == 0 ){
1021  // Neighbor in dim 1 should continue the algorithm.
1022  send_Back(0, -1);
1023  send_Back(1, i);
1024  } else {
1025  // Neighbors do not continue the search.
1026  send_Back(0, -1);
1027  send_Back(1, -1);
1028  }
1029 
1030 }
1031 
1032 
1033 
1034 void phase_comp_back(){
1035 
1036 
1037  // reset to max point to do alignment
1038  // We start at the maximum point
1039  int i=lxMax;
1040  int j=lyMax;
1041 
1042  // Position of the output proteins.
1043  int pos=lsize[0]+lsize[1]+1;
1044  int outsize=0;
1045 
1046  // Backtracking algorithm
1047  while (i>0 && j>0 && h(i,j) > 0){
1048 
1049  int movex = (i>xback(i,j));
1050  int movey = (j>yback(i,j));
1051 
1052  pos--;
1053  outsize++;
1054 
1055  if (movex && movey) {
1056  lout1[pos]=lp1[i];
1057  lout2[pos]=lp2[j];
1058  } else if (movey) {
1059  lout1[pos]=GAP_AA;
1060  lout2[pos]=lp2[j];
1061  } else if (movex) {
1062  lout1[pos]=lp1[i];
1063  lout2[pos]=GAP_AA;
1064  }
1065 
1066  // Move to the next point.
1067  int tempi=i;
1068  int tempj=j;
1069  i=xback(tempi,tempj);
1070  j=yback(tempi,tempj);
1071  }
1072 
1073  pos--;
1074 
1075  lmatch_pos = pos;
1076 
1077  lmatch_length = outsize;
1078  lmatch_length1 = lxMax - i;
1079  lmatch_length2 = lyMax - j;
1080 
1081  lxBegin = i;
1082  lyBegin = j;
1083 
1084 }
1085 
1086 
1087 
1089 
1090 #ifdef DEBUG
1091 
1092  // Local result print
1093 ALL_SEQUENTIAL_BEGIN(seq5);
1094  printf("Processor: %d/%d, ",rank,nProcs);
1095  printf("Coords %d x %d\n",lcoords[0],lcoords[1]);
1096 
1097  printf(" First point (%d,%d) -> Last point (%d,%d)\n",xMax,yMax,lxBegin,lyBegin);
1098 
1099  printf(" Local Match length of A: %d\n", lmatch_length1);
1100  printf(" Local Match length of B: %d\n", lmatch_length2);
1101  printf(" Local Alignment length: %d\n", lmatch_length);
1102 
1103  printf(" MCH1: "); printProteinMatch(lp1,lsize[0],lxBegin,lxMax);
1104  printf(" MCH2: "); printProteinMatch(lp2,lsize[1],lyBegin,lyMax);
1105 
1106  printf(" LOCAL OUT1: "); printProtein(&lout1[lmatch_pos],lmatch_length);
1107  printf(" LOCAL OUT2: "); printProtein(&lout2[lmatch_pos],lmatch_length);
1108  printf("\n");
1109 ALL_SEQUENTIAL_END(seq5);
1110 #endif
1111 
1112 
1113 }
1114 
1115 
1116 void distributeProtein(aa_t * protein, aa_t *lprotein, int dim){
1117 
1118  // IRecv
1119  MPI_Request request;
1120  MPI_Irecv(&lprotein[1], lsize[dim], MPI_AA, 0, 0, comm, &request);
1121 
1122  // Non root processor can exit
1123  if(rank != 0){
1124  MPI_Wait(&request,MPI_STATUS_IGNORE);
1125  return;
1126  }
1127 
1128  // Send from root
1129  int p;
1130  for(p=0; p<nProcs; p++){
1131 
1132  // Get the coordinates to calculate the begin and size
1133  int pcoords[2];
1134  MPI_Cart_coords(comm,p,2,pcoords);
1135  int pbegin = BEGIN(pcoords[dim],dims[dim],size[dim]);
1136  int psize = SIZE(pcoords[dim],dims[dim],size[dim]);
1137 
1138  // Pointer to the beginning of the local part
1139  void * buff = &protein[pbegin + 1];
1140 
1141  MPI_Send(buff, psize, MPI_AA, p, 0, comm);
1142 
1143  }
1144 
1145  // Root waits for its own part
1146  MPI_Wait(&request,MPI_STATUS_IGNORE);
1147 
1148 }
1149 
1150 
1151 
1152 
1153 
1154 void send_Back(int dir, int pos){
1155 
1156  int pcoords[2];
1157 
1158  pcoords[0] = lcoords[0];
1159  pcoords[1] = lcoords[1];
1160 
1161  if(dir == 0)
1162  pcoords[0]--;
1163  else
1164  pcoords[1]--;
1165 
1166  // Exit if we have no one to send.
1167  if(pcoords[0] < 0) return;
1168  if(pcoords[1] < 0) return;
1169 
1170 
1171  // Get the rank of the neighbor
1172  int prank;
1173  MPI_Cart_rank(comm,pcoords,&prank);
1174  MPI_Send(&pos, 1, MPI_INT, prank, 0, comm);
1175 
1176 }
1177 
1178 
1179 void recv_Back(int dir, int * pos){
1180 
1181  int pcoords[2];
1182 
1183  pcoords[0] = lcoords[0];
1184  pcoords[1] = lcoords[1];
1185 
1186  if(dir == 0)
1187  pcoords[0]++;
1188  else
1189  pcoords[1]++;
1190 
1191  if(pcoords[0] >= dims[0]) return;
1192  if(pcoords[1] >= dims[1]) return;
1193 
1194  // Get the rank of the neighbor
1195  int prank;
1196  MPI_Cart_rank(comm,pcoords,&prank);
1197  MPI_Recv(pos, 1, MPI_INT, prank, 0, comm, MPI_STATUS_IGNORE);
1198 
1199 }
1200 
1201 
1202 
1203 
1204 void phase_compose_sequence(){
1205 
1206  // Calculate the local part of the proteins that is being composed
1207  xBegin = lbegin[0] + lxBegin;
1208  xMax = lbegin[0] + lxMax;
1209  yBegin = lbegin[1] + lyBegin;
1210  yMax = lbegin[1] + lyMax;
1211 
1212 
1213  // Compose the global proteins
1216 
1217 
1218 }
1219 
1220 
1221 
1222 
1223 #define MIN(a,b) (((a)<(b))?(a):(b))
1224 #define MAX(a,b) (((a)>(b))?(a):(b))
1225 
1226 
1227 
1228 void composeOutput(aa_t * out, int * outpos, int * outsize, aa_t * lout, int loutsize, int * gBegin, int * gEnd){
1229 
1230 
1231  // MPI_Request for all the Sends.
1232  MPI_Request request[4];
1233 
1234  // Each processor sends to root its size.
1235  MPI_Isend(&loutsize, 1, MPI_INT, 0, 1, comm, &request[0]);
1236 
1237  // If the size is not 0, they also send the begin and end position of
1238  // the alignment and the array with the local alignment.
1239  if(loutsize > 0){
1240  MPI_Isend(gBegin, 1, MPI_INT, 0, 1, comm, &request[1]);
1241  MPI_Isend(gEnd, 1, MPI_INT, 0, 1, comm, &request[2]);
1242  MPI_Isend(lout, loutsize, MPI_AA, 0, 1, comm, &request[3]);
1243  }
1244 
1245  // Non-root processors exit after waiting the isends.
1246  if(rank != 0){
1247  MPI_Wait(&request[0],MPI_STATUS_IGNORE);
1248  if(loutsize > 0){
1249  MPI_Wait(&request[1],MPI_STATUS_IGNORE);
1250  MPI_Wait(&request[2],MPI_STATUS_IGNORE);
1251  MPI_Wait(&request[3],MPI_STATUS_IGNORE);
1252  }
1253  return;
1254  }
1255 
1256  // Reset the first index of the alignment (output)
1257  // and the size of the combined output (outsize).
1258  *outpos = size[0]+size[1];
1259  *outsize = 0;
1260 
1261  // Variables to keep the global indexes of the alignment.
1262  int first_flag = 1;
1263  int minBegin = 0;
1264  int maxEnd = 0;
1265 
1266 
1267  // Root receives from every processor
1268  int p;
1269  for(p=nProcs-1; p>=0; p--){
1270 
1271  // Receive the output local size of the p processor.
1272  int poutsize;
1273  MPI_Recv(&poutsize,1,MPI_INT,p,1,comm,MPI_STATUS_IGNORE);
1274 
1275  // Skip processors without alignment.
1276  if(!(poutsize > 0)) continue;
1277 
1278  // Get the begin and end
1279  int pgBegin;
1280  int pgEnd;
1281  MPI_Recv(&pgBegin,1,MPI_INT,p,1,comm,MPI_STATUS_IGNORE);
1282  MPI_Recv(&pgEnd,1,MPI_INT,p,1,comm,MPI_STATUS_IGNORE);
1283 
1284  // Update the begin and end indexes.
1285  if(first_flag){
1286  minBegin = pgBegin;
1287  maxEnd = pgEnd;
1288  first_flag = 0;
1289  } else {
1290  minBegin = MIN(minBegin,pgBegin);
1291  maxEnd = MAX(maxEnd,pgEnd);
1292  }
1293 
1294  // Update the outpos and outsize
1295  (*outpos) -= poutsize;
1296  (*outsize) += poutsize;
1297 
1298  // Receive the alignment.
1299  MPI_Recv(&out[(*outpos)],poutsize,MPI_AA,p,1,comm,MPI_STATUS_IGNORE);
1300 
1301  }
1302 
1303  // This is because the proteins start at 1.
1304  (*outpos)--;
1305 
1306 
1307  // Wait for the isends.
1308  // This could be not necessary since root uses blocking recv calls.
1309  MPI_Wait(&request[0],MPI_STATUS_IGNORE);
1310  if(loutsize > 0){
1311  MPI_Wait(&request[1],MPI_STATUS_IGNORE);
1312  MPI_Wait(&request[2],MPI_STATUS_IGNORE);
1313  MPI_Wait(&request[3],MPI_STATUS_IGNORE);
1314 
1315  }
1316 
1317  // Return the global Begin and End.
1318  *gBegin = minBegin;
1319  *gEnd = maxEnd;
1320 
1321 }
1322 
1323 
void phase_distribute_sequences()
Definition: SWpar.c:637
#define MAX(a, b)
Definition: SWpar_ref.c:1224
HitTile_h_t H
Definition: SWpar.c:243
HitTile_aa_t p1
Definition: SWpar.c:206
int lbegin[2]
Definition: SWpar_ref.c:191
#define MIN(a, b)
Definition: SWpar_ref.c:1223
double h_t
Definition: SWcommon.h:61
int match_pos
Definition: SWpar.c:227
int lyBegin
Definition: SWpar.c:262
#define T_SYNH
int iter
Definition: mmult_bit.c:68
h_t lMax
Definition: SWpar.c:254
void printProteinMatch(HitTile_aa_t p, int begin, int end)
Definition: SWcommon_hit.c:90
HitTile_aa_t out1
Definition: SWpar.c:210
void readProtein(ProteinFile *fprotein, HitTile_aa_t *protein, int psize)
Definition: SWcommon_hit.c:48
void debug_show_matrices()
Definition: SWpar.c:782
#define T_DIST
#define CHECK_NULL(check)
Definition: SWcommon.h:275
#define pam(a, b)
Definition: SWcommon.h:181
void phase_comp_back()
Definition: SWpar.c:955
void openProtein(ProteinFile *fprotein, char *filename)
Definition: SWcommon.c:241
HitTile_aa_t lout2
Definition: SWpar.c:240
#define BEGIN(coord, dims, size)
Definition: SWpar_ref.c:72
char * prot_name1
Definition: SWcommon.c:49
int yMax
Definition: SWpar.c:215
int lyMax
Definition: SWpar.c:258
int lcoords[2]
Definition: SWpar_ref.c:189
#define T_READ
int lxMax
Definition: SWpar.c:257
double gapPenalty
Definition: SWcommon.c:54
int rank
Definition: SWpar_ref.c:181
void recv_Back(int dir, int *pos)
Definition: SWpar.c:1100
#define T_HMAT
int nProcs
Definition: SWpar_ref.c:183
void phase_send_back()
Definition: SWpar.c:927
int lmatch_length2
Definition: SWpar.c:267
MPI_Datatype * types
Definition: refMPIluBack.c:109
int dims[2]
Definition: SWpar_ref.c:187
MPI_Comm comm
Definition: SWpar_ref.c:193
void input_parameters(int argc, char *argv[])
Definition: SWcommon.c:75
void phase_compose_sequence()
Definition: SWpar.c:1036
#define T_BACK
ProteinFile pfile1
Definition: SWcommon.c:63
#define xback(i, j)
Definition: SWpar_ref.c:248
#define T_SYNB
HitTile_trace_t yTraceback
Definition: SWpar.c:248
int lxBegin
Definition: SWpar.c:261
#define x
HitTile_aa_t lout1
Definition: SWpar.c:239
void initPAM(char *filename)
Definition: SWcommon.c:128
int size[2]
Definition: SWcommon.c:57
int match_length
Definition: SWpar.c:222
void printProtein(HitTile_aa_t p, int offset, int psize)
Definition: SWcommon_hit.c:78
#define h(i, j)
Definition: SWpar_ref.c:243
int match_length2
Definition: SWpar.c:224
void phase_recv_back()
Definition: SWpar.c:862
#define T_COMP
int lmatch_length1
Definition: SWpar.c:266
void debug_show_backtracking()
Definition: SWpar.c:1009
void composeOutput(HitTile_aa_t out, int *outpos, int *outsize, HitTile_aa_t lout, int loutpos, int loutsize, int *gBegin, int *gEnd)
Definition: SWpar.c:1123
void distributeProtein(HitTile_aa_t protein, HitTile_aa_t lprotein, int dim)
Definition: SWpar.c:1051
void phase_send_hmatrix()
Definition: SWpar.c:684
ProteinFile pfile2
Definition: SWcommon.c:64
HitTile_aa_t lp2
Definition: SWpar.c:236
#define T_LOCAL
short aa_t
Definition: SWcommon.h:59
HitTile_trace_t xTraceback
Definition: SWpar.c:247
int main(int argc, char *argv[])
Definition: cannonAsync.c:62
#define iterations
Definition: mg.c:119
void send_Back(int dir, int pos)
Definition: SWpar.c:1084
HitTile_aa_t p2
Definition: SWpar.c:207
char * pam_name
Definition: SWcommon.c:51
#define is_root()
Definition: SWcommon.h:123
int match_length1
Definition: SWpar.c:223
#define allTimers(timer)
HitTile_aa_t lp1
Definition: SWpar.c:235
void phase_recv_hmatrix()
Definition: SWpar.c:662
#define T_AVG
int lsize[2]
Definition: SWpar.c:188
int xMax
Definition: SWpar.c:214
short trace_t
Definition: SWcommon.h:63
int xBegin
Definition: SWpar.c:218
#define GAP_AA
Definition: SWcommon.h:243
#define printfroot(...)
Definition: SWcommon.h:135
HitTile_aa_t out2
Definition: SWpar.c:211
double currentTime()
Definition: SWcommon_ref.c:106
#define max(a, b)
Definition: cannonAsync.c:47
int lmatch_pos
Definition: SWpar.c:270
char * prot_name2
Definition: SWcommon.c:50
int yBegin
Definition: SWpar.c:219
void phase_read_sequences()
Definition: SWpar.c:619
#define SIZE(coord, dims, size)
Definition: SWpar_ref.c:82
int lmatch_length
Definition: SWpar.c:265
MPI_Op opHMaxLoc
Definition: SWpar.c:191
void phase_comp_hmatrix()
Definition: SWpar.c:708
#define yback(i, j)
Definition: SWpar_ref.c:249