binaries/src/fasta34/compacc.c

   1
   2 /* copyright (c) 1996, 1997, 1998, 1999 William R. Pearson and the
   3    U. of Virginia */
   4
   5 /* $Name: fa_34_26_5 $ - $Id: compacc.c,v 1.61 2007/04/26 18:37:18 wrp Exp $ */
   6
   7 /* Concurrent read version */
   8
   9 #include <stdio.h>
  10 #include <stdlib.h>
  11 #if defined(UNIX) || defined(WIN32)
  12 #include <sys/types.h>
  13 #endif
  14
  15 #include <limits.h>
  16 #include <float.h>
  17
  18 #include <string.h>
  19 #include <time.h>
  20 #include <math.h>
  21
  22 #include "defs.h"
  23 #include "param.h"
  24 #include "structs.h"
  25
  26 #ifndef PCOMPLIB
  27 #include "mw.h"
  28 #else
  29 #include "p_mw.h"
  30 #endif
  31
  32 #define XTERNAL
  33 #include "uascii.h"
  34 #include "upam.h"
  35 #undef XTERNAL
  36
  37 #ifdef PCOMPLIB
  38 #include "msg.h"
  39 extern int nnodes;
  40 #ifdef PVM_SRC
  41 #include "pvm3.h"
  42 extern int pinums[],hosttid;
  43 #endif
  44 #ifdef MPI_SRC
  45 #include "mpi.h"
  46 #endif
  47 #endif
  48
  49 extern time_t tdone, tstart;            /* Timing */
  50 extern void abort ();
  51 extern void ptime ();
  52
  53 /* because it is used to pre-allocate space, maxn has various
  54    constraints.  For "simple" comparisons, it is simply the length of
  55    the longest library sequence.  But for translated comparisons, it
  56    must be 3 or 6X the length of the query sequence.
  57
  58    In addition, however, it can be reduced to make certain that
  59    sequences are read in smaller chunks.  And, maxn affect how large
  60    overlaps must be when sequences are read in chunks.
  61 */
  62
  63 int
  64 reset_maxn(struct mngmsg *m_msg, int maxn) {
  65
  66   /* reduce maxn if requested */
  67   if (m_msg->maxn > 0 && m_msg->maxn < maxn) maxn = m_msg->maxn;
  68
  69   if (m_msg->qdnaseq==m_msg->ldnaseq || m_msg->qdnaseq==SEQT_DNA ||
  70       m_msg->qdnaseq == SEQT_RNA) {/* !TFAST - either FASTA or FASTX*/
  71
  72    if (m_msg->n0> m_msg->max_tot/3) {
  73       fprintf(stderr," query sequence is too long %d > %d %s\n",
  74               m_msg->n0,
  75               m_msg->max_tot/3,
  76               m_msg->sqnam);
  77       exit(1);
  78     }
  79     m_msg->loff = m_msg->n0;
  80     m_msg->maxt3 = maxn-m_msg->loff;
  81   }
  82   else {        /* is TFAST */
  83     if (m_msg->n0 > MAXTST) {
  84       fprintf(stderr," query sequence is too long %d %s\n",m_msg->n0,m_msg->sqnam);
  85       exit(1);
  86     }
  87
  88     if (m_msg->n0*3 > maxn ) {  /* n0*3 for the three frames - this
  89                                    will only happen if maxn has been
  90                                    set low manually */
  91
  92       if (m_msg->n0*4+2 < m_msg->max_tot) { /* m_msg0*3 + m_msg0 */
  93         fprintf(stderr,
  94                 " query sequence too long for library segment: %d - resetting to %d\n",
  95               maxn,m_msg->n0*3);
  96         maxn = m_msg->maxn = m_msg->n0*3;
  97       }
  98       else {
  99         fprintf(stderr," query sequence too long for translated search: %d * 4 > %d %s\n",
 100               m_msg->n0,maxn, m_msg->sqnam);
 101         exit(1);
 102       }
 103     }
 104
 105     /* set up some constants for overlaps */
 106     m_msg->loff = 3*m_msg->n0;
 107     m_msg->maxt3 = maxn-m_msg->loff-3;
 108     m_msg->maxt3 -= m_msg->maxt3%3;
 109     m_msg->maxt3++;
 110
 111     maxn = maxn - 3; maxn -= maxn%3; maxn++;
 112   }
 113   return maxn;
 114 }
 115
 116
 117 int
 118 scanseq(unsigned char *seq, int n, char *str) {
 119   int tot,i;
 120   char aaray[128];              /* this must be set > nsq */
 121
 122   for (i=0; i<128; i++)  aaray[i]=0;
 123   for (i=0; (size_t)i < strlen(str); i++) aaray[qascii[str[i]]]=1;
 124   for (i=tot=0; i<n; i++) tot += aaray[seq[i]];
 125   return tot;
 126 }
 127
 128 /* subs_env takes a string, possibly with ${ENV}, and looks up all the
 129    potential environment variables and substitutes them into the
 130    string */
 131
 132 void subs_env(char *dest, char *src, int dest_size) {
 133   char *last_src, *bp, *bp1;
 134
 135   last_src = src;
 136
 137   if ((bp = strchr(src,'$'))==NULL) {
 138     strncpy(dest, src, dest_size);
 139     dest[dest_size-1] = '\0';
 140   }
 141   else {
 142     *dest = '\0';
 143     while (strlen(dest) < dest_size-1 && bp != NULL ) {
 144       /* copy stuff before ${*/
 145       *bp = '\0';
 146       strncpy(dest, last_src, dest_size);
 147       *bp = '$';
 148
 149       /* copy ENV */
 150       if (*(bp+1) != '{') {
 151         strncat(dest, "$", dest_size - strlen(dest) -1);
 152         dest[dest_size-1] = '\0';
 153         bp += 1;
 154       }
 155       else {    /* have  ${ENV} - put it in */
 156         if ((bp1 = strchr(bp+2,'}'))==NULL) {
 157           fprintf(stderr, "Unterminated ENV: %s\n",src);
 158           break;
 159         }
 160         else {
 161           *bp1 = '\0';
 162           if (getenv(bp+2)!=NULL) {
 163             strncat(dest, getenv(bp+2), dest_size - strlen(dest) - 1);
 164             dest[dest_size-1] = '\0';
 165             *bp1 = '}';
 166           }
 167           bp = bp1+1;   /* bump bp even if getenv == NULL */
 168         }
 169       }
 170       last_src = bp;
 171
 172       /* now get the next ${ENV} if present */
 173       bp = strchr(last_src,'$');
 174     }
 175     /* now copy the last stuff */
 176     strncat(dest, last_src, dest_size - strlen(dest) - 1);
 177     dest[dest_size-1]='\0';
 178   }
 179 }
 180
 181
 182 void selectbest(bptr,k,n)       /* k is rank in array */
 183      struct beststr **bptr;
 184      int k,n;
 185 {
 186   int v, i, j, l, r;
 187   struct beststr *tmptr;
 188
 189   l=0; r=n-1;
 190
 191   while ( r > l ) {
 192     v = bptr[r]->score[0];
 193     i = l-1;
 194     j = r;
 195     do {
 196       while (bptr[++i]->score[0] > v) ;
 197       while (bptr[--j]->score[0] < v) ;
 198       tmptr = bptr[i]; bptr[i]=bptr[j]; bptr[j]=tmptr;
 199     } while (j > i);
 200     bptr[j]=bptr[i]; bptr[i]=bptr[r]; bptr[r]=tmptr;
 201     if (i>=k) r = i-1;
 202     if (i<=k) l = i+1;
 203   }
 204 }
 205
 206 void selectbestz(bptr,k,n)      /* k is rank in array */
 207      struct beststr **bptr;
 208      int k,n;
 209 {
 210   int i, j, l, r;
 211   struct beststr *tmptr;
 212   double v;
 213
 214   l=0; r=n-1;
 215
 216   while ( r > l ) {
 217     v = bptr[r]->zscore;
 218     i = l-1;
 219     j = r;
 220     do {
 221       while (bptr[++i]->zscore > v) ;
 222       while (bptr[--j]->zscore < v) ;
 223       tmptr = bptr[i]; bptr[i]=bptr[j]; bptr[j]=tmptr;
 224     } while (j > i);
 225     bptr[j]=bptr[i]; bptr[i]=bptr[r]; bptr[r]=tmptr;
 226     if (i>=k) r = i-1;
 227     if (i<=k) l = i+1;
 228   }
 229 }
 230
 231 /* improved shellsort with high-performance increments */
 232 /*
 233 shellsort(itemType a[], int l, int r)
 234 { int i, j, k, h; itemType v;
 235  int incs[16] = { 1391376, 463792, 198768, 86961, 33936,
 236                   13776, 4592, 1968, 861, 336,
 237                   112, 48, 21, 7, 3, 1 };
 238  for ( k = 0; k < 16; k++)
 239    for (h = incs[k], i = l+h; i <= r; i++)
 240      {
 241        v = a[i]; j = i;
 242        while (j > h && a[j-h] > v)
 243          { a[j] = a[j-h]; j -= h; }
 244        a[j] = v;
 245      }
 246 }
 247 */
 248
 249 /* ?improved? version of sortbestz using optimal increments and fewer
 250    exchanges */
 251 void sortbestz(struct beststr **bptr, int nbest)
 252 {
 253   int gap, i, j, k;
 254   struct beststr *tmp;
 255   double v;
 256   int incs[16] = { 1391376, 463792, 198768, 86961, 33936,
 257                    13776, 4592, 1968, 861, 336,
 258                    112, 48, 21, 7, 3, 1 };
 259
 260   for ( k = 0; k < 16; k++) {
 261     gap = incs[k];
 262     for (i=gap; i < nbest; i++) {
 263       tmp = bptr[i];
 264       j = i;
 265       v = bptr[i]->zscore;
 266       while ( j >= gap && bptr[j-gap]->zscore < v) {
 267         bptr[j] = bptr[j - gap];
 268         j -= gap;
 269       }
 270       bptr[j] = tmp;
 271     }
 272   }
 273 }
 274
 275
 276 void sortbeste(struct beststr **bptr, int nbest)
 277 {
 278   int gap, i, j, k;
 279   struct beststr *tmp;
 280   double v;
 281   int incs[16] = { 1391376, 463792, 198768, 86961, 33936,
 282                    13776, 4592, 1968, 861, 336,
 283                    112, 48, 21, 7, 3, 1 };
 284
 285   for ( k = 0; k < 16; k++) {
 286     gap = incs[k];
 287     for (i=gap; i < nbest; i++) {
 288       j = i;
 289       tmp = bptr[i];
 290       v = tmp->escore;
 291       while ( j >= gap && bptr[j-gap]->escore > v) {
 292         bptr[j] = bptr[j - gap];
 293         j -= gap;
 294       }
 295       bptr[j] = tmp;
 296     }
 297   }
 298
 299   /* sometimes there are many high scores with E()==0.0, sort
 300      those by z() score */
 301
 302   j = 0;
 303   while (j < nbest && bptr[j]->escore <= 2.0*DBL_MIN ) {j++;}
 304   if (j > 1) sortbestz(bptr,j);
 305 }
 306
 307 extern double zs_to_Ec(double zs, long entries);
 308
 309 /*
 310 extern double ks_dev;
 311 extern int ks_df; */
 312 extern char hstring1[];
 313
 314 void
 315 prhist(FILE *fd, struct mngmsg m_msg,
 316        struct pstruct pst,
 317        struct hist_str hist,
 318        int nstats,
 319        struct db_str ntt,
 320        char *gstring2)
 321 {
 322   int i,j,hl,hll, el, ell, ev;
 323   char hline[80], pch, *bp;
 324   int mh1, mht;
 325   int maxval, maxvalt, dotsiz, ddotsiz,doinset;
 326   double cur_e, prev_e, f_int;
 327   double max_dev, x_tmp;
 328   double db_tt;
 329   int n_chi_sq, cum_hl=0, max_i;
 330
 331
 332   fprintf(fd,"\n");
 333
 334   if (pst.zsflag_f < 0) {
 335     fprintf(fd, "%7ld residues in %5ld sequences\n", ntt.length,ntt.entries);
 336     fprintf(fd,"\n%s\n",gstring2);
 337     return;
 338   }
 339
 340   if (nstats > 20) {
 341     max_dev = 0.0;
 342     mh1 = hist.maxh-1;
 343     mht = (3*hist.maxh-3)/4 - 1;
 344
 345     if (!m_msg.nohist && mh1 > 0) {
 346       for (i=0,maxval=0,maxvalt=0; i<hist.maxh; i++) {
 347         if (hist.hist_a[i] > maxval) maxval = hist.hist_a[i];
 348         if (i >= mht &&  hist.hist_a[i]>maxvalt) maxvalt = hist.hist_a[i];
 349       }
 350       n_chi_sq = 0;
 351       cum_hl = -hist.hist_a[0];
 352       dotsiz = (maxval-1)/60+1;
 353       ddotsiz = (maxvalt-1)/50+1;
 354       doinset = (ddotsiz < dotsiz && dotsiz > 2);
 355
 356       if (pst.zsflag_f>=0)
 357         fprintf(fd,"       opt      E()\n");
 358       else
 359         fprintf(fd,"     opt\n");
 360
 361       prev_e =  zs_to_Ec((double)(hist.min_hist-hist.histint/2),hist.entries);
 362       for (i=0; i<=mh1; i++) {
 363         pch = (i==mh1) ? '>' : ' ';
 364         pch = (i==0) ? '<' : pch;
 365         hll = hl = hist.hist_a[i];
 366         if (pst.zsflag_f>=0) {
 367           cum_hl += hl;
 368           f_int = (double)(i*hist.histint+hist.min_hist)+(double)hist.histint/2.0;
 369           cur_e = zs_to_Ec(f_int,hist.entries);
 370           ev = el = ell = (int)(cur_e - prev_e + 0.5);
 371           if (hl > 0  && i > 5 && i < (90-hist.min_hist)/hist.histint) {
 372             x_tmp  = fabs(cum_hl - cur_e);
 373             if ( x_tmp > max_dev) {
 374               max_dev = x_tmp;
 375               max_i = i;
 376             }
 377             n_chi_sq++;
 378           }
 379           if ((el=(el+dotsiz-1)/dotsiz) > 60) el = 60;
 380           if ((ell=(ell+ddotsiz-1)/ddotsiz) > 40) ell = 40;
 381           fprintf(fd,"%c%3d %5d %5d:",
 382                   pch,(i<mh1)?(i)*hist.histint+hist.min_hist :
 383                   mh1*hist.histint+hist.min_hist,hl,ev);
 384         }
 385         else fprintf(fd,"%c%3d %5d :",
 386                      pch,(i<mh1)?(i)*hist.histint+hist.min_hist :
 387                      mh1*hist.histint+hist.min_hist,hl);
 388
 389         if ((hl=(hl+dotsiz-1)/dotsiz) > 60) hl = 60;
 390         if ((hll=(hll+ddotsiz-1)/ddotsiz) > 40) hll = 40;
 391         for (j=0; j<hl; j++) hline[j]='=';
 392         if (pst.zsflag_f>=0) {
 393           if (el <= hl ) {
 394             if (el > 0) hline[el-1]='*';
 395             hline[hl]='\0';
 396           }
 397           else {
 398             for (j = hl; j < el; j++) hline[j]=' ';
 399             hline[el-1]='*';
 400             hline[hl=el]='\0';
 401           }
 402         }
 403         else hline[hl] = 0;
 404         if (i==1) {
 405           for (j=hl; j<10; j++) hline[j]=' ';
 406           sprintf(&hline[10]," one = represents %d library sequences",dotsiz);
 407         }
 408         if (doinset && i == mht-2) {
 409           for (j = hl; j < 10; j++) hline[j]=' ';
 410           sprintf(&hline[10]," inset = represents %d library sequences",ddotsiz);
 411         }
 412         if (i >= mht&& doinset ) {
 413           for (j = hl; j < 10; j++) hline[j]=' ';
 414           hline[10]=':';
 415           for (j = 11; j<11+hll; j++) hline[j]='=';
 416           hline[11+hll]='\0';
 417           if (pst.zsflag_f>=0) {
 418             if (ell <= hll) hline[10+ell]='*';
 419             else {
 420               for (j = 11+hll; j < 10+ell; j++) hline[j]=' ';
 421               hline[10+ell] = '*';
 422               hline[11+ell] = '\0';
 423             }
 424           }
 425         }
 426
 427         fprintf(fd,"%s\n",hline);
 428         prev_e = cur_e;
 429       }
 430     }
 431   }
 432
 433   if (ntt.carry==0) {
 434     fprintf(fd, "%7ld residues in %5ld sequences\n", ntt.length, ntt.entries);
 435   }
 436   else {
 437     db_tt = (double)ntt.carry*(double)LONG_MAX + (double)ntt.length;
 438     fprintf(fd, "%.0f residues in %5ld library sequences\n", db_tt, ntt.entries);
 439   }
 440
 441   if (pst.zsflag_f>=0) {
 442     if (MAXSTATS < hist.entries)
 443 #ifdef SAMP_STATS
 444       fprintf(fd," statistics sampled from %d to %ld sequences\n",
 445               MAXSTATS,hist.entries);
 446 #else
 447       fprintf(fd," statistics extrapolated from %d to %ld sequences\n",
 448               MAXSTATS,hist.entries);
 449 #endif
 450     /*    summ_stats(stat_info); */
 451     fprintf(fd," %s\n",hist.stat_info);
 452     if (!m_msg.nohist && cum_hl > 0)
 453       fprintf(fd," Kolmogorov-Smirnov  statistic: %6.4f (N=%d) at %3d\n",
 454               max_dev/(float)cum_hl, n_chi_sq,max_i*hist.histint+hist.min_hist);
 455     if (m_msg.markx & MX_M10FORM) {
 456       while ((bp=strchr(hist.stat_info,'\n'))!=NULL) *bp=' ';
 457       if (cum_hl <= 0) cum_hl = -1;
 458       sprintf(hstring1,"; mp_extrap: %d %ld\n; mp_stats: %s\n; mp_KS: %6.4f (N=%d) at %3d\n",
 459               MAXSTATS,hist.entries,hist.stat_info,max_dev/(float)cum_hl, n_chi_sq,max_i*hist.histint+hist.min_hist);
 460     }
 461   }
 462   fprintf(fd,"\n%s\n",gstring2);
 463   fflush(fd);
 464 }
 465
 466 extern char prog_name[], *verstr;
 467
 468 void s_abort (char *p,  char *p1)
 469 {
 470   int i;
 471
 472   fprintf (stderr, "\n***[%s] %s%s***\n", prog_name, p, p1);
 473 #ifdef PCOMPLIB
 474 #ifdef PVM_SRC
 475   for (i=FIRSTNODE; i< nnodes; i++) pvm_kill(pinums[i]);
 476   pvm_exit();
 477 #endif
 478 #ifdef MPI_SRC
 479   MPI_Abort(MPI_COMM_WORLD,1);
 480   MPI_Finalize();
 481 #endif
 482 #endif
 483   exit (1);
 484 }
 485
 486 #ifndef MPI_SRC
 487 void w_abort (char *p, char *p1)
 488 {
 489   fprintf (stderr, "\n***[%s] %s%s***\n\n", prog_name, p, p1);
 490   exit (1);
 491 }
 492 #endif
 493
 494 #ifndef PCOMPLIB
 495 /* copies from from to to shuffling */
 496
 497 extern int nrand(int);
 498
 499 void
 500 shuffle(unsigned char *from, unsigned char *to, int n)
 501 {
 502   int i,j; unsigned char tmp;
 503
 504   if (from != to) memcpy((void *)to,(void *)from,n);
 505
 506   for (i=n; i>0; i--) {
 507     j = nrand(i);
 508     tmp = to[j];
 509     to[j] = to[i-1];
 510     to[i-1] = tmp;
 511   }
 512   to[n] = 0;
 513 }
 514
 515 /* copies from from to from shuffling, ieven changed for threads */
 516 void
 517 wshuffle(unsigned char *from, unsigned char *to, int n, int wsiz, int *ieven)
 518 {
 519   int i,j, k, mm;
 520   unsigned char tmp, *top;
 521
 522   memcpy((void *)to,(void *)from,n);
 523
 524   mm = n%wsiz;
 525
 526   if (*ieven) {
 527     for (k=0; k<(n-wsiz); k += wsiz) {
 528       top = &to[k];
 529       for (i=wsiz; i>0; i--) {
 530         j = nrand(i);
 531         tmp = top[j];
 532         top[j] = top[i-1];
 533         top[i-1] = tmp;
 534       }
 535     }
 536     top = &to[n-mm];
 537     for (i=mm; i>0; i--) {
 538       j = nrand(i);
 539       tmp = top[j];
 540       top[j] = top[i-1];
 541       top[i-1] = tmp;
 542     }
 543     *ieven = 0;
 544   }
 545   else {
 546     for (k=n; k>=wsiz; k -= wsiz) {
 547       top = &to[k-wsiz];
 548       for (i=wsiz; i>0; i--) {
 549         j = nrand(i);
 550         tmp = top[j];
 551         top[j] = top[i-1];
 552         top[i-1] = tmp;
 553       }
 554     }
 555     top = &to[0];
 556     for (i=mm; i>0; i--) {
 557       j = nrand(i);
 558       tmp = top[j];
 559       top[j] = top[i-1];
 560       top[i-1] = tmp;
 561     }
 562     *ieven = 1;
 563   }
 564   to[n] = 0;
 565 }
 566
 567 #endif
 568
 569 int
 570 sfn_cmp(int *q, int *s)
 571 {
 572   if (*q == *s) return *q;
 573   while (*q && *s) {
 574     if (*q == *s) return *q;
 575     else if (*q < *s) q++;
 576     else if (*q > *s) s++;
 577   }
 578   return 0;
 579 }
 580
 581 #ifndef MPI_SRC
 582
 583 #define ESS 49
 584
 585 void
 586 revcomp(unsigned char *seq, int n, int *c_nt)
 587 {
 588   unsigned char tmp;
 589   int i, ni;
 590
 591   for (i=0, ni = n-1; i< n/2; i++,ni--) {
 592     tmp = c_nt[seq[i]];
 593     seq[i] = c_nt[seq[ni]];
 594     seq[ni] = tmp;
 595   }
 596   if ((n%2)==1) {
 597     i = n/2;
 598     seq[i] = c_nt[seq[i]];
 599   }
 600   seq[n]=0;
 601 }
 602 #endif
 603
 604 #ifdef PCOMPLIB
 605
 606 /* init_stage2 sets up the data structures necessary to send a subset
 607    of sequences to the nodes, and then collects the results
 608 */
 609
 610 /* wstage2[] FIRSTNODE .. nnodes has the next sequence to be do_opt()/do_walign()ed */
 611 /* wstage2p[] is a list of sequence numbers/frames, to be sent to workers */
 612 /* wstage2b[] is a list of bptr's that shares the index with wstage2p[] */
 613
 614 static int  wstage2[MAXWRKR +1];        /* count of second stage scores */
 615 static struct stage2_str  *wstage2p[MAXWRKR+1]; /* list of second stage sequences */
 616 static int  wstage2i[MAXWRKR+1];        /* index into second stage sequences */
 617 static struct beststr *bbptr,
 618      **wstage2b[MAXWRKR+1];     /* reverse pointers to bestr */
 619
 620 void
 621 do_stage2(struct beststr **bptr, int nbest, struct mngmsg m_msg0,
 622           int s_func, struct qmng_str *qm_msp) {
 623
 624   int i, is, ib, iw, nres;
 625   int node, snode, node_done;
 626   int bufid, numt, tid;
 627   char errstr[120];
 628   struct comstr2 bestr2[BFR2+1];        /* temporary structure array */
 629   char *seqc_buff, *seqc;
 630   int seqc_buff_len, aln_code_n;
 631 #ifdef MPI_SRC
 632   MPI_Status mpi_status;
 633 #endif
 634
 635   /* initialize the counter for each worker to 0 */
 636   for (iw = FIRSTNODE; iw < nnodes; iw++) wstage2[iw] = 0;
 637
 638   /* for each result, bump the counter for the worker that has
 639      the sequence */
 640   for (ib = 0; ib < nbest; ib++ ) { wstage2[bptr[ib]->wrkr]++;  }
 641
 642   /* now allocate enough space to send each worker a
 643      list of its sequences stage2_str {seqnm, frame} */
 644   for (iw = FIRSTNODE; iw < nnodes; iw++) {
 645     if (wstage2[iw]>0) {
 646       if ((wstage2p[iw]=
 647            (struct stage2_str *)
 648            calloc(wstage2[iw],sizeof(struct stage2_str)))==NULL) {
 649         sprintf(errstr," cannot allocate sequence listp %d %d",
 650                 iw,wstage2[iw]);
 651         s_abort(errstr,"");
 652       }
 653
 654       /* allocate space to remember the bptr's for each result */
 655       if ((wstage2b[iw]=(struct beststr **)
 656            calloc(wstage2[iw],sizeof(struct beststr *)))==NULL) {
 657         sprintf(errstr," cannot allocate sequence listb %d %d",
 658                 iw,wstage2[iw]);
 659         s_abort(errstr,"");
 660       }
 661       wstage2i[iw]=0;
 662     }
 663     else {
 664       wstage2p[iw] = NULL;
 665       wstage2b[iw] = NULL;
 666     }
 667   }
 668
 669   /* for each result, set wstage2p[worker][result_index_in_worker] */
 670   for (is = 0; is < nbest; is++) {
 671     iw=bptr[is]->wrkr;
 672     wstage2p[iw][wstage2i[iw]].seqnm = bptr[is]->seqnm;
 673     wstage2p[iw][wstage2i[iw]].frame = bptr[is]->frame;
 674     wstage2b[iw][wstage2i[iw]] = bptr[is];
 675     wstage2i[iw]++;
 676   }
 677
 678
 679   /* at this point, wstage2i[iw] should equal wstage2[iw] */
 680   node_done = 0;
 681   for (node = FIRSTNODE; node < nnodes; node++) {
 682
 683     /*    fprintf(stderr,"node: %d stage2: %d\n",node,wstage2[node]); */
 684
 685     /* if a worker has no results, move on */
 686     if (wstage2[node]<=0) { node_done++; continue;}
 687
 688     qm_msp->slist = wstage2[node];      /* set number of results to return */
 689     qm_msp->s_func = s_func;            /* set s_funct for do_opt/do_walign */
 690 #ifdef PVM_SRC
 691     pvm_initsend(PvmDataRaw);
 692     pvm_pkbyte((char *)qm_msp,sizeof(struct qmng_str),1);
 693     pvm_send(pinums[node],MSEQTYPE);    /* send qm_msp */
 694     pvm_initsend(PvmDataRaw);   /* send the list of seqnm/frame */
 695     pvm_pkbyte((char *)wstage2p[node],wstage2[node]*sizeof(struct stage2_str),1);
 696     pvm_send(pinums[node],LISTTYPE);
 697 #endif
 698 #ifdef MPI_SRC
 699     MPI_Send(qm_msp,sizeof(struct qmng_str),MPI_BYTE,node,MSEQTYPE,
 700              MPI_COMM_WORLD);
 701     MPI_Send((char *)wstage2p[node],wstage2[node]*
 702              sizeof(struct stage2_str),MPI_BYTE,node,LISTTYPE,
 703              MPI_COMM_WORLD);
 704 #endif
 705   }
 706
 707   /* all the workers have their list of sequences */
 708   /* reset the index of results to obtain */
 709   for (iw = 0; iw < nnodes; iw++) wstage2i[iw]=0;
 710
 711   while (node_done < nnodes-FIRSTNODE) {
 712 #ifdef PVM_SRC
 713     bufid = pvm_recv(-1,LISTRTYPE);     /* wait for results */
 714     pvm_bufinfo(bufid,NULL,NULL,&tid);
 715     /* get a chunk of comstr2 results */
 716     pvm_upkbyte((char *)&bestr2[0],sizeof(struct comstr2)*(BFR2+1),1);
 717     snode = (iw=tidtonode(tid));
 718     pvm_freebuf(bufid);
 719 #endif
 720 #ifdef MPI_SRC
 721     MPI_Recv((char *)&bestr2[0],sizeof(struct comstr2)*(BFR2+1),
 722              MPI_BYTE,MPI_ANY_SOURCE,LISTRTYPE,MPI_COMM_WORLD,
 723              &mpi_status);
 724     snode = mpi_status.MPI_SOURCE;
 725     iw = snode;
 726 #endif
 727
 728     seqc_buff = NULL;
 729     if (s_func == DO_OPT_FLG && m_msg0.show_code==SHOW_CODE_ALIGN) {
 730 #ifdef PVM_SRC
 731       bufid = pvm_recv(tid,CODERTYPE);
 732       pvm_upkint(&seqc_buff_len,1,1);   /* get the code string length */
 733 #endif
 734 #ifdef MPI_SRC
 735       MPI_Recv((char *)&seqc_buff_len,1,MPI_INT, snode,
 736                CODERTYPE,MPI_COMM_WORLD, &mpi_status);
 737 #endif
 738
 739       seqc=seqc_buff = NULL;
 740       if (seqc_buff_len > 0) {          /* allocate space for it */
 741         if ((seqc=seqc_buff=calloc(seqc_buff_len,sizeof(char)))==NULL) {
 742           fprintf(stderr,"Cannot allocate seqc_buff: %d\n",seqc_buff_len);
 743           seqc_buff_len=0;
 744         }
 745         else {
 746 #ifdef PVM_SRC
 747           pvm_upkbyte(seqc_buff,seqc_buff_len*sizeof(char),1);
 748 #endif
 749 #ifdef MPI_SRC
 750           MPI_Recv((char *)seqc_buff,seqc_buff_len*sizeof(char),
 751                    MPI_BYTE,snode,CODERTYPE,MPI_COMM_WORLD, &mpi_status);
 752 #endif
 753         }
 754       }
 755 #ifdef PVM_SRC
 756       pvm_freebuf(bufid);
 757 #endif
 758     }
 759
 760     /* get number of results in this message */
 761     nres = bestr2[BFR2].seqnm & ~FINISHED;
 762     /* check to see if finished */
 763     if (bestr2[BFR2].seqnm&FINISHED) {node_done++;}
 764
 765     seqc = seqc_buff;
 766
 767     /* count through results from a specific worker */
 768     for (i=0,is=wstage2i[iw]; i < nres; i++,is++) {
 769
 770       /* get the (saved) bptr for this result */
 771       bbptr=wstage2b[iw][is];
 772       /* consistency check seqnm's must agree */
 773       if (wstage2p[iw][is].seqnm ==  bbptr->seqnm) {
 774         if (s_func == DO_CALC_FLG && m_msg0.last_calc_flg) {
 775           bbptr->score[0] = bestr2[i].score[0];
 776           bbptr->score[1] = bestr2[i].score[1];
 777           bbptr->score[2] = bestr2[i].score[2];
 778           bbptr->escore = bestr2[i].escore;
 779           bbptr->segnum = bestr2[i].segnum;
 780           bbptr->seglen = bestr2[i].seglen;
 781         }
 782         else if (m_msg0.stages > 1) {
 783           bbptr->score[0] = bestr2[i].score[0];
 784           bbptr->score[1] = bestr2[i].score[1];
 785           bbptr->score[2] = bestr2[i].score[2];
 786         }
 787
 788         if (s_func == DO_OPT_FLG && m_msg0.markx & MX_M9SUMM) {
 789           /* get score, alignment information, percents */
 790           bbptr->sw_score = bestr2[i].sw_score;
 791           memcpy(bbptr->aln_d,&bestr2[i].aln_d,sizeof(struct a_struct));
 792           bbptr->percent = bestr2[i].percent;
 793           bbptr->gpercent = bestr2[i].gpercent;
 794
 795           if (m_msg0.show_code == 2) {  /* if show code */
 796             /* length of encoding */
 797             aln_code_n = bbptr->aln_code_n = bestr2[i].aln_code_n;
 798             if (aln_code_n > 0) {
 799               if ((bbptr->aln_code =
 800                    (char *)calloc(aln_code_n+1,sizeof(char)))==NULL) {
 801                 fprintf(stderr,"cannot allocate seq_code[%d:%d]: %d\n",
 802                         bbptr->wrkr,bbptr->seqnm,aln_code_n);
 803                 seqc += aln_code_n+1;
 804                 bbptr->aln_code_n = 0;
 805               }
 806               else {
 807                 strncpy(bbptr->aln_code,seqc,aln_code_n);
 808                 bbptr->aln_code[aln_code_n]='\0';
 809                 seqc += aln_code_n+1;
 810               }
 811             }
 812             else {
 813               fprintf(stderr," aln_code_n <=0: %d\n",aln_code_n);
 814             }
 815           }
 816         }
 817       }
 818       else fprintf(stderr,"phase error in phase II return %d %d", iw,i);
 819     }
 820     if (seqc_buff != NULL) {
 821       free(seqc_buff);
 822       seqc_buff = NULL;
 823     }
 824     wstage2i[iw] += nres;
 825   }
 826
 827   for (iw=FIRSTNODE; iw < nnodes; iw++) {
 828     if ((void *)wstage2p[iw]!=NULL) free((void *)wstage2p[iw]);
 829     if ((void *)wstage2b[iw]!=NULL) free((void *)wstage2b[iw]);
 830   }
 831 }
 832
 833 #endif