source: branches/help/GDE/TREEPUZZLE/src/ppuzzle.c

Last change on this file was 19518, checked in by westram, 3 months ago
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 70.7 KB
Line 
1/*
2 *   ppuzzle.c
3 *
4 *
5 * Part of TREE-PUZZLE 5.0 (June 2000)
6 *
7 * (c) 1999-2000 by Heiko A. Schmidt, Korbinian Strimmer,
8 *                  M. Vingron, and Arndt von Haeseler
9 * (c) 1995-1999 by Korbinian Strimmer and Arndt von Haeseler
10 *
11 * All parts of the source except where indicated are distributed under
12 * the GNU public licence.  See http://www.opensource.org for details.
13 */
14
15#define EXTERN extern
16 
17#include <mpi.h>
18#include <time.h>
19#include "ppuzzle.h"
20 
21
22int PP_IamMaster;
23int PP_IamSlave;
24int PP_Myid;
25int PP_MyMaster;
26int PP_NumProcs;
27MPI_Comm PP_Comm;
28
29int *freeslaves;           /* Queue of free slaves */
30int firstslave,            /* headpointer of queue */
31    lastslave;             /* tailpointer of queue */
32
33int *permutsent,
34    *permutrecved,
35    *quartsent,
36    *quartrecved,
37    *doquartsent,
38    *doquartrecved,
39    *splitsent,
40    *splitrecved,
41    *permutsentn,
42    *permutrecvedn,
43    *quartsentn,
44    *quartrecvedn,
45    *doquartsentn,
46    *doquartrecvedn,
47    *splitsentn,
48    *splitrecvedn;
49
50double *walltimes,
51       *cputimes;
52double *fullwalltimes,
53       *fullcputimes;
54double *altwalltimes,
55       *altcputimes;
56
57int PP_permutsent = 0;          /* # of  */
58int PP_permutrecved = 0;        /* # of  */
59int PP_quartsent = 0;           /* # of  */
60int PP_quartrecved = 0;         /* # of  */
61int PP_doquartsent = 0;         /* # of  */
62int PP_doquartrecved = 0;       /* # of  */
63int PP_splitsent = 0;           /* # of  */
64int PP_splitrecved = 0;         /* # of  */
65int PP_permutsentn = 0;          /* # of  */
66int PP_permutrecvedn = 0;        /* # of  */
67int PP_quartsentn = 0;           /* # of  */
68int PP_quartrecvedn = 0;         /* # of  */
69int PP_doquartsentn = 0;         /* # of  */
70int PP_doquartrecvedn = 0;       /* # of  */
71int PP_splitsentn = 0;           /* # of  */
72int PP_splitrecvedn = 0;         /* # of  */
73
74double PP_starttime     = 0,
75 PP_stoptime      = 0,
76 PP_inittime      = 0,
77 PP_paramcomptime = 0,
78 PP_paramsendtime = 0,
79 PP_quartcomptime = 0,
80 PP_quartsendtime = 0,
81 PP_puzzletime    = 0,
82 PP_treetime      = 0,
83 PP_lasttime      = 0;
84
85int PP_MaxSlave = 0;
86 
87
88/*********************************************************************
89*  miscellaneous utilities                                           *
90*********************************************************************/
91
92int dcmp(const void *a, const void *b)
93{
94        if (*(double *)a > *(double *)b) return (-1);
95        else if (*(double *)a < *(double *)b) return 1;
96        else return 0;
97}
98
99/******************/
100
101void PP_cmpd(int rank, double a, double b)
102{
103  if (a != b)
104     FPRINTF(STDOUTFILE "(%2d) *** %.3f != %.3f\n", rank, a, b);
105}
106
107/******************/
108
109void PP_cmpi(int rank, int a, int b)
110{
111  if (a != b)
112     FPRINTF(STDOUTFILE "(%2d) *** %d != %d\n", rank, a, b);
113}
114
115/******************/
116
117double PP_timer()
118{
119  double tmptime;
120  if (PP_lasttime == 0) {
121     PP_lasttime = MPI_Wtime();
122     return(0);
123     }
124  else {
125     tmptime = PP_lasttime;
126     PP_lasttime = MPI_Wtime();
127     return(PP_lasttime - tmptime);
128  }
129}
130
131/******************/
132 
133void PP_Printerror(FILE *of, int id, int err)
134{
135        char  errstr[MPI_MAX_ERROR_STRING];
136        int   errstrlen;
137
138  if ((err > MPI_SUCCESS) && (err <= MPI_ERR_LASTCODE)) {
139    MPI_Error_string(err, errstr, &errstrlen);
140    fprintf(of, "(%2d) MPI ERROR %d : %s\n", id, err, errstr);
141    }
142  else {
143    if (err == MPI_SUCCESS) 
144      fprintf(of, "(%2d) MPI ERROR %d : No error\n", id, err);
145    else
146      fprintf(of, "(%2d) MPI ERROR %d : unknown error number\n", id, err);
147  }
148} /* PP_Printerror */
149
150/******************/
151
152void PP_Printbiparts(cmatrix biparts)
153{ int n1, n2;
154    for (n1=0; n1<(Maxspc-3); n1++) {
155       if (n1==0) FPRINTF(STDOUTFILE "(%2d) bipartition : ", PP_Myid);
156       else       FPRINTF(STDOUTFILE "(%2d)             : ", PP_Myid);
157       for (n2=0; n2<Maxspc; n2++)
158          FPRINTF(STDOUTFILE "%c", biparts[n1][n2]);
159    FPRINTF(STDOUTFILE "\n");
160    }
161}  /* PP_Printbiparts */
162
163/******************/
164
165#if 0
166void num2quart(uli qnum, int *a, int *b, int *c, int *d)
167{
168        double temp;
169        uli aa, bb, cc, dd;
170        uli lowval=0, highval=0;
171
172        aa=0; bb=1; cc=2; dd=3;
173
174        temp = (double)(24 * qnum);
175        temp = sqrt(temp);
176        temp = sqrt(temp);
177        /* temp = pow(temp, (double)(1/4)); */
178        dd = (uli) floor(temp) + 1;
179        if (dd < 3) dd = 3;
180        lowval =  (uli) dd*(dd-1)*(dd-2)*(dd-3)/24;
181        highval = (uli) (dd+1)*dd*(dd-1)*(dd-2)/24;
182        if (lowval >= qnum)
183            while ((lowval > qnum)) {
184                dd -= 1; lowval = (uli) dd*(dd-1)*(dd-2)*(dd-3)/24;
185            }
186        else {
187            while (highval <= qnum) {
188                dd += 1; highval = (uli) (dd+1)*dd*(dd-1)*(dd-2)/24;
189            }
190            lowval = (uli) dd*(dd-1)*(dd-2)*(dd-3)/24;
191        }
192        qnum -= lowval;
193        if (qnum > 0) {
194            temp = (double)(6 * qnum);
195            temp = pow(temp, (double)(1/3));
196            cc = (uli) floor(temp);
197            if (cc < 2) cc= 2;
198            lowval =  (uli) cc*(cc-1)*(cc-2)/6;
199            highval = (uli) (cc+1)*cc*(cc-1)/6;
200            if (lowval >= qnum)
201                while ((lowval > qnum)) {
202                   cc -= 1; lowval = (uli) cc*(cc-1)*(cc-2)/6;
203                }
204            else {
205                while (highval <= qnum) {
206                   cc += 1; highval = (uli) (cc+1)*cc*(cc-1)/6;
207                }
208                lowval = (uli) cc*(cc-1)*(cc-2)/6;
209            }
210            qnum -= lowval;
211            if (qnum > 0) {
212                temp = (double)(2 * qnum);
213                temp = sqrt(temp);
214                bb = (uli) floor(temp);
215                if (bb < 1) bb= 1;
216                lowval =  (uli) bb*(bb-1)/2;
217                highval = (uli) (bb+1)*bb/2;
218                if (lowval >= qnum)
219                    while ((lowval > qnum)) {
220                        bb -= 1; lowval = (uli) bb*(bb-1)/2;
221                    }
222                else {
223                    while (highval <= qnum) {
224                       bb += 1; highval = (uli) (bb+1)*bb/2;
225                    }
226                    lowval = (uli) bb*(bb-1)/2;
227                }
228                qnum -= lowval;
229                if (qnum > 0) {
230                   aa = (uli) qnum;
231                   if (aa < 0) aa= 0;
232                }
233            }
234        }
235        *d = (int)dd;
236        *c = (int)cc;
237        *b = (int)bb;
238        *a = (int)aa;
239}  /* num2quart */
240
241/******************/
242
243uli numquarts(int maxspc)
244{
245   uli tmp;
246   int a, b, c, d;
247
248   if (maxspc < 4)
249     return (uli)0;
250   else {
251      maxspc--;
252      a = maxspc-3;
253      b = maxspc-2;
254      c = maxspc-1;
255      d = maxspc;
256
257      tmp = (uli) 1 + a +
258            (uli) b * (b-1) / 2 +
259            (uli) c * (c-1) * (c-2) / 6 +
260            (uli) d * (d-1) * (d-2) * (d-3) / 24;
261      return (tmp);
262   }
263}  /* numquarts */
264
265/******************/
266
267uli quart2num (int a, int b, int c, int d)
268{
269      uli tmp;
270      if ((a>b) || (b>c) || (c>d)) {
271         fprintf(stderr, "Error PP5 not (%d <= %d <= %d <= %d) !!!\n", a, b, c, d);
272         exit (1);
273      }
274      tmp = (uli) a +
275            (uli) b * (b-1) / 2 +
276            (uli) c * (c-1) * (c-2) / 6 +
277            (uli) d * (d-1) * (d-2) * (d-3) / 24;
278      return (tmp);
279}  /* quart2num */
280#endif
281/******************/
282
283
284/*********************************************************************
285*  queue for storing the ranks of slaves waiting for work            *
286*********************************************************************/
287
288void PP_initslavequeue()
289{
290  int n;
291  freeslaves = new_ivector(PP_NumProcs);
292    firstslave = 0;
293    PP_MaxSlave = PP_NumProcs-1;
294    lastslave  = PP_MaxSlave-1;
295    freeslaves[PP_MaxSlave] = PP_MaxSlave;
296    for (n=0; n<PP_MaxSlave; n++){
297      freeslaves[n] = n+1;
298    }
299}
300
301/******************/
302
303int PP_fullslave()
304{
305return (freeslaves[PP_MaxSlave] == PP_MaxSlave);
306}
307
308/******************/
309
310int PP_emptyslave()
311{
312return (freeslaves[PP_MaxSlave] == 0);
313}
314
315/******************/
316
317void PP_putslave(int sl)
318{
319  if (freeslaves[PP_MaxSlave] == PP_MaxSlave) {
320    FPRINTF(STDOUTFILE "\n\n\nHALT: PLEASE REPORT ERROR PP1 TO DEVELOPERS\n\n\n");
321    MPI_Finalize();
322    exit(1);
323  }
324  (freeslaves[PP_MaxSlave])++;
325  lastslave = (lastslave + 1) % PP_MaxSlave;
326  freeslaves[lastslave] = sl;
327}
328
329/******************/
330
331int PP_getslave()
332{
333  int tmp;
334  if (freeslaves[PP_MaxSlave] == 0)
335    return MPI_PROC_NULL;
336  else {
337    tmp = freeslaves[firstslave];
338    firstslave = (firstslave + 1) % PP_MaxSlave;
339    (freeslaves[PP_MaxSlave])--;
340    return tmp;
341  }
342}
343
344/*********************************************************************
345*  procedures to parallelize the puzzling step                       *
346*********************************************************************/
347
348void PP_slave_do_puzzling(ivector trueID)
349{
350int    i, a, b, c;
351uli    nq;
352#if SEQUENTIAL
353       double tc2, mintogo, minutes, hours;
354#endif
355
356     /* initialize tree */
357     inittree();
358
359     /* adding all other leafs */
360     for (i = 3; i < Maxspc; i++) { 
361 
362         /* clear all edgeinfos */
363         resetedgeinfo();
364 
365         /* clear counter of quartets */
366         nq = 0;
367 
368         /*
369           * core of quartet puzzling algorithm
370          */
371 
372          for (a = 0; a < nextleaf - 2; a++)
373              for (b = a + 1; b < nextleaf - 1; b++)
374                  for (c = b + 1; c < nextleaf; c++) {
375 
376                      /*      check which two _leaves_ out of a, b, c
377                          are closer related to each other than
378                          to leaf i according to a least squares
379                          fit of the continuous Baysian weights to the
380                          seven trivial "attractive regions". We assign
381                          a score of 1 to all edges between these two leaves
382                          chooseA and chooseB */
383 
384                      checkquartet(a, b, c, i);
385                      incrementedgeinfo(chooseA, chooseB);
386 
387                      nq++;
388 
389#                     if SEQUENTIAL
390                         /* generate message every 15 minutes */
391   
392                         /* check timer */
393                         time(&time2);
394                         if ( (time2 - time1) > 900) {
395                              /* every 900 seconds */
396                              /* percentage of completed trees */
397                              if (mflag == 0) {
398                                      FPRINTF(STDOUTFILE "\n");
399                                      mflag = 1;
400                              }
401                              tc2 = 100.0*Currtrial/Numtrial + 
402                                    100.0*nq/Numquartets/Numtrial;
403                              mintogo = (100.0-tc2) *
404                                        (double) (time2-time0)/60.0/tc2;
405                              hours = floor(mintogo/60.0);
406                              minutes = mintogo - 60.0*hours;
407                              FPRINTF(STDOUTFILE "%2.2f%%", tc2);
408                              FPRINTF(STDOUTFILE " completed (remaining");
409                              FPRINTF(STDOUTFILE " time: %.0f", hours);
410                              FPRINTF(STDOUTFILE " hours %.0f", minutes);
411                              FPRINTF(STDOUTFILE " minutes)\n");
412                              time1 = time2;
413                         }
414#                     endif /* SEQUENTIAL */
415              }
416
417          /* find out which edge has the lowest edgeinfo */
418          minimumedgeinfo();
419 
420          /* add the next leaf on minedge */
421          addnextleaf(minedge);
422     }
423 
424     /* compute bipartitions of current tree */
425     computebiparts();
426
427#if  PARALLEL
428       if (PP_IamMaster) makenewsplitentries();
429#    else
430       makenewsplitentries();
431#    endif
432
433        {
434                int *ctree, startnode;
435                char *trstr;
436                ctree = initctree();
437                copytree(ctree);
438                startnode = sortctree(ctree);
439                trstr=sprintfctree(ctree, psteptreestrlen);
440                (void) addtree2list(&trstr, 1, &psteptreelist, &psteptreenum, &psteptreesum);
441#               ifdef PVERBOSE2
442                        /* fprintf(STDOUT, "%s\n", trstr); */
443                        printfpstrees(psteptreelist);
444#               endif
445                freectree(&ctree);
446        }
447
448
449     /* free tree before building the next tree */
450     freetree();
451
452} /* PP_slave_do_puzzling */
453
454/******************/
455
456void PP_do_puzzling(ivector trueID)
457{
458int    dest;
459
460#    if PARALLEL
461       dest = PP_getslave();
462       PP_SendPermut(dest, Maxspc, trueID);
463#    endif
464
465     /* initialize tree */
466     inittree();
467
468     PP_RecvSplits(Maxspc, biparts); 
469
470#    ifdef PVERBOSE3
471       PP_Printbiparts(biparts);
472#    endif /* PVERBOSE3 */
473
474     makenewsplitentries();
475
476     /* free tree before building the next tree */
477     freetree();
478
479} /* PP_do_puzzling */
480 
481/******************/
482
483
484void PP_do_write_quart(int e,
485                       int f,
486                       int g,
487                       int h,
488                       double d1,
489                       double d2,
490                       double d3,
491                       uli   *numbq,
492                       uli   *bqarr)
493{
494        double lhs[3],
495               temp,
496               wlist[6],
497               plist[6];
498        unsigned char qpbranching;
499        int badquartet;
500
501        lhs[0] = d1;
502        lhs[1] = d2;
503        lhs[2] = d3;
504
505         badquartet = FALSE;
506
507         /* compute Bayesian weights */
508         temp = (lhs[0] + lhs[1] + lhs[2])/3.0;
509         lhs[0] = exp(lhs[0] - temp);
510         lhs[1] = exp(lhs[1] - temp);
511         lhs[2] = exp(lhs[2] - temp);
512         temp = lhs[0] + lhs[1] + lhs[2];
513         wlist[0] = lhs[0] / temp;
514         wlist[1] = 1.0;
515         wlist[2] = lhs[1] / temp;
516         wlist[3] = 2.0;
517         wlist[4] = lhs[2] / temp;
518         wlist[5] = 4.0;
519
520         /* sort in descending order */
521         qsort(wlist, 3, 2*sizeof(double), dcmp);
522
523         /* check out the three possibilities */
524
525         /* 100 distribution */
526         plist[0] = (1.0 - wlist[0])*(1.0 - wlist[0]) +
527                    (0.0 - wlist[2])*(0.0 - wlist[2]) +
528                    (0.0 - wlist[4])*(0.0 - wlist[4]);
529         plist[1] = wlist[1];
530
531         /* 110 distribution */
532         plist[2] = (0.5 - wlist[0])*(0.5 - wlist[0]) +
533                    (0.5 - wlist[2])*(0.5 - wlist[2]) +
534                    (0.0 - wlist[4])*(0.0 - wlist[4]);
535         plist[3] = wlist[1] + wlist[3];
536
537         /* 111 distribution */
538         temp = 1.0/3.0;
539         plist[4] = (temp - wlist[0])*(temp - wlist[0]) +
540                    (temp - wlist[2])*(temp - wlist[2]) +
541                    (temp - wlist[4])*(temp - wlist[4]);
542         plist[5] = wlist[1] + wlist[3] + wlist[5];
543
544         /* sort in descending order */
545         qsort(plist, 3, 2*sizeof(double), dcmp);
546
547         qpbranching = (unsigned char) plist[5];
548         writequartet(e, f, g, h, qpbranching);
549
550         /* a bad quartet is a quartet that shows
551            equal weights for all three possible topologies */
552         if (qpbranching == 7) badquartet = TRUE;
553
554         if (badquartet) {
555               bqarr[(*numbq)++] = quart2num(e, f, g, h);
556#              ifdef PVERBOSE3
557                  FPRINTF(STDOUTFILE "(%2d) bad quartet: %d  %d  %d %d -> %ld\n", 
558                          PP_Myid, e, f, g, h, quart2num(e, f, g, h));
559#              endif /* PVERBOSE3 */
560               badqs++;
561               badtaxon[e]++;
562               badtaxon[f]++;
563               badtaxon[g]++;
564               badtaxon[h]++;
565         } /* if badquartet */
566} /* PP_do_write_quart */
567
568/*********************************************************************
569*  sending/receiving the important sizes and parameter  (M->S)       *
570*********************************************************************/
571 
572void PP_SendSizes(int    mspc, 
573                  int    msite,
574                  int    ncats,
575                  int    nptrn,
576                  int    rad,
577                  int    outgr,
578                  double frconst,
579                  int    rseed)
580{
581# define NUMINT 7
582# define NUMDBL 1
583  int    ints[NUMINT];
584  double doubles[NUMDBL];
585  MPI_Datatype Dtypes[2] =    {MPI_INT, MPI_DOUBLE};
586  int          Dtypelens[2] = {NUMINT , NUMDBL};
587  MPI_Aint     Dtypeaddr[2];
588  MPI_Datatype PP_Sizes;
589  int          dest;
590  int          error;
591
592# ifdef PVERBOSE2
593    FPRINTF(STDOUTFILE "(%2d) Sending: Maxspc=%d  Maxsite=%d  numcats=%d\n", PP_Myid, mspc, msite, ncats);
594    FPRINTF(STDOUTFILE "(%2d)          Numprtn=%d  tpmradix=%d  fracconst=%.3f\n", PP_Myid, nptrn, rad, frconst);
595# endif /* PVERBOSE2 */
596
597  ints[0] = mspc;
598  ints[1] = msite;
599  ints[2] = ncats;
600  ints[3] = nptrn;
601  ints[4] = rad;
602  ints[5] = outgr;
603  ints[6] = rseed;
604  doubles[0] = frconst;
605
606  MPI_Address(ints,     Dtypeaddr);
607  MPI_Address(doubles, (Dtypeaddr+1));
608 
609  MPI_Type_struct(2, Dtypelens, Dtypeaddr, Dtypes, &PP_Sizes);
610  MPI_Type_commit(&PP_Sizes);
611
612  for (dest=1; dest<PP_NumProcs; dest++) {
613
614    error = MPI_Ssend(MPI_BOTTOM, 1, PP_Sizes, dest, PP_SIZES, PP_Comm);
615    if (error != MPI_SUCCESS)
616      PP_Printerror(STDOUT, 600+PP_Myid, error);
617
618#   ifdef PVERBOSE3
619       FPRINTF(STDOUTFILE "(%2d) -> (%2d) Sent Sizes\n", PP_Myid, dest);
620#   endif /* PVERBOSE3 */
621
622  } /* for each slave */
623
624  MPI_Type_free(&PP_Sizes);
625
626# ifdef PVERBOSE3
627    FPRINTF(STDOUTFILE "(%2d) ... Sent Sizes\n", PP_Myid);
628# endif /* PVERBOSE3 */
629
630# undef NUMINT
631# undef NUMDBL
632} /* PP_SendSizes */
633
634
635/******************/
636
637void PP_RecvSizes(int    *mspc, 
638                  int    *msite,
639                  int    *ncats,
640                  int    *nptrn,
641                  int    *rad,
642                  int    *outgr,
643                  double *frconst,
644                  int    *rseed)
645{
646# define NUMINT 7
647# define NUMDBL 1
648  int    ints[NUMINT];
649  double doubles[NUMDBL];
650  MPI_Datatype Dtypes[2] =    {MPI_INT, MPI_DOUBLE};
651  int          Dtypelens[2] = {NUMINT , NUMDBL};
652  MPI_Aint     Dtypeaddr[2];
653  MPI_Datatype PP_Sizes;
654  MPI_Status   stat;
655  int          error;
656 
657# ifdef PVERBOSE3
658    FPRINTF(STDOUTFILE "(%2d) Receiving Sizes ...\n", PP_Myid);
659# endif /* PVERBOSE3 */
660
661  MPI_Address(ints,     Dtypeaddr);
662  MPI_Address(doubles, (Dtypeaddr+1));
663 
664  MPI_Type_struct(2, Dtypelens, Dtypeaddr, Dtypes, &PP_Sizes);
665  MPI_Type_commit(&PP_Sizes);
666 
667  error = MPI_Probe(PP_MyMaster, MPI_ANY_TAG, PP_Comm, &stat);
668  if (error != MPI_SUCCESS)
669    PP_Printerror(STDOUT, 700+PP_Myid, error);
670  if (stat.MPI_TAG != PP_SIZES) {
671        if (stat.MPI_TAG == PP_DONE) {
672                PP_RecvDone();
673#               ifdef PVERBOSE1
674                        FPRINTF(STDOUTFILE "(%2d) Finishing...\n", PP_Myid);
675#               endif /* PVERBOSE1 */
676                MPI_Finalize();
677                exit(1);
678        } else {
679                FPRINTF(STDOUTFILE "(%2d) Error: unexpected TAG received...\n", PP_Myid);
680                MPI_Finalize();
681                exit(1);
682        }
683  }
684
685  error = MPI_Recv(MPI_BOTTOM, 1, PP_Sizes, PP_MyMaster, MPI_ANY_TAG, PP_Comm, &stat);
686  if (error != MPI_SUCCESS)
687    PP_Printerror(STDOUT, 700+PP_Myid, error);
688  if (stat.MPI_TAG != PP_SIZES) {
689        FPRINTF(STDOUTFILE "(%2d) Error: unexpected TAG received...\n", PP_Myid);
690        MPI_Finalize();
691        exit(1);
692  }
693
694  *mspc    = ints[0];
695  *msite   = ints[1];
696  *ncats   = ints[2];
697  *nptrn   = ints[3];
698  *rad     = ints[4];
699  *outgr   = ints[5];
700  *rseed   = ints[6];
701  *frconst = doubles[0];
702 
703  MPI_Type_free(&PP_Sizes);
704 
705# ifdef PVERBOSE2
706    FPRINTF(STDOUTFILE "(%2d) <- (%2d) Received: Maxspec=%d  Maxsite=%d  numcats=%d\n", PP_Myid, PP_MyMaster, *mspc, *msite, *ncats);
707    FPRINTF(STDOUTFILE "(%2d)                   Numprtn=%d  tpmradix=%d  fracconst=%.3f\n", PP_Myid, *nptrn, *rad, *frconst);
708# endif /* PVERBOSE2 */
709
710# undef NUMINT
711# undef NUMDBL
712} /* PP_RecvSizes */
713 
714
715
716/*********************************************************************
717*  sending/receiving the data matrices  (M->S)                       *
718*********************************************************************/
719
720void PP_RecvData(
721        cmatrix Seqpat,           /* cmatrix (Maxspc x Numptrn)             */
722        ivector Alias,            /* ivector (Maxsite)                      */
723        ivector Weight,           /* ivector (Numptrn)                      */
724        ivector constpat,
725        dvector Rates,            /* dvector (numcats)                      */
726        dvector Eval,             /* dvector (tpmradix)                     */
727        dvector Freqtpm,
728        dmatrix Evec,             /* dmatrix (tpmradix x tpmradix)          */
729        dmatrix Ievc,
730        dmatrix iexp,
731        dmatrix Distanmat,        /* dmatrix (Maxspc x Maxspc)              */
732        dcube   ltprobr)          /* dcube (numcats x tpmradix x tpmradix)  */
733{
734  MPI_Datatype Dtypes[12];
735  int          Dtypelens[12];
736  MPI_Aint     Dtypeaddr[12];
737  MPI_Datatype PP_Data;
738  MPI_Status   stat;
739  int          error;
740 
741# ifdef PVERBOSE3
742    FPRINTF(STDOUTFILE "(%2d) Receiving Sizes ...\n", PP_Myid);
743# endif /* PVERBOSE2 */
744
745  Dtypes  [0] = MPI_CHAR; Dtypelens [0] = Maxspc * Numptrn;
746  MPI_Address(&(Seqpat[0][0]), &(Dtypeaddr[0]));
747  Dtypes  [1] = MPI_INT; Dtypelens  [1] = Maxsite ;
748  MPI_Address(&(Alias[0]), &(Dtypeaddr[1]));
749  Dtypes  [2] = MPI_INT; Dtypelens  [2] = Numptrn ;
750  MPI_Address(&(Weight[0]), &(Dtypeaddr[2]));
751  Dtypes  [3] = MPI_INT; Dtypelens  [3] = Numptrn ;
752  MPI_Address(&(constpat[0]), &(Dtypeaddr[3]));
753  Dtypes  [4] = MPI_DOUBLE; Dtypelens  [4] = numcats ;
754  MPI_Address(&(Rates[0]), &(Dtypeaddr[4]));
755  Dtypes  [5] = MPI_DOUBLE; Dtypelens  [5] = tpmradix ;
756  MPI_Address(&(Eval[0]), &(Dtypeaddr[5]));
757  Dtypes  [6] = MPI_DOUBLE; Dtypelens  [6] = tpmradix ;
758  MPI_Address(&(Freqtpm[0]), &(Dtypeaddr[6]));
759  Dtypes  [7] = MPI_DOUBLE; Dtypelens  [7] = tpmradix * tpmradix ;
760  MPI_Address(&(Evec[0][0]), &(Dtypeaddr[7]));
761  Dtypes  [8] = MPI_DOUBLE; Dtypelens  [8] = tpmradix * tpmradix ;
762  MPI_Address(&(Ievc[0][0]), &(Dtypeaddr[8]));
763  Dtypes  [9] = MPI_DOUBLE; Dtypelens [9] = tpmradix * tpmradix ;
764  MPI_Address(&(iexp[0][0]), &(Dtypeaddr[9]));
765  Dtypes [10] = MPI_DOUBLE; Dtypelens [10] = Maxspc * Maxspc ;
766  MPI_Address(&(Distanmat[0][0]), &(Dtypeaddr[10]));
767  Dtypes [11] = MPI_DOUBLE; Dtypelens [11] = numcats * tpmradix * tpmradix ;
768  MPI_Address(&(ltprobr[0][0][0]), &(Dtypeaddr[11]));
769 
770  MPI_Type_struct(12, Dtypelens, Dtypeaddr, Dtypes, &PP_Data);
771  MPI_Type_commit(&PP_Data);
772 
773
774  error = MPI_Probe(PP_MyMaster, MPI_ANY_TAG, PP_Comm, &stat);
775  if (error != MPI_SUCCESS)
776    PP_Printerror(STDOUT, 700+PP_Myid, error);
777  if (stat.MPI_TAG != PP_DATA) {
778        if (stat.MPI_TAG == PP_DONE) {
779                PP_RecvDone();
780#               ifdef PVERBOSE1
781                        FPRINTF(STDOUTFILE "(%2d) Finishing...\n", PP_Myid);
782#               endif /* PVERBOSE1 */
783                MPI_Finalize();
784                exit(1);
785        } else {
786                FPRINTF(STDOUTFILE "(%2d) Error: unexpected TAG received...\n", PP_Myid);
787                MPI_Finalize();
788                exit(1);
789        }
790  }
791
792
793  error = MPI_Recv(MPI_BOTTOM, 1, PP_Data, PP_MyMaster, PP_DATA, PP_Comm, &stat);
794  if (error != MPI_SUCCESS)
795    PP_Printerror(STDOUT, 900+PP_Myid, error);
796
797# ifdef PVERBOSE2
798    FPRINTF(STDOUTFILE "(%2d) <- (%2d) Received : Alias(0)=%d - Weight(0)=%d - constpat(0)=%d\n", PP_Myid, PP_MyMaster, Alias[0], Weight[0], constpat[0]);
799    FPRINTF(STDOUTFILE "(%2d)                    Rates(0)=%.3f - Eval(0)=%.3f - Freqtpm(0)=%.3f\n", PP_Myid, Rates[0], Eval[0], Freqtpm[0]);
800    FPRINTF(STDOUTFILE "(%2d)                    Evec(0,0)=%.3f - Ievc(0,0)=%.3f - iexp(0,0)=%.3f - Distanmat(0,1)=%.3f\n", PP_Myid, Evec[0][0], Ievc[0][0], iexp[0][0], Distanmat[0][1]);
801    FPRINTF(STDOUTFILE "(%2d)                    Distanmat(0,1)=%.3f\n", PP_Myid, Distanmat[0][1]);
802    FPRINTF(STDOUTFILE "(%2d)                    ltprobr(0,0,0)=%.3f\n", PP_Myid, ltprobr[0][0][0]);
803# endif /* PVERBOSE2 */
804 
805  MPI_Type_free(&PP_Data);
806 
807} /* PP_RecvData */
808
809
810/******************/
811
812void PP_SendData(
813        cmatrix Seqpat,           /* cmatrix (Maxspc x Numptrn)             */
814        ivector Alias,            /* ivector (Maxsite)                      */
815        ivector Weight,           /* ivector (Numptrn)                      */
816        ivector constpat,
817        dvector Rates,            /* dvector (numcats)                      */
818        dvector Eval,             /* dvector (tpmradix)                     */
819        dvector Freqtpm,
820        dmatrix Evec,             /* dmatrix (tpmradix x tpmradix)          */
821        dmatrix Ievc,
822        dmatrix iexp,
823        dmatrix Distanmat,        /* dmatrix (Maxspc x Maxspc)              */
824        dcube   ltprobr)          /* dcube (numcats x tpmradix x tpmradix)  */
825{
826  MPI_Datatype Dtypes[12];
827  int          Dtypelens[12];
828  MPI_Aint     Dtypeaddr[12];
829  MPI_Datatype PP_Data;
830  int          dest;
831  int          error;
832 
833# ifdef PVERBOSE2
834    FPRINTF(STDOUTFILE "(%2d) Sending: Alias(0)=%d - Weight(0)=%d - constpat(0)=%d\n", PP_Myid, Alias[0], Weight[0], constpat[0]);
835    FPRINTF(STDOUTFILE "(%2d)          Rates(0)=%.3f - Eval(0)=%.3f - Freqtpm(0)=%.3f\n", PP_Myid, Rates[0], Eval[0], Freqtpm[0]);
836    FPRINTF(STDOUTFILE "(%2d)          Evec(0,0)=%.3f - Ievc(0,0)=%.3f - iexp(0,0)=%.3f - Distanmat(0,1)=%.3f\n", PP_Myid, Evec[0][0], Ievc[0][0], iexp[0][0], Distanmat[0][1]);
837    FPRINTF(STDOUTFILE "(%2d)          ltprobr(0,0,0)=%.3f\n", PP_Myid, ltprobr[0][0][0]);
838# endif /* PVERBOSE2 */
839 
840  Dtypes  [0] = MPI_CHAR; Dtypelens [0] = Maxspc * Numptrn;
841  MPI_Address(&(Seqpat[0][0]), &(Dtypeaddr[0]));
842  Dtypes  [1] = MPI_INT; Dtypelens  [1] = Maxsite ;
843  MPI_Address(&(Alias[0]), &(Dtypeaddr[1]));
844  Dtypes  [2] = MPI_INT; Dtypelens  [2] = Numptrn ;
845  MPI_Address(&(Weight[0]), &(Dtypeaddr[2]));
846  Dtypes  [3] = MPI_INT; Dtypelens  [3] = Numptrn ;
847  MPI_Address(&(constpat[0]), &(Dtypeaddr[3]));
848  Dtypes  [4] = MPI_DOUBLE; Dtypelens  [4] = numcats ;
849  MPI_Address(&(Rates[0]), &(Dtypeaddr[4]));
850  Dtypes  [5] = MPI_DOUBLE; Dtypelens  [5] = tpmradix ;
851  MPI_Address(&(Eval[0]), &(Dtypeaddr[5]));
852  Dtypes  [6] = MPI_DOUBLE; Dtypelens  [6] = tpmradix ;
853  MPI_Address(&(Freqtpm[0]), &(Dtypeaddr[6]));
854  Dtypes  [7] = MPI_DOUBLE; Dtypelens  [7] = tpmradix * tpmradix ;
855  MPI_Address(&(Evec[0][0]), &(Dtypeaddr[7]));
856  Dtypes  [8] = MPI_DOUBLE; Dtypelens  [8] = tpmradix * tpmradix ;
857  MPI_Address(&(Ievc[0][0]), &(Dtypeaddr[8]));
858  Dtypes  [9] = MPI_DOUBLE; Dtypelens  [9] = tpmradix * tpmradix ;
859  MPI_Address(&(iexp[0][0]), &(Dtypeaddr [9]));
860  Dtypes [10] = MPI_DOUBLE; Dtypelens [10] = Maxspc * Maxspc ;
861  MPI_Address(&(Distanmat[0][0]), &(Dtypeaddr[10]));
862  Dtypes [11] = MPI_DOUBLE; Dtypelens [11] = numcats * tpmradix * tpmradix ;
863  MPI_Address(&(ltprobr[0][0][0]), &(Dtypeaddr[11]));
864 
865  MPI_Type_struct(12, Dtypelens, Dtypeaddr, Dtypes, &PP_Data);
866  MPI_Type_commit(&PP_Data);
867 
868  for (dest=1; dest<PP_NumProcs; dest++) {
869 
870    error = MPI_Ssend(MPI_BOTTOM, 1, PP_Data, dest, PP_DATA, PP_Comm);
871    if (error != MPI_SUCCESS)
872      PP_Printerror(STDOUT, 1100+PP_Myid, error);
873 
874#     ifdef PVERBOSE3
875         FPRINTF(STDOUTFILE "(%2d) -> (%2d) Sent Data\n", PP_Myid, dest);
876#     endif /* PVERBOSE2 */
877 
878  } /* for each slave */
879 
880  MPI_Type_free(&PP_Data);
881
882# ifdef PVERBOSE3
883    FPRINTF(STDOUTFILE "(%2d) ... Sent Data\n", PP_Myid);
884# endif /* PVERBOSE2 */
885 
886} /* PP_SendData */
887
888
889/**************************************************************************
890*  procedures to send the request to compute a single quartet  (M->S)     *
891**************************************************************************/
892 
893void PP_SendDoQuart(int    dest,
894                    int    a, 
895                    int    b,
896                    int    c,
897                    int    d,
898                    int    approx)
899{
900# define NUMINT 5
901  int    ints[NUMINT];
902  int    error;
903 
904  ints[0] = a;
905  ints[1] = b;
906  ints[2] = c;
907  ints[3] = d;
908  ints[4] = approx;
909
910  PP_doquartsent++;
911  PP_doquartsentn++;
912
913# ifdef PVERBOSE2
914    FPRINTF(STDOUTFILE "(%2d) Sending  -> (%2d): Quart(%d,%d,%d,%d)\n", PP_Myid, dest, a, b, c, d);
915# endif /* PVERBOSE2 */
916 
917  error = MPI_Ssend(ints, NUMINT, MPI_INT, dest, PP_DOQUART, PP_Comm);
918  if (error != MPI_SUCCESS)
919     PP_Printerror(STDOUT, PP_Myid, error);
920 
921# ifdef PVERBOSE3
922    FPRINTF(STDOUTFILE "(%2d) ... Sent \n", PP_Myid);
923# endif /* PVERBOSE3 */
924# undef NUMINT
925 
926} /* PP_SendDoQuart */
927
928
929
930/******************/
931
932void PP_RecvDoQuart(int    *a, 
933                    int    *b,
934                    int    *c,
935                    int    *d,
936                    int    *approx)
937{
938# define NUMINT 5
939  int    ints[NUMINT];
940  int    error;
941  MPI_Status   stat;
942  PP_doquartrecved++;
943  PP_doquartrecvedn++;
944
945# ifdef PVERBOSE3
946    FPRINTF(STDOUTFILE "(%2d) Receiving: Quart\n", PP_Myid);
947# endif /* PVERBOSE3 */
948 
949  error = MPI_Recv(ints, NUMINT, MPI_INT, PP_MyMaster, PP_DOQUART, PP_Comm, &stat);
950  if (error != MPI_SUCCESS)
951    PP_Printerror(STDOUT, 200+PP_Myid, error);
952 
953  *a      = ints[0];
954  *b      = ints[1];
955  *c      = ints[2];
956  *d      = ints[3];
957  *approx = ints[4];
958 
959# ifdef PVERBOSE2
960    FPRINTF(STDOUTFILE "(%2d) Received: Quart(%d,%d,%d,%d,%c)\n", PP_Myid, *a, *b, *c, *d, (approx ? 'A' : 'E'));
961# endif /* PVERBOSE2 */
962# undef NUMINT
963 
964} /* PP_RecvDoQuart */
965 
966 
967/**************************************************************************
968*  procedures to send the result of a single quartet  (S->M)              *
969**************************************************************************/
970 
971void PP_SendQuart(int    a, 
972                  int    b,
973                  int    c,
974                  int    d,
975                  double d1,
976                  double d2,
977                  double d3,
978                  int    approx)
979{
980# define NUMINT 5
981# define NUMDBL 3
982  int    ints[NUMINT];
983  double doubles[NUMDBL];
984  MPI_Datatype Dtypes[2] =    {MPI_INT, MPI_DOUBLE};
985  int          Dtypelens[2] = {NUMINT , NUMDBL};
986  MPI_Aint     Dtypeaddr[2];
987  MPI_Datatype PP_Quart;
988  int          error;
989 
990  PP_quartsent++;
991  PP_quartsentn++;
992  ints[0] = a;
993  ints[1] = b;
994  ints[2] = c;
995  ints[3] = d;
996  ints[4] = approx;
997  doubles[0] = d1;
998  doubles[1] = d2;
999  doubles[2] = d3;
1000 
1001  MPI_Address(ints,     Dtypeaddr);
1002  MPI_Address(doubles, (Dtypeaddr+1));
1003 
1004  MPI_Type_struct(2, Dtypelens, Dtypeaddr, Dtypes, &PP_Quart);
1005  MPI_Type_commit(&PP_Quart);
1006 
1007# ifdef PVERBOSE2
1008    FPRINTF(STDOUTFILE "(%2d) Sending: Quart(%d,%d,%d,%d) = (%.3f, %.3f, %.3f)\n", PP_Myid, a, b, c, d, d1, d2, d3);
1009# endif /* PVERBOSE2 */
1010 
1011  error = MPI_Ssend(MPI_BOTTOM, 1, PP_Quart, PP_MyMaster, PP_QUART, PP_Comm);
1012  if (error != MPI_SUCCESS)
1013    PP_Printerror(STDOUT, 400+PP_Myid, error);
1014 
1015  MPI_Type_free(&PP_Quart);
1016
1017# ifdef PVERBOSE3
1018    FPRINTF(STDOUTFILE "(%2d) ... Sent \n", PP_Myid);
1019# endif /* PVERBOSE3 */
1020# undef NUMINT
1021# undef NUMDBL
1022
1023} /* PP_SendQuart */
1024 
1025 
1026 
1027/******************/
1028 
1029void PP_RecvQuart(int    *a, 
1030                  int    *b,
1031                  int    *c,
1032                  int    *d,
1033                  double *d1,
1034                  double *d2,
1035                  double *d3,
1036                  int    *approx)
1037{
1038# define NUMINT 5
1039# define NUMDBL 3
1040  int          ints[NUMINT];
1041  double       doubles[NUMDBL];
1042  MPI_Datatype Dtypes[2] =    {MPI_INT, MPI_DOUBLE};
1043  int          Dtypelens[2] = {NUMINT , NUMDBL};
1044  MPI_Aint     Dtypeaddr[2];
1045  MPI_Datatype PP_Quart;
1046  int          error;
1047  MPI_Status   stat;
1048 
1049  PP_quartrecved++;
1050  PP_quartrecvedn++;
1051  MPI_Address(ints,     Dtypeaddr);
1052  MPI_Address(doubles, (Dtypeaddr+1));
1053 
1054  MPI_Type_struct(2, Dtypelens, Dtypeaddr, Dtypes, &PP_Quart);
1055  MPI_Type_commit(&PP_Quart);
1056 
1057  error = MPI_Recv(MPI_BOTTOM, 1, PP_Quart, MPI_ANY_SOURCE, PP_QUART, PP_Comm, &stat);
1058
1059  if (error != MPI_SUCCESS)
1060    PP_Printerror(STDOUT, 500+PP_Myid, error);
1061
1062  PP_putslave(stat.MPI_SOURCE);
1063 
1064  *a      = ints[0];
1065  *b      = ints[1];
1066  *c      = ints[2];
1067  *d      = ints[3];
1068  *d1     = doubles[0];
1069  *d2     = doubles[1];
1070  *d3     = doubles[2];
1071  *approx = ints[4];
1072
1073  MPI_Type_free(&PP_Quart);
1074 
1075# ifdef PVERBOSE2
1076    FPRINTF(STDOUTFILE "(%2d) Received <- (%2d): Quart(%d,%d,%d,%d)=(%.3f, %.3f, %.3f)\n", PP_Myid, stat.MPI_SOURCE, *a, *b, *c, *d, *d1, *d2, *d3);
1077# endif /* PVERBOSE2 */
1078# undef NUMINT
1079# undef NUMDBL
1080 
1081} /* PP_RecvQuart */
1082
1083
1084 
1085/**************************************************************************
1086*  procedures to send the request to compute a block of quartets  (M->S)  *
1087**************************************************************************/
1088
1089void PP_SendDoQuartBlock(int dest, uli firstq, uli amount, int approx)
1090{
1091#       define NUMULI 3
1092        uli           ulongs[NUMULI];
1093        int           error;
1094 
1095        PP_doquartsent += amount;
1096        PP_doquartsentn++;
1097#       ifdef PVERBOSE2
1098           FPRINTF(STDOUTFILE "(%2d) Sending: DOQuartBlock Signal\n", PP_Myid);
1099#       endif /* PVERBOSE2 */
1100 
1101        ulongs[0] = firstq;
1102        ulongs[1] = amount;
1103        ulongs[2] = (uli)approx;
1104
1105        error = MPI_Ssend(ulongs, NUMULI, MPI_UNSIGNED_LONG, dest, PP_DOQUARTBLOCK, PP_Comm);
1106        if (error != MPI_SUCCESS) PP_Printerror(STDOUT, 2100+PP_Myid, error);
1107 
1108#       ifdef PVERBOSE3
1109           FPRINTF(STDOUTFILE "(%2d) ... Sent DOQuartBlock Signal (addr:%ld, num:%ld)\n", PP_Myid, firstq, amount);
1110#       endif /* PVERBOSE3 */
1111#       undef NUMULI
1112 
1113} /* PP_SendDoQuartBlock */
1114
1115/******************/
1116
1117void PP_RecvDoQuartBlock(uli *firstq, uli *amount, uli **bq, int *approx)
1118{
1119#       define NUMULI 3
1120        uli           ulongs[NUMULI];
1121        MPI_Status    stat;
1122        int           error;
1123 
1124#       ifdef PVERBOSE2
1125          FPRINTF(STDOUTFILE "(%2d) Receiving: DOQuartBlock Signal\n", PP_Myid);
1126#       endif /* PVERBOSE2 */
1127 
1128        error = MPI_Recv(&ulongs, NUMULI, MPI_UNSIGNED_LONG, PP_MyMaster, PP_DOQUARTBLOCK, PP_Comm, &stat);
1129        if (error != MPI_SUCCESS)
1130           PP_Printerror(STDOUT, 2100+PP_Myid, error);
1131
1132        *firstq=ulongs[0];
1133        *amount=ulongs[1];
1134        *approx= (int)ulongs[2];
1135
1136        *bq = malloc((unsigned)*amount * sizeof(uli));
1137
1138        PP_doquartrecved += *amount;
1139        PP_doquartrecvedn++;
1140
1141#       ifdef PVERBOSE3
1142           FPRINTF(STDOUTFILE "(%2d) ... DOQuartBlock (addr:%ld, num:%ld)\n", 
1143                              PP_Myid, *firstq, *amount);
1144#       endif /* PVERBOSE3 */
1145 
1146#       undef NUMULI
1147} /* PP_RecvDoQuartBlock */
1148
1149/*********************************************************************
1150*  procedures to send the results of a block of quartets  (S->M)     *
1151*********************************************************************/
1152
1153void PP_SendQuartBlock(uli startq,
1154                       uli numofq,
1155                       unsigned char *quartetinfo,
1156                       uli  numofbq,
1157                       uli *bq,
1158                       int approx)
1159{
1160# define NUMULI 3
1161# define NUMINT 1
1162  unsigned char *trueaddr;
1163  uli            truenum;
1164  int            error;
1165  int    ints[NUMINT];
1166  uli    ulis[NUMULI];
1167  MPI_Datatype Dtypes[2] =    {MPI_UNSIGNED_LONG, MPI_INT};
1168  int          Dtypelens[2] = {NUMULI,            NUMINT};
1169  MPI_Aint     Dtypeaddr[2];
1170  MPI_Datatype PP_QBlockSpecs;
1171  MPI_Datatype DtypesRes[2] = {MPI_UNSIGNED_CHAR, MPI_UNSIGNED_LONG};
1172  int          DtypelensRes[2];
1173  MPI_Aint     DtypeaddrRes[2];
1174  MPI_Datatype PP_QBlockRes;
1175
1176/*
1177  uli *bq;
1178  uli  numofbq;
1179*/
1180
1181  PP_quartsent += numofq;
1182  PP_quartsentn++;
1183
1184  truenum = (uli)((numofq+1)/2);
1185  trueaddr = (unsigned char *)(quartetinfo + (uli)(startq/2));
1186
1187# ifdef PVERBOSE2
1188    FPRINTF(STDOUTFILE "(%2d) Sending: startq=%lud numofq=%lud\n", PP_Myid, startq, numofq);
1189    FPRINTF(STDOUTFILE "(%2d)          approx=%c\n", PP_Myid, (approx ? 'A' : 'E'));
1190# endif /* PVERBOSE2 */
1191
1192  ints[0] = approx;
1193  ulis[0] = startq;
1194  ulis[1] = numofq;
1195  ulis[2] = numofbq;
1196
1197  MPI_Address(ulis,  Dtypeaddr);
1198  MPI_Address(ints, (Dtypeaddr+1));
1199 
1200  MPI_Type_struct(2, Dtypelens, Dtypeaddr, Dtypes, &PP_QBlockSpecs);
1201  MPI_Type_commit(&PP_QBlockSpecs);
1202 
1203# ifdef PVERBOSE2
1204    FPRINTF(STDOUTFILE "(%2d) Sending: xxPP_QuartBlockSpecs(0,%lu)=%d,%d\n", PP_Myid, truenum-1, trueaddr[0], trueaddr[truenum-1]);
1205# endif /* PVERBOSE2 */
1206 
1207
1208  error = MPI_Ssend(MPI_BOTTOM, 1, PP_QBlockSpecs, PP_MyMaster, PP_QUARTBLOCKSPECS, PP_Comm);
1209# ifdef PVERBOSE3
1210    FPRINTF(STDOUTFILE "(%2d) ... Sent QuartBlockSpecs (%ld, %ld, %ld, %d)\n", PP_Myid, ulis[0], ulis[1], ulis[2], ints[0]);
1211# endif /* PVERBOSE3 */
1212
1213  MPI_Address(trueaddr, DtypeaddrRes);
1214  DtypelensRes[0] = truenum;
1215
1216  MPI_Address(bq, (DtypeaddrRes + 1));
1217  DtypelensRes[1] = numofbq;
1218  MPI_Type_struct(2, DtypelensRes, DtypeaddrRes, DtypesRes, &PP_QBlockRes);
1219  MPI_Type_commit(&PP_QBlockRes);
1220
1221  error = MPI_Ssend(MPI_BOTTOM, 1, PP_QBlockRes, PP_MyMaster, PP_QUARTBLOCK, PP_Comm);
1222  if (error != MPI_SUCCESS)
1223    PP_Printerror(STDOUT, PP_Myid, error);
1224 
1225  MPI_Type_free(&PP_QBlockSpecs);
1226  MPI_Type_free(&PP_QBlockRes);
1227# ifdef PVERBOSE3
1228    FPRINTF(STDOUTFILE "(%2d) ... Sent xxPP_QuartBlock(0,%lu)=%d,%d\n", PP_Myid, truenum-1, trueaddr[0], trueaddr[truenum-1]);
1229# endif /* PVERBOSE3 */
1230 
1231# undef NUMULI
1232# undef NUMINT
1233} /* PP_SendQuartBlock */
1234 
1235 
1236
1237/******************/
1238
1239void PP_RecvQuartBlock(int slave,
1240                       uli *startq,
1241                       uli *numofq,
1242                       unsigned char *quartetinfo,
1243                       int *approx)
1244{
1245#       define NUMULI 3
1246#       define NUMINT 1
1247        unsigned char *trueaddr;
1248        uli            truenum;
1249        int            error;
1250        int            dest;
1251        int    ints[NUMINT];
1252        uli    ulis[NUMULI];
1253        MPI_Datatype Dtypes[2] =    {MPI_UNSIGNED_LONG, MPI_INT};
1254        int          Dtypelens[2] = {NUMULI,             NUMINT};
1255        MPI_Aint     Dtypeaddr[2];
1256        MPI_Datatype PP_QBlockSpecs;
1257        MPI_Datatype DtypesRes[2] =    {MPI_UNSIGNED_CHAR, MPI_UNSIGNED_LONG};
1258        int          DtypelensRes[2];
1259        MPI_Aint     DtypeaddrRes[2];
1260        MPI_Datatype PP_QBlockRes;
1261        MPI_Status   stat;
1262        uli          count;
1263uli num;
1264uli *numofbq;
1265uli *bq;
1266numofbq=&num;
1267#       ifdef PVERBOSE3
1268            FPRINTF(STDOUTFILE "(%2d) Receiving QuartBlock ...\n", PP_Myid);
1269#       endif /* PVERBOSE3 */
1270        MPI_Address(ulis,  Dtypeaddr);
1271        MPI_Address(ints, (Dtypeaddr+1));
1272 
1273        MPI_Type_struct(2, Dtypelens, Dtypeaddr, Dtypes, &PP_QBlockSpecs);
1274        MPI_Type_commit(&PP_QBlockSpecs);
1275 
1276        MPI_Probe(MPI_ANY_SOURCE, PP_QUARTBLOCKSPECS, PP_Comm, &stat);
1277        dest = stat.MPI_SOURCE;
1278        error = MPI_Recv(MPI_BOTTOM, 1, PP_QBlockSpecs, dest, PP_QUARTBLOCKSPECS, PP_Comm, &stat);
1279        if (error != MPI_SUCCESS)
1280              PP_Printerror(STDOUT, PP_Myid, error);
1281
1282        *approx  = ints[0];
1283        *startq  = ulis[0];
1284        *numofq  = ulis[1];
1285        *numofbq = ulis[2];
1286
1287        PP_quartrecved += *numofq;
1288        PP_quartrecvedn++;
1289        truenum = (uli)((*numofq+1)/2);
1290        trueaddr = (unsigned char *)(quartetinfo + (uli)(*startq/2));
1291#       ifdef PVERBOSE3
1292           FPRINTF(STDOUTFILE "(%2d) ... Recv QuartBlockSpecs (%ld, %ld, %ld, %d)\n", PP_Myid, ulis[0], ulis[1], ulis[2], ints[0]);
1293#       endif /* PVERBOSE3 */
1294
1295        DtypelensRes[0] =  truenum;
1296        MPI_Address(trueaddr,  DtypeaddrRes);
1297
1298        bq = malloc((unsigned) *numofbq * sizeof(uli));
1299
1300        DtypelensRes[1] = *numofbq;
1301        MPI_Address(bq, (DtypeaddrRes+1));
1302        MPI_Type_struct(2, DtypelensRes, DtypeaddrRes, DtypesRes, &PP_QBlockRes);
1303        MPI_Type_commit(&PP_QBlockRes);
1304 
1305        error = MPI_Recv(MPI_BOTTOM, 1, PP_QBlockRes, dest, PP_QUARTBLOCK, PP_Comm, &stat);
1306        if (error != MPI_SUCCESS)
1307            PP_Printerror(STDOUT, PP_Myid, error);
1308#       ifdef PVERBOSE3
1309           FPRINTF(STDOUTFILE "(%2d) ... Recv QuartBlock \n", PP_Myid);
1310#       endif /* PVERBOSE3 */
1311
1312        PP_putslave(dest);
1313
1314        for(count = 0; count < *numofbq; count++){
1315          int a, b, c, d;
1316          num2quart(bq[count], &a, &b, &c, &d);
1317#         ifdef PVERBOSE2
1318             FPRINTF(STDOUTFILE "(%2d) %ld. bad quarted (%d, %d, %d, %d) = %ld\n", PP_Myid, count, a, b, c, d, bq[count]);
1319#         endif /* PVERBOSE2 */
1320
1321          badqs++;
1322          badtaxon[a]++;
1323          badtaxon[b]++;
1324          badtaxon[c]++;
1325          badtaxon[d]++;
1326          if (show_optn) {
1327             fputid10(unresfp, a);
1328             fprintf(unresfp, "  ");
1329             fputid10(unresfp, b);
1330             fprintf(unresfp, "  ");
1331             fputid10(unresfp, c);
1332             fprintf(unresfp, "  ");
1333             fputid(unresfp, d);
1334             fprintf(unresfp, "\n");
1335          }
1336        }
1337        free(bq);
1338        MPI_Type_free(&PP_QBlockSpecs);
1339        MPI_Type_free(&PP_QBlockRes);
1340#       ifdef PVERBOSE2
1341           FPRINTF(STDOUTFILE "(%2d) <- (%2d) ... Recv xxPP_QuartBlock(0,%lu)=%d,%d\n", PP_Myid, dest, truenum-1, trueaddr[0], trueaddr[truenum-1]);
1342#       endif /* PVERBOSE2 */
1343 
1344#       undef NUMULI
1345#       undef NUMINT
1346} /* PP_RecvQuartBlock */
1347 
1348
1349/*********************************************************************
1350*  send/receive array with all quartets  (M->S)                      *
1351*********************************************************************/
1352
1353void PP_SendAllQuarts(unsigned long  Numquartets,
1354                      unsigned char *local_quartetinfo)
1355{
1356  MPI_Datatype  Dtypes[1] =    {MPI_UNSIGNED_CHAR};
1357  int           Dtypelens[1];
1358  MPI_Aint      Dtypeaddr[1];
1359  MPI_Datatype  PP_AllQuarts;
1360  int           dest;
1361  int           error;
1362 
1363# ifdef PVERBOSE2
1364    FPRINTF(STDOUTFILE "(%2d) Sending: PP_AllQuart(0)=%d\n", PP_Myid, local_quartetinfo[0]);
1365# endif /* PVERBOSE2 */
1366 
1367  /* compute number of quartets */
1368  if (Numquartets % 2 == 0) { /* even number */
1369     Dtypelens[0] = (Numquartets)/2;
1370  } else { /* odd number */
1371     Dtypelens[0] = (Numquartets + 1)/2;
1372  }
1373
1374  MPI_Address(&(local_quartetinfo[0]), Dtypeaddr);
1375  MPI_Type_struct(1, Dtypelens, Dtypeaddr, Dtypes, &PP_AllQuarts);
1376  MPI_Type_commit(&PP_AllQuarts);
1377 
1378  for (dest=1; dest<PP_NumProcs; dest++) {
1379 
1380    error = MPI_Ssend(MPI_BOTTOM, 1, PP_AllQuarts, dest, PP_ALLQUARTS, PP_Comm);
1381    if (error != MPI_SUCCESS)
1382      PP_Printerror(STDOUT, 1200+PP_Myid, error);
1383 
1384#   ifdef PVERBOSE3
1385       FPRINTF(STDOUTFILE "(%2d) -> (%2d) ... Sent xxAllQuart(0,%d)=%d,%d (%luq -> %db)\n", 
1386                           PP_Myid, dest, Dtypelens[0]-1, local_quartetinfo[0], local_quartetinfo[Dtypelens[0]-1],
1387                           Numquartets, Dtypelens[0]-1);
1388#   endif /* PVERBOSE3 */
1389  } /* for each slave */
1390 
1391  MPI_Type_free(&PP_AllQuarts);
1392 
1393 
1394} /* PP_SendAllQuarts */
1395 
1396 
1397
1398/******************/
1399
1400void PP_RecvAllQuarts(int            taxa,
1401                      unsigned long *Numquartets,
1402                      unsigned char *local_quartetinfo)
1403{
1404  MPI_Datatype  Dtypes[1] =    {MPI_UNSIGNED_CHAR};
1405  int           Dtypelens[1];
1406  MPI_Aint      Dtypeaddr[1];
1407  MPI_Datatype  PP_AllQuarts;
1408  MPI_Status    stat;
1409  int           error;
1410 
1411# ifdef PVERBOSE3
1412    FPRINTF(STDOUTFILE "(%2d) Receiving AllQuarts ...\n", PP_Myid);
1413# endif /* PVERBOSE3 */
1414 
1415  /* compute number of quartets */
1416  *Numquartets = (uli) taxa*(taxa-1)*(taxa-2)*(taxa-3)/24;
1417  if (*Numquartets % 2 == 0) { /* even number */
1418    Dtypelens[0] = (*Numquartets)/2;
1419  } else { /* odd number */
1420    Dtypelens[0] = (*Numquartets + 1)/2;
1421  }
1422 
1423  MPI_Address(&(local_quartetinfo[0]), Dtypeaddr);
1424  MPI_Type_struct(1, Dtypelens, Dtypeaddr, Dtypes, &PP_AllQuarts);
1425  MPI_Type_commit(&PP_AllQuarts);
1426 
1427  error = MPI_Recv(MPI_BOTTOM, 1, PP_AllQuarts, PP_MyMaster, PP_ALLQUARTS, PP_Comm, &stat);
1428  if (error != MPI_SUCCESS)
1429    PP_Printerror(STDOUT, 1300+PP_Myid, error);
1430
1431  MPI_Type_free(&PP_AllQuarts);
1432 
1433# ifdef PVERBOSE2
1434     FPRINTF(STDOUTFILE "(%2d) <- (%2d) ... Recv xxAllQuart(0,%d)=%d,%d (%luq -> %db)\n", 
1435                        PP_Myid, PP_MyMaster, Dtypelens[0]-1, local_quartetinfo[0], local_quartetinfo[Dtypelens[0]-1],
1436                        *Numquartets, Dtypelens[0]-1);
1437# endif /* PVERBOSE2 */
1438 
1439} /* PP_RecvAllQuarts */
1440 
1441
1442
1443/*********************************************************************
1444*  procedures to send request for a single puzzle tree               *
1445*********************************************************************/
1446
1447void PP_SendPermut(int      dest,
1448                   int      taxa, 
1449                   ivector  permut)
1450{
1451  MPI_Datatype  Dtypes[1] =    {MPI_INT};
1452  int           Dtypelens[1];
1453  MPI_Aint      Dtypeaddr[1];
1454  MPI_Datatype  PP_Permut;
1455  int           error;
1456 
1457  PP_permutsent++;
1458  PP_permutsentn++;
1459  Dtypelens[0] = taxa;
1460 
1461  MPI_Address(&(permut[0]), Dtypeaddr);
1462  MPI_Type_struct(1, Dtypelens, Dtypeaddr, Dtypes, &PP_Permut);
1463  MPI_Type_commit(&PP_Permut);
1464 
1465# ifdef PVERBOSE2
1466    FPRINTF(STDOUTFILE "(%2d) Sending  -> (%2d): PP_Permut(0)=%d\n", PP_Myid, dest, permut[0]);
1467# endif /* PVERBOSE2 */
1468 
1469  error = MPI_Ssend(MPI_BOTTOM, 1, PP_Permut, dest, PP_DOPUZZLE, PP_Comm);
1470  if (error != MPI_SUCCESS)
1471    PP_Printerror(STDOUT, 1500+PP_Myid, error);
1472 
1473  MPI_Type_free(&PP_Permut);
1474 
1475# ifdef PVERBOSE3
1476    FPRINTF(STDOUTFILE "(%2d) ... Sent PP_Permut\n", PP_Myid);
1477# endif /* PVERBOSE3 */
1478 
1479} /* PP_SendPermut */
1480 
1481/******************/
1482 
1483void PP_RecvPermut(int      taxa,
1484                   ivector  permut)
1485{
1486  MPI_Datatype  Dtypes[1] =    {MPI_INT};
1487  int           Dtypelens[1];
1488  MPI_Aint      Dtypeaddr[1];
1489  MPI_Datatype  PP_Permut;
1490  MPI_Status    stat;
1491  int           error;
1492 
1493  PP_permutrecved++;
1494  PP_permutrecvedn++;
1495# ifdef PVERBOSE3
1496    FPRINTF(STDOUTFILE "(%2d) Receiving: PP_Permut\n", PP_Myid);
1497# endif /* PVERBOSE3 */
1498 
1499  Dtypelens[0] = taxa;
1500 
1501  MPI_Address(&(permut[0]), Dtypeaddr);
1502  MPI_Type_struct(1, Dtypelens, Dtypeaddr, Dtypes, &PP_Permut);
1503  MPI_Type_commit(&PP_Permut);
1504 
1505  error = MPI_Recv(MPI_BOTTOM, 1, PP_Permut, PP_MyMaster, PP_DOPUZZLE, PP_Comm, &stat);
1506  if (error != MPI_SUCCESS)
1507    PP_Printerror(STDOUT, 1700+PP_Myid, error);
1508 
1509  MPI_Type_free(&PP_Permut);
1510 
1511# ifdef PVERBOSE2
1512    FPRINTF(STDOUTFILE "(%2d) Received: PP_Permut(0)=%d\n", PP_Myid, permut[0]);
1513# endif /* PVERBOSE2 */
1514 
1515} /* PP_RecvPermut */
1516
1517/*********************************************************************
1518*  procedures to send the splits of a puzzle tree to the master      *
1519*********************************************************************/
1520
1521void PP_SendSplitsBlock(int               taxa, 
1522                        uli               blocksize,
1523                        cmatrix          *biparts,
1524                        int               pstnum,
1525                        treelistitemtype *pstlist)
1526{
1527  MPI_Datatype     *Dtypes;
1528  int              *Dtypelens;
1529  MPI_Aint         *Dtypeaddr;
1530  MPI_Datatype      PP_Biparts;
1531  int               error;
1532  int               n;
1533  int               ints[3];
1534  int              *pstnumarr;
1535  treelistitemtype *pstptr;
1536
1537  PP_splitsent+=blocksize;
1538  PP_splitsentn++;
1539 
1540  ints[0] = taxa;
1541  ints[1] = (int) blocksize;
1542  ints[2] = pstnum;
1543  error = MPI_Ssend(ints, 3, MPI_INT, PP_MyMaster, PP_PUZZLEBLOCKSPECS, PP_Comm);
1544  if (error != MPI_SUCCESS)
1545    PP_Printerror(STDOUT, 1800+PP_Myid, error);
1546 
1547  Dtypes    = malloc((blocksize + pstnum + 1) * sizeof(MPI_Datatype));
1548  Dtypelens = malloc((blocksize + pstnum + 1) * sizeof(int));
1549  Dtypeaddr = malloc((blocksize + pstnum + 1) * sizeof(MPI_Aint));
1550  pstnumarr = malloc(pstnum * sizeof(int));
1551
1552# ifdef PVERBOSE2
1553    FPRINTF(STDOUTFILE "(%2d) Sending: PP_bipartsblock(0..%lu,0,0)8=\"%c\"\n", PP_Myid, blocksize, biparts[0][0][0]);
1554# endif /* PVERBOSE2 */
1555 
1556  for (n=0; n<(int)blocksize; n++) {
1557    Dtypes[n]    = MPI_CHAR;
1558    Dtypelens[n] = (taxa - 3) * taxa;
1559    MPI_Address(&(biparts[n][0][0]), &(Dtypeaddr[n]));
1560  }
1561  pstptr = pstlist;
1562  for (n=0; n<pstnum; n++) {
1563    Dtypes[(int)blocksize + n]    = MPI_CHAR;
1564    Dtypelens[(int)blocksize + n] = psteptreestrlen;
1565    MPI_Address((*pstptr).tree, &(Dtypeaddr[(int)blocksize + n]));
1566    pstnumarr[n] = (*pstptr).count;
1567#   ifdef PVERBOSE3
1568       FPRINTF(STDOUTFILE "(%2d) Sent tree item ->%d: [%d/%d] #=%d  \"%s\"\n",
1569               PP_Myid, PP_MyMaster, n, pstnum, pstnumarr[n], (*pstptr).tree);
1570#   endif /* PVERBOSE3 */
1571    pstptr = (*pstptr).succ;
1572  }
1573  Dtypes[((int)blocksize + pstnum)]    = MPI_INT;
1574  Dtypelens[((int)blocksize + pstnum)] = pstnum;
1575  MPI_Address(&(pstnumarr[0]), &(Dtypeaddr[((int)blocksize + pstnum)]));
1576
1577  MPI_Type_struct(((int)blocksize + pstnum + 1), Dtypelens, Dtypeaddr, Dtypes, &PP_Biparts);
1578  MPI_Type_commit(&PP_Biparts);
1579 
1580  error = MPI_Ssend(MPI_BOTTOM, 1, PP_Biparts, PP_MyMaster, PP_PUZZLEBLOCK, PP_Comm);
1581  if (error != MPI_SUCCESS)
1582    PP_Printerror(STDOUT, 1800+PP_Myid, error);
1583 
1584  MPI_Type_free(&PP_Biparts);
1585  free(Dtypes);
1586  free(Dtypelens);
1587  free(Dtypeaddr);
1588  free(pstnumarr);
1589 
1590# ifdef PVERBOSE3
1591    FPRINTF(STDOUTFILE "(%2d) ... Sent PP_bipartsblock\n", PP_Myid);
1592# endif /* PVERBOSE3 */
1593 
1594} /* PP_SendSplitsBlock */
1595 
1596/******************/
1597
1598void PP_RecvSplitsBlock(int               *taxa, 
1599                        uli               *blocksize,
1600                        cmatrix          **bip,
1601                        treelistitemtype **pstlist,
1602                        int               *pstnum,
1603                        int               *pstsum)
1604/* bp -> (*bip) */
1605{
1606  MPI_Datatype      *Dtypes;
1607  int               *Dtypelens;
1608  MPI_Aint          *Dtypeaddr;
1609  MPI_Datatype       PP_Biparts;
1610  MPI_Status         stat;
1611  int                error;
1612  int                n;
1613  int                dest;
1614  int                ints[3];
1615  int                pstlistnum;
1616  int                tmpnum;
1617  int                tmpsum;
1618  int               *pstnumarr;
1619  char             **pstarr;
1620  treelistitemtype  *treeitem;
1621
1622  error = MPI_Recv(ints, 3, MPI_INT, MPI_ANY_SOURCE, PP_PUZZLEBLOCKSPECS, PP_Comm, &stat);
1623  if (error != MPI_SUCCESS)
1624    PP_Printerror(STDOUT, 1900+PP_Myid, error);
1625
1626    dest       =       stat.MPI_SOURCE; 
1627   *taxa       =       ints[0];
1628   *blocksize  = (uli) ints[1];
1629    pstlistnum =       ints[2];
1630 
1631# ifdef PVERBOSE2
1632    FPRINTF(STDOUTFILE "(%2d) Received<-%d: PP_bipartsblockspec(t=%d,b=%ld,p=%d)\n", PP_Myid, dest, *taxa, *blocksize, pstlistnum);
1633# endif /* PVERBOSE2 */
1634
1635  PP_splitrecved += *blocksize;
1636  PP_splitrecvedn++;
1637 
1638  Dtypes    = malloc((*blocksize + pstlistnum + 1) * sizeof(MPI_Datatype));
1639  Dtypelens = malloc((*blocksize + pstlistnum + 1) * sizeof(int));
1640  Dtypeaddr = malloc((*blocksize + pstlistnum + 1) * sizeof(MPI_Aint));
1641  (*bip)    = (cmatrix *) malloc(*blocksize * sizeof(void *));
1642  pstnumarr = (int *) malloc(pstlistnum * sizeof(int));
1643  pstarr    = (char **) malloc(pstlistnum * sizeof(char *));
1644 
1645/*  pstarr[0] = (char *) malloc(psteptreestrlen * pstlistnum * sizeof(char));
1646    for(n=1; n<pstlistnum; n++)
1647    pstarr[n] = &(pstarr[0][n * psteptreestrlen]);
1648*/
1649
1650  for (n=0; n<(int)*blocksize; n++) {
1651    (*bip)[n]        = new_cmatrix(*taxa - 3, *taxa);
1652    Dtypes[n]    = MPI_CHAR;
1653    Dtypelens[n] = (*taxa - 3) * *taxa;
1654    MPI_Address(&((*bip)[n][0][0]), &(Dtypeaddr[n]));
1655  }
1656  for (n=0; n<pstlistnum; n++) {
1657    pstarr[n]        = (char *)malloc(psteptreestrlen * sizeof(char)); 
1658    Dtypes[(int)*blocksize + n]    = MPI_CHAR;
1659    Dtypelens[(int)*blocksize + n] = psteptreestrlen;
1660    MPI_Address(&(pstarr[n][0]), &(Dtypeaddr[(int)*blocksize + n]));
1661  }
1662 
1663  Dtypes[(int)*blocksize + pstlistnum]    = MPI_INT;
1664  Dtypelens[(int)*blocksize + pstlistnum] = pstlistnum;
1665  MPI_Address(&(pstnumarr[0]), &(Dtypeaddr[(int)*blocksize + pstlistnum]));
1666
1667  MPI_Type_struct(((int)*blocksize + pstlistnum + 1), Dtypelens, Dtypeaddr, Dtypes, &PP_Biparts);
1668  MPI_Type_commit(&PP_Biparts);
1669 
1670  error = MPI_Recv(MPI_BOTTOM, 1, PP_Biparts, dest, PP_PUZZLEBLOCK, PP_Comm, &stat);
1671  if (error != MPI_SUCCESS)
1672    PP_Printerror(STDOUT, 1910+PP_Myid, error);
1673
1674  tmpsum = *pstsum;
1675  for (n=0; n<pstlistnum; n++) tmpsum += pstnumarr[n];
1676  tmpnum = *pstnum;
1677
1678  for (n=0; n<pstlistnum; n++) {
1679#   ifdef PVERBOSE3
1680       FPRINTF(STDOUTFILE "(%2d) Received tree item <-%d: [%d/%d] #=%d  \"%s\"\n",
1681               PP_Myid, dest, n, pstlistnum, pstnumarr[n], pstarr[n]);
1682#   endif /* PVERBOSE3 */
1683
1684    treeitem = addtree2list(&(pstarr[n]), pstnumarr[n], pstlist, pstnum, pstsum);
1685    if((listqptrees == PSTOUT_LIST) 
1686      || (listqptrees == PSTOUT_LISTORDER)) {
1687       /* print: order no/# topol per this id/tree id/sum of unique topologies/sum of trees so far */
1688       fprintf(qptlist, "%d.\t%d\t%d\t%d\t%d\t%d\n", 
1689              PP_splitrecvedn, pstnumarr[n], (*treeitem).count, (*treeitem).id, *pstnum, tmpsum);
1690    }
1691
1692
1693  }
1694
1695  MPI_Type_free(&PP_Biparts);
1696  free(Dtypes);
1697  free(Dtypelens);
1698  free(Dtypeaddr);
1699  free(pstnumarr);
1700  free(pstarr);
1701 
1702# ifdef PVERBOSE2
1703    FPRINTF(STDOUTFILE "(%2d) Received<-%d: PP_bipartsblock(0..%lu,0,0):=%c\n", PP_Myid, dest, *blocksize, (*bip)[0][0][0]);
1704# endif /* PVERBOSE2 */
1705 
1706  PP_putslave(dest); 
1707 
1708  /* MPI_Type_free(&PP_Biparts); */
1709
1710} /* PP_RecvSplitsBlock */
1711 
1712/******************/
1713
1714void PP_SendSplits(int      taxa, 
1715                   cmatrix  biparts)
1716{
1717  MPI_Datatype  Dtypes[1] =    {MPI_CHAR};
1718  int           Dtypelens[1];
1719  MPI_Aint      Dtypeaddr[1];
1720  MPI_Datatype  PP_Biparts;
1721  int           error;
1722
1723  PP_splitsent++;
1724  PP_splitsentn++;
1725 
1726# ifdef PVERBOSE2
1727    FPRINTF(STDOUTFILE "(%2d) Sending: PP_biparts(0,0)=%c\n", PP_Myid, biparts[0][0]);
1728# endif /* PVERBOSE2 */
1729 
1730  Dtypelens[0] = (taxa - 3) * taxa;
1731
1732  MPI_Address(&(biparts[0][0]), Dtypeaddr);
1733  MPI_Type_struct(1, Dtypelens, Dtypeaddr, Dtypes, &PP_Biparts);
1734  MPI_Type_commit(&PP_Biparts);
1735 
1736  error = MPI_Ssend(MPI_BOTTOM, 1, PP_Biparts, PP_MyMaster, PP_PUZZLE, PP_Comm);
1737  if (error != MPI_SUCCESS)
1738    PP_Printerror(STDOUT, 1800+PP_Myid, error);
1739 
1740  MPI_Type_free(&PP_Biparts);
1741 
1742# ifdef PVERBOSE3
1743    FPRINTF(STDOUTFILE "(%2d) ... Sent PP_biparts\n", PP_Myid);
1744# endif /* PVERBOSE3 */
1745 
1746} /* PP_SendSplits */
1747 
1748/******************/
1749
1750void PP_RecvSplits(int      taxa,
1751                   cmatrix  biparts)
1752{
1753  MPI_Datatype  Dtypes[1] =    {MPI_CHAR};
1754  int           Dtypelens[1];
1755  MPI_Aint      Dtypeaddr[1];
1756  MPI_Datatype  PP_Biparts;
1757  MPI_Status    stat;
1758  int           error;
1759
1760  PP_splitrecved++;
1761  PP_splitrecvedn++;
1762 
1763# ifdef PVERBOSE3
1764    FPRINTF(STDOUTFILE "(%2d) Receiving: PP_biparts\n", PP_Myid);
1765# endif /* PVERBOSE3 */
1766 
1767  Dtypelens[0] = (taxa - 3) * taxa;
1768 
1769  MPI_Address(&(biparts[0][0]), Dtypeaddr);
1770  MPI_Type_struct(1, Dtypelens, Dtypeaddr, Dtypes, &PP_Biparts);
1771  MPI_Type_commit(&PP_Biparts);
1772 
1773  error = MPI_Recv(MPI_BOTTOM, 1, PP_Biparts, MPI_ANY_SOURCE, PP_PUZZLE, PP_Comm, &stat);
1774  if (error != MPI_SUCCESS)
1775    PP_Printerror(STDOUT, 1920+PP_Myid, error);
1776
1777  PP_putslave(stat.MPI_SOURCE); 
1778 
1779  MPI_Type_free(&PP_Biparts);
1780 
1781# ifdef PVERBOSE2
1782    FPRINTF(STDOUTFILE "(%2d) Received <- (%2d): PP_biparts(0,0)=%c\n", PP_Myid, stat.MPI_SOURCE, biparts[0][0]);
1783# endif /* PVERBOSE2 */
1784 
1785} /* PP_RecvSplits */
1786 
1787/*********************************************************************
1788*  procedures to send the request to compute puzzle tree             *
1789*********************************************************************/
1790
1791void PP_SendDoPermutBlock(uli puzzlings)
1792{
1793  int           dest;
1794  int           error;
1795  uli           numpuzzlings;
1796  uli           puzzlingssent = 0;
1797  schedtype     sched; 
1798  int           bipnum;
1799  int           tx;
1800  uli           bs;
1801  cmatrix      *bp;
1802
1803
1804 
1805  initsched(&sched, puzzlings, PP_NumProcs-1, 4);
1806  /* numpuzzlings = sc(&sched); */
1807  numpuzzlings = sgss(&sched);
1808  while (numpuzzlings > 0) {
1809     if (PP_emptyslave()) {
1810        PP_RecvSplitsBlock(&tx, &bs, &bp, &psteptreelist, &psteptreenum, &psteptreesum);
1811        for (bipnum=0; bipnum<bs; bipnum++) {
1812               inittree();
1813               biparts = bp[bipnum];
1814               makenewsplitentries();
1815               freetree();
1816               free_cmatrix(bp[bipnum]);
1817        }
1818        free(bp);
1819        puzzlingssent -= bs;
1820     }
1821
1822     dest = PP_getslave();
1823
1824#    ifdef PVERBOSE2
1825        FPRINTF(STDOUTFILE "(%2d) Sending: DOPuzzleBlock Signal\n", PP_Myid);
1826#    endif /* PVERBOSE2 */
1827
1828     error = MPI_Ssend(&numpuzzlings, 1, MPI_UNSIGNED_LONG, dest, PP_DOPUZZLEBLOCK, PP_Comm);
1829     if (error != MPI_SUCCESS)
1830        PP_Printerror(STDOUT, 2100+PP_Myid, error);
1831
1832#    ifdef PVERBOSE3
1833        FPRINTF(STDOUTFILE "(%2d) ... Sent DOPuzzleBlock Signal\n", PP_Myid);
1834#    endif /* PVERBOSE3 */
1835
1836     puzzlingssent += numpuzzlings;
1837     PP_permutsent += numpuzzlings;
1838     PP_permutsentn ++;
1839 
1840     numpuzzlings = sgss(&sched);
1841
1842  } /* while */
1843 
1844  while (puzzlingssent > 0) {
1845        PP_RecvSplitsBlock(&tx, &bs, &bp, &psteptreelist, &psteptreenum, &psteptreesum);
1846        for (bipnum=0; bipnum<bs; bipnum++) {
1847               inittree();
1848               biparts = bp[bipnum];
1849               makenewsplitentries();
1850               freetree();
1851               free_cmatrix(bp[bipnum]);
1852        }
1853        free(bp);
1854        puzzlingssent -= bs;
1855  }
1856 
1857} /* PP_SendDoPermutBlock */
1858
1859/******************/
1860
1861
1862void PP_RecvDoPermutBlock(uli *taxa)
1863{
1864  MPI_Status    stat;
1865  int           error;
1866 
1867# ifdef PVERBOSE2
1868    FPRINTF(STDOUTFILE "(%2d) Receiving: DOPuzzleBlock Signal\n", PP_Myid);
1869# endif /* PVERBOSE2 */
1870 
1871  error = MPI_Recv(taxa, 1, MPI_UNSIGNED_LONG, PP_MyMaster, PP_DOPUZZLEBLOCK, PP_Comm, &stat);
1872  if (error != MPI_SUCCESS)
1873    PP_Printerror(STDOUT, 2100+PP_Myid, error);
1874
1875  PP_permutrecved += *taxa;
1876  PP_permutrecvedn ++;
1877
1878# ifdef PVERBOSE3
1879    FPRINTF(STDOUTFILE "(%2d) ... DOPuzzleBlock Signal\n", PP_Myid);
1880# endif /* PVERBOSE3 */
1881 
1882} /* PP_RecvDoPermutBlock */
1883
1884/*********************************************************************
1885*  procedures to make the slaves to finalize                         *
1886*********************************************************************/
1887
1888void PP_SendDone()
1889{
1890# define NUMINT 8
1891# define NUMDBL 6
1892  int           dest;
1893  int           error;
1894  MPI_Status    stat;
1895  int    ints[NUMINT];
1896  double doubles[NUMDBL];
1897  MPI_Datatype Dtypes[2] =    {MPI_INT, MPI_DOUBLE};
1898  int          Dtypelens[2] = {NUMINT , NUMDBL};
1899  MPI_Aint     Dtypeaddr[2];
1900  MPI_Datatype PP_Stats;
1901 
1902# ifdef PVERBOSE2
1903    FPRINTF(STDOUTFILE "(%2d) Sending: DONE Signal\n", PP_Myid);
1904# endif /* PVERBOSE2 */
1905 
1906  for (dest=1; dest<PP_NumProcs; dest++) {
1907 
1908    error = MPI_Ssend(&dest, 0, MPI_INT, dest, PP_DONE, PP_Comm);
1909    if (error != MPI_SUCCESS)
1910      PP_Printerror(STDOUT, 2100+PP_Myid, error);
1911 
1912  } /* for each slave */
1913 
1914# ifdef PVERBOSE3
1915    FPRINTF(STDOUTFILE "(%2d) ... Sent DONE Signal\n", PP_Myid);
1916# endif /* PVERBOSE3 */
1917 
1918  MPI_Address(ints,     Dtypeaddr);
1919  MPI_Address(doubles, (Dtypeaddr+1));
1920
1921  MPI_Type_struct(2, Dtypelens, Dtypeaddr, Dtypes, &PP_Stats);
1922  MPI_Type_commit(&PP_Stats);
1923
1924  doquartrecved[0]  = 0;
1925  doquartrecvedn[0] = 0;
1926  quartsent[0]      = 0;
1927  quartsentn[0]     = 0;
1928  permutrecved[0]   = 0;
1929  permutrecvedn[0]  = 0;
1930  splitsent[0]      = 0;
1931  splitsentn[0]     = 0;
1932  cputimes[0]       = 0.0;
1933  walltimes[0]      = 0.0;
1934  fullcputimes[0]   = 0.0;
1935  fullwalltimes[0]  = 0.0;
1936  altcputimes[0]    = 0.0;
1937  altwalltimes[0]   = 0.0;
1938
1939  for (dest=1; dest<PP_NumProcs; dest++) {
1940
1941     error = MPI_Recv(MPI_BOTTOM, 1, PP_Stats, dest, PP_STATS, PP_Comm, &stat);
1942     if (error != MPI_SUCCESS)
1943       PP_Printerror(STDOUT, 2100+PP_Myid, error);
1944
1945     doquartrecved[dest]  = ints[0];
1946     doquartrecvedn[dest] = ints[1];
1947     quartsent[dest]      = ints[2];
1948     quartsentn[dest]     = ints[3];
1949     permutrecved[dest]   = ints[4];
1950     permutrecvedn[dest]  = ints[5];
1951     splitsent[dest]      = ints[6];
1952     splitsentn[dest]     = ints[7];
1953     cputimes[dest]       = doubles[0];
1954     walltimes[dest]      = doubles[1];
1955     fullcputimes[dest]   = doubles[2];
1956     fullwalltimes[dest]  = doubles[3];
1957     altcputimes[dest]    = doubles[4];
1958     altwalltimes[dest]   = doubles[5];
1959
1960#    ifdef PVERBOSE1
1961        FPRINTF(STDOUTFILE "(%2d) ... Stats received (%d/%d, %d/%d, %d/%d, %d/%d, %f, %f, %f, %f, %f, %f)\n", 
1962        PP_Myid, doquartrecved[dest], doquartrecvedn[dest], quartsent[dest], quartsentn[dest],
1963                 permutrecved[dest], permutrecvedn[dest], splitsent[dest], splitsentn[dest],
1964                 cputimes[dest], walltimes[dest], fullcputimes[dest], fullwalltimes[dest], 
1965                 altcputimes[dest], altwalltimes[dest]);
1966#    endif /* PVERBOSE1 */
1967  } /* for each slave */
1968 
1969  MPI_Type_free(&PP_Stats);
1970
1971# undef NUMINT
1972# undef NUMDBL
1973} /* PP_SendDone */
1974
1975/******************/
1976
1977
1978
1979void PP_RecvDone()
1980{
1981# define NUMINT 8
1982# define NUMDBL 6
1983
1984  int           dummy;
1985  MPI_Status    stat;
1986  int           error;
1987  int    ints[NUMINT];
1988  double doubles[NUMDBL];
1989  MPI_Datatype Dtypes[2] =    {MPI_INT, MPI_DOUBLE};
1990  int          Dtypelens[2] = {NUMINT , NUMDBL};
1991  MPI_Aint     Dtypeaddr[2];
1992  MPI_Datatype PP_Stats;
1993 
1994# ifdef PVERBOSE2
1995    FPRINTF(STDOUTFILE "(%2d) Receiving: DONE Signal\n", PP_Myid);
1996# endif /* PVERBOSE2 */
1997 
1998  error = MPI_Recv(&dummy, 0, MPI_INT, PP_MyMaster, PP_DONE, PP_Comm, &stat);
1999  if (error != MPI_SUCCESS)
2000    PP_Printerror(STDOUT, 2100+PP_Myid, error);
2001
2002# ifdef PVERBOSE3
2003    FPRINTF(STDOUTFILE "(%2d) ... DONE Signal\n", PP_Myid);
2004# endif /* PVERBOSE3 */
2005 
2006  addtimes(OVERALL, &tarr);
2007# ifdef TIMEDEBUG
2008  printtimearr(&tarr);
2009# endif /* TIMEDEBUG */
2010
2011  time(&walltimestop);
2012  walltime = difftime(walltimestop, walltimestart);
2013  cputimestop = clock();
2014  cputime = ((double) (cputimestop - cputimestart) / CLOCKS_PER_SEC);
2015
2016# ifdef PVERBOSE1
2017    FPRINTF(STDOUTFILE "(%2d) total wallclock time: %0.2f sec (= %0.2f min = %0.2f hours)\n",
2018            PP_Myid, walltime, walltime / 60, walltime / 3600);
2019    FPRINTF(STDOUTFILE "(%2d) total CPU time: %0.2f sec (= %0.2f min = %0.2f hours)\n",
2020            PP_Myid, cputime, cputime / 60, cputime / 3600);
2021# endif /* PVERBOSE1 */
2022
2023  ints[0] = PP_doquartrecved;
2024  ints[1] = PP_doquartrecvedn;
2025  ints[2] = PP_quartsent;
2026  ints[3] = PP_quartsentn;
2027  ints[4] = PP_permutrecved;
2028  ints[5] = PP_permutrecvedn;
2029  ints[6] = PP_splitsent;
2030  ints[7] = PP_splitsentn;
2031  doubles[0] = cputime;
2032  doubles[1] = walltime;
2033  doubles[2] = tarr.fullcpu;
2034  doubles[3] = tarr.fulltime;
2035  doubles[4] = tarr.cpu;
2036  doubles[5] = tarr.time;
2037
2038  MPI_Address(ints,     Dtypeaddr);
2039  MPI_Address(doubles, (Dtypeaddr+1));
2040
2041  MPI_Type_struct(2, Dtypelens, Dtypeaddr, Dtypes, &PP_Stats);
2042  MPI_Type_commit(&PP_Stats);
2043
2044  error = MPI_Ssend(MPI_BOTTOM, 1, PP_Stats, PP_MyMaster, PP_STATS, PP_Comm);
2045  if (error != MPI_SUCCESS)
2046    PP_Printerror(STDOUT, 2100+PP_Myid, error);
2047
2048  MPI_Type_free(&PP_Stats);
2049
2050# ifdef PVERBOSE3
2051        FPRINTF(STDOUTFILE "(%2d) ... Stats sent (%d/%d, %d/%d, %d/%d, %d/%d, %f, %f)\n", 
2052        PP_Myid, PP_doquartrecved, PP_doquartrecvedn, PP_quartsent, PP_quartsentn,
2053                 PP_permutrecved, PP_permutrecvedn, PP_splitsent, PP_splitsentn,
2054                 cputime, walltime, fullcputime, fullwalltime, altcputime, altwalltime);
2055# endif /* PVERBOSE3 */
2056
2057# undef NUMINT
2058# undef NUMDBL
2059} /* PP_RecvDone */
2060
2061/*********************************************************************
2062*  procedures to initialize and cleanup                              *
2063*********************************************************************/
2064
2065void PP_Init(int *argc, char **argv[])
2066{
2067  MPI_Init(argc, argv);
2068  PP_Comm = MPI_COMM_WORLD;
2069  MPI_Comm_rank(PP_Comm, &PP_Myid);
2070  MPI_Comm_size(PP_Comm, &PP_NumProcs);
2071
2072  if (PP_NumProcs <= 1) {
2073        FPRINTF(STDOUTFILE "\nHalt: TREE-PUZZLE needs at least 2 processes for a parallel run!!!\n\n");
2074        MPI_Finalize();
2075        exit(1);
2076  }
2077
2078  PP_MyMaster = 0;
2079
2080  if (PP_Myid == PP_MyMaster) {
2081    PP_IamMaster=1;
2082    PP_IamSlave =0;
2083    PP_initslavequeue();
2084
2085    permutsent     = new_ivector(PP_NumProcs);
2086    permutrecved   = new_ivector(PP_NumProcs);
2087    quartsent      = new_ivector(PP_NumProcs);
2088    quartrecved    = new_ivector(PP_NumProcs);
2089    doquartsent    = new_ivector(PP_NumProcs);
2090    doquartrecved  = new_ivector(PP_NumProcs);
2091    splitsent      = new_ivector(PP_NumProcs);
2092    splitrecved    = new_ivector(PP_NumProcs);
2093    permutsentn    = new_ivector(PP_NumProcs);
2094    permutrecvedn  = new_ivector(PP_NumProcs);
2095    quartsentn     = new_ivector(PP_NumProcs);
2096    quartrecvedn   = new_ivector(PP_NumProcs);
2097    doquartsentn   = new_ivector(PP_NumProcs);
2098    doquartrecvedn = new_ivector(PP_NumProcs);
2099    splitsentn     = new_ivector(PP_NumProcs);
2100    splitrecvedn   = new_ivector(PP_NumProcs);
2101
2102    walltimes      = new_dvector(PP_NumProcs);
2103    cputimes       = new_dvector(PP_NumProcs);
2104    fullwalltimes  = new_dvector(PP_NumProcs);
2105    fullcputimes   = new_dvector(PP_NumProcs);
2106    altwalltimes   = new_dvector(PP_NumProcs);
2107    altcputimes    = new_dvector(PP_NumProcs);
2108    { int n;
2109      for (n = 0; n<PP_NumProcs; n++){
2110        permutsent     [n] = -1; 
2111        permutrecved   [n] = -1; 
2112        quartsent      [n] = -1; 
2113        quartrecved    [n] = -1; 
2114        doquartsent    [n] = -1; 
2115        doquartrecved  [n] = -1; 
2116        splitsent      [n] = -1; 
2117        splitrecved    [n] = -1; 
2118        permutsentn    [n] = -1; 
2119        permutrecvedn  [n] = -1; 
2120        quartsentn     [n] = -1; 
2121        quartrecvedn   [n] = -1; 
2122        doquartsentn   [n] = -1; 
2123        doquartrecvedn [n] = -1; 
2124        splitsentn     [n] = -1; 
2125        splitrecvedn   [n] = -1; 
2126        walltimes      [n] = -1.0; 
2127        cputimes       [n] = -1.0; 
2128        fullwalltimes  [n] = -1.0; 
2129        fullcputimes   [n] = -1.0; 
2130        altwalltimes   [n] = -1.0; 
2131        altcputimes    [n] = -1.0; 
2132      }
2133    }
2134  } else {
2135    PP_IamMaster=0;
2136    PP_IamSlave =1;
2137  }
2138
2139# ifdef PVERBOSE1
2140    { char ma[7] = "master";
2141      char sl[7] = "slave ";
2142      char *st;
2143      char PP_ProcName[MPI_MAX_PROCESSOR_NAME] = "empty";
2144      int  flag;
2145      int  PP_ProcNameLen = 6;
2146      int  PP_Clocksync = 0;
2147      int  PP_Maxtag    = 0;
2148      int  PP_IO        = 0;
2149      int  PP_Hostexist = 0;
2150    if (PP_IamMaster)
2151      st = ma;
2152    else
2153      st = sl;
2154    FPRINTF(STDOUTFILE "(%2d) Init: ID %d, %d processes running \n", PP_Myid, PP_Myid, PP_NumProcs); 
2155    FPRINTF(STDOUTFILE "(%2d)       I am %s \n", PP_Myid, st); 
2156    MPI_Get_processor_name(PP_ProcName, &PP_ProcNameLen);
2157    FPRINTF(STDOUTFILE "(%2d)       I am running on \"%s\"\n", PP_Myid, PP_ProcName); 
2158    MPI_Attr_get(PP_Comm, MPI_IO, &PP_IO, &flag);
2159    if (flag) FPRINTF(STDOUTFILE "(%2d)       next IO Proc=%d - ", PP_Myid, PP_IO); 
2160    MPI_Attr_get(PP_Comm, MPI_TAG_UB, &PP_Maxtag, &flag);
2161    if (flag) FPRINTF(STDOUTFILE "MaxTag=%d - ", PP_Maxtag); 
2162    MPI_Attr_get(PP_Comm, MPI_HOST, &PP_Hostexist, &flag);
2163    if (flag) {
2164      if (PP_Hostexist == MPI_PROC_NULL)
2165        FPRINTF(STDOUTFILE "No HostProc - "); 
2166      else
2167        FPRINTF(STDOUTFILE "HostProc=%d - ", PP_Hostexist); 
2168    }
2169    MPI_Attr_get(PP_Comm, MPI_WTIME_IS_GLOBAL, &PP_Clocksync, &flag);
2170    if (PP_Clocksync)
2171      FPRINTF(STDOUTFILE "global time"); 
2172    else
2173      FPRINTF(STDOUTFILE "local time"); 
2174    FPRINTF(STDOUTFILE "\n");
2175    }
2176# endif /* PVERBOSE1 */
2177} /* PP_Init */
2178
2179/******************/
2180
2181void PP_Finalize()
2182{
2183  if (PP_IamMaster) {
2184    free_ivector(freeslaves);
2185
2186    free_ivector(permutsent);
2187    free_ivector(permutrecved);
2188    free_ivector(quartsent);
2189    free_ivector(quartrecved);
2190    free_ivector(doquartsent);
2191    free_ivector(doquartrecved);
2192    free_ivector(splitsent);
2193    free_ivector(splitrecved);
2194    free_ivector(permutsentn);
2195    free_ivector(permutrecvedn);
2196    free_ivector(quartsentn);
2197    free_ivector(quartrecvedn);
2198    free_ivector(doquartsentn);
2199    free_ivector(doquartrecvedn);
2200    free_ivector(splitsentn);
2201    free_ivector(splitrecvedn);
2202
2203    free_dvector(walltimes);
2204    free_dvector(cputimes);
2205    free_dvector(fullwalltimes);
2206    free_dvector(fullcputimes);
2207    free_dvector(altwalltimes);
2208    free_dvector(altcputimes);
2209  }
2210
2211# ifdef PVERBOSE1
2212    FPRINTF(STDOUTFILE "(%2d) Finished ...\n", PP_Myid);
2213
2214    FPRINTF(STDOUTFILE "(%2d) doqu(s%6d/%6d,%6d/%6d) qu(s%6d/%6d,%6d/%6d)\n",
2215            PP_Myid, PP_doquartsent, PP_doquartsentn, PP_doquartrecved, PP_doquartrecvedn, 
2216                     PP_quartsent, PP_quartsentn, PP_quartrecved, PP_quartrecvedn);
2217    FPRINTF(STDOUTFILE "(%2d) perm(s%6d/%6d,%6d/%6d) spli(s%6d/%6d,%6d/%6d)\n",
2218            PP_Myid, PP_permutsent, PP_permutsentn, PP_permutrecved, PP_permutrecvedn, 
2219                     PP_splitsent, PP_splitsentn, PP_splitrecved, PP_splitrecvedn);
2220    if (PP_IamMaster) {
2221       FPRINTF(STDOUTFILE "(%2d)  Init: %2.4f   Param(Comp: %2.4f  Send: %2.4f)\n", PP_Myid, PP_inittime, PP_paramcomptime, PP_paramsendtime);
2222       FPRINTF(STDOUTFILE "(%2d)  Quart(Comp: %2.4f  Send: %2.4f)   Puzzle: %2.4f   Tree: %2.4f\n", PP_Myid, PP_quartcomptime, PP_quartsendtime, PP_puzzletime, PP_treetime);
2223    }
2224# endif /* PVERBOSE1 */
2225
2226  MPI_Finalize();
2227}
2228
2229/***********************************************************
2230*  main part of the slave process                          *
2231***********************************************************/
2232
2233int slave_main(int argc, char *argv[])
2234{       
2235        int i, a, b, c, d, approx;
2236        double d1, d2, d3;
2237
2238        int        notdone;
2239        MPI_Status stat;
2240        int        PP_AllQuartsReceived = 0;
2241       
2242        PP_RecvSizes(&Maxspc, &Maxsite, &numcats, &Numptrn, &tpmradix, &outgroup, &fracconst, &randseed);
2243        psteptreestrlen = (Maxspc * (int)(1 + log10(Maxspc))) +
2244                          (Maxspc * 3);
2245
2246
2247        Maxbrnch = 2*Maxspc - 3;
2248
2249        initrandom(randseed);
2250
2251        /* initialise ML (from PP_mlstart) */
2252        Seqpat    = new_cmatrix(Maxspc, Numptrn);
2253        Alias     = new_ivector(Maxsite);
2254        Weight    = new_ivector(Numptrn);
2255        constpat  = new_ivector(Numptrn);
2256        Rates     = new_dvector(numcats);
2257        Eval      = new_dvector(tpmradix);
2258        Freqtpm   = new_dvector(tpmradix);
2259        Evec      = new_dmatrix(tpmradix,tpmradix);
2260        Ievc      = new_dmatrix(tpmradix,tpmradix);
2261        iexp      = new_dmatrix(tpmradix,tpmradix);
2262        Distanmat = new_dmatrix(Maxspc, Maxspc);
2263        ltprobr   = new_dcube(numcats, tpmradix,tpmradix);
2264        PP_RecvData(Seqpat,                       /* cmatrix */
2265                    Alias, Weight, constpat,      /* ivector */
2266                    Rates, Eval, Freqtpm,         /* dvector */
2267                    Evec, Ievc, iexp, Distanmat,  /* dmatrix */
2268                    ltprobr);                     /* dcube   */
2269
2270        Ctree = new_quartet(Numptrn, Seqpat);
2271        Numbrnch = NUMQBRNCH;
2272        Numibrnch = NUMQIBRNCH;
2273        Numspc = NUMQSPC;
2274
2275        /* allocate variable used for randomizing input order */
2276        trueID = new_ivector(Maxspc);
2277
2278        /* allocate and initialize badtaxonvector */
2279        badtaxon = new_ulivector(Maxspc);
2280        for (i = 0; i < Maxspc; i++) badtaxon[i] = 0;
2281
2282        /* allocate memory for quartets */
2283        quartetinfo = mallocquartets(Maxspc);
2284
2285        /* prepare for consensus tree analysis */
2286        initconsensus();
2287
2288        MPI_Probe(PP_MyMaster, MPI_ANY_TAG, PP_Comm, &stat);
2289       
2290        notdone = (stat.MPI_TAG != PP_DONE);
2291        if (!notdone) {
2292          PP_RecvDone();
2293#         ifdef TIMEDEBUG
2294             printtimearr(&tarr);
2295#         endif /* TIMEDEBUG */
2296
2297        }
2298
2299        while (notdone) {
2300          switch (stat.MPI_TAG) {
2301            case PP_ALLQUARTS: {
2302                PP_RecvAllQuarts(Maxspc, &Numquartets, quartetinfo);
2303                PP_AllQuartsReceived = 1;
2304                break;
2305                }
2306            case PP_DOQUARTBLOCK:
2307            case PP_DOQUARTBLOCKSPECS:
2308                {
2309                uli qtodo, qstart, qend, idx, nofbq, *bqarr;
2310                int a, b, c, i;
2311                double d1, d2, d3;
2312
2313                nofbq=0;
2314                PP_RecvDoQuartBlock(&qstart, &qtodo, &bqarr, &approx);
2315                qend   = (qstart + qtodo) - 1;
2316#               ifdef PVERBOSE1
2317                        FPRINTF(STDOUTFILE "(%2d) Quartets %4ld->%4ld (%dx%ld)\n", PP_Myid, qstart, qend, PP_NumProcs-1, qtodo);
2318#               endif
2319
2320                addtimes(GENERAL, &tarr);
2321                for (i = 3; i < Maxspc; i++)
2322                   for (c = 2; c < i; c++)
2323                      for (b = 1; b < c; b++)
2324                         for (a = 0; a < b; a++) {
2325
2326                            idx = (uli) a + 
2327                                  (uli) b*(b-1)/2 +
2328                                  (uli) c*(c-1)*(c-2)/6 +
2329                                  (uli) i*(i-1)*(i-2)*(i-3)/24;
2330                            if ((idx >= qstart) && (idx <= qend)) {
2331#                           ifdef PVERBOSE4
2332                               FPRINTF(STDOUTFILE "(%2d)  %4ld <---> (%d,%d,%d,%d)\n",PP_Myid, idx, a,b,c,i);
2333#                           endif
2334                               compute_quartlklhds(a,b,c,i,&d1,&d2,&d3,approx);
2335                               PP_do_write_quart(a,b,c,i,d1,d2,d3,&nofbq,bqarr);
2336                               addtimes(QUARTETS, &tarr);
2337                            } /* if idx */
2338                         } /* for for for for */   
2339                PP_SendQuartBlock(qstart, qtodo, quartetinfo, nofbq, bqarr, approx);
2340
2341                free(bqarr); bqarr=NULL;
2342
2343                break;
2344                }
2345
2346            case PP_DOPUZZLEBLOCK: {
2347                if (PP_AllQuartsReceived){
2348                   uli Numtrial, ptodo;
2349                   cmatrix *bp;
2350                   int n;
2351
2352                   PP_RecvDoPermutBlock(&Numtrial);
2353                   ptodo = Numtrial;
2354
2355                   bp = (cmatrix *) malloc(Numtrial * sizeof(void *));
2356                   for(n=0; n<Numtrial; n++) {
2357                      bp[n] = new_cmatrix(Maxspc - 3, Maxspc);
2358                   }
2359
2360                   addtimes(GENERAL, &tarr);
2361                   for (Currtrial = 0; Currtrial < ptodo; Currtrial++) {
2362                      /* randomize input order */
2363                      biparts=bp[Currtrial];
2364                      chooser(Maxspc, Maxspc, trueID);
2365
2366                      PP_slave_do_puzzling(trueID);
2367                      addtimes(PUZZLING, &tarr);
2368                   }
2369                   PP_SendSplitsBlock(Maxspc, Numtrial, bp, psteptreenum, psteptreelist);
2370                   for (Currtrial = 0; Currtrial < ptodo; Currtrial++) {
2371                      free_cmatrix(bp[Currtrial]);
2372                   }
2373                   free(bp);
2374                } else {
2375                   FPRINTF(STDOUTFILE "ERROR: Requested Puzzle Tree without receiving Quartets!!!\n");
2376                   notdone = 0;
2377                }
2378                freetreelist(&psteptreelist, &psteptreenum, &psteptreesum);
2379                break;
2380                }
2381            case PP_DOQUART: {
2382                PP_RecvDoQuart(&a,&b,&c,&d, &approx);
2383                addtimes(GENERAL, &tarr);
2384                compute_quartlklhds(a,b,c,d,&d1,&d2,&d3,approx);
2385                addtimes(QUARTETS, &tarr);
2386                PP_SendQuart(a,b,c,d,d1,d2,d3, approx);
2387                break;
2388                }
2389            case PP_DOPUZZLE: {
2390                if (PP_AllQuartsReceived){
2391                   PP_RecvPermut(Maxspc, trueID);
2392                   addtimes(GENERAL, &tarr);
2393                   PP_slave_do_puzzling(trueID);
2394                   addtimes(PUZZLING, &tarr);
2395                   }
2396                else {
2397                   FPRINTF(STDOUTFILE "ERROR: Requested Puzzle Tree without receiving Quartets!!!\n");
2398                   notdone = 0;
2399                }
2400                break;
2401                }
2402            default: {
2403                FPRINTF(STDOUTFILE "ERROR: Unknown Message Tag received\n");
2404                MPI_Abort(PP_Comm, 1);
2405                }
2406          } /* switch stat.MPI_TAG */ 
2407
2408          MPI_Probe(PP_MyMaster, MPI_ANY_TAG, PP_Comm, &stat);
2409          notdone = (stat.MPI_TAG != PP_DONE);
2410          if (!notdone)
2411            PP_RecvDone();
2412
2413        } /* while notdone */
2414
2415        return (0);
2416
2417} /* slave_main */
2418
Note: See TracBrowser for help on using the repository browser.