Hitmap 1.3
 All Data Structures Namespaces Files Functions Variables Typedefs Friends Macros Groups Pages
SWpar.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 // Hitmap Include
60 #include <hitmap.h>
61 
62 // Include the common variables and utilities for all the versions
63 #include "SWcommon.h"
64 
65 // Include the common variables and utilities for the Hitmap versions
66 #include "SWcommon_hit.h"
67 
68 
69 
70 // Hitmap tile integer type
71 hit_tileNewType(int);
72 
73 
75 #define aa_hit_t MPI_SHORT
76 
77 #define h_hit_t MPI_DOUBLE
78 
79 
83 void distributeProtein(HitTile_aa_t protein, HitTile_aa_t lprotein, int dim);
84 
85 
90 hit_tileNewType(double_3int);
91 
92 
99 void recv_Back(int dir, int * pos);
100 
105 void send_Back(int dir, int pos);
106 
107 
111 void composeOutput(HitTile_aa_t out, int * outpos, int * outsize, HitTile_aa_t lout, int loutpos, int loutsize, int * gBegin, int * gEnd);
112 
113 
114 
115 
116 
117 // **********************************************************************
118 // Declarations of the different phases
119 // **********************************************************************
120 
124 void phase_read_sequences();
125 
130 
134 void phase_recv_hmatrix();
135 
139 void phase_comp_hmatrix();
140 
144 void phase_send_hmatrix();
145 
149 void phase_recv_back();
150 
154 void phase_comp_back();
155 
159 void phase_send_back();
160 
165 
166 
170 void debug_show_matrices();
171 
176 
177 
178 
179 
180 
181 
182 
183 // **********************************************************************
184 // Parallel related global variables
185 // **********************************************************************
186 
188 int lsize[2];
189 
191 MPI_Op opHMaxLoc;
192 
195 
198 
199 
200 
201 // **********************************************************************
202 // Global root variables
203 // **********************************************************************
204 
205 // Protein sequences
206 HitTile_aa_t p1;
207 HitTile_aa_t p2;
208 
209 // Output proteins with the alignment
210 HitTile_aa_t out1;
211 HitTile_aa_t out2;
212 
213 // X and Y coordinates of the maximum value
214 int xMax;
215 int yMax;
216 
217 // Begin coordinates of the match
218 int xBegin;
219 int yBegin;
220 
221 // Match lengths (total, sequence 1, sequence 2)
225 
226 // Output match begin
228 
229 
230 // **********************************************************************
231 // Global (C) local (parallel) variables
232 // **********************************************************************
233 
234 // Protein sequences
235 HitTile_aa_t lp1;
236 HitTile_aa_t lp2;
237 
238 // Local output proteins with the alignment
239 HitTile_aa_t lout1;
240 HitTile_aa_t lout2;
241 
242 // Similarity matrix
243 HitTile_h_t H;
244 #define h(i,j) (hit_tileElemAtNoStride(H,2,i,j))
245 
246 // Traceback matrices
247 HitTile_trace_t xTraceback;
248 HitTile_trace_t yTraceback;
249 #define xback(i,j) (hit_tileElemAtNoStride(xTraceback,2,i,j))
250 #define yback(i,j) (hit_tileElemAtNoStride(yTraceback,2,i,j))
251 
252 
253 // Local Max value
255 
256 // X and Y coordinates of the maximum value
257 int lxMax;
258 int lyMax;
259 
260 // Begin coordinates of the match
263 
264 // Match lengths (total, sequence 1, sequence 2)
268 
269 // Output match begin
271 
272 
274 
275 
276 
280 int main(int argc, char *argv[]) {
281 
282  // **********************************************************************
283  // Hitmap init
284  // **********************************************************************
285  hit_comInit( &argc, &argv );
286 
287 
288  // **********************************************************************
289  // HIT datatypes
290  // **********************************************************************
291 
292  hit_comTypeStruct(&MPI_DOUBLE_3INT, double_3int, 4,
293  val, 1, HIT_DOUBLE,
294  rank, 1, HIT_INT,
295  i, 1, HIT_INT,
296  j, 1, HIT_INT);
297 
298 
299  // Create the operation to reduce the double_3int values.
300  MPI_Op_create( HMaxLoc, 1, &opHMaxLoc );
301 
302  // **********************************************************************
303  // Check input parameters
304  // **********************************************************************
305  input_parameters(argc, argv);
306 
307  // **********************************************************************
308  // Layout
309  // **********************************************************************
310 
311  // Topology
312  HitTopology topo = hit_topology( plug_topArray2DComplete );
313 
314  // Print processor info
315  printfroot("=== Parallel ===\n"
316  " Processors: %d\n"
317  " Processors: %d x %d\n\n",hit_NProcs,hit_topDimCard(topo,0), hit_topDimCard(topo,1));
318 
319  // Layout
320  HitShape shape = hit_shape(2,hit_sig(1,size[0],1),hit_sig(1,size[1],1));
321  layout = hit_layout( plug_layBlocks, topo, shape );
322 
323  // Shape from the layout
324  layshape = hit_layShape(layout);
325  lsize[0] = hit_shapeSigCard(layshape,0);
326  lsize[1] = hit_shapeSigCard(layshape,1);
327 
328  // Local shape to alloc the H matriz
329  lshape = hit_shapeDimExpand(layshape,0,HIT_SHAPE_BEGIN,-1);
330  lshape = hit_shapeDimExpand(lshape ,1,HIT_SHAPE_BEGIN,-1);
331 
332  // Shapes for the proteins
333  HitShape protshape[2];
334  protshape[0] = hit_shape(1,hit_shapeSig(lshape,0));
335  protshape[1] = hit_shape(1,hit_shapeSig(lshape,1));
336 
337  // Print info of each processor
338 #ifdef DEBUG
339  ALL_SEQUENTIAL_BEGIN(seq1)
340  printf("Processor: %d/%d, ",hit_Rank,hit_NProcs);
341  printf("coords %d x %d\n", hit_laySelfRanksDim(layout,0), hit_laySelfRanksDim(layout,1));
342  printf("Dim1: Block %2d- (%2d)\n", hit_shapeSig(layshape,0).begin , hit_shapeSigCard(layshape,0));
343  printf("Dim2: Block %2d- (%2d)\n", hit_shapeSig(layshape,1).begin , hit_shapeSigCard(layshape,1));
344  printf("\n");
345  ALL_SEQUENTIAL_END(seq1)
346 #endif
347 
348 
349  // **********************************************************************
350  // Timers
351  // **********************************************************************
352 
353  HitClock time_init;
354  HitClock time_comp;
356  HitClock time_comp_read;
357  HitClock time_comp_hmatrix;
358  HitClock time_comp_back;
360  HitClock time_comp_distr;
361  HitClock time_comp_syncH;
362  HitClock time_comp_syncB;
363  HitClock time_comp_compose;
365  hit_clockReset(time_comp_read);
366  hit_clockReset(time_comp_hmatrix);
367  hit_clockReset(time_comp_back);
368 
369  hit_clockReset(time_comp_distr);
370  hit_clockReset(time_comp_syncH);
371  hit_clockReset(time_comp_syncB);
372  hit_clockReset(time_comp_compose);
373 
374 
375  // **********************************************************************
376  // INIT PHASE: Structure allocation and PAM matrix load
377  // **********************************************************************
378 
379  // Init time
381  hit_clockStart(time_init);
382 
383  // read the PAM matrix
384  // Every one do this, the pam matrix is small and all procs
385  // need it, we don't bother to read it and distribute it.
386  initPAM(pam_name);
387 
388  // Open the proteins files
391 
392  // Allocate memory for the proteins (Only root has them)
393  if(is_root()){
394  hit_tileDomainAlloc(&p1,aa_t,1,size[0]+1);
395  hit_tileDomainAlloc(&p2,aa_t,1,size[1]+1);
396  }
397 
398  // Local proteins
399  hit_tileDomainShapeAlloc(&lp1,aa_t,protshape[0]);
400  hit_tileDomainShapeAlloc(&lp2,aa_t,protshape[1]);
401 
402  // Local Similarity matrix
403  hit_tileDomainShapeAlloc(&H,h_t,lshape);
404 
405  // Traceback matrices
408 
409  // Allocate size of the output (local)
410  hit_tileDomainAlloc(&lout1, aa_t, 1 ,lsize[0]+lsize[1]+1 );
411  hit_tileDomainAlloc(&lout2, aa_t, 1 ,lsize[0]+lsize[1]+1 );
412 
413 
414  // Allocate size of the output (global, only root)
415  if(is_root()){
416  hit_tileDomainAlloc(&out1, aa_t, 1 ,size[0]+size[1]+1 );
417  hit_tileDomainAlloc(&out2, aa_t, 1 ,size[0]+size[1]+1 );
418  }
419 
420  // Init time
421  hit_clockStop(time_init);
422 
423 
424 
425  // **********************************************************************
426  // WHOLE COMPUTATION PHASE: iteration of: Read, H matrix, Backtracking
427  // **********************************************************************
428 
429  // Computation time
430  hit_clockStart(time_comp);
431 
432  int iter;
433  for(iter=0;iter<iterations;iter++){
434 
435 
436 #ifdef DEBUG
437  printfroot("=== Iteration: %d/%d ===\n\n",iter,iterations);
438 #endif
439 
440  // **********************************************************************
441  // READ PHASE
442  // **********************************************************************
443  hit_clockContinue(time_comp_read);
445  hit_clockStop(time_comp_read);
446 
447  // **********************************************************************
448  // DISTRIBUTION PHASE
449  // **********************************************************************
450  hit_clockContinue(time_comp_distr);
452  hit_clockStop(time_comp_distr);
453 
454 
455  // **********************************************************************
456  // RECV H PHASE
457  // **********************************************************************
458  hit_clockContinue(time_comp_syncH);
459  // Recv the needed H values from the neighbors
461  hit_clockStop(time_comp_syncH);
462 
463 
464  // **********************************************************************
465  // H MATRIX PHASE
466  // **********************************************************************
467  hit_clockContinue(time_comp_hmatrix);
469  hit_clockStop(time_comp_hmatrix);
470 
471 
472  // **********************************************************************
473  // SEND H PHASE
474  // **********************************************************************
475  hit_clockContinue(time_comp_syncH);
476  // Send the needed H values to the neighbors
478  hit_clockStop(time_comp_syncH);
479 
480 #ifdef DEBUG
482 #endif
483 
484 
485  // **********************************************************************
486  // RECV BACKTRAING PHASE
487  // **********************************************************************
488  hit_clockContinue(time_comp_syncB);
489  phase_recv_back();
490  hit_clockStop(time_comp_syncB);
491 
492 
493  // **********************************************************************
494  // BACKTRACKING PHASE
495  // **********************************************************************
496  hit_clockContinue(time_comp_back);
497  phase_comp_back();
498  hit_clockStop(time_comp_back);
499 
500 
501  // **********************************************************************
502  // SEND BACKTRAING PHASE
503  // **********************************************************************
504  hit_clockContinue(time_comp_syncB);
505  phase_send_back();
506  hit_clockStop(time_comp_syncB);
507 
508 
509  // **********************************************************************
510  // COMPOSE FINAL PROTEIN PHASE
511  // **********************************************************************
512  hit_clockContinue(time_comp_compose);
514  hit_clockStop(time_comp_compose);
515 
516 
517 #ifdef DEBUG
519 #endif
520 
521  } // Iterations Loop
522 
523  // Computation time
525  hit_clockStop(time_comp);
526 
527 
528  // **********************************************************************
529  // Show results
530  // **********************************************************************
531 
532  // Reduce the timers
533  hit_clockReduce(layout,time_comp_read);
534  hit_clockReduce(layout,time_comp_hmatrix);
535  hit_clockReduce(layout,time_comp_back);
536 
537  hit_clockReduce(layout,time_comp_distr);
538  hit_clockReduce(layout,time_comp_syncH);
539  hit_clockReduce(layout,time_comp_syncB);
540  hit_clockReduce(layout,time_comp_compose);
541 
542 
543 #define allTimers(timer) hit_clockGetMinSeconds(timer),hit_clockGetAvgSeconds(timer),hit_clockGetMaxSeconds(timer)
544 
545 
546 
547  // Final result.
548  if(is_root()){
549 
552 
553  printf("=== Result ===\n");
554  printf("Init time: %f\n",hit_clockGetSeconds(time_init));
555  printf("Comp time: %f\n",hit_clockGetSeconds(time_comp));
556  printf(" Read protein: min: %f, avg: %f, max: %f \n",allTimers(time_comp_read));
557  printf(" Distribute: min: %f, avg: %f, max: %f \n",allTimers(time_comp_distr));
558  printf(" H matrix: min: %f, avg: %f, max: %f \n",allTimers(time_comp_hmatrix));
559  printf(" sync H mat: min: %f, avg: %f, max: %f \n",allTimers(time_comp_syncH));
560  printf(" Backtracking: min: %f, avg: %f, max: %f \n",allTimers(time_comp_back));
561  printf(" Sync Back: min: %f, avg: %f, max: %f \n",allTimers(time_comp_syncB));
562  printf(" Compose: min: %f, avg: %f, max: %f \n",allTimers(time_comp_compose));
563 
564  printf("Last Match length of A: %d\n", match_length1);
565  printf("Last Match length of B: %d\n", match_length2);
566  printf("Last Alignment length: %d\n", match_length);
567 
568 
569  // Print the result for small inputsets
570  if(match_length < 50){
571  printf("\n");
572  printf("MCH1: "); printProteinMatch(p1,xBegin,xMax);
573  printf("MCH2: "); printProteinMatch(p2,yBegin,yMax);
574 
575  printf("\n");
576  printf("OUT1: "); printProtein(out1,match_pos,match_length);
577  printf("OUT2: "); printProtein(out2,match_pos,match_length);
578  } else {
579  printf("\nInputset too big to print to stdout\n");
580  }
581  }
582 
583 
584  // @todo
585  // Free all the resources
586  hit_tileFree(lp1);
587  hit_tileFree(lp2);
590  hit_tileFree(H);
593  MPI_Op_free(&opHMaxLoc);
594  MPI_Type_free(&MPI_DOUBLE_3INT);
595 
596  hit_layFree(layout);
597  hit_topFree(topo);
598 
599  if(is_root()){
600  hit_tileFree(p1);
601  hit_tileFree(p2);
604  }
605 
606 
607  // Finalize Hitmap
608  hit_comFinalize();
609 
610  return EXIT_SUCCESS;
611 
612 
613 }
614 
615 
616 
617 
618 
620 
621  if(! is_root() ) return;
622 
623  // Read the proteins
624  readProtein(&pfile1, &p1, size[0]);
625  readProtein(&pfile2, &p2, size[1]);
626 
627 #ifdef DEBUG
628  printf("PRT1: "); printProtein(p1,0,size[0]);
629  printf("PRT2: "); printProtein(p2,0,size[1]);
630  printf("\n");
631 #endif
632 
633 
634 }
635 
636 
638 
639  // Send and Recv the proteins from the root proc.
642 
643  // Print info of each processor
644 #ifdef DEBUG
645  ALL_SEQUENTIAL_BEGIN(seq2)
646  printf("Processor: %d/%d, ",hit_Rank,hit_NProcs);
647  printf("Coords %d x %d\n", hit_laySelfRanksDim(layout,0), hit_laySelfRanksDim(layout,1));
648  printf("Local protein 1: ");
649  printProtein(lp1,0,lsize[0]);
650  printf("Local protein 2: ");
651  printProtein(lp2,0,lsize[1]);
652  printf("\n");
653  ALL_SEQUENTIAL_END(seq2)
654 #endif
655 
656 }
657 
658 
659 
660 
661 
663 
664  HitPattern pattern = hit_pattern( HIT_PAT_UNORDERED );
665 
666  HitShape selectX = hit_shape(2, hit_sigIndex(0), hit_sig(1,hit_tileDimCard(H,1)-1,1) );
667  HitRanks neighX = hit_layNeighbor(layout, 0, -1);
668  HitCom commX = hit_comRecvSelect(layout,neighX,&H,selectX,HIT_COM_TILECOORDS, h_hit_t);
669  hit_patternAdd(&pattern,commX);
670 
671 
672  HitShape selectY = hit_shape(2,hit_sig(0,hit_tileDimCard(H,0)-1,1),hit_sigIndex(0));
673  HitRanks neighY = hit_layNeighbor(layout, 1, -1);
674  HitCom commY = hit_comRecvSelect(layout,neighY,&H,selectY,HIT_COM_TILECOORDS, h_hit_t);
675  hit_patternAdd(&pattern,commY);
676 
677  hit_patternDo(pattern);
678  hit_patternFree(&pattern);
679 
680 }
681 
682 
683 
685 
686 
687  HitPattern pattern = hit_pattern( HIT_PAT_UNORDERED );
688 
689  HitShape selectX = hit_shape(2, hit_sigIndex(hit_tileDimCard(H,0)-1), hit_sig(1,hit_tileDimCard(H,1)-1,1 ) );
690  HitRanks neighX = hit_layNeighbor(layout, 0, 1);
691  HitCom commX = hit_comSendSelect(layout,neighX,&H,selectX,HIT_COM_TILECOORDS, h_hit_t);
692  hit_patternAdd(&pattern,commX);
693 
694 
695  HitShape selectY = hit_shape(2,hit_sig(0,hit_tileDimCard(H,0)-1,1 ), hit_sigIndex(hit_tileDimCard(H,1)-1));
696  HitRanks neighY = hit_layNeighbor(layout, 1, 1);
697  HitCom commY = hit_comSendSelect(layout,neighY,&H,selectY,HIT_COM_TILECOORDS, h_hit_t);
698  hit_patternAdd(&pattern,commY);
699 
700  hit_patternDo(pattern);
701  hit_patternFree(&pattern);
702 
703 
704 }
705 
706 
707 
709 
710  // Compute the Similarity Matrix (H)
711 
712  // Initialize values
713  int i,j;
714 
715  if(hit_laySelfRanks(layout).rank[1] == 0){
716  for (i=0;i<=lsize[0];i++){
717  xback(i,0) = -1;
718  yback(i,0) = -1;
719  h(i,0)=0;
720  }
721  }
722 
723  if(hit_laySelfRanks(layout).rank[0] == 0){
724  for (j=0;j<=lsize[1];j++){
725  xback(0,j) = -1;
726  yback(0,j) = -1;
727  h(0,j)=0;
728  }
729  }
730 
731  // Maximum values to retrieve the alignment
732  lMax = -1;
733  xMax = -1;
734  yMax = -1;
735 
736  // Perform wavefront algorithm
737  for(i=1; i<lsize[0]+1; i++){
738  for(j=1; j<lsize[1]+1; j++){
739 
740  // Get the scores
741  h_t diag = h(i-1,j-1) + pam(hit_tileElemAt(lp1,1,i),hit_tileElemAt(lp2,1,j));
742  h_t down = h(i-1,j ) + gapPenalty;
743  h_t right = h(i ,j-1) + gapPenalty;
744 
745  // Calculate the maximum
746  int idx;
747  h_t max = MAX4(diag,down,right,0,&idx);
748  h(i,j) = max;
749 
750  // Set the back trace variables.
751  if ( idx == 0 ){ // max == diag
752  xback(i,j) = (trace_t) (i-1);
753  yback(i,j) = (trace_t) (j-1);
754  } else if ( idx == 1 ) { // max == down
755  xback(i,j) = (trace_t) (i-1);
756  yback(i,j) = (trace_t) j;
757  } else if ( idx == 2 ) { // max == right
758  xback(i,j) = (trace_t) i;
759  yback(i,j) = (trace_t) (j-1);
760  } else { // max == 0
761  xback(i,j) = -1;
762  yback(i,j) = -1;
763  }
764 
765  // Maximum
766  if(max >= lMax){
767  lxMax = i;
768  lyMax = j;
769  lMax = max;
770  }
771 
772  } // End loop j;
773  } // End loop i;
774 
775 
776 }
777 
778 
779 
780 
781 
783 
784 
785 
786  // Print info of each processor
787  #ifdef DEBUG
788  int i,j;
789 
790  ALL_SEQUENTIAL_BEGIN(seq3)
791 
792  printf("Processor: %d/%d, ",hit_Rank,hit_NProcs);
793  printf("coords %d x %d\n", hit_laySelfRanksDim(layout,0), hit_laySelfRanksDim(layout,1));
794 
795  // Print
796  printf(" ");
797  for (j=0;j<=lsize[1];j++){
798  if(j==0)
799  printf(" -");
800  else
801  printf("%6c",AA2char(hit_tileElemAt(lp2,1,j)));
802  }
803  printf("\n");
804  for (i=0;i<=lsize[0];i++){
805 
806  if(i==0)
807  printf("- |");
808  else
809  printf("%c |",AA2char(hit_tileElemAt(lp1,1,i)));
810 
811  for (j=0;j<=lsize[1];j++){
812  printf("%6.1f",(double)h(i,j));
813  }
814  printf("\n");
815  }
816  printf("\n");
817 
818  ALL_SEQUENTIAL_END(seq3)
819  #endif
820 
821 
822 
823  #ifdef DEBUG
824 
825  ALL_SEQUENTIAL_BEGIN(seq4)
826 
827  printf("Processor: %d/%d, ",hit_Rank,hit_NProcs);
828  printf("coords %d x %d\n", hit_laySelfRanksDim(layout,0), hit_laySelfRanksDim(layout,1));
829 
830  // Print
831  printf(" ");
832  for (j=0;j<=lsize[1];j++){
833  if(j==0)
834  printf(" -");
835  else
836  printf("%6c",AA2char(hit_tileElemAt(lp2,1,j)));
837  }
838  printf("\n");
839  for (i=0;i<=lsize[0];i++){
840 
841  if(i==0)
842  printf("- |");
843  else
844  printf("%c |",AA2char(hit_tileElemAt(lp1,1,i)));
845 
846  for (j=0;j<=lsize[1];j++){
847  printf("(%2d,%2d) ",xback(i,j),yback(i,j));
848  }
849  printf("\n");
850  }
851  printf("\n");
852 
853  ALL_SEQUENTIAL_END(seq4)
854 
855  #endif
856 
857 
858 }
859 
860 
861 
863 
864  // Determine who has the maximum value
865  HitTile_double_3int Max, outMax;
866  hit_tileDomainAlloc(&Max,double_3int,1,1);
867  hit_tileDomainAlloc(&outMax,double_3int,1,1);
868 
869  // Fill Max with our values.
870  hit_tileElemAt(Max,1,0).val = lMax;
871  hit_tileElemAt(Max,1,0).rank = hit_Rank;
872  hit_tileElemAt(Max,1,0).i = hit_shapeSig(layshape,0).begin + xMax;
873  hit_tileElemAt(Max,1,0).j = hit_shapeSig(layshape,1).begin + yMax;
874 
875  // Reduce the values.
876  HitCom comm_reduce = hit_comReduce(layout,HIT_RANKS_NULL,&Max,&outMax,MPI_DOUBLE_3INT,opHMaxLoc);
877  hit_comDo(&comm_reduce);
878  hit_comFree(comm_reduce);
879 
880 #ifdef DEBUG
881  // Print the maximum
882  printfroot("Max %f in proc: %d\n\n",hit_tileElemAt(outMax,1,0).val, hit_tileElemAt(outMax,1,0).rank);
883 #endif
884 
885  // Receive the point to continue the backtracking algorithm.
886  if(hit_Rank == hit_tileElemAt(outMax,1,0).rank){
887 
888  // Avoid possible deadlock, if we are the processor with
889  // the maximum value, we receive the messages but we ignore them.
890  int pos;
891  recv_Back(0, &pos);
892  recv_Back(1, &pos);
893 
894  // Get the position from a neighbor
895  } else {
896 
897  // By default we do not have the path.
898  lxMax=-1;
899  lyMax=-1;
900 
901  // Receive the path from the neigh from dim 0
902  int pos = -1;
903  recv_Back(0, &pos);
904 
905  if(pos != -1){
906  lxMax = lsize[0];
907  lyMax = pos;
908  }
909 
910  // Receive the path from the neigh from dim 1
911  pos = -1;
912  recv_Back(1, &pos);
913 
914  if(pos != -1){
915  lxMax = pos;
916  lyMax = lsize[1];
917  }
918  }
919 
920 
921  hit_tileFree(Max);
922  hit_tileFree(outMax);
923 
924 }
925 
926 
928 
929  int i = lxBegin;
930  int j = lyBegin;
931 
932  // Send a message to the neighbors to continue or stop
933  // backtracking algorithm
934  if( i == -1 || j == -1 ){
935  // The current processors do not have any point of the path of the alignment
936  send_Back(0, -1);
937  send_Back(1, -1);
938  } else if(h(i,j) > 0 && i == 0 ){
939  // Neighbor in dim 0 should continue the algorithm.
940  send_Back(0, j);
941  send_Back(1, -1);
942  } else if (h(i,j) > 0 && j == 0 ){
943  // Neighbor in dim 1 should continue the algorithm.
944  send_Back(0, -1);
945  send_Back(1, i);
946  } else {
947  // Neighbors do not continue the search.
948  send_Back(0, -1);
949  send_Back(1, -1);
950  }
951 
952 }
953 
954 
956 
957  // reset to max point to do alignment
958  // We start at the maximum point
959  int i=lxMax;
960  int j=lyMax;
961 
962  // Position of the output proteins.
963  int pos=lsize[0]+lsize[1]+1;
964  int outsize=0;
965 
966  // Backtracking algorithm
967  while (i>0 && j>0 && h(i,j) > 0){
968 
969  int movex = (i>xback(i,j));
970  int movey = (j>yback(i,j));
971 
972  pos--;
973  outsize++;
974 
975  if (movex && movey) {
978  } else if (movey) {
979  hit_tileElemAt(lout1,1,pos)=GAP_AA;
981  } else if (movex) {
983  hit_tileElemAt(lout2,1,pos)=GAP_AA;
984  }
985 
986  // Move to the next point.
987  int tempi=i;
988  int tempj=j;
989  i=xback(tempi,tempj);
990  j=yback(tempi,tempj);
991  }
992 
993  pos--;
994 
995  lmatch_pos = pos;
996 
997  lmatch_length = outsize;
998  lmatch_length1 = lxMax - i;
999  lmatch_length2 = lyMax - j;
1000 
1001  lxBegin = i;
1002  lyBegin = j;
1003 
1004 }
1005 
1006 
1007 
1008 
1010 
1011 #ifdef DEBUG
1012 
1013  // Local result print
1014 ALL_SEQUENTIAL_BEGIN(seq5);
1015  printf("Processor: %d/%d, ",hit_Rank,hit_NProcs);
1016  printf("Coords %d x %d\n", hit_laySelfRanksDim(layout,0), hit_laySelfRanksDim(layout,1));
1017 
1018  printf(" First point (%d,%d) -> Last point (%d,%d)\n",lxMax,lyMax,lxBegin,lyBegin);
1019 
1020  printf(" Local Match length of A: %d\n", lmatch_length1);
1021  printf(" Local Match length of B: %d\n", lmatch_length2);
1022  printf(" Local Alignment length: %d\n", lmatch_length);
1023 
1024  printf(" MCH1: "); printProteinMatch(lp1,lxBegin,lxMax);
1025  printf(" MCH2: "); printProteinMatch(lp2,lyBegin,lyMax);
1026 
1027  printf(" LOCAL OUT1: "); printProtein(lout1,lmatch_pos,lmatch_length);
1028  printf(" LOCAL OUT2: "); printProtein(lout2,lmatch_pos,lmatch_length);
1029  printf("\n");
1030 ALL_SEQUENTIAL_END(seq5);
1031 #endif
1032 
1033 
1034 }
1035 
1037 
1038  // Calculate the local part of the proteins that is being composed
1039  xBegin =hit_shapeSig(lshape,0).begin + lxBegin;
1040  xMax = hit_shapeSig(lshape,0).begin + lxMax;
1041  yBegin = hit_shapeSig(lshape,1).begin + lyBegin;
1042  yMax = hit_shapeSig(lshape,1).begin + lyMax;
1043 
1044  // Compose the global proteins
1047 
1048 }
1049 
1050 
1051 void distributeProtein(HitTile_aa_t protein, HitTile_aa_t lprotein, int dim){
1052 
1053 
1054  HitPattern pattern = hit_pattern( HIT_PAT_UNORDERED );
1055 
1056  if(hit_Rank == 0){
1057  int p;
1058  for(p=0; p<hit_layNumActives(layout); p++){
1059 
1060  HitRanks ranks = hit_layActiveIdRanks(layout, p);
1061  HitShape matrixshape = hit_layShapeOther(layout,ranks);
1062  HitShape shape = hit_shape(1,hit_shapeSig(matrixshape,dim));
1063 
1064  HitCom comm = hit_comSendSelect(layout,ranks,&protein,shape,HIT_COM_ARRAYCOORDS, aa_hit_t);
1065  hit_patternAdd(&pattern,comm);
1066  }
1067  }
1068 
1069  // @todo Use the one precalcualted
1070  HitShape lprotein_shape = hit_shape(1,hit_shapeSig(hit_layShape(layout),dim));
1071 
1072 
1073  HitCom comm = hit_comRecvSelect(layout, hit_ranks2(0,0), &lprotein, lprotein_shape, HIT_COM_ARRAYCOORDS, aa_hit_t);
1074  hit_patternAdd(&pattern,comm);
1075  hit_patternDo(pattern);
1076 
1077  hit_patternFree(&pattern);
1078 
1079 }
1080 
1081 
1082 
1083 
1084 void send_Back(int dir, int pos){
1085 
1086  HitTile_int hitpos;
1087  hit_tileDomainAlloc(&hitpos,int,1,1);
1088  hit_tileElemAt(hitpos,1,0) = pos;
1089  HitRanks neigh = hit_layNeighbor(layout, dir, -1);
1090  HitCom comm = hit_comSend(layout,neigh,&hitpos,HIT_INT);
1091  hit_comDo(&comm);
1092 
1093  hit_tileFree(hitpos);
1094  hit_comFree(comm);
1095 
1096 }
1097 
1098 
1099 
1100 void recv_Back(int dir, int * pos){
1101 
1102  HitTile_int hitpos;
1103  hit_tileDomainAlloc(&hitpos,int,1,1);
1104  hit_tileElemAt(hitpos,1,0) = -1;
1105  HitRanks neigh = hit_layNeighbor(layout, dir, 1);
1106  HitCom comm = hit_comRecv(layout,neigh,&hitpos,HIT_INT);
1107  hit_comDo(&comm);
1108  *pos = hit_tileElemAt(hitpos,1,0);
1109 
1110  hit_tileFree(hitpos);
1111  hit_comFree(comm);
1112 
1113 }
1114 
1115 
1116 
1117 
1118 #define MIN(a,b) (((a)<(b))?(a):(b))
1119 #define MAX(a,b) (((a)>(b))?(a):(b))
1120 
1121 
1122 
1123 void composeOutput(HitTile_aa_t out, int * outpos, int * outsize, HitTile_aa_t lout, int loutpos, int loutsize, int * gBegin, int * gEnd){
1124 
1125  // Send the local out size
1126  HitTile_int hitdata;
1127  hit_tileDomainAlloc(&hitdata,int,1,3);
1128  hit_tileElemAt(hitdata,1,0) = loutsize;
1129  hit_tileElemAt(hitdata,1,1) = *gBegin;
1130  hit_tileElemAt(hitdata,1,2) = *gEnd;
1131 
1132  // Each processor sends to root its size.
1133  HitCom comm1 = hit_comSend(layout,hit_ranks2(0,0),&hitdata,HIT_INT);
1134  hit_comStartSend(&comm1);
1135 
1136  HitCom comm2;
1137 
1138 
1139  // If the size is not 0, they also send the begin and end position of
1140  // the alignment and the array with the local alignment.
1141  if(loutsize > 0){
1142 
1143  HitShape selection = hit_shape(1,hit_sig(loutpos,loutpos+loutsize-1,1));
1144 
1145  comm2 = hit_comSendSelect(layout,hit_ranks2(0,0),&lout,selection,HIT_COM_TILECOORDS,aa_hit_t);
1146  hit_comStartSend(&comm2);
1147 
1148  }
1149 
1150 
1151  // Non-root processors exit after waiting the isends.
1152  if(hit_Rank != 0){
1153  hit_comEndSend(&comm1);
1154  if(loutsize > 0){
1155  hit_comEndSend(&comm2);
1156  hit_comFree(comm2);
1157  }
1158  hit_comFree(comm1);
1159  hit_tileFree(hitdata);
1160  return;
1161  }
1162 
1163 
1164  // Reset the first index of the alignment (output)
1165  // and the size of the combined output (outsize).
1166  *outpos = size[0]+size[1];
1167  *outsize = 0;
1168 
1169 
1170  // Variables to keep the global indexes of the alignment.
1171  int first_flag = 1;
1172  int minBegin = 0;
1173  int maxEnd = 0;
1174 
1175  // Recv data
1176  HitTile_int hitdataRecv;
1177  hit_tileDomainAlloc(&hitdataRecv,int,1,3);
1178 
1179  // Root receives from every processor
1180  int p;
1181  for(p=hit_NProcs-1; p>=0; p--){
1182 
1183  // Receive the output local size of the p processor.
1184  HitCom comm = hit_comRecv(layout,hit_layActiveIdRanks(layout,p),&hitdataRecv,HIT_INT);
1185  hit_comDo(&comm);
1186  hit_comFree(comm);
1187 
1188  // Each processor sends to root its size.
1189  int poutsize = hit_tileElemAt(hitdataRecv,1,0);
1190 
1191  // Get the begin and end
1192  int pgBegin = hit_tileElemAt(hitdataRecv,1,1);
1193  int pgEnd = hit_tileElemAt(hitdataRecv,1,2);
1194 
1195 
1196 
1197  // Skip processors without alignment.
1198  if(!(poutsize > 0)) continue;
1199 
1200  // Update the begin and end indexes.
1201  if(first_flag){
1202  minBegin = pgBegin;
1203  maxEnd = pgEnd;
1204  first_flag = 0;
1205  } else {
1206  minBegin = MIN(minBegin,pgBegin);
1207  maxEnd = MAX(maxEnd,pgEnd);
1208  }
1209 
1210  // Update the outpos and outsize
1211  (*outpos) -= poutsize;
1212  (*outsize) += poutsize;
1213 
1214  // Receive the alignment.
1215  HitShape selection = hit_shape(1,hit_sig((*outpos),(*outpos)+poutsize-1,1));
1216 
1217  HitCom commB = hit_comRecvSelect(layout,hit_layActiveIdRanks(layout,p),&out,selection,HIT_COM_TILECOORDS,aa_hit_t);
1218  hit_comDo(&commB);
1219  hit_comFree(commB);
1220 
1221  }
1222 
1223  hit_tileFree(hitdataRecv);
1224 
1225  // This is because the proteins start at 1.
1226  (*outpos)--;
1227 
1228  // Return the global Begin and End.
1229  *gBegin = minBegin;
1230  *gEnd = maxEnd;
1231 
1232  hit_comEndSend(&comm1);
1233  hit_comFree(comm1);
1234 
1235  if(loutsize > 0){
1236  hit_comEndSend(&comm2);
1237  hit_comFree(comm2);
1238 
1239  }
1240 
1241  hit_tileFree(hitdata);
1242 
1243 }
1244 
1245 
1246 
1247 
void phase_distribute_sequences()
Definition: SWpar.c:637
HitTile_h_t H
Definition: SWpar.c:243
#define hit_shape(nd,...)
Definition: hit_sshape.h:175
HitTile_aa_t p1
Definition: SWpar.c:206
#define hit_layShape(lay)
Definition: hit_layout.h:650
#define hit_tileNewType(baseType)
Definition: hit_tile.h:163
HitShape layshape
Definition: SWpar.c:197
double h_t
Definition: SWcommon.h:61
int match_pos
Definition: SWpar.c:227
int lyBegin
Definition: SWpar.c:262
#define hit_clockSynchronizeAll()
Definition: hit_utils.h:55
#define HIT_INT
Definition: hit_com.h:74
#define hit_comSendSelect(lay, sendTo, tileP, selection, modeSelect, baseType)
Definition: hit_com.h:442
#define hit_Rank
Definition: hit_com.h:140
int iter
Definition: mmult_bit.c:68
#define hit_topDimCard(topo, i)
Definition: hit_topology.h:413
#define hit_tileElemAt(var, ndims,...)
Definition: hit_tile.h:519
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
#define hit_laySelfRanks(lay)
Definition: hit_layout.h:848
void readProtein(ProteinFile *fprotein, HitTile_aa_t *protein, int psize)
Definition: SWcommon_hit.c:48
void debug_show_matrices()
Definition: SWpar.c:782
void hit_comFree(HitCom issue)
Definition: hit_com.c:1995
void hit_comDo(HitCom *issue)
Definition: hit_com.c:2408
#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 hit_clockContinue(c)
Definition: hit_utils.h:121
char * prot_name1
Definition: SWcommon.c:49
#define xback(i, j)
Definition: SWpar.c:249
int yMax
Definition: SWpar.c:215
#define hit_tileDomainAlloc(newVarP, baseType, numDims,...)
Definition: hit_tile.h:336
HitLayout layout[k_num]
Definition: mg.c:150
int lyMax
Definition: SWpar.c:258
#define hit_topology(name,...)
Definition: hit_topology.h:308
#define hit_layShapeOther(lay, ranks)
Definition: hit_layout.h:732
HitRanks HIT_RANKS_NULL
Definition: hit_topology.c:68
int lxMax
Definition: SWpar.c:257
double gapPenalty
Definition: SWcommon.c:54
void hit_comEndSend(HitCom *issue)
Definition: hit_com.c:2062
int rank
Definition: SWpar_ref.c:181
void recv_Back(int dir, int *pos)
Definition: SWpar.c:1100
#define MAX(a, b)
Definition: SWpar.c:1119
void phase_send_back()
Definition: SWpar.c:927
#define HIT_COM_ARRAYCOORDS
Definition: hit_com.h:345
#define hit_laySelfRanksDim(lay, dim)
Definition: hit_layout.h:857
int lmatch_length2
Definition: SWpar.c:267
#define HIT_PAT_UNORDERED
Definition: hit_pattern.h:81
void hit_patternAdd(HitPattern *pattern, HitCom comm)
Definition: hit_pattern.c:61
MPI_Comm comm
Definition: SWpar_ref.c:193
#define hit_NProcs
Definition: hit_com.h:142
void input_parameters(int argc, char *argv[])
Definition: SWcommon.c:75
void phase_compose_sequence()
Definition: SWpar.c:1036
#define hit_tileDimCard(var, dim)
Definition: hit_tile.h:750
ProteinFile pfile1
Definition: SWcommon.c:63
#define hit_tileDomainShapeAlloc(newVarP, baseType, shape)
Definition: hit_tile.h:354
HitTile_trace_t yTraceback
Definition: SWpar.c:248
int lxBegin
Definition: SWpar.c:261
#define yback(i, j)
Definition: SWpar.c:250
#define aa_hit_t
Definition: SWpar.c:75
HitTile_aa_t lout1
Definition: SWpar.c:239
#define hit_clockReset(c)
Definition: hit_utils.h:78
void initPAM(char *filename)
Definition: SWcommon.c:128
void hit_topFree(HitTopology topo)
Definition: hit_topology.c:129
#define allTimers(timer)
int size[2]
Definition: SWcommon.c:57
HitRanks hit_layActiveIdRanks(HitLayout lay, int id)
Definition: hit_layout.c:1935
int match_length
Definition: SWpar.c:222
void printProtein(HitTile_aa_t p, int offset, int psize)
Definition: SWcommon_hit.c:78
int match_length2
Definition: SWpar.c:224
void hit_patternFree(HitPattern *pattern)
Definition: hit_pattern.c:94
void phase_recv_back()
Definition: SWpar.c:862
int lmatch_length1
Definition: SWpar.c:266
void debug_show_backtracking()
Definition: SWpar.c:1009
void hit_comInit(int *pargc, char **pargv[])
Definition: hit_com.c:111
HitShape shape
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
HitShape lshape
Definition: SWpar.c:273
short aa_t
Definition: SWcommon.h:59
#define hit_shapeSigCard(shape, dim)
Definition: hit_sshape.h:412
HitTile_trace_t xTraceback
Definition: SWpar.c:247
void hit_layFree(HitLayout lay)
Definition: hit_layout.c:2251
#define hit_clockStop(c)
Definition: hit_utils.h:109
int main(int argc, char *argv[])
Definition: cannonAsync.c:62
#define iterations
Definition: mg.c:119
#define hit_comRecvSelect(lay, receiveFrom, tileP, selection, modeSelect, baseType)
Definition: hit_com.h:489
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 hit_layNeighbor(lay, dim, shift)
Definition: hit_layout.h:722
#define hit_comRecv(lay, receiveFrom, tileP, baseType)
Definition: hit_com.h:499
#define hit_comReduce(lay, root, tilePSend, tilePRecv, baseType, operation)
Definition: hit_com.h:725
#define hit_comTypeStruct(new_type, Nstruct, n,...)
Definition: hit_com.h:208
#define is_root()
Definition: SWcommon.h:123
HitShape hit_shapeDimExpand(HitShape shape, int dim, int position, int offset)
Definition: hit_shape.c:226
int match_length1
Definition: SWpar.c:223
int hit_layNumActives(HitLayout lay)
Definition: hit_layout.c:2213
HitTile_aa_t lp1
Definition: SWpar.c:235
void phase_recv_hmatrix()
Definition: SWpar.c:662
int lsize[2]
Definition: SWpar.c:188
#define hit_layout(name, topo,...)
Definition: hit_layout.h:415
#define hit_shapeSig(shape, dim)
Definition: hit_sshape.h:400
#define HIT_COM_TILECOORDS
Definition: hit_com.h:343
int xMax
Definition: SWpar.c:214
#define h_hit_t
Definition: SWpar.c:77
short trace_t
Definition: SWcommon.h:63
int xBegin
Definition: SWpar.c:218
#define hit_tileFree(var)
Definition: hit_tile.h:369
#define GAP_AA
Definition: SWcommon.h:243
#define printfroot(...)
Definition: SWcommon.h:135
HitTile_aa_t out2
Definition: SWpar.c:211
#define max(a, b)
Definition: cannonAsync.c:47
#define HIT_DOUBLE
Definition: hit_com.h:78
int lmatch_pos
Definition: SWpar.c:270
char * prot_name2
Definition: SWcommon.c:50
int yBegin
Definition: SWpar.c:219
#define MIN(a, b)
Definition: SWpar.c:1118
#define hit_patternDo(pattern)
Definition: hit_pattern.h:157
void hit_comFinalize()
Definition: hit_com.c:159
#define hit_clockGetSeconds(c)
Definition: hit_utils.h:190
void phase_read_sequences()
Definition: SWpar.c:619
#define h(i, j)
Definition: SWpar.c:244
#define hit_clockReduce(lay, c)
Definition: hit_utils.h:139
int lmatch_length
Definition: SWpar.c:265
MPI_Op opHMaxLoc
Definition: SWpar.c:191
#define hit_comSend(lay, sendTo, tileP, baseType)
Definition: hit_com.h:452
#define hit_clockStart(c)
Definition: hit_utils.h:87
void phase_comp_hmatrix()
Definition: SWpar.c:708
void hit_comStartSend(HitCom *issue)
Definition: hit_com.c:2051
#define HIT_SHAPE_BEGIN
Definition: hit_sshape.h:515