forester/archive/RIO/others/puzzle_mod/src/ppuzzle.c

   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
  22 int PP_IamMaster;
  23 int PP_IamSlave;
  24 int PP_Myid;
  25 int PP_MyMaster;
  26 int PP_NumProcs;
  27 MPI_Comm PP_Comm;
  28
  29 int *freeslaves;           /* Queue of free slaves */
  30 int firstslave,            /* headpointer of queue */
  31     lastslave;             /* tailpointer of queue */
  32
  33 int *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
  50 double *walltimes,
  51        *cputimes;
  52 double *fullwalltimes,
  53        *fullcputimes;
  54 double *altwalltimes,
  55        *altcputimes;
  56
  57 int PP_permutsent = 0;          /* # of  */
  58 int PP_permutrecved = 0;        /* # of  */
  59 int PP_quartsent = 0;           /* # of  */
  60 int PP_quartrecved = 0;         /* # of  */
  61 int PP_doquartsent = 0;         /* # of  */
  62 int PP_doquartrecved = 0;       /* # of  */
  63 int PP_splitsent = 0;           /* # of  */
  64 int PP_splitrecved = 0;         /* # of  */
  65 int PP_permutsentn = 0;          /* # of  */
  66 int PP_permutrecvedn = 0;        /* # of  */
  67 int PP_quartsentn = 0;           /* # of  */
  68 int PP_quartrecvedn = 0;         /* # of  */
  69 int PP_doquartsentn = 0;         /* # of  */
  70 int PP_doquartrecvedn = 0;       /* # of  */
  71 int PP_splitsentn = 0;           /* # of  */
  72 int PP_splitrecvedn = 0;         /* # of  */
  73
  74 double 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
  85 int PP_MaxSlave = 0;
  86
  87
  88 /*********************************************************************
  89 *  miscellaneous utilities                                           *
  90 *********************************************************************/
  91
  92 int 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
 101 void 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
 109 void 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
 117 double 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
 133 void 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
 152 void 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
 166 void 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
 243 uli 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
 267 uli 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
 288 void 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
 303 int PP_fullslave()
 304 {
 305 return (freeslaves[PP_MaxSlave] == PP_MaxSlave);
 306 }
 307
 308 /******************/
 309
 310 int PP_emptyslave()
 311 {
 312 return (freeslaves[PP_MaxSlave] == 0);
 313 }
 314
 315 /******************/
 316
 317 void 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
 331 int 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
 348 void PP_slave_do_puzzling(ivector trueID)
 349 {
 350 int    i, a, b, c;
 351 uli    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 continous 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
 456 void PP_do_puzzling(ivector trueID)
 457 {
 458 int    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
 484 void 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
 572 void 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
 637 void 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 matrizes  (M->S)                       *
 718 *********************************************************************/
 719
 720 void 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
 812 void 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
 893 void 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
 932 void 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
 971 void 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
1029 void 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
1089 void 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
1117 void 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
1153 void 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
1239 void 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;
1263 uli num;
1264 uli *numofbq;
1265 uli *bq;
1266 numofbq=&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
1353 void PP_SendAllQuarts(unsigned long  Numquartets,
1354                       unsigned char *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, 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(&(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, quartetinfo[0], 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
1400 void PP_RecvAllQuarts(int            taxa,
1401                       unsigned long *Numquartets,
1402                       unsigned char *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(&(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, quartetinfo[0], 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
1447 void 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
1483 void 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
1521 void 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
1598 void 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
1714 void 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
1750 void 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
1791 void 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
1862 void 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
1888 void 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
1979 void 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
2065 void 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
2181 void 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
2233 int 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