source: branches/nameserver/GDE/RAxML/makenewzGenericSpecial.c

Last change on this file was 13845, checked in by westram, 10 years ago
  • remove executable flag
File size: 111.0 KB
Line 
1/*  RAxML-VI-HPC (version 2.2) a program for sequential and parallel estimation of phylogenetic trees
2 *  Copyright August 2006 by Alexandros Stamatakis
3 *
4 *  Partially derived from
5 *  fastDNAml, a program for estimation of phylogenetic trees from sequences by Gary J. Olsen
6 *
7 *  and
8 *
9 *  Programs of the PHYLIP package by Joe Felsenstein.
10 *
11 *  This program is free software; you may redistribute it and/or modify its
12 *  under the terms of the GNU General Public License as published by the Free
13 *  Software Foundation; either version 2 of the License, or (at your option)
14 *  any later version.
15 *
16 *  This program is distributed in the hope that it will be useful, but
17 *  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 *  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
19 *  for more details.
20 *
21 *
22 *  For any other enquiries send an Email to Alexandros Stamatakis
23 *  Alexandros.Stamatakis@epfl.ch
24 *
25 *  When publishing work that is based on the results from RAxML-VI-HPC please cite:
26 *
27 *  Alexandros Stamatakis:"RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with
28 *  thousands of taxa and mixed models".
29 *  Bioinformatics 2006; doi: 10.1093/bioinformatics/btl446
30 */
31
32#ifndef WIN32
33#include <unistd.h>
34#endif
35
36
37
38#include <math.h>
39#include <time.h>
40#include <stdlib.h>
41#include <stdio.h>
42#include <ctype.h>
43#include <string.h>
44#include "axml.h"
45
46#ifdef __SIM_SSE3
47#include <xmmintrin.h>
48#include <pmmintrin.h>
49/*#include <tmmintrin.h>*/
50#endif
51
52#ifdef _USE_PTHREADS
53extern volatile double *reductionBuffer;
54extern volatile double *reductionBufferTwo;
55extern volatile int NumberOfThreads;
56#endif
57
58extern const unsigned int mask32[32];
59
60
61/*******************/
62
63static void sumCAT_BINARY(int tipCase, double *sum, double *x1_start, double *x2_start, double *tipVector,
64                          unsigned char *tipX1, unsigned char *tipX2, int n)
65{
66  int i;
67 
68#ifndef __SIM_SSE3
69  int j;
70#endif
71  double *x1, *x2;
72
73  switch(tipCase)
74    {
75    case TIP_TIP:
76      for (i = 0; i < n; i++)
77        {
78          x1 = &(tipVector[2 * tipX1[i]]);
79          x2 = &(tipVector[2 * tipX2[i]]);
80
81#ifndef __SIM_SSE3
82          for(j = 0; j < 2; j++)
83            sum[i * 2 + j]     = x1[j] * x2[j];
84#else
85          _mm_store_pd(&sum[i * 2], _mm_mul_pd( _mm_load_pd(x1), _mm_load_pd(x2)));
86#endif
87        }
88      break;
89    case TIP_INNER:
90      for (i = 0; i < n; i++)
91        {
92          x1 = &(tipVector[2 * tipX1[i]]);
93          x2 = &x2_start[2 * i];
94
95#ifndef __SIM_SSE3
96          for(j = 0; j < 2; j++)
97            sum[i * 2 + j]     = x1[j] * x2[j];
98#else
99          _mm_store_pd(&sum[i * 2], _mm_mul_pd( _mm_load_pd(x1), _mm_load_pd(x2))); 
100#endif
101        }
102      break;
103    case INNER_INNER:
104      for (i = 0; i < n; i++)
105        {
106          x1 = &x1_start[2 * i];
107          x2 = &x2_start[2 * i];
108#ifndef __SIM_SSE3
109          for(j = 0; j < 2; j++)
110            sum[i * 2 + j]     = x1[j] * x2[j];
111#else
112          _mm_store_pd(&sum[i * 2], _mm_mul_pd( _mm_load_pd(x1), _mm_load_pd(x2)));   
113#endif
114        }
115      break;
116    default:
117      assert(0);
118    }
119}
120#ifndef __SIM_SSE3
121static void sumCAT(int tipCase, double *sum, double *x1_start, double *x2_start, double *tipVector,
122                   unsigned char *tipX1, unsigned char *tipX2, int n)
123{
124  int i, j;
125  double *x1, *x2;
126
127  switch(tipCase)
128    {
129    case TIP_TIP:
130      for (i = 0; i < n; i++)
131        {
132          x1 = &(tipVector[4 * tipX1[i]]);
133          x2 = &(tipVector[4 * tipX2[i]]);
134
135          for(j = 0; j < 4; j++)
136            sum[i * 4 + j]     = x1[j] * x2[j];
137
138        }
139      break;
140    case TIP_INNER:
141      for (i = 0; i < n; i++)
142        {
143          x1 = &(tipVector[4 * tipX1[i]]);
144          x2 = &x2_start[4 * i];
145
146          for(j = 0; j < 4; j++)
147            sum[i * 4 + j]     = x1[j] * x2[j];
148
149        }
150      break;
151    case INNER_INNER:
152      for (i = 0; i < n; i++)
153        {
154          x1 = &x1_start[4 * i];
155          x2 = &x2_start[4 * i];
156
157
158          for(j = 0; j < 4; j++)
159            sum[i * 4 + j]     = x1[j] * x2[j];
160
161        }
162      break;
163    default:
164      assert(0);
165    }
166}
167
168#else
169
170static void sumCAT_SAVE(int tipCase, double *sum, double *x1_start, double *x2_start, double *tipVector,
171    unsigned char *tipX1, unsigned char *tipX2, int n, double *x1_gapColumn, double *x2_gapColumn, unsigned int *x1_gap, unsigned int *x2_gap)
172{
173  int i;
174  double 
175    *x1, 
176    *x2,   
177    *x1_ptr = x1_start,
178    *x2_ptr = x2_start;
179
180  switch(tipCase)
181  {
182    case TIP_TIP:
183      for (i = 0; i < n; i++)
184      {
185        x1 = &(tipVector[4 * tipX1[i]]);
186        x2 = &(tipVector[4 * tipX2[i]]);
187
188        _mm_store_pd( &sum[i*4 + 0], _mm_mul_pd( _mm_load_pd( &x1[0] ), _mm_load_pd( &x2[0] )));
189        _mm_store_pd( &sum[i*4 + 2], _mm_mul_pd( _mm_load_pd( &x1[2] ), _mm_load_pd( &x2[2] )));
190      }
191      break;
192    case TIP_INNER:
193      for (i = 0; i < n; i++)
194      {
195        x1 = &(tipVector[4 * tipX1[i]]);
196        if(isGap(x2_gap, i))
197          x2 = x2_gapColumn;
198        else
199        {
200          x2 = x2_ptr;
201          x2_ptr += 4;
202        }
203
204        _mm_store_pd( &sum[i*4 + 0], _mm_mul_pd( _mm_load_pd( &x1[0] ), _mm_load_pd( &x2[0] )));
205        _mm_store_pd( &sum[i*4 + 2], _mm_mul_pd( _mm_load_pd( &x1[2] ), _mm_load_pd( &x2[2] )));
206      }
207      break;
208    case INNER_INNER:
209      for (i = 0; i < n; i++)
210      {
211        if(isGap(x1_gap, i))
212          x1 = x1_gapColumn;
213        else
214        {
215          x1 = x1_ptr;
216          x1_ptr += 4;
217        }
218
219        if(isGap(x2_gap, i))
220          x2 = x2_gapColumn;
221        else
222        {
223          x2 = x2_ptr;
224          x2_ptr += 4;
225        }
226
227        _mm_store_pd( &sum[i*4 + 0], _mm_mul_pd( _mm_load_pd( &x1[0] ), _mm_load_pd( &x2[0] )));
228        _mm_store_pd( &sum[i*4 + 2], _mm_mul_pd( _mm_load_pd( &x1[2] ), _mm_load_pd( &x2[2] )));
229
230      }   
231      break;
232    default:
233      assert(0);
234  }
235}
236
237
238static void sumCAT(int tipCase, double *sum, double *x1_start, double *x2_start, double *tipVector,
239                   unsigned char *tipX1, unsigned char *tipX2, int n)
240{
241  int i;
242  double 
243    *x1, 
244    *x2;
245
246  switch(tipCase)
247    {
248    case TIP_TIP:
249      for (i = 0; i < n; i++)
250        {
251          x1 = &(tipVector[4 * tipX1[i]]);
252          x2 = &(tipVector[4 * tipX2[i]]);
253
254          _mm_store_pd( &sum[i*4 + 0], _mm_mul_pd( _mm_load_pd( &x1[0] ), _mm_load_pd( &x2[0] )));
255          _mm_store_pd( &sum[i*4 + 2], _mm_mul_pd( _mm_load_pd( &x1[2] ), _mm_load_pd( &x2[2] )));
256        }
257      break;
258    case TIP_INNER:
259      for (i = 0; i < n; i++)
260        {
261          x1 = &(tipVector[4 * tipX1[i]]);
262          x2 = &x2_start[4 * i];
263
264          _mm_store_pd( &sum[i*4 + 0], _mm_mul_pd( _mm_load_pd( &x1[0] ), _mm_load_pd( &x2[0] )));
265          _mm_store_pd( &sum[i*4 + 2], _mm_mul_pd( _mm_load_pd( &x1[2] ), _mm_load_pd( &x2[2] )));
266        }
267      break;
268    case INNER_INNER:
269      for (i = 0; i < n; i++)
270        {
271          x1 = &x1_start[4 * i];
272          x2 = &x2_start[4 * i];
273
274          _mm_store_pd( &sum[i*4 + 0], _mm_mul_pd( _mm_load_pd( &x1[0] ), _mm_load_pd( &x2[0] )));
275          _mm_store_pd( &sum[i*4 + 2], _mm_mul_pd( _mm_load_pd( &x1[2] ), _mm_load_pd( &x2[2] )));
276
277        }   
278      break;
279    default:
280      assert(0);
281    }
282}
283
284#endif
285
286
287
288
289static void coreGTRCAT_BINARY(int upper, int numberOfCategories, double *sum,
290                              volatile double *d1, volatile double *d2, 
291                              double *rptr, double *EIGN, int *cptr, double lz, int *wgt)
292{
293  int i;
294  double
295    *d, *d_start,
296    tmp_0, inv_Li, dlnLidlz, d2lnLidlz2,
297    dlnLdlz = 0.0,
298    d2lnLdlz2 = 0.0;
299  double e[2];
300  double dd1;
301
302  e[0] = EIGN[0];
303  e[1] = EIGN[0] * EIGN[0];
304
305
306  d = d_start = (double *)rax_malloc(numberOfCategories * sizeof(double));
307
308  dd1 = e[0] * lz;
309
310  for(i = 0; i < numberOfCategories; i++)
311    d[i] = EXP(dd1 * rptr[i]);
312
313  for (i = 0; i < upper; i++)
314    {
315      double
316        r = rptr[cptr[i]],
317        wr1 = r * wgt[i],
318        wr2 = r * r * wgt[i];
319     
320      d = &d_start[cptr[i]];
321
322      inv_Li = sum[2 * i];
323      inv_Li += (tmp_0 = d[0] * sum[2 * i + 1]);
324
325      inv_Li = 1.0/inv_Li;
326
327      dlnLidlz   = tmp_0 * e[0];
328      d2lnLidlz2 = tmp_0 * e[1];
329
330      dlnLidlz   *= inv_Li;
331      d2lnLidlz2 *= inv_Li;
332
333      dlnLdlz   += wr1 * dlnLidlz;
334      d2lnLdlz2 += wr2 * (d2lnLidlz2 - dlnLidlz * dlnLidlz);
335    }
336
337  *d1 = dlnLdlz;
338  *d2 = d2lnLdlz2;
339
340  rax_free(d_start);
341}
342
343#ifdef __SIM_SSE3
344
345static void coreGTRCAT(int upper, int numberOfCategories, double *sum,
346                           volatile double *d1, volatile double *d2,
347                       double *rptr, double *EIGN, int *cptr, double lz, int *wgt)
348{
349  int i;
350  double
351    *d, *d_start,
352    inv_Li, dlnLidlz, d2lnLidlz2,
353    dlnLdlz = 0.0,
354    d2lnLdlz2 = 0.0;
355  double e1[4] __attribute__ ((aligned (BYTE_ALIGNMENT)));
356  double e2[4] __attribute__ ((aligned (BYTE_ALIGNMENT)));
357  double dd1, dd2, dd3;
358
359  __m128d
360    e1v[2],
361    e2v[2];
362
363  e1[0] = 0.0;
364  e2[0] = 0.0;
365  e1[1] = EIGN[0];
366  e2[1] = EIGN[0] * EIGN[0];
367  e1[2] = EIGN[1];
368  e2[2] = EIGN[1] * EIGN[1];
369  e1[3] = EIGN[2];
370  e2[3] = EIGN[2] * EIGN[2];
371
372  e1v[0]= _mm_load_pd(&e1[0]);
373  e1v[1]= _mm_load_pd(&e1[2]);
374
375  e2v[0]= _mm_load_pd(&e2[0]);
376  e2v[1]= _mm_load_pd(&e2[2]);
377
378  d = d_start = (double *)rax_malloc(numberOfCategories * 4 * sizeof(double));
379
380  dd1 = EIGN[0] * lz;
381  dd2 = EIGN[1] * lz;
382  dd3 = EIGN[2] * lz;
383
384  for(i = 0; i < numberOfCategories; i++)
385    {
386      d[i * 4 + 0] = 1.0;
387      d[i * 4 + 1] = EXP(dd1 * rptr[i]);
388      d[i * 4 + 2] = EXP(dd2 * rptr[i]);
389      d[i * 4 + 3] = EXP(dd3 * rptr[i]);
390    }
391
392  for (i = 0; i < upper; i++)
393    {
394      double 
395        *s = &sum[4 * i],
396        r = rptr[cptr[i]],
397        wr1 = r * wgt[i],
398        wr2 = r * r * wgt[i];
399     
400      d = &d_start[4 * cptr[i]]; 
401     
402      __m128d tmp_0v =_mm_mul_pd(_mm_load_pd(&d[0]),_mm_load_pd(&s[0]));
403      __m128d tmp_1v =_mm_mul_pd(_mm_load_pd(&d[2]),_mm_load_pd(&s[2]));
404
405      __m128d inv_Liv    = _mm_add_pd(tmp_0v, tmp_1v);     
406                 
407      __m128d dlnLidlzv   = _mm_add_pd(_mm_mul_pd(tmp_0v, e1v[0]), _mm_mul_pd(tmp_1v, e1v[1]));   
408      __m128d d2lnLidlz2v = _mm_add_pd(_mm_mul_pd(tmp_0v, e2v[0]), _mm_mul_pd(tmp_1v, e2v[1]));
409
410
411      inv_Liv   = _mm_hadd_pd(inv_Liv, inv_Liv);
412      dlnLidlzv = _mm_hadd_pd(dlnLidlzv, dlnLidlzv);
413      d2lnLidlz2v = _mm_hadd_pd(d2lnLidlz2v, d2lnLidlz2v);                 
414 
415      _mm_storel_pd(&inv_Li, inv_Liv);     
416      _mm_storel_pd(&dlnLidlz, dlnLidlzv);                 
417      _mm_storel_pd(&d2lnLidlz2, d2lnLidlz2v);     
418
419      inv_Li = 1.0/inv_Li;
420
421      dlnLidlz   *= inv_Li;
422      d2lnLidlz2 *= inv_Li;
423
424      dlnLdlz   += wr1 * dlnLidlz;
425      d2lnLdlz2 += wr2 * (d2lnLidlz2 - dlnLidlz * dlnLidlz);
426    }
427
428  *d1 = dlnLdlz;
429  *d2 = d2lnLdlz2;
430
431  rax_free(d_start);
432}
433
434
435#else
436
437static void coreGTRCAT(int upper, int numberOfCategories, double *sum,
438                           volatile double *d1, volatile double *d2, 
439                       double *rptr, double *EIGN, int *cptr, double lz, int *wgt)
440{
441  int i;
442  double
443    *d, *d_start,
444    tmp_0, tmp_1, tmp_2, inv_Li, dlnLidlz, d2lnLidlz2,
445    dlnLdlz = 0.0,
446    d2lnLdlz2 = 0.0;
447  double e[6];
448  double dd1, dd2, dd3;
449
450  e[0] = EIGN[0];
451  e[1] = EIGN[0] * EIGN[0];
452  e[2] = EIGN[1];
453  e[3] = EIGN[1] * EIGN[1];
454  e[4] = EIGN[2];
455  e[5] = EIGN[2] * EIGN[2];
456
457  d = d_start = (double *)rax_malloc(numberOfCategories * 4 * sizeof(double));
458
459  dd1 = e[0] * lz;
460  dd2 = e[2] * lz;
461  dd3 = e[4] * lz;
462
463  for(i = 0; i < numberOfCategories; i++)
464    {
465      d[i * 4] = EXP(dd1 * rptr[i]);
466      d[i * 4 + 1] = EXP(dd2 * rptr[i]);
467      d[i * 4 + 2] = EXP(dd3 * rptr[i]);
468    }
469
470  for (i = 0; i < upper; i++)
471    {
472      double
473        r = rptr[cptr[i]],
474        wr1 = r * wgt[i],
475        wr2 = r * r * wgt[i];
476
477      d = &d_start[4 * cptr[i]];
478
479      inv_Li = sum[4 * i];
480      inv_Li += (tmp_0 = d[0] * sum[4 * i + 1]);
481      inv_Li += (tmp_1 = d[1] * sum[4 * i + 2]);
482      inv_Li += (tmp_2 = d[2] * sum[4 * i + 3]);
483
484      inv_Li = 1.0/inv_Li;
485
486      dlnLidlz   = tmp_0 * e[0];
487      d2lnLidlz2 = tmp_0 * e[1];
488
489      dlnLidlz   += tmp_1 * e[2];
490      d2lnLidlz2 += tmp_1 * e[3];
491
492      dlnLidlz   += tmp_2 * e[4];
493      d2lnLidlz2 += tmp_2 * e[5];
494
495      dlnLidlz   *= inv_Li;
496      d2lnLidlz2 *= inv_Li;
497
498
499      dlnLdlz   += wr1 * dlnLidlz;
500      d2lnLdlz2 += wr2 * (d2lnLidlz2 - dlnLidlz * dlnLidlz);
501    }
502
503  *d1 = dlnLdlz;
504  *d2 = d2lnLdlz2;
505
506  rax_free(d_start);
507}
508
509#endif
510
511
512
513
514
515#ifdef __SIM_SSE3
516static void sumGTRCATPROT_SAVE(int tipCase, double *sumtable, double *x1, double *x2, double *tipVector,
517    unsigned char *tipX1, unsigned char *tipX2, int n, 
518    double *x1_gapColumn, double *x2_gapColumn, unsigned int *x1_gap, unsigned int *x2_gap)
519{
520  int 
521    i, 
522  l;
523
524  double 
525    *sum, 
526    *left, 
527    *right,
528    *left_ptr = x1,
529    *right_ptr = x2;
530
531  switch(tipCase)
532  {
533    case TIP_TIP:
534      for (i = 0; i < n; i++)
535      {
536        left  = &(tipVector[20 * tipX1[i]]);
537        right = &(tipVector[20 * tipX2[i]]);
538        sum = &sumtable[20 * i];
539
540        for(l = 0; l < 20; l+=2)
541        {
542          __m128d sumv = _mm_mul_pd(_mm_load_pd(&left[l]), _mm_load_pd(&right[l]));
543
544          _mm_store_pd(&sum[l], sumv);           
545        }
546
547      }
548      break;
549    case TIP_INNER:
550      for (i = 0; i < n; i++)
551      {
552        left = &(tipVector[20 * tipX1[i]]);       
553
554        if(isGap(x2_gap, i))
555          right = x2_gapColumn;
556        else
557        {
558          right = right_ptr;
559          right_ptr += 20;
560        }
561
562        sum = &sumtable[20 * i];
563
564        for(l = 0; l < 20; l+=2)
565        {
566          __m128d sumv = _mm_mul_pd(_mm_load_pd(&left[l]), _mm_load_pd(&right[l]));
567
568          _mm_store_pd(&sum[l], sumv);           
569        }
570
571      }
572      break;
573    case INNER_INNER:
574      for (i = 0; i < n; i++)
575      { 
576        if(isGap(x1_gap, i))
577          left = x1_gapColumn;
578        else
579        {
580          left = left_ptr;
581          left_ptr += 20;
582        }
583
584        if(isGap(x2_gap, i))
585          right = x2_gapColumn;
586        else
587        {
588          right = right_ptr;
589          right_ptr += 20;
590        }
591
592        sum = &sumtable[20 * i];
593
594        for(l = 0; l < 20; l+=2)
595        {
596          __m128d sumv = _mm_mul_pd(_mm_load_pd(&left[l]), _mm_load_pd(&right[l]));
597
598          _mm_store_pd(&sum[l], sumv);           
599        }
600      }
601      break;
602    default:
603      assert(0);
604  }
605}
606
607
608#endif
609
610
611static void sumGTRCATPROT(int tipCase, double *sumtable, double *x1, double *x2, double *tipVector,
612                          unsigned char *tipX1, unsigned char *tipX2, int n)
613{
614  int i, l;
615  double *sum, *left, *right;
616
617  switch(tipCase)
618    {
619    case TIP_TIP:
620      for (i = 0; i < n; i++)
621        {
622          left  = &(tipVector[20 * tipX1[i]]);
623          right = &(tipVector[20 * tipX2[i]]);
624          sum = &sumtable[20 * i];
625#ifdef __SIM_SSE3
626          for(l = 0; l < 20; l+=2)
627            {
628              __m128d sumv = _mm_mul_pd(_mm_load_pd(&left[l]), _mm_load_pd(&right[l]));
629             
630              _mm_store_pd(&sum[l], sumv);               
631            }
632#else
633          for(l = 0; l < 20; l++)
634            sum[l] = left[l] * right[l];
635#endif
636        }
637      break;
638    case TIP_INNER:
639      for (i = 0; i < n; i++)
640        {
641          left = &(tipVector[20 * tipX1[i]]);
642          right = &x2[20 * i];
643          sum = &sumtable[20 * i];
644#ifdef __SIM_SSE3
645          for(l = 0; l < 20; l+=2)
646            {
647              __m128d sumv = _mm_mul_pd(_mm_load_pd(&left[l]), _mm_load_pd(&right[l]));
648             
649              _mm_store_pd(&sum[l], sumv);               
650            }
651#else
652          for(l = 0; l < 20; l++)
653            sum[l] = left[l] * right[l];
654#endif
655        }
656      break;
657    case INNER_INNER:
658      for (i = 0; i < n; i++)
659        {
660          left  = &x1[20 * i];
661          right = &x2[20 * i];
662          sum = &sumtable[20 * i];
663#ifdef __SIM_SSE3
664          for(l = 0; l < 20; l+=2)
665            {
666              __m128d sumv = _mm_mul_pd(_mm_load_pd(&left[l]), _mm_load_pd(&right[l]));
667             
668              _mm_store_pd(&sum[l], sumv);               
669            }
670#else
671          for(l = 0; l < 20; l++)
672            sum[l] = left[l] * right[l];
673#endif
674        }
675      break;
676    default:
677      assert(0);
678    }
679}
680
681
682
683static void sumGTRCATSECONDARY(int tipCase, double *sumtable, double *x1, double *x2, double *tipVector,
684                               unsigned char *tipX1, unsigned char *tipX2, int n)
685{
686  int i, l;
687  double *sum, *left, *right;
688
689  switch(tipCase)
690    {
691    case TIP_TIP:
692      for (i = 0; i < n; i++)
693        {
694          left  = &(tipVector[16 * tipX1[i]]);
695          right = &(tipVector[16 * tipX2[i]]);
696          sum = &sumtable[16 * i];
697
698          for(l = 0; l < 16; l++)
699            sum[l] = left[l] * right[l];
700        }
701      break;
702    case TIP_INNER:
703      for (i = 0; i < n; i++)
704        {
705          left = &(tipVector[16 * tipX1[i]]);
706          right = &x2[16 * i];
707          sum = &sumtable[16 * i];
708
709          for(l = 0; l < 16; l++)
710            sum[l] = left[l] * right[l];
711        }
712      break;
713    case INNER_INNER:
714      for (i = 0; i < n; i++)
715        {
716          left  = &x1[16 * i];
717          right = &x2[16 * i];
718          sum = &sumtable[16 * i];
719
720          for(l = 0; l < 16; l++)
721            sum[l] = left[l] * right[l];
722        }
723      break;
724    default:
725      assert(0);
726    }
727}
728
729
730static void sumCatFlex(int tipCase, double *sumtable, double *x1, double *x2, double *tipVector,
731                       unsigned char *tipX1, unsigned char *tipX2, int n, const int numStates)
732{
733  int i, l;
734  double *sum, *left, *right;
735
736  switch(tipCase)
737    {
738    case TIP_TIP:
739      for (i = 0; i < n; i++)
740        {
741          left  = &(tipVector[numStates * tipX1[i]]);
742          right = &(tipVector[numStates * tipX2[i]]);
743          sum = &sumtable[numStates * i];
744
745          for(l = 0; l < numStates; l++)
746            sum[l] = left[l] * right[l];
747        }
748      break;
749    case TIP_INNER:
750      for (i = 0; i < n; i++)
751        {
752          left = &(tipVector[numStates * tipX1[i]]);
753          right = &x2[numStates * i];
754          sum = &sumtable[numStates * i];
755
756          for(l = 0; l < numStates; l++)
757            sum[l] = left[l] * right[l];
758        }
759      break;
760    case INNER_INNER:
761      for (i = 0; i < n; i++)
762        {
763          left  = &x1[numStates * i];
764          right = &x2[numStates * i];
765          sum = &sumtable[numStates * i];
766
767          for(l = 0; l < numStates; l++)
768            sum[l] = left[l] * right[l];
769        }
770      break;
771    default:
772      assert(0);
773    }
774}
775
776
777static void sumGTRCATSECONDARY_6(int tipCase, double *sumtable, double *x1, double *x2, double *tipVector,
778                                 unsigned char *tipX1, unsigned char *tipX2, int n)
779{
780  int i, l;
781  double *sum, *left, *right;
782
783  switch(tipCase)
784    {
785    case TIP_TIP:
786      for (i = 0; i < n; i++)
787        {
788          left  = &(tipVector[6 * tipX1[i]]);
789          right = &(tipVector[6 * tipX2[i]]);
790          sum = &sumtable[6 * i];
791
792          for(l = 0; l < 6; l++)
793            sum[l] = left[l] * right[l];
794        }
795      break;
796    case TIP_INNER:
797      for (i = 0; i < n; i++)
798        {
799          left = &(tipVector[6 * tipX1[i]]);
800          right = &x2[6 * i];
801          sum = &sumtable[6 * i];
802
803          for(l = 0; l < 6; l++)
804            sum[l] = left[l] * right[l];
805        }
806      break;
807    case INNER_INNER:
808      for (i = 0; i < n; i++)
809        {
810          left  = &x1[6 * i];
811          right = &x2[6 * i];
812          sum = &sumtable[6 * i];
813
814          for(l = 0; l < 6; l++)
815            sum[l] = left[l] * right[l];
816        }
817      break;
818    default:
819      assert(0);
820    }
821}
822
823static void sumGTRCATSECONDARY_7(int tipCase, double *sumtable, double *x1, double *x2, double *tipVector,
824                                 unsigned char *tipX1, unsigned char *tipX2, int n)
825{
826  int i, l;
827  double *sum, *left, *right;
828
829  switch(tipCase)
830    {
831    case TIP_TIP:
832      for (i = 0; i < n; i++)
833        {
834          left  = &(tipVector[7 * tipX1[i]]);
835          right = &(tipVector[7 * tipX2[i]]);
836          sum = &sumtable[7 * i];
837
838          for(l = 0; l < 7; l++)
839            sum[l] = left[l] * right[l];
840        }
841      break;
842    case TIP_INNER:
843      for (i = 0; i < n; i++)
844        {
845          left = &(tipVector[7 * tipX1[i]]);
846          right = &x2[7 * i];
847          sum = &sumtable[7 * i];
848
849          for(l = 0; l < 7; l++)
850            sum[l] = left[l] * right[l];
851        }
852      break;
853    case INNER_INNER:
854      for (i = 0; i < n; i++)
855        {
856          left  = &x1[7 * i];
857          right = &x2[7 * i];
858          sum = &sumtable[7 * i];
859
860          for(l = 0; l < 7; l++)
861            sum[l] = left[l] * right[l];
862        }
863      break;
864    default:
865      assert(0);
866    }
867}
868
869
870#ifdef __SIM_SSE3
871
872static void coreGTRCATPROT(double *EIGN, double lz, int numberOfCategories, double *rptr, int *cptr, int upper,
873                           volatile double *ext_dlnLdlz,  volatile double *ext_d2lnLdlz2, double *sumtable, int *wgt)
874{
875  int i, l;
876  double *d1, *d_start, *sum;
877  double 
878    e[20] __attribute__ ((aligned (BYTE_ALIGNMENT))), 
879    s[20] __attribute__ ((aligned (BYTE_ALIGNMENT))), 
880    dd[20] __attribute__ ((aligned (BYTE_ALIGNMENT)));
881  double inv_Li, dlnLidlz, d2lnLidlz2;
882  double  dlnLdlz = 0.0;
883  double  d2lnLdlz2 = 0.0;
884
885  d1 = d_start = (double *)rax_malloc(numberOfCategories * 20 * sizeof(double));
886
887  e[0] = 0.0;
888  s[0] = 0.0; 
889
890  for(l = 1; l < 20; l++)
891    {
892      e[l]  = EIGN[l-1] * EIGN[l-1];
893      s[l]  = EIGN[l-1];
894      dd[l] = s[l] * lz;
895    }
896
897  for(i = 0; i < numberOfCategories; i++)
898    {     
899      d1[20 * i] = 1.0;
900      for(l = 1; l < 20; l++)
901        d1[20 * i + l] = EXP(dd[l] * rptr[i]);
902    }
903
904  for (i = 0; i < upper; i++)
905    {
906      double
907        r = rptr[cptr[i]],
908        wr1 = r * wgt[i],
909        wr2 = r * r * wgt[i];
910
911      __m128d a0 = _mm_setzero_pd();
912      __m128d a1 = _mm_setzero_pd();
913      __m128d a2 = _mm_setzero_pd();
914
915      d1 = &d_start[20 * cptr[i]];
916      sum = &sumtable[20 * i];
917         
918      for(l = 0; l < 20; l+=2)
919        {         
920          __m128d tmpv = _mm_mul_pd(_mm_load_pd(&d1[l]), _mm_load_pd(&sum[l]));
921         
922          a0 = _mm_add_pd(a0, tmpv);
923          __m128d sv = _mm_load_pd(&s[l]);       
924         
925          a1 = _mm_add_pd(a1, _mm_mul_pd(tmpv, sv));
926          __m128d ev = _mm_load_pd(&e[l]);       
927
928          a2 = _mm_add_pd(a2, _mm_mul_pd(tmpv, ev));
929        }
930
931      a0 = _mm_hadd_pd(a0, a0);
932      a1 = _mm_hadd_pd(a1, a1);
933      a2 = _mm_hadd_pd(a2, a2);
934
935      _mm_storel_pd(&inv_Li, a0);     
936      _mm_storel_pd(&dlnLidlz, a1);                 
937      _mm_storel_pd(&d2lnLidlz2, a2);
938     
939      inv_Li = 1.0/inv_Li;
940
941      dlnLidlz   *= inv_Li;
942      d2lnLidlz2 *= inv_Li;
943
944      dlnLdlz  += wr1 * dlnLidlz;
945      d2lnLdlz2 += wr2 * (d2lnLidlz2 - dlnLidlz * dlnLidlz);
946    }
947
948  *ext_dlnLdlz   = dlnLdlz;
949  *ext_d2lnLdlz2 = d2lnLdlz2;
950
951  rax_free(d_start);
952}
953
954#else
955
956static void coreGTRCATPROT(double *EIGN, double lz, int numberOfCategories, double *rptr, int *cptr, int upper,
957                           volatile double *ext_dlnLdlz,  volatile double *ext_d2lnLdlz2, double *sumtable, int *wgt)
958{
959  int i, l;
960  double *d1, *d_start, *sum;
961  double e[20], s[20], dd[20];
962  double tmp;
963  double inv_Li, dlnLidlz, d2lnLidlz2;
964  double  dlnLdlz = 0.0;
965  double  d2lnLdlz2 = 0.0;
966
967  d1 = d_start = (double *)rax_malloc(numberOfCategories * 20 * sizeof(double));
968
969  for(l = 1; l < 20; l++)
970    {
971      e[l]  = EIGN[l-1] * EIGN[l-1];
972      s[l]  = EIGN[l-1];
973      dd[l] = s[l] * lz;
974    }
975
976  for(i = 0; i < numberOfCategories; i++)
977    {
978      for(l = 1; l < 20; l++)
979        d1[20 * i + l] = EXP(dd[l] * rptr[i]);
980    }
981
982  for (i = 0; i < upper; i++)
983    {
984      double
985        r = rptr[cptr[i]],
986        wr1 = r * wgt[i],
987        wr2 = r * r * wgt[i];
988     
989      d1 = &d_start[20 * cptr[i]];
990      sum = &sumtable[20 * i];
991
992      inv_Li     = sum[0];
993      dlnLidlz   = 0.0;
994      d2lnLidlz2 = 0.0;
995
996      for(l = 1; l < 20; l++)
997        {
998          inv_Li     += (tmp = d1[l] * sum[l]);
999          dlnLidlz   += tmp *  s[l];
1000          d2lnLidlz2 += tmp *  e[l];
1001        }
1002
1003      inv_Li = 1.0/inv_Li;
1004
1005      dlnLidlz   *= inv_Li;
1006      d2lnLidlz2 *= inv_Li;
1007
1008      dlnLdlz  += wr1 * dlnLidlz;
1009      d2lnLdlz2 += wr2 * (d2lnLidlz2 - dlnLidlz * dlnLidlz);
1010    }
1011
1012  *ext_dlnLdlz   = dlnLdlz;
1013  *ext_d2lnLdlz2 = d2lnLdlz2;
1014
1015  rax_free(d_start);
1016}
1017
1018#endif
1019
1020
1021
1022static void coreCatFlex(double *EIGN, double lz, int numberOfCategories, double *rptr, int *cptr, int upper,
1023                        volatile double *ext_dlnLdlz,  volatile double *ext_d2lnLdlz2, double *sumtable,
1024                        const int numStates, int *wgt)
1025{
1026  int i, l;
1027  double 
1028    *d1, 
1029    *d_start, 
1030    *sum,
1031    e[64], 
1032    s[64], 
1033    dd[64],
1034    tmp,
1035    inv_Li, 
1036    dlnLidlz, 
1037    d2lnLidlz2,
1038    dlnLdlz = 0.0,
1039    d2lnLdlz2 = 0.0;
1040
1041  d1 = d_start = (double *)rax_malloc(numberOfCategories * numStates * sizeof(double));
1042
1043  for(l = 1; l < numStates; l++)
1044    {
1045      e[l]  = EIGN[l-1] * EIGN[l-1];
1046      s[l]  = EIGN[l-1];
1047      dd[l] = s[l] * lz;
1048    }
1049
1050  for(i = 0; i < numberOfCategories; i++)
1051    {
1052      for(l = 1; l < numStates; l++)
1053        d1[numStates * i + l] = EXP(dd[l] * rptr[i]);
1054    }
1055
1056  for (i = 0; i < upper; i++)
1057    {
1058      double
1059        r = rptr[cptr[i]],
1060        wr1 = r * wgt[i],
1061        wr2 = r * r * wgt[i];
1062
1063      d1 = &d_start[numStates * cptr[i]];
1064      sum = &sumtable[numStates * i];
1065
1066      inv_Li     = sum[0];
1067      dlnLidlz   = 0.0;
1068      d2lnLidlz2 = 0.0;
1069
1070      for(l = 1; l < numStates; l++)
1071        {
1072          inv_Li     += (tmp = d1[l] * sum[l]);
1073          dlnLidlz   += tmp *  s[l];
1074          d2lnLidlz2 += tmp *  e[l];
1075        }
1076
1077      inv_Li = 1.0/inv_Li;
1078
1079      dlnLidlz   *= inv_Li;
1080      d2lnLidlz2 *= inv_Li;
1081
1082      dlnLdlz  += wr1 * dlnLidlz;
1083      d2lnLdlz2 += wr2 * (d2lnLidlz2 - dlnLidlz * dlnLidlz);
1084    }
1085
1086  *ext_dlnLdlz   = dlnLdlz;
1087  *ext_d2lnLdlz2 = d2lnLdlz2;
1088
1089  rax_free(d_start);
1090}
1091
1092static void coreGammaFlex(double *gammaRates, double *EIGN, double *sumtable, int upper, int *wrptr,
1093                          volatile double *ext_dlnLdlz,  volatile double *ext_d2lnLdlz2, double lz, const int numStates)
1094{
1095  double 
1096    *sum, 
1097    diagptable[1024],
1098    dlnLdlz = 0.0,
1099    d2lnLdlz2 = 0.0,
1100    ki, 
1101    kisqr,
1102    tmp,
1103    inv_Li, 
1104    dlnLidlz, 
1105    d2lnLidlz2;
1106
1107  int     
1108    i, 
1109    j, 
1110    l; 
1111
1112  const int 
1113    gammaStates = 4 * numStates;
1114
1115  for(i = 0; i < 4; i++)
1116    {
1117      ki = gammaRates[i];
1118      kisqr = ki * ki;
1119
1120      for(l = 1; l < numStates; l++)
1121        {
1122          diagptable[i * gammaStates + l * 4]     = EXP(EIGN[l-1] * ki * lz);
1123          diagptable[i * gammaStates + l * 4 + 1] = EIGN[l-1] * ki;
1124          diagptable[i * gammaStates + l * 4 + 2] = EIGN[l-1] * EIGN[l-1] * kisqr;
1125        }
1126    }
1127
1128  for (i = 0; i < upper; i++)
1129    {
1130      sum = &sumtable[i * gammaStates];
1131      inv_Li   = 0.0;
1132      dlnLidlz = 0.0;
1133      d2lnLidlz2 = 0.0;
1134
1135      for(j = 0; j < 4; j++)
1136        {
1137          inv_Li += sum[j * numStates];
1138
1139          for(l = 1; l < numStates; l++)
1140            {
1141              inv_Li     += (tmp = diagptable[j * gammaStates + l * 4] * sum[j * numStates + l]);
1142              dlnLidlz   +=  tmp * diagptable[j * gammaStates + l * 4 + 1];
1143              d2lnLidlz2 +=  tmp * diagptable[j * gammaStates + l * 4 + 2];
1144            }
1145        }
1146
1147      inv_Li = 1.0 / inv_Li;
1148
1149      dlnLidlz   *= inv_Li;
1150      d2lnLidlz2 *= inv_Li;
1151
1152      dlnLdlz   += wrptr[i] * dlnLidlz;
1153      d2lnLdlz2 += wrptr[i] * (d2lnLidlz2 - dlnLidlz * dlnLidlz);
1154    }
1155
1156  *ext_dlnLdlz   = dlnLdlz;
1157  *ext_d2lnLdlz2 = d2lnLdlz2;
1158}
1159
1160
1161
1162
1163static void coreGammaInvarFlex(double *gammaRates, double *EIGN, double *sumtable, int upper, int *wrptr,
1164                               volatile double *ext_dlnLdlz,  volatile double *ext_d2lnLdlz2, double lz, double *frequencies,
1165                               double propInvar, int *iptr, const int numStates)
1166{
1167  double 
1168    *sum, 
1169    diagptable[1024],
1170    dlnLdlz = 0.0,
1171    d2lnLdlz2 = 0.0,
1172    ki, 
1173    kisqr,
1174    freqs[64],
1175    scaler =  0.25 * (1.0 - propInvar),
1176    tmp,
1177    inv_Li, 
1178    dlnLidlz, 
1179    d2lnLidlz2;
1180 
1181  int     
1182    i, 
1183    l, 
1184    j;
1185 
1186  const int 
1187    gammaStates = 4 * numStates;
1188 
1189  for(i = 0; i < numStates; i++)
1190    freqs[i] = frequencies[i] * propInvar;
1191
1192  for(i = 0; i < 4; i++)
1193    {
1194      ki = gammaRates[i];
1195      kisqr = ki * ki;
1196
1197      for(l = 1; l < numStates; l++)
1198        {
1199          diagptable[i * gammaStates + l * 4]     = EXP(EIGN[l-1] * ki * lz);
1200          diagptable[i * gammaStates + l * 4 + 1] = EIGN[l-1] * ki;
1201          diagptable[i * gammaStates + l * 4 + 2] = EIGN[l-1] * EIGN[l-1] * kisqr;
1202        }
1203    }
1204
1205
1206   for(i = 0; i < upper; i++)
1207     {
1208       sum = &sumtable[i * gammaStates];
1209       inv_Li   = 0.0;
1210       dlnLidlz = 0.0;
1211       d2lnLidlz2 = 0.0;
1212
1213       for(j = 0; j < 4; j++)
1214        {
1215          inv_Li += sum[j * numStates];
1216
1217          for(l = 1; l < numStates; l++)
1218            {
1219              inv_Li     += (tmp = diagptable[j * gammaStates + l * 4] * sum[j * numStates + l]);
1220              dlnLidlz   +=  tmp * diagptable[j * gammaStates + l * 4 + 1];
1221              d2lnLidlz2 +=  tmp * diagptable[j * gammaStates + l * 4 + 2];
1222            }
1223        }
1224
1225       inv_Li *= scaler;
1226
1227       if(iptr[i] < numStates)
1228         inv_Li += freqs[iptr[i]];
1229
1230       inv_Li = 1.0 / inv_Li;
1231
1232       dlnLidlz   *= inv_Li;
1233       d2lnLidlz2 *= inv_Li;
1234
1235       dlnLidlz *= scaler;
1236       d2lnLidlz2 *= scaler;
1237
1238       dlnLdlz  += wrptr[i] * dlnLidlz;
1239       d2lnLdlz2 += wrptr[i] * (d2lnLidlz2 - dlnLidlz * dlnLidlz);
1240     }
1241
1242   *ext_dlnLdlz   = dlnLdlz;
1243   *ext_d2lnLdlz2 = d2lnLdlz2;
1244}
1245
1246
1247static void coreGTRCATSECONDARY(double *EIGN, double lz, int numberOfCategories, double *rptr, int *cptr, int upper,
1248                                volatile double *ext_dlnLdlz,  volatile double *ext_d2lnLdlz2, double *sumtable, int *wgt)
1249{
1250  int i, l;
1251  double *d1, *d_start, *sum;
1252  double e[16], s[16], dd[16];
1253  double tmp;
1254  double inv_Li, dlnLidlz, d2lnLidlz2;
1255  double  dlnLdlz = 0.0;
1256  double  d2lnLdlz2 = 0.0;
1257
1258  d1 = d_start = (double *)rax_malloc(numberOfCategories * 16 * sizeof(double));
1259
1260  for(l = 1; l < 16; l++)
1261    {
1262      e[l]  = EIGN[l-1] * EIGN[l-1];
1263      s[l]  = EIGN[l-1];
1264      dd[l] = s[l] * lz;
1265    }
1266
1267  for(i = 0; i < numberOfCategories; i++)
1268    {
1269      for(l = 1; l < 16; l++)
1270        d1[16 * i + l] = EXP(dd[l] * rptr[i]);
1271    }
1272
1273  for (i = 0; i < upper; i++)
1274    {
1275      double
1276        r = rptr[cptr[i]],
1277        wr1 = r * wgt[i],
1278        wr2 = r * r * wgt[i];
1279     
1280      d1 = &d_start[16 * cptr[i]];
1281      sum = &sumtable[16 * i];
1282
1283      inv_Li     = sum[0];
1284      dlnLidlz   = 0.0;
1285      d2lnLidlz2 = 0.0;
1286
1287      for(l = 1; l < 16; l++)
1288        {
1289          inv_Li     += (tmp = d1[l] * sum[l]);
1290          dlnLidlz   += tmp *  s[l];
1291          d2lnLidlz2 += tmp *  e[l];
1292        }
1293
1294      inv_Li = 1.0/inv_Li;
1295
1296      dlnLidlz   *= inv_Li;
1297      d2lnLidlz2 *= inv_Li;
1298
1299      dlnLdlz  += wr1 * dlnLidlz;
1300      d2lnLdlz2 += wr2 * (d2lnLidlz2 - dlnLidlz * dlnLidlz);
1301    }
1302
1303  *ext_dlnLdlz   = dlnLdlz;
1304  *ext_d2lnLdlz2 = d2lnLdlz2;
1305
1306  rax_free(d_start);
1307}
1308
1309static void coreGTRCATSECONDARY_6(double *EIGN, double lz, int numberOfCategories, double *rptr, int *cptr, int upper,
1310                                  volatile double *ext_dlnLdlz,  volatile double *ext_d2lnLdlz2, double *sumtable, int *wgt)
1311{
1312  int i, l;
1313  double *d1, *d_start, *sum;
1314  double e[6], s[6], dd[6];
1315  double tmp;
1316  double inv_Li, dlnLidlz, d2lnLidlz2;
1317  double  dlnLdlz = 0.0;
1318  double  d2lnLdlz2 = 0.0;
1319
1320  d1 = d_start = (double *)rax_malloc(numberOfCategories * 6 * sizeof(double));
1321
1322  for(l = 1; l < 6; l++)
1323    {
1324      e[l]  = EIGN[l-1] * EIGN[l-1];
1325      s[l]  = EIGN[l-1];
1326      dd[l] = s[l] * lz;
1327    }
1328
1329  for(i = 0; i < numberOfCategories; i++)
1330    {
1331      for(l = 1; l < 6; l++)
1332        d1[6 * i + l] = EXP(dd[l] * rptr[i]);
1333    }
1334
1335  for (i = 0; i < upper; i++)
1336    {
1337      double
1338        r = rptr[cptr[i]],
1339        wr1 = r * wgt[i],
1340        wr2 = r * r * wgt[i];
1341     
1342      d1 = &d_start[6 * cptr[i]];
1343      sum = &sumtable[6 * i];
1344
1345      inv_Li     = sum[0];
1346      dlnLidlz   = 0.0;
1347      d2lnLidlz2 = 0.0;
1348
1349      for(l = 1; l < 6; l++)
1350        {
1351          inv_Li     += (tmp = d1[l] * sum[l]);
1352          dlnLidlz   += tmp *  s[l];
1353          d2lnLidlz2 += tmp *  e[l];
1354        }
1355
1356      inv_Li = 1.0/inv_Li;
1357
1358      dlnLidlz   *= inv_Li;
1359      d2lnLidlz2 *= inv_Li;
1360
1361      dlnLdlz  += wr1 * dlnLidlz;
1362      d2lnLdlz2 += wr2 * (d2lnLidlz2 - dlnLidlz * dlnLidlz);
1363    }
1364
1365  *ext_dlnLdlz   = dlnLdlz;
1366  *ext_d2lnLdlz2 = d2lnLdlz2;
1367
1368  rax_free(d_start);
1369}
1370
1371static void coreGTRCATSECONDARY_7(double *EIGN, double lz, int numberOfCategories, double *rptr, int *cptr, int upper,
1372                                  volatile double *ext_dlnLdlz,  volatile double *ext_d2lnLdlz2, double *sumtable, int *wgt)
1373{
1374  int i, l;
1375  double *d1, *d_start, *sum;
1376  double e[7], s[7], dd[7];
1377  double tmp;
1378  double inv_Li, dlnLidlz, d2lnLidlz2;
1379  double  dlnLdlz = 0.0;
1380  double  d2lnLdlz2 = 0.0;
1381
1382  d1 = d_start = (double *)rax_malloc(numberOfCategories * 7 * sizeof(double));
1383
1384  for(l = 1; l < 7; l++)
1385    {
1386      e[l]  = EIGN[l-1] * EIGN[l-1];
1387      s[l]  = EIGN[l-1];
1388      dd[l] = s[l] * lz;
1389    }
1390
1391  for(i = 0; i < numberOfCategories; i++)
1392    {
1393      for(l = 1; l < 7; l++)
1394        d1[7 * i + l] = EXP(dd[l] * rptr[i]);
1395    }
1396
1397  for (i = 0; i < upper; i++)
1398    {
1399      double
1400        r = rptr[cptr[i]],
1401        wr1 = r * wgt[i],
1402        wr2 = r * r * wgt[i];
1403
1404      d1 = &d_start[7 * cptr[i]];
1405      sum = &sumtable[7 * i];
1406
1407      inv_Li     = sum[0];
1408      dlnLidlz   = 0.0;
1409      d2lnLidlz2 = 0.0;
1410
1411      for(l = 1; l < 7; l++)
1412        {
1413          inv_Li     += (tmp = d1[l] * sum[l]);
1414          dlnLidlz   += tmp *  s[l];
1415          d2lnLidlz2 += tmp *  e[l];
1416        }
1417
1418      inv_Li = 1.0/inv_Li;
1419
1420      dlnLidlz   *= inv_Li;
1421      d2lnLidlz2 *= inv_Li;
1422
1423      dlnLdlz  += wr1 * dlnLidlz;
1424      d2lnLdlz2 += wr2 * (d2lnLidlz2 - dlnLidlz * dlnLidlz);
1425    }
1426
1427  *ext_dlnLdlz   = dlnLdlz;
1428  *ext_d2lnLdlz2 = d2lnLdlz2;
1429
1430  rax_free(d_start);
1431}
1432
1433
1434
1435static void coreGTRGAMMAINVAR_BINARY(double propInvar, double *frequencies, double *gammaRates, double *EIGN,
1436                                     double *sumtable,  volatile double *ext_dlnLdlz,  volatile double *ext_d2lnLdlz2,
1437                                     int *iptr, int *wrptr, int upper, double lz)
1438{
1439  double  *sum, diagptable[12];
1440  int     i, j;
1441    double  dlnLdlz = 0;
1442  double d2lnLdlz2 = 0;
1443  double ki, kisqr;
1444  double freqs[4];
1445  double scaler =  0.25 * (1.0 - propInvar);
1446  double tmp_1;
1447  double inv_Li, dlnLidlz, d2lnLidlz2;
1448
1449  freqs[0] = frequencies[0] * propInvar;
1450  freqs[1] = frequencies[1] * propInvar;
1451
1452  for(i = 0; i < 4; i++)
1453    {
1454      ki = gammaRates[i];
1455      kisqr = ki * ki;
1456
1457      diagptable[i * 3]     = EXP (EIGN[0] * ki * lz);
1458      diagptable[i * 3 + 1] = EIGN[0] * ki;
1459      diagptable[i * 3 + 2] = EIGN[0] * EIGN[0] * kisqr;
1460    }
1461
1462  for (i = 0; i < upper; i++)
1463    {
1464      sum = &(sumtable[i * 8]);
1465
1466      inv_Li      = 0.0;
1467      dlnLidlz    = 0.0;
1468      d2lnLidlz2  = 0.0;
1469
1470      for(j = 0; j < 4; j++)
1471        {
1472          inv_Li += sum[2 * j];
1473
1474
1475          tmp_1      =  diagptable[3 * j] * sum[2 * j + 1];
1476          inv_Li     += tmp_1;
1477          dlnLidlz   += tmp_1 * diagptable[3 * j + 1];
1478          d2lnLidlz2 += tmp_1 * diagptable[3 * j + 2];
1479         }
1480
1481      inv_Li *= scaler;
1482
1483      if(iptr[i] < 2)
1484        inv_Li += freqs[iptr[i]];
1485
1486      inv_Li = 1.0 / inv_Li;
1487
1488      dlnLidlz   *= inv_Li;
1489      d2lnLidlz2 *= inv_Li;
1490
1491      dlnLidlz *= scaler;
1492      d2lnLidlz2 *= scaler;
1493
1494      dlnLdlz  += wrptr[i] * dlnLidlz;
1495      d2lnLdlz2 += wrptr[i] * (d2lnLidlz2 - dlnLidlz * dlnLidlz);
1496    }
1497
1498  *ext_dlnLdlz   = dlnLdlz;
1499  *ext_d2lnLdlz2 = d2lnLdlz2;
1500}
1501
1502
1503
1504static void coreGTRGAMMAINVAR(double propInvar, double *frequencies, double *gammaRates, double *EIGN,
1505                              double *sumtable,  volatile double *ext_dlnLdlz,  volatile double *ext_d2lnLdlz2,
1506                              int *iptr, int *wrptr, int upper, double lz)
1507{
1508  double  *sum, diagptable[64];
1509  int     i, j, k;
1510    double  dlnLdlz = 0;
1511  double d2lnLdlz2 = 0;
1512  double ki, kisqr;
1513  double freqs[4];
1514  double scaler =  0.25 * (1.0 - propInvar);
1515  double tmp_1;
1516  double inv_Li, dlnLidlz, d2lnLidlz2;
1517
1518  freqs[0] = frequencies[0] * propInvar;
1519  freqs[1] = frequencies[1] * propInvar;
1520  freqs[2] = frequencies[2] * propInvar;
1521  freqs[3] = frequencies[3] * propInvar;
1522
1523  for(i = 0; i < 4; i++)
1524    {
1525      ki = gammaRates[i];
1526      kisqr = ki * ki;
1527
1528      diagptable[i * 16]     = EXP (EIGN[0] * ki * lz);
1529      diagptable[i * 16 + 1] = EXP (EIGN[1] * ki * lz);
1530      diagptable[i * 16 + 2] = EXP (EIGN[2] * ki * lz);
1531
1532      diagptable[i * 16 + 3] = EIGN[0] * ki;
1533      diagptable[i * 16 + 4] = EIGN[0] * EIGN[0] * kisqr;
1534
1535      diagptable[i * 16 + 5] = EIGN[1] * ki;
1536      diagptable[i * 16 + 6] = EIGN[1] * EIGN[1] * kisqr;
1537
1538      diagptable[i * 16 + 7] = EIGN[2] * ki;
1539      diagptable[i * 16 + 8] = EIGN[2] * EIGN[2] * kisqr;
1540    }
1541
1542  for (i = 0; i < upper; i++)
1543    {
1544      sum = &(sumtable[i * 16]);
1545
1546       inv_Li      = 0.0;
1547       dlnLidlz    = 0.0;
1548       d2lnLidlz2  = 0.0;
1549
1550       for(j = 0; j < 4; j++)
1551         {
1552           inv_Li += sum[4 * j];
1553
1554           for(k = 0; k < 3; k++)
1555             {
1556               tmp_1      =  diagptable[16 * j + k] * sum[4 * j + k + 1];
1557               inv_Li     += tmp_1;
1558               dlnLidlz   += tmp_1 * diagptable[16 * j + k * 2 + 3];
1559               d2lnLidlz2 += tmp_1 * diagptable[16 * j + k * 2 + 4];
1560             }
1561         }
1562
1563      inv_Li *= scaler;
1564
1565      if(iptr[i] < 4)
1566        inv_Li += freqs[iptr[i]];
1567
1568      inv_Li = 1.0 / inv_Li;
1569
1570      dlnLidlz   *= inv_Li;
1571      d2lnLidlz2 *= inv_Li;
1572
1573      dlnLidlz *= scaler;
1574      d2lnLidlz2 *= scaler;
1575
1576      dlnLdlz  += wrptr[i] * dlnLidlz;
1577      d2lnLdlz2 += wrptr[i] * (d2lnLidlz2 - dlnLidlz * dlnLidlz);
1578    }
1579
1580  *ext_dlnLdlz   = dlnLdlz;
1581  *ext_d2lnLdlz2 = d2lnLdlz2;
1582}
1583
1584
1585/* C-OPT the two functions below are used to optimize the branch lengths of the tree.
1586   together, they account for approx 25-30% of overall run time sumGAMMA requires less
1587   overall time than coreGAMMA though. */
1588
1589static void sumGAMMA_BINARY(int tipCase, double *sumtable, double *x1_start, double *x2_start, double *tipVector,
1590                            unsigned char *tipX1, unsigned char *tipX2, int n)
1591{
1592  double *x1, *x2, *sum;
1593  int i, j;
1594#ifndef __SIM_SSE3
1595  int k;
1596#endif
1597
1598  /* C-OPT once again switch over possible configurations at inner node */
1599
1600  switch(tipCase)
1601    {
1602    case TIP_TIP:
1603      /* C-OPT main for loop overt alignment length */
1604      for (i = 0; i < n; i++)
1605        {
1606          x1 = &(tipVector[2 * tipX1[i]]);
1607          x2 = &(tipVector[2 * tipX2[i]]);
1608          sum = &sumtable[i * 8];
1609#ifndef __SIM_SSE3       
1610          for(j = 0; j < 4; j++)
1611            for(k = 0; k < 2; k++)
1612              sum[j * 2 + k] = x1[k] * x2[k];
1613#else
1614          for(j = 0; j < 4; j++)
1615            _mm_store_pd( &sum[j*2], _mm_mul_pd( _mm_load_pd( &x1[0] ), _mm_load_pd( &x2[0] )));         
1616#endif
1617        }
1618      break;
1619    case TIP_INNER:
1620      for (i = 0; i < n; i++)
1621        {
1622          x1  = &(tipVector[2 * tipX1[i]]);
1623          x2  = &x2_start[8 * i];
1624          sum = &sumtable[8 * i];
1625
1626#ifndef __SIM_SSE3
1627          for(j = 0; j < 4; j++)
1628            for(k = 0; k < 2; k++)
1629              sum[j * 2 + k] = x1[k] * x2[j * 2 + k];
1630#else
1631          for(j = 0; j < 4; j++)
1632            _mm_store_pd( &sum[j*2], _mm_mul_pd( _mm_load_pd( &x1[0] ), _mm_load_pd( &x2[j * 2] )));
1633#endif
1634        }
1635      break;
1636    case INNER_INNER:
1637      for (i = 0; i < n; i++)
1638        {
1639          x1  = &x1_start[8 * i];
1640          x2  = &x2_start[8 * i];
1641          sum = &sumtable[8 * i];
1642#ifndef __SIM_SSE3
1643          for(j = 0; j < 4; j++)
1644            for(k = 0; k < 2; k++)
1645              sum[j * 2 + k] = x1[j * 2 + k] * x2[j * 2 + k];
1646#else
1647          for(j = 0; j < 4; j++)
1648            _mm_store_pd( &sum[j*2], _mm_mul_pd( _mm_load_pd( &x1[j * 2] ), _mm_load_pd( &x2[j * 2] )));
1649#endif
1650        }
1651      break;
1652    default:
1653      assert(0);
1654    }
1655}
1656
1657
1658
1659static void sumGAMMA_GAPPED_SAVE(int tipCase, double *sumtable, double *x1_start, double *x2_start, double *tipVector,
1660                                 unsigned char *tipX1, unsigned char *tipX2, int n, 
1661                                 double *x1_gapColumn, double *x2_gapColumn, unsigned int *x1_gap, unsigned int *x2_gap)
1662{
1663  double 
1664    *x1, 
1665    *x2, 
1666    *sum,
1667    *x1_ptr = x1_start,
1668    *x2_ptr = x2_start;
1669 
1670  int i, j, k; 
1671
1672  switch(tipCase)
1673    {
1674    case TIP_TIP:     
1675      for (i = 0; i < n; i++)
1676        {
1677          x1 = &(tipVector[4 * tipX1[i]]);
1678          x2 = &(tipVector[4 * tipX2[i]]);
1679          sum = &sumtable[i * 16];
1680#ifndef __SIM_SSE3
1681          for(j = 0; j < 4; j++)           
1682            for(k = 0; k < 4; k++)
1683              sum[j * 4 + k] = x1[k] * x2[k];
1684#else
1685          for(j = 0; j < 4; j++)           
1686            for(k = 0; k < 4; k+=2)
1687              _mm_store_pd( &sum[j*4 + k], _mm_mul_pd( _mm_load_pd( &x1[k] ), _mm_load_pd( &x2[k] )));
1688#endif
1689        }
1690      break;
1691    case TIP_INNER:
1692      for (i = 0; i < n; i++)
1693        {
1694          x1  = &(tipVector[4 * tipX1[i]]);
1695         
1696          if(x2_gap[i / 32] & mask32[i % 32])
1697            x2 = x2_gapColumn;
1698          else
1699            {
1700              x2  = x2_ptr;
1701              x2_ptr += 16;
1702            }
1703         
1704          sum = &sumtable[16 * i];
1705#ifndef __SIM_SSE3
1706          for(j = 0; j < 4; j++)
1707            for(k = 0; k < 4; k++)
1708              sum[j * 4 + k] = x1[k] * x2[j * 4 + k];
1709#else
1710          for(j = 0; j < 4; j++)           
1711            for(k = 0; k < 4; k+=2)
1712              _mm_store_pd( &sum[j*4 + k], _mm_mul_pd( _mm_load_pd( &x1[k] ), _mm_load_pd( &x2[j * 4 + k] )));
1713#endif
1714        }
1715      break;
1716    case INNER_INNER:
1717      for (i = 0; i < n; i++)
1718        {
1719          if(x1_gap[i / 32] & mask32[i % 32])
1720            x1 = x1_gapColumn;
1721          else
1722            {
1723              x1  = x1_ptr;
1724              x1_ptr += 16;
1725            }
1726         
1727          if(x2_gap[i / 32] & mask32[i % 32])
1728            x2 = x2_gapColumn;
1729          else
1730            {
1731              x2  = x2_ptr;
1732              x2_ptr += 16;
1733            }
1734
1735          sum = &sumtable[16 * i];
1736         
1737#ifndef __SIM_SSE3
1738          for(j = 0; j < 4; j++)
1739            for(k = 0; k < 4; k++)
1740              sum[j * 4 + k] = x1[j * 4 + k] * x2[j * 4 + k];
1741#else
1742           for(j = 0; j < 4; j++)           
1743            for(k = 0; k < 4; k+=2)
1744              _mm_store_pd( &sum[j*4 + k], _mm_mul_pd( _mm_load_pd( &x1[j * 4 + k] ), _mm_load_pd( &x2[j * 4 + k] )));
1745#endif
1746        }
1747      break;
1748    default:
1749      assert(0);
1750    }
1751}
1752
1753
1754
1755
1756
1757static void sumGAMMA(int tipCase, double *sumtable, double *x1_start, double *x2_start, double *tipVector,
1758                     unsigned char *tipX1, unsigned char *tipX2, int n)
1759{
1760  double *x1, *x2, *sum;
1761  int i, j, k;
1762
1763  /* C-OPT once again switch over possible configurations at inner node */
1764
1765  switch(tipCase)
1766    {
1767    case TIP_TIP:
1768      /* C-OPT main for loop overt alignment length */
1769      for (i = 0; i < n; i++)
1770        {
1771          x1 = &(tipVector[4 * tipX1[i]]);
1772          x2 = &(tipVector[4 * tipX2[i]]);
1773          sum = &sumtable[i * 16];
1774#ifndef __SIM_SSE3
1775          for(j = 0; j < 4; j++)           
1776            for(k = 0; k < 4; k++)
1777              sum[j * 4 + k] = x1[k] * x2[k];
1778#else
1779          for(j = 0; j < 4; j++)           
1780            for(k = 0; k < 4; k+=2)
1781              _mm_store_pd( &sum[j*4 + k], _mm_mul_pd( _mm_load_pd( &x1[k] ), _mm_load_pd( &x2[k] )));
1782#endif
1783        }
1784      break;
1785    case TIP_INNER:
1786      for (i = 0; i < n; i++)
1787        {
1788          x1  = &(tipVector[4 * tipX1[i]]);
1789          x2  = &x2_start[16 * i];
1790          sum = &sumtable[16 * i];
1791#ifndef __SIM_SSE3
1792          for(j = 0; j < 4; j++)
1793            for(k = 0; k < 4; k++)
1794              sum[j * 4 + k] = x1[k] * x2[j * 4 + k];
1795#else
1796          for(j = 0; j < 4; j++)           
1797            for(k = 0; k < 4; k+=2)
1798              _mm_store_pd( &sum[j*4 + k], _mm_mul_pd( _mm_load_pd( &x1[k] ), _mm_load_pd( &x2[j * 4 + k] )));
1799#endif
1800        }
1801      break;
1802    case INNER_INNER:
1803      for (i = 0; i < n; i++)
1804        {
1805          x1  = &x1_start[16 * i];
1806          x2  = &x2_start[16 * i];
1807          sum = &sumtable[16 * i];
1808#ifndef __SIM_SSE3
1809          for(j = 0; j < 4; j++)
1810            for(k = 0; k < 4; k++)
1811              sum[j * 4 + k] = x1[j * 4 + k] * x2[j * 4 + k];
1812#else
1813           for(j = 0; j < 4; j++)           
1814            for(k = 0; k < 4; k+=2)
1815              _mm_store_pd( &sum[j*4 + k], _mm_mul_pd( _mm_load_pd( &x1[j * 4 + k] ), _mm_load_pd( &x2[j * 4 + k] )));
1816#endif
1817        }
1818      break;
1819    default:
1820      assert(0);
1821    }
1822}
1823
1824
1825
1826
1827#ifndef __SIM_SSE3
1828static void coreGTRGAMMA_BINARY(const int upper, double *sumtable,
1829                                volatile double *d1,   volatile double *d2, double *EIGN, double *gammaRates, double lz, int *wrptr)
1830{
1831  int i, j;
1832  double
1833    *diagptable, *diagptable_start, *sum,
1834    tmp_1, inv_Li, dlnLidlz, d2lnLidlz2, ki, kisqr,
1835    dlnLdlz = 0.0,
1836    d2lnLdlz2 = 0.0;
1837
1838  diagptable = diagptable_start = (double *)rax_malloc(sizeof(double) * 12);
1839
1840  for(i = 0; i < 4; i++)
1841    {
1842      ki = gammaRates[i];
1843      kisqr = ki * ki;
1844
1845      diagptable[i * 3]     = EXP (EIGN[0] * ki * lz);
1846      diagptable[i * 3 + 1] = EIGN[0] * ki;
1847      diagptable[i * 3 + 2] = EIGN[0] * EIGN[0] * kisqr;
1848    }
1849
1850  for (i = 0; i < upper; i++)
1851    {
1852      diagptable = diagptable_start;
1853      sum = &(sumtable[i * 8]);
1854
1855      inv_Li      = 0.0;
1856      dlnLidlz    = 0.0;
1857      d2lnLidlz2  = 0.0;
1858
1859      for(j = 0; j < 4; j++)
1860        {
1861          inv_Li += sum[2 * j];
1862
1863          tmp_1      =  diagptable[3 * j] * sum[2 * j + 1];
1864          inv_Li     += tmp_1;
1865          dlnLidlz   += tmp_1 * diagptable[3 * j + 1];
1866          d2lnLidlz2 += tmp_1 * diagptable[3 * j + 2];
1867        }
1868
1869      inv_Li = 1.0 / inv_Li;
1870
1871      dlnLidlz   *= inv_Li;
1872      d2lnLidlz2 *= inv_Li;
1873
1874
1875      dlnLdlz  += wrptr[i] * dlnLidlz;
1876      d2lnLdlz2 += wrptr[i] * (d2lnLidlz2 - dlnLidlz * dlnLidlz);
1877    }
1878
1879  *d1 = dlnLdlz;
1880  *d2 = d2lnLdlz2;
1881
1882  rax_free(diagptable_start);
1883}
1884#else
1885static void coreGTRGAMMA_BINARY(const int upper, double *sumtable,
1886                                volatile double *d1,   volatile double *d2, double *EIGN, double *gammaRates, double lz, int *wrptr)
1887{
1888  double 
1889    dlnLdlz = 0.0,
1890    d2lnLdlz2 = 0.0,
1891    ki, 
1892    kisqr, 
1893    inv_Li, 
1894    dlnLidlz, 
1895    d2lnLidlz2, 
1896    *sum, 
1897    diagptable0[8] __attribute__ ((aligned (BYTE_ALIGNMENT))),
1898    diagptable1[8] __attribute__ ((aligned (BYTE_ALIGNMENT))),
1899    diagptable2[8] __attribute__ ((aligned (BYTE_ALIGNMENT)));   
1900   
1901  int     
1902    i, 
1903    j;
1904 
1905  for(i = 0; i < 4; i++)
1906    {
1907      ki = gammaRates[i];
1908      kisqr = ki * ki;
1909     
1910      diagptable0[i * 2] = 1.0;
1911      diagptable1[i * 2] = 0.0;
1912      diagptable2[i * 2] = 0.0;
1913     
1914      diagptable0[i * 2 + 1] = EXP(EIGN[0] * ki * lz);
1915      diagptable1[i * 2 + 1] = EIGN[0] * ki;
1916      diagptable2[i * 2 + 1] = EIGN[0] * EIGN[0] * kisqr;   
1917    }
1918
1919  for (i = 0; i < upper; i++)
1920    { 
1921      __m128d a0 = _mm_setzero_pd();
1922      __m128d a1 = _mm_setzero_pd();
1923      __m128d a2 = _mm_setzero_pd();
1924
1925      sum = &sumtable[i * 8];         
1926
1927      for(j = 0; j < 4; j++)
1928        {                       
1929          double           
1930            *d0 = &diagptable0[j * 2],
1931            *d1 = &diagptable1[j * 2],
1932            *d2 = &diagptable2[j * 2];
1933                         
1934          __m128d tmpv = _mm_mul_pd(_mm_load_pd(d0), _mm_load_pd(&sum[j * 2]));
1935          a0 = _mm_add_pd(a0, tmpv);
1936          a1 = _mm_add_pd(a1, _mm_mul_pd(tmpv, _mm_load_pd(d1)));
1937          a2 = _mm_add_pd(a2, _mm_mul_pd(tmpv, _mm_load_pd(d2)));
1938                         
1939        }
1940
1941      a0 = _mm_hadd_pd(a0, a0);
1942      a1 = _mm_hadd_pd(a1, a1);
1943      a2 = _mm_hadd_pd(a2, a2);
1944
1945      _mm_storel_pd(&inv_Li, a0);     
1946      _mm_storel_pd(&dlnLidlz, a1);
1947      _mm_storel_pd(&d2lnLidlz2, a2); 
1948
1949      inv_Li = 1.0 / inv_Li;
1950     
1951      dlnLidlz   *= inv_Li;
1952      d2lnLidlz2 *= inv_Li;     
1953
1954      dlnLdlz   += wrptr[i] * dlnLidlz;
1955      d2lnLdlz2 += wrptr[i] * (d2lnLidlz2 - dlnLidlz * dlnLidlz);
1956    }
1957
1958 
1959  *d1   = dlnLdlz;
1960  *d2 = d2lnLdlz2; 
1961}
1962
1963
1964#endif
1965
1966#ifndef __SIM_SSE3
1967static void coreGTRGAMMA(const int upper, double *sumtable,
1968                         volatile double *d1,   volatile double *d2, double *EIGN, double *gammaRates, double lz, int *wrptr)
1969{
1970  int i, j, k;
1971  double
1972    *diagptable, *diagptable_start, *sum,
1973    tmp_1, inv_Li, dlnLidlz, d2lnLidlz2, ki, kisqr,
1974    dlnLdlz = 0.0,
1975    d2lnLdlz2 = 0.0;
1976
1977  diagptable = diagptable_start = (double *)rax_malloc(sizeof(double) * 64);
1978
1979
1980
1981  for(i = 0; i < 4; i++)
1982    {
1983      ki = gammaRates[i];
1984      kisqr = ki * ki;
1985
1986      diagptable[i * 16]     = EXP (EIGN[0] * ki * lz);
1987      diagptable[i * 16 + 1] = EXP (EIGN[1] * ki * lz);
1988      diagptable[i * 16 + 2] = EXP (EIGN[2] * ki * lz);
1989
1990      diagptable[i * 16 + 3] = EIGN[0] * ki;
1991      diagptable[i * 16 + 4] = EIGN[0] * EIGN[0] * kisqr;
1992
1993      diagptable[i * 16 + 5] = EIGN[1] * ki;
1994      diagptable[i * 16 + 6] = EIGN[1] * EIGN[1] * kisqr;
1995
1996      diagptable[i * 16 + 7] = EIGN[2] * ki;
1997      diagptable[i * 16 + 8] = EIGN[2] * EIGN[2] * kisqr;
1998    }
1999
2000  for (i = 0; i < upper; i++)
2001    {
2002      sum = &(sumtable[i * 16]);
2003
2004      inv_Li      = 0.0;
2005      dlnLidlz    = 0.0;
2006      d2lnLidlz2  = 0.0;
2007
2008      for(j = 0; j < 4; j++)
2009        {
2010          inv_Li += sum[4 * j];
2011
2012          for(k = 0; k < 3; k++)
2013            {
2014              tmp_1      =  diagptable[16 * j + k] * sum[4 * j + k + 1];
2015              inv_Li     += tmp_1;
2016              dlnLidlz   += tmp_1 * diagptable[16 * j + k * 2 + 3];
2017              d2lnLidlz2 += tmp_1 * diagptable[16 * j + k * 2 + 4];
2018            }
2019        }
2020
2021      inv_Li = 1.0 / inv_Li;
2022
2023      dlnLidlz   *= inv_Li;
2024      d2lnLidlz2 *= inv_Li;
2025
2026
2027
2028      dlnLdlz  += wrptr[i] * dlnLidlz;
2029      d2lnLdlz2 += wrptr[i] * (d2lnLidlz2 - dlnLidlz * dlnLidlz);
2030    }
2031
2032  *d1 = dlnLdlz;
2033  *d2 = d2lnLdlz2;
2034
2035  rax_free(diagptable_start);
2036}
2037#else
2038
2039static void coreGTRGAMMA(const int upper, double *sumtable,
2040                         volatile double *ext_dlnLdlz,  volatile double *ext_d2lnLdlz2, double *EIGN, double *gammaRates, double lz, int *wrptr)
2041{
2042  double 
2043    dlnLdlz = 0.0,
2044    d2lnLdlz2 = 0.0,
2045    ki, 
2046    kisqr, 
2047    inv_Li, 
2048    dlnLidlz, 
2049    d2lnLidlz2, 
2050    *sum, 
2051    diagptable0[16] __attribute__ ((aligned (BYTE_ALIGNMENT))),
2052    diagptable1[16] __attribute__ ((aligned (BYTE_ALIGNMENT))),
2053    diagptable2[16] __attribute__ ((aligned (BYTE_ALIGNMENT)));   
2054   
2055  int     
2056    i, 
2057    j, 
2058    l;
2059 
2060  for(i = 0; i < 4; i++)
2061    {
2062      ki = gammaRates[i];
2063      kisqr = ki * ki;
2064     
2065      diagptable0[i * 4] = 1.0;
2066      diagptable1[i * 4] = 0.0;
2067      diagptable2[i * 4] = 0.0;
2068
2069      for(l = 1; l < 4; l++)
2070        {
2071          diagptable0[i * 4 + l] = EXP(EIGN[l-1] * ki * lz);
2072          diagptable1[i * 4 + l] = EIGN[l-1] * ki;
2073          diagptable2[i * 4 + l] = EIGN[l-1] * EIGN[l-1] * kisqr;
2074        }
2075    }
2076
2077  for (i = 0; i < upper; i++)
2078    { 
2079      __m128d a0 = _mm_setzero_pd();
2080      __m128d a1 = _mm_setzero_pd();
2081      __m128d a2 = _mm_setzero_pd();
2082
2083      sum = &sumtable[i * 16];         
2084
2085      for(j = 0; j < 4; j++)
2086        {                       
2087          double           
2088            *d0 = &diagptable0[j * 4],
2089            *d1 = &diagptable1[j * 4],
2090            *d2 = &diagptable2[j * 4];
2091                 
2092          for(l = 0; l < 4; l+=2)
2093            {
2094              __m128d tmpv = _mm_mul_pd(_mm_load_pd(&d0[l]), _mm_load_pd(&sum[j * 4 + l]));
2095              a0 = _mm_add_pd(a0, tmpv);
2096              a1 = _mm_add_pd(a1, _mm_mul_pd(tmpv, _mm_load_pd(&d1[l])));
2097              a2 = _mm_add_pd(a2, _mm_mul_pd(tmpv, _mm_load_pd(&d2[l])));
2098            }             
2099        }
2100
2101      a0 = _mm_hadd_pd(a0, a0);
2102      a1 = _mm_hadd_pd(a1, a1);
2103      a2 = _mm_hadd_pd(a2, a2);
2104
2105      _mm_storel_pd(&inv_Li, a0);     
2106      _mm_storel_pd(&dlnLidlz, a1);
2107      _mm_storel_pd(&d2lnLidlz2, a2); 
2108
2109      inv_Li = 1.0 / inv_Li;
2110     
2111      dlnLidlz   *= inv_Li;
2112      d2lnLidlz2 *= inv_Li;     
2113
2114      dlnLdlz   += wrptr[i] * dlnLidlz;
2115      d2lnLdlz2 += wrptr[i] * (d2lnLidlz2 - dlnLidlz * dlnLidlz);
2116    }
2117
2118 
2119  *ext_dlnLdlz   = dlnLdlz;
2120  *ext_d2lnLdlz2 = d2lnLdlz2; 
2121}
2122
2123
2124#endif
2125
2126
2127
2128
2129
2130
2131
2132
2133
2134static void sumGAMMAPROT(int tipCase, double *sumtable, double *x1, double *x2, double *tipVector,
2135                         unsigned char *tipX1, unsigned char *tipX2, int n)
2136{
2137  int i, l, k;
2138  double *left, *right, *sum;
2139
2140  switch(tipCase)
2141    {
2142    case TIP_TIP:
2143      for(i = 0; i < n; i++)
2144        {
2145          left  = &(tipVector[20 * tipX1[i]]);
2146          right = &(tipVector[20 * tipX2[i]]);
2147
2148          for(l = 0; l < 4; l++)
2149            {
2150              sum = &sumtable[i * 80 + l * 20];
2151#ifdef __SIM_SSE3
2152              for(k = 0; k < 20; k+=2)
2153                {
2154                  __m128d sumv = _mm_mul_pd(_mm_load_pd(&left[k]), _mm_load_pd(&right[k]));
2155                 
2156                  _mm_store_pd(&sum[k], sumv);           
2157                }
2158#else
2159              for(k = 0; k < 20; k++)
2160                sum[k] = left[k] * right[k];
2161#endif
2162            }
2163        }
2164      break;
2165    case TIP_INNER:
2166      for(i = 0; i < n; i++)
2167        {
2168          left = &(tipVector[20 * tipX1[i]]);
2169
2170          for(l = 0; l < 4; l++)
2171            {
2172              right = &(x2[80 * i + l * 20]);
2173              sum = &sumtable[i * 80 + l * 20];
2174#ifdef __SIM_SSE3
2175              for(k = 0; k < 20; k+=2)
2176                {
2177                  __m128d sumv = _mm_mul_pd(_mm_load_pd(&left[k]), _mm_load_pd(&right[k]));
2178                 
2179                  _mm_store_pd(&sum[k], sumv);           
2180                }
2181#else
2182              for(k = 0; k < 20; k++)
2183                sum[k] = left[k] * right[k];
2184#endif
2185            }
2186        }
2187      break;
2188    case INNER_INNER:
2189      for(i = 0; i < n; i++)
2190        {
2191          for(l = 0; l < 4; l++)
2192            {
2193              left  = &(x1[80 * i + l * 20]);
2194              right = &(x2[80 * i + l * 20]);
2195              sum   = &(sumtable[i * 80 + l * 20]);
2196
2197#ifdef __SIM_SSE3
2198              for(k = 0; k < 20; k+=2)
2199                {
2200                  __m128d sumv = _mm_mul_pd(_mm_load_pd(&left[k]), _mm_load_pd(&right[k]));
2201                 
2202                  _mm_store_pd(&sum[k], sumv);           
2203                }
2204#else
2205              for(k = 0; k < 20; k++)
2206                sum[k] = left[k] * right[k];
2207#endif
2208            }
2209        }
2210      break;
2211    default:
2212      assert(0);
2213    }
2214}
2215
2216static void sumGAMMAPROT_LG4(int tipCase, double *sumtable, double *x1, double *x2, double *tipVector[4],
2217                             unsigned char *tipX1, unsigned char *tipX2, int n)
2218{
2219  int i, l, k;
2220  double *left, *right, *sum;
2221
2222  switch(tipCase)
2223    {
2224    case TIP_TIP:
2225      for(i = 0; i < n; i++)
2226        {         
2227          for(l = 0; l < 4; l++)
2228            {
2229              left  = &(tipVector[l][20 * tipX1[i]]);
2230              right = &(tipVector[l][20 * tipX2[i]]);
2231
2232              sum = &sumtable[i * 80 + l * 20];
2233#ifdef __SIM_SSE3
2234              for(k = 0; k < 20; k+=2)
2235                {
2236                  __m128d sumv = _mm_mul_pd(_mm_load_pd(&left[k]), _mm_load_pd(&right[k]));
2237                 
2238                  _mm_store_pd(&sum[k], sumv);           
2239                }
2240#else
2241              for(k = 0; k < 20; k++)
2242                sum[k] = left[k] * right[k];
2243#endif
2244            }
2245        }
2246      break;
2247    case TIP_INNER:
2248      for(i = 0; i < n; i++)
2249        {
2250         
2251
2252          for(l = 0; l < 4; l++)
2253            { 
2254              left = &(tipVector[l][20 * tipX1[i]]);
2255              right = &(x2[80 * i + l * 20]);
2256              sum = &sumtable[i * 80 + l * 20];
2257#ifdef __SIM_SSE3
2258              for(k = 0; k < 20; k+=2)
2259                {
2260                  __m128d sumv = _mm_mul_pd(_mm_load_pd(&left[k]), _mm_load_pd(&right[k]));
2261                 
2262                  _mm_store_pd(&sum[k], sumv);           
2263                }
2264#else
2265              for(k = 0; k < 20; k++)
2266                sum[k] = left[k] * right[k];
2267#endif
2268            }
2269        }
2270      break;
2271    case INNER_INNER:
2272      for(i = 0; i < n; i++)
2273        {
2274          for(l = 0; l < 4; l++)
2275            {
2276              left  = &(x1[80 * i + l * 20]);
2277              right = &(x2[80 * i + l * 20]);
2278              sum   = &(sumtable[i * 80 + l * 20]);
2279
2280#ifdef __SIM_SSE3
2281              for(k = 0; k < 20; k+=2)
2282                {
2283                  __m128d sumv = _mm_mul_pd(_mm_load_pd(&left[k]), _mm_load_pd(&right[k]));
2284                 
2285                  _mm_store_pd(&sum[k], sumv);           
2286                }
2287#else
2288              for(k = 0; k < 20; k++)
2289                sum[k] = left[k] * right[k];
2290#endif
2291            }
2292        }
2293      break;
2294    default:
2295      assert(0);
2296    }
2297}
2298
2299
2300#ifdef __SIM_SSE3
2301static void sumGAMMAPROT_GAPPED_SAVE(int tipCase, double *sumtable, double *x1, double *x2, double *tipVector,
2302    unsigned char *tipX1, unsigned char *tipX2, int n, 
2303    double *x1_gapColumn, double *x2_gapColumn, unsigned int *x1_gap, unsigned int *x2_gap)
2304{
2305  int i, l, k;
2306  double 
2307    *left, 
2308    *right, 
2309    *sum,
2310    *x1_ptr = x1,
2311    *x2_ptr = x2,
2312    *x1v,
2313    *x2v;
2314
2315  switch(tipCase)
2316  {
2317    case TIP_TIP:
2318      for(i = 0; i < n; i++)
2319      {
2320        left  = &(tipVector[20 * tipX1[i]]);
2321        right = &(tipVector[20 * tipX2[i]]);
2322
2323        for(l = 0; l < 4; l++)
2324        {
2325          sum = &sumtable[i * 80 + l * 20];
2326
2327          for(k = 0; k < 20; k+=2)
2328          {
2329            __m128d sumv = _mm_mul_pd(_mm_load_pd(&left[k]), _mm_load_pd(&right[k]));
2330
2331            _mm_store_pd(&sum[k], sumv);                 
2332          }
2333
2334        }
2335      }
2336      break;
2337    case TIP_INNER:
2338      for(i = 0; i < n; i++)
2339      {
2340        left = &(tipVector[20 * tipX1[i]]);
2341
2342        if(x2_gap[i / 32] & mask32[i % 32])
2343          x2v = x2_gapColumn;
2344        else
2345        {
2346          x2v = x2_ptr;
2347          x2_ptr += 80;
2348        }
2349
2350        for(l = 0; l < 4; l++)
2351        {
2352          right = &(x2v[l * 20]);
2353          sum = &sumtable[i * 80 + l * 20];
2354
2355          for(k = 0; k < 20; k+=2)
2356          {
2357            __m128d sumv = _mm_mul_pd(_mm_load_pd(&left[k]), _mm_load_pd(&right[k]));
2358
2359            _mm_store_pd(&sum[k], sumv);                 
2360          }
2361        }
2362      }
2363      break;
2364    case INNER_INNER:
2365      for(i = 0; i < n; i++)
2366      {
2367        if(x1_gap[i / 32] & mask32[i % 32])
2368          x1v = x1_gapColumn;
2369        else
2370        {
2371          x1v  = x1_ptr;
2372          x1_ptr += 80;
2373        }
2374
2375        if(x2_gap[i / 32] & mask32[i % 32])
2376          x2v = x2_gapColumn;
2377        else
2378        {
2379          x2v  = x2_ptr;
2380          x2_ptr += 80;
2381        }
2382
2383        for(l = 0; l < 4; l++)
2384        {
2385          left  = &(x1v[l * 20]);
2386          right = &(x2v[l * 20]);
2387          sum   = &(sumtable[i * 80 + l * 20]);
2388
2389          for(k = 0; k < 20; k+=2)
2390          {
2391            __m128d sumv = _mm_mul_pd(_mm_load_pd(&left[k]), _mm_load_pd(&right[k]));
2392
2393            _mm_store_pd(&sum[k], sumv);                 
2394          }
2395        }
2396      }
2397      break;
2398    default:
2399      assert(0);
2400  }
2401}
2402
2403#endif
2404
2405
2406static void sumGammaFlex(int tipCase, double *sumtable, double *x1, double *x2, double *tipVector,
2407                         unsigned char *tipX1, unsigned char *tipX2, int n, const int numStates)
2408{
2409  int i, l, k;
2410  double *left, *right, *sum;
2411
2412  const int gammaStates = numStates * 4;
2413
2414  switch(tipCase)
2415    {
2416    case TIP_TIP:
2417      for(i = 0; i < n; i++)
2418        {
2419          left  = &(tipVector[numStates * tipX1[i]]);
2420          right = &(tipVector[numStates * tipX2[i]]);
2421
2422          for(l = 0; l < 4; l++)
2423            {
2424              sum = &sumtable[i * gammaStates + l * numStates];
2425              for(k = 0; k < numStates; k++)
2426                sum[k] = left[k] * right[k];
2427            }
2428        }
2429      break;
2430    case TIP_INNER:
2431      for(i = 0; i < n; i++)
2432        {
2433          left = &(tipVector[numStates * tipX1[i]]);
2434
2435          for(l = 0; l < 4; l++)
2436            {
2437              right = &(x2[gammaStates * i + l * numStates]);
2438              sum = &sumtable[i * gammaStates + l * numStates];
2439
2440              for(k = 0; k < numStates; k++)
2441                sum[k] = left[k] * right[k];
2442            }
2443        }
2444      break;
2445    case INNER_INNER:
2446      for(i = 0; i < n; i++)
2447        {
2448          for(l = 0; l < 4; l++)
2449            {
2450              left  = &(x1[gammaStates * i + l * numStates]);
2451              right = &(x2[gammaStates * i + l * numStates]);
2452              sum   = &(sumtable[i * gammaStates + l * numStates]);
2453
2454              for(k = 0; k < numStates; k++)
2455                sum[k] = left[k] * right[k];
2456            }
2457        }
2458      break;
2459    default:
2460      assert(0);
2461    }
2462}
2463
2464static void sumGAMMASECONDARY(int tipCase, double *sumtable, double *x1, double *x2, double *tipVector,
2465                              unsigned char *tipX1, unsigned char *tipX2, int n)
2466{
2467  int i, l, k;
2468  double *left, *right, *sum;
2469
2470  switch(tipCase)
2471    {
2472    case TIP_TIP:
2473      for(i = 0; i < n; i++)
2474        {
2475          left  = &(tipVector[16 * tipX1[i]]);
2476          right = &(tipVector[16 * tipX2[i]]);
2477
2478          for(l = 0; l < 4; l++)
2479            {
2480              sum = &sumtable[i * 64 + l * 16];
2481              for(k = 0; k < 16; k++)
2482                sum[k] = left[k] * right[k];
2483            }
2484        }
2485      break;
2486    case TIP_INNER:
2487      for(i = 0; i < n; i++)
2488        {
2489          left = &(tipVector[16 * tipX1[i]]);
2490
2491          for(l = 0; l < 4; l++)
2492            {
2493              right = &(x2[64 * i + l * 16]);
2494              sum = &sumtable[i * 64 + l * 16];
2495
2496              for(k = 0; k < 16; k++)
2497                sum[k] = left[k] * right[k];
2498            }
2499        }
2500      break;
2501    case INNER_INNER:
2502      for(i = 0; i < n; i++)
2503        {
2504          for(l = 0; l < 4; l++)
2505            {
2506              left  = &(x1[64 * i + l * 16]);
2507              right = &(x2[64 * i + l * 16]);
2508              sum   = &(sumtable[i * 64 + l * 16]);
2509
2510              for(k = 0; k < 16; k++)
2511                sum[k] = left[k] * right[k];
2512            }
2513        }
2514      break;
2515    default:
2516      assert(0);
2517    }
2518}
2519
2520static void sumGAMMASECONDARY_6(int tipCase, double *sumtable, double *x1, double *x2, double *tipVector,
2521                                unsigned char *tipX1, unsigned char *tipX2, int n)
2522{
2523  int i, l, k;
2524  double *left, *right, *sum;
2525
2526  switch(tipCase)
2527    {
2528    case TIP_TIP:
2529      for(i = 0; i < n; i++)
2530        {
2531          left  = &(tipVector[6 * tipX1[i]]);
2532          right = &(tipVector[6 * tipX2[i]]);
2533
2534          for(l = 0; l < 4; l++)
2535            {
2536              sum = &sumtable[i * 24 + l * 6];
2537              for(k = 0; k < 6; k++)
2538                sum[k] = left[k] * right[k];
2539            }
2540        }
2541      break;
2542    case TIP_INNER:
2543      for(i = 0; i < n; i++)
2544        {
2545          left = &(tipVector[6 * tipX1[i]]);
2546
2547          for(l = 0; l < 4; l++)
2548            {
2549              right = &(x2[24 * i + l * 6]);
2550              sum = &sumtable[i * 24 + l * 6];
2551
2552              for(k = 0; k < 6; k++)
2553                sum[k] = left[k] * right[k];
2554            }
2555        }
2556      break;
2557    case INNER_INNER:
2558      for(i = 0; i < n; i++)
2559        {
2560          for(l = 0; l < 4; l++)
2561            {
2562              left  = &(x1[24 * i + l * 6]);
2563              right = &(x2[24 * i + l * 6]);
2564              sum   = &(sumtable[i * 24 + l * 6]);
2565
2566              for(k = 0; k < 6; k++)
2567                sum[k] = left[k] * right[k];
2568            }
2569        }
2570      break;
2571    default:
2572      assert(0);
2573    }
2574}
2575
2576static void sumGAMMASECONDARY_7(int tipCase, double *sumtable, double *x1, double *x2, double *tipVector,
2577                                unsigned char *tipX1, unsigned char *tipX2, int n)
2578{
2579  int i, l, k;
2580  double *left, *right, *sum;
2581
2582  switch(tipCase)
2583    {
2584    case TIP_TIP:
2585      for(i = 0; i < n; i++)
2586        {
2587          left  = &(tipVector[7 * tipX1[i]]);
2588          right = &(tipVector[7 * tipX2[i]]);
2589
2590          for(l = 0; l < 4; l++)
2591            {
2592              sum = &sumtable[i * 28 + l * 7];
2593              for(k = 0; k < 7; k++)
2594                sum[k] = left[k] * right[k];
2595            }
2596        }
2597      break;
2598    case TIP_INNER:
2599      for(i = 0; i < n; i++)
2600        {
2601          left = &(tipVector[7 * tipX1[i]]);
2602
2603          for(l = 0; l < 4; l++)
2604            {
2605              right = &(x2[28 * i + l * 7]);
2606              sum = &sumtable[i * 28 + l * 7];
2607
2608              for(k = 0; k < 7; k++)
2609                sum[k] = left[k] * right[k];
2610            }
2611        }
2612      break;
2613    case INNER_INNER:
2614      for(i = 0; i < n; i++)
2615        {
2616          for(l = 0; l < 4; l++)
2617            {
2618              left  = &(x1[28 * i + l * 7]);
2619              right = &(x2[28 * i + l * 7]);
2620              sum   = &(sumtable[i * 28 + l * 7]);
2621
2622              for(k = 0; k < 7; k++)
2623                sum[k] = left[k] * right[k];
2624            }
2625        }
2626      break;
2627    default:
2628      assert(0);
2629    }
2630}
2631
2632#ifdef __SIM_SSE3
2633
2634static void coreGTRGAMMAPROT(double *gammaRates, double *EIGN, double *sumtable, int upper, int *wrptr,
2635                              volatile double *ext_dlnLdlz,  volatile double *ext_d2lnLdlz2, double lz)
2636{
2637  double  *sum, 
2638    diagptable0[80] __attribute__ ((aligned (BYTE_ALIGNMENT))),
2639    diagptable1[80] __attribute__ ((aligned (BYTE_ALIGNMENT))),
2640    diagptable2[80] __attribute__ ((aligned (BYTE_ALIGNMENT)));   
2641  int     i, j, l;
2642  double  dlnLdlz = 0;
2643  double d2lnLdlz2 = 0;
2644  double ki, kisqr; 
2645  double inv_Li, dlnLidlz, d2lnLidlz2;
2646
2647  for(i = 0; i < 4; i++)
2648    {
2649      ki = gammaRates[i];
2650      kisqr = ki * ki;
2651     
2652      diagptable0[i * 20] = 1.0;
2653      diagptable1[i * 20] = 0.0;
2654      diagptable2[i * 20] = 0.0;
2655
2656      for(l = 1; l < 20; l++)
2657        {
2658          diagptable0[i * 20 + l] = EXP(EIGN[l-1] * ki * lz);
2659          diagptable1[i * 20 + l] = EIGN[l-1] * ki;
2660          diagptable2[i * 20 + l] = EIGN[l-1] * EIGN[l-1] * kisqr;
2661        }
2662    }
2663
2664  for (i = 0; i < upper; i++)
2665    { 
2666      __m128d a0 = _mm_setzero_pd();
2667      __m128d a1 = _mm_setzero_pd();
2668      __m128d a2 = _mm_setzero_pd();
2669
2670      sum = &sumtable[i * 80];         
2671
2672      for(j = 0; j < 4; j++)
2673        {                       
2674          double           
2675            *d0 = &diagptable0[j * 20],
2676            *d1 = &diagptable1[j * 20],
2677            *d2 = &diagptable2[j * 20];
2678                 
2679          for(l = 0; l < 20; l+=2)
2680            {
2681              __m128d tmpv = _mm_mul_pd(_mm_load_pd(&d0[l]), _mm_load_pd(&sum[j * 20 +l]));
2682              a0 = _mm_add_pd(a0, tmpv);
2683              a1 = _mm_add_pd(a1, _mm_mul_pd(tmpv, _mm_load_pd(&d1[l])));
2684              a2 = _mm_add_pd(a2, _mm_mul_pd(tmpv, _mm_load_pd(&d2[l])));
2685            }             
2686        }
2687
2688      a0 = _mm_hadd_pd(a0, a0);
2689      a1 = _mm_hadd_pd(a1, a1);
2690      a2 = _mm_hadd_pd(a2, a2);
2691
2692      _mm_storel_pd(&inv_Li, a0);
2693      _mm_storel_pd(&dlnLidlz, a1);
2694      _mm_storel_pd(&d2lnLidlz2, a2);
2695
2696      inv_Li = 1.0 / inv_Li;
2697
2698      dlnLidlz   *= inv_Li;
2699      d2lnLidlz2 *= inv_Li;
2700
2701      dlnLdlz   += wrptr[i] * dlnLidlz;
2702      d2lnLdlz2 += wrptr[i] * (d2lnLidlz2 - dlnLidlz * dlnLidlz);
2703    }
2704
2705  *ext_dlnLdlz   = dlnLdlz;
2706  *ext_d2lnLdlz2 = d2lnLdlz2;
2707}
2708
2709static void coreGTRGAMMAPROT_LG4(double *gammaRates, double *EIGN[4], double *sumtable, int upper, int *wrptr,
2710                                 volatile double *ext_dlnLdlz,  volatile double *ext_d2lnLdlz2, double lz, double *weights)
2711{
2712  double  *sum, 
2713    diagptable0[80] __attribute__ ((aligned (BYTE_ALIGNMENT))),
2714    diagptable1[80] __attribute__ ((aligned (BYTE_ALIGNMENT))),
2715    diagptable2[80] __attribute__ ((aligned (BYTE_ALIGNMENT)));   
2716  int     i, j, l;
2717  double  dlnLdlz = 0;
2718  double d2lnLdlz2 = 0;
2719  double ki, kisqr;   
2720
2721  for(i = 0; i < 4; i++)
2722    {
2723      ki = gammaRates[i];
2724      kisqr = ki * ki;
2725     
2726      diagptable0[i * 20] = 1.0;
2727      diagptable1[i * 20] = 0.0;
2728      diagptable2[i * 20] = 0.0;
2729
2730      for(l = 1; l < 20; l++)
2731        {
2732          diagptable0[i * 20 + l] = EXP(EIGN[i][l-1] * ki * lz);
2733          diagptable1[i * 20 + l] = EIGN[i][l-1] * ki;
2734          diagptable2[i * 20 + l] = EIGN[i][l-1] * EIGN[i][l-1] * kisqr;
2735        }
2736    }
2737
2738  for (i = 0; i < upper; i++)
2739    {       
2740      double   
2741        inv_Li = 0.0, 
2742        dlnLidlz = 0.0, 
2743        d2lnLidlz2 = 0.0;
2744
2745      sum = &sumtable[i * 80];         
2746
2747      for(j = 0; j < 4; j++)
2748        {                       
2749          double
2750            l0,
2751            l1,
2752            l2,           
2753            *d0 = &diagptable0[j * 20],
2754            *d1 = &diagptable1[j * 20],
2755            *d2 = &diagptable2[j * 20];
2756         
2757          __m128d 
2758            a0 = _mm_setzero_pd(),
2759            a1 = _mm_setzero_pd(),
2760            a2 = _mm_setzero_pd();
2761                 
2762          for(l = 0; l < 20; l+=2)
2763            {
2764              __m128d 
2765                tmpv = _mm_mul_pd(_mm_load_pd(&d0[l]), _mm_load_pd(&sum[j * 20 +l]));
2766             
2767              a0 = _mm_add_pd(a0, tmpv);
2768              a1 = _mm_add_pd(a1, _mm_mul_pd(tmpv, _mm_load_pd(&d1[l])));
2769              a2 = _mm_add_pd(a2, _mm_mul_pd(tmpv, _mm_load_pd(&d2[l])));
2770            }
2771         
2772          a0 = _mm_hadd_pd(a0, a0);
2773          a1 = _mm_hadd_pd(a1, a1);
2774          a2 = _mm_hadd_pd(a2, a2);
2775
2776          _mm_storel_pd(&l0, a0);
2777          _mm_storel_pd(&l1, a1);
2778          _mm_storel_pd(&l2, a2);
2779         
2780          inv_Li     += weights[j] * l0;
2781          dlnLidlz   += weights[j] * l1;
2782          d2lnLidlz2 += weights[j] * l2; 
2783        }
2784
2785      inv_Li = 1.0 / inv_Li;
2786
2787      dlnLidlz   *= inv_Li;
2788      d2lnLidlz2 *= inv_Li;
2789
2790      dlnLdlz   += wrptr[i] * dlnLidlz;
2791      d2lnLdlz2 += wrptr[i] * (d2lnLidlz2 - dlnLidlz * dlnLidlz);
2792    }
2793
2794  *ext_dlnLdlz   = dlnLdlz;
2795  *ext_d2lnLdlz2 = d2lnLdlz2;
2796}
2797
2798
2799#else
2800
2801static void coreGTRGAMMAPROT(double *gammaRates, double *EIGN, double *sumtable, int upper, int *wrptr,
2802                              volatile double *ext_dlnLdlz,  volatile double *ext_d2lnLdlz2, double lz)
2803{
2804  double  *sum, diagptable[512];
2805  int     i, j, l;
2806  double  dlnLdlz = 0;
2807  double d2lnLdlz2 = 0;
2808  double ki, kisqr;
2809  double tmp;
2810  double inv_Li, dlnLidlz, d2lnLidlz2;
2811
2812  for(i = 0; i < 4; i++)
2813    {
2814      ki = gammaRates[i];
2815      kisqr = ki * ki;
2816
2817      for(l = 1; l < 20; l++)
2818        {
2819          diagptable[i * 128 + l * 4]     = EXP(EIGN[l-1] * ki * lz);
2820          diagptable[i * 128 + l * 4 + 1] = EIGN[l-1] * ki;
2821          diagptable[i * 128 + l * 4 + 2] = EIGN[l-1] * EIGN[l-1] * kisqr;
2822        }
2823    }
2824
2825  for (i = 0; i < upper; i++)
2826    {
2827      sum = &sumtable[i * 80];
2828      inv_Li   = 0.0;
2829      dlnLidlz = 0.0;
2830      d2lnLidlz2 = 0.0;
2831
2832      for(j = 0; j < 4; j++)
2833        {
2834          inv_Li += sum[j * 20];
2835
2836          for(l = 1; l < 20; l++)
2837            {
2838              inv_Li     += (tmp = diagptable[j * 128 + l * 4] * sum[j * 20 + l]);
2839              dlnLidlz   +=  tmp * diagptable[j * 128 + l * 4 + 1];
2840              d2lnLidlz2 +=  tmp * diagptable[j * 128 + l * 4 + 2];
2841            }
2842        }
2843
2844      inv_Li = 1.0 / inv_Li;
2845
2846      dlnLidlz   *= inv_Li;
2847      d2lnLidlz2 *= inv_Li;
2848
2849      dlnLdlz   += wrptr[i] * dlnLidlz;
2850      d2lnLdlz2 += wrptr[i] * (d2lnLidlz2 - dlnLidlz * dlnLidlz);
2851    }
2852
2853  *ext_dlnLdlz   = dlnLdlz;
2854  *ext_d2lnLdlz2 = d2lnLdlz2;
2855}
2856
2857static void coreGTRGAMMAPROT_LG4(double *gammaRates, double *EIGN[4], double *sumtable, int upper, int *wrptr,
2858                                 volatile double *ext_dlnLdlz,  volatile double *ext_d2lnLdlz2, double lz, double *weights)
2859{
2860  double  *sum, diagptable[512];
2861  int     i, j, l;
2862  double  dlnLdlz = 0;
2863  double d2lnLdlz2 = 0;
2864  double ki, kisqr;
2865  double tmp;
2866  double inv_Li, dlnLidlz, d2lnLidlz2;
2867
2868  for(i = 0; i < 4; i++)
2869    {
2870      ki = gammaRates[i];
2871      kisqr = ki * ki;
2872
2873      for(l = 1; l < 20; l++)
2874        {
2875          diagptable[i * 128 + l * 4]     = EXP(EIGN[i][l-1] * ki * lz);
2876          diagptable[i * 128 + l * 4 + 1] = EIGN[i][l-1] * ki;
2877          diagptable[i * 128 + l * 4 + 2] = EIGN[i][l-1] * EIGN[i][l-1] * kisqr;
2878        }
2879    }
2880
2881  for (i = 0; i < upper; i++)
2882    {
2883      sum = &sumtable[i * 80];
2884      inv_Li   = 0.0;
2885      dlnLidlz = 0.0;
2886      d2lnLidlz2 = 0.0;
2887
2888      for(j = 0; j < 4; j++)
2889        {
2890          double 
2891            a1 = 0.0,
2892            a2 = 0.0,
2893            a3 = 0.0;
2894
2895          a1 += sum[j * 20];
2896
2897          for(l = 1; l < 20; l++)
2898            {
2899              a1 += (tmp = diagptable[j * 128 + l * 4] * sum[j * 20 + l]);
2900              a2 +=  tmp * diagptable[j * 128 + l * 4 + 1];
2901              a3 +=  tmp * diagptable[j * 128 + l * 4 + 2];
2902            }
2903
2904          inv_Li     += weights[j] * a1;
2905          dlnLidlz   += weights[j] * a2;
2906          d2lnLidlz2 += weights[j] * a3;
2907        }
2908
2909      inv_Li = 1.0 / inv_Li;
2910
2911      dlnLidlz   *= inv_Li;
2912      d2lnLidlz2 *= inv_Li;
2913
2914      dlnLdlz   += wrptr[i] * dlnLidlz;
2915      d2lnLdlz2 += wrptr[i] * (d2lnLidlz2 - dlnLidlz * dlnLidlz);
2916    }
2917
2918  *ext_dlnLdlz   = dlnLdlz;
2919  *ext_d2lnLdlz2 = d2lnLdlz2;
2920}
2921
2922
2923#endif
2924
2925
2926
2927
2928static void coreGTRGAMMASECONDARY(double *gammaRates, double *EIGN, double *sumtable, int upper, int *wrptr,
2929                                  volatile double *ext_dlnLdlz,  volatile double *ext_d2lnLdlz2, double lz)
2930{
2931  double  *sum, diagptable[256];
2932  int     i, j, l;
2933  double  dlnLdlz = 0;
2934  double d2lnLdlz2 = 0;
2935  double ki, kisqr;
2936  double tmp;
2937  double inv_Li, dlnLidlz, d2lnLidlz2;
2938
2939  for(i = 0; i < 4; i++)
2940    {
2941      ki = gammaRates[i];
2942      kisqr = ki * ki;
2943
2944      for(l = 1; l < 16; l++)
2945        {
2946          diagptable[i * 64 + l * 4]     = EXP(EIGN[l-1] * ki * lz);
2947          diagptable[i * 64 + l * 4 + 1] = EIGN[l-1] * ki;
2948          diagptable[i * 64 + l * 4 + 2] = EIGN[l-1] * EIGN[l-1] * kisqr;
2949        }
2950    }
2951
2952  for (i = 0; i < upper; i++)
2953    {
2954      sum = &sumtable[i * 64];
2955      inv_Li   = 0.0;
2956      dlnLidlz = 0.0;
2957      d2lnLidlz2 = 0.0;
2958
2959      for(j = 0; j < 4; j++)
2960        {
2961          inv_Li += sum[j * 16];
2962
2963          for(l = 1; l < 16; l++)
2964            {
2965              inv_Li     += (tmp = diagptable[j * 64 + l * 4] * sum[j * 16 + l]);
2966              dlnLidlz   +=  tmp * diagptable[j * 64 + l * 4 + 1];
2967              d2lnLidlz2 +=  tmp * diagptable[j * 64 + l * 4 + 2];
2968            }
2969        }
2970
2971      inv_Li = 1.0 / inv_Li;
2972
2973      dlnLidlz   *= inv_Li;
2974      d2lnLidlz2 *= inv_Li;
2975
2976      dlnLdlz   += wrptr[i] * dlnLidlz;
2977      d2lnLdlz2 += wrptr[i] * (d2lnLidlz2 - dlnLidlz * dlnLidlz);
2978    }
2979
2980  *ext_dlnLdlz   = dlnLdlz;
2981  *ext_d2lnLdlz2 = d2lnLdlz2;
2982}
2983
2984static void coreGTRGAMMASECONDARY_6(double *gammaRates, double *EIGN, double *sumtable, int upper, int *wrptr,
2985                                    volatile double *ext_dlnLdlz,  volatile double *ext_d2lnLdlz2, double lz)
2986{
2987  double  *sum, diagptable[96];
2988  int     i, j, l;
2989  double  dlnLdlz = 0;
2990  double d2lnLdlz2 = 0;
2991  double ki, kisqr;
2992  double tmp;
2993  double inv_Li, dlnLidlz, d2lnLidlz2;
2994
2995  for(i = 0; i < 4; i++)
2996    {
2997      ki = gammaRates[i];
2998      kisqr = ki * ki;
2999
3000      for(l = 1; l < 6; l++)
3001        {
3002          diagptable[i * 24 + l * 4]     = EXP(EIGN[l-1] * ki * lz);
3003          diagptable[i * 24 + l * 4 + 1] = EIGN[l-1] * ki;
3004          diagptable[i * 24 + l * 4 + 2] = EIGN[l-1] * EIGN[l-1] * kisqr;
3005        }
3006    }
3007
3008  for (i = 0; i < upper; i++)
3009    {
3010      sum = &sumtable[i * 24];
3011      inv_Li   = 0.0;
3012      dlnLidlz = 0.0;
3013      d2lnLidlz2 = 0.0;
3014
3015      for(j = 0; j < 4; j++)
3016        {
3017          inv_Li += sum[j * 6];
3018
3019          for(l = 1; l < 6; l++)
3020            {
3021              inv_Li     += (tmp = diagptable[j * 24 + l * 4] * sum[j * 6 + l]);
3022              dlnLidlz   +=  tmp * diagptable[j * 24 + l * 4 + 1];
3023              d2lnLidlz2 +=  tmp * diagptable[j * 24 + l * 4 + 2];
3024            }
3025        }
3026
3027      inv_Li = 1.0 / inv_Li;
3028
3029      dlnLidlz   *= inv_Li;
3030      d2lnLidlz2 *= inv_Li;
3031
3032      dlnLdlz   += wrptr[i] * dlnLidlz;
3033      d2lnLdlz2 += wrptr[i] * (d2lnLidlz2 - dlnLidlz * dlnLidlz);
3034    }
3035
3036  *ext_dlnLdlz   = dlnLdlz;
3037  *ext_d2lnLdlz2 = d2lnLdlz2;
3038}
3039
3040static void coreGTRGAMMASECONDARY_7(double *gammaRates, double *EIGN, double *sumtable, int upper, int *wrptr,
3041                                    volatile double *ext_dlnLdlz,  volatile double *ext_d2lnLdlz2, double lz)
3042{
3043  double  *sum, diagptable[112];
3044  int     i, j, l;
3045  double  dlnLdlz = 0;
3046  double d2lnLdlz2 = 0;
3047  double ki, kisqr;
3048  double tmp;
3049  double inv_Li, dlnLidlz, d2lnLidlz2;
3050
3051  for(i = 0; i < 4; i++)
3052    {
3053      ki = gammaRates[i];
3054      kisqr = ki * ki;
3055
3056      for(l = 1; l < 7; l++)
3057        {
3058          diagptable[i * 28 + l * 4]     = EXP(EIGN[l-1] * ki * lz);
3059          diagptable[i * 28 + l * 4 + 1] = EIGN[l-1] * ki;
3060          diagptable[i * 28 + l * 4 + 2] = EIGN[l-1] * EIGN[l-1] * kisqr;
3061        }
3062    }
3063
3064  for (i = 0; i < upper; i++)
3065    {
3066      sum = &sumtable[i * 28];
3067      inv_Li   = 0.0;
3068      dlnLidlz = 0.0;
3069      d2lnLidlz2 = 0.0;
3070
3071      for(j = 0; j < 4; j++)
3072        {
3073          inv_Li += sum[j * 7];
3074
3075          for(l = 1; l < 7; l++)
3076            {
3077              inv_Li     += (tmp = diagptable[j * 28 + l * 4] * sum[j * 7 + l]);
3078              dlnLidlz   +=  tmp * diagptable[j * 28 + l * 4 + 1];
3079              d2lnLidlz2 +=  tmp * diagptable[j * 28 + l * 4 + 2];
3080            }
3081        }
3082
3083      inv_Li = 1.0 / inv_Li;
3084
3085      dlnLidlz   *= inv_Li;
3086      d2lnLidlz2 *= inv_Li;
3087
3088      dlnLdlz   += wrptr[i] * dlnLidlz;
3089      d2lnLdlz2 += wrptr[i] * (d2lnLidlz2 - dlnLidlz * dlnLidlz);
3090    }
3091
3092  *ext_dlnLdlz   = dlnLdlz;
3093  *ext_d2lnLdlz2 = d2lnLdlz2;
3094}
3095
3096static void coreGTRGAMMAPROTINVAR(double *gammaRates, double *EIGN, double *sumtable, int upper, int *wrptr,
3097                                   volatile double *ext_dlnLdlz,  volatile double *ext_d2lnLdlz2, double lz, double *frequencies,
3098                                  double propInvar, int *iptr)
3099{
3100  double  *sum, diagptable[512];
3101  int     i, l, j;
3102  double  dlnLdlz = 0;
3103  double d2lnLdlz2 = 0;
3104  double ki, kisqr;
3105  double freqs[20];
3106  double scaler =  0.25 * (1.0 - propInvar);
3107  double tmp;
3108  double inv_Li, dlnLidlz, d2lnLidlz2;
3109
3110  for(i = 0; i < 20; i++)
3111    freqs[i] = frequencies[i] * propInvar;
3112
3113  for(i = 0; i < 4; i++)
3114    {
3115      ki = gammaRates[i];
3116      kisqr = ki * ki;
3117
3118      for(l = 1; l < 20; l++)
3119        {
3120          diagptable[i * 128 + l * 4]     = EXP(EIGN[l-1] * ki * lz);
3121          diagptable[i * 128 + l * 4 + 1] = EIGN[l-1] * ki;
3122          diagptable[i * 128 + l * 4 + 2] = EIGN[l-1] * EIGN[l-1] * kisqr;
3123        }
3124    }
3125
3126
3127   for(i = 0; i < upper; i++)
3128     {
3129       sum = &sumtable[i * 80];
3130       inv_Li   = 0.0;
3131       dlnLidlz = 0.0;
3132       d2lnLidlz2 = 0.0;
3133
3134       for(j = 0; j < 4; j++)
3135        {
3136          inv_Li += sum[j * 20];
3137
3138          for(l = 1; l < 20; l++)
3139            {
3140              inv_Li     += (tmp = diagptable[j * 128 + l * 4] * sum[j * 20 + l]);
3141              dlnLidlz   +=  tmp * diagptable[j * 128 + l * 4 + 1];
3142              d2lnLidlz2 +=  tmp * diagptable[j * 128 + l * 4 + 2];
3143            }
3144        }
3145
3146       inv_Li *= scaler;
3147
3148       if(iptr[i] < 20)
3149         inv_Li += freqs[iptr[i]];
3150
3151       inv_Li = 1.0 / inv_Li;
3152
3153       dlnLidlz   *= inv_Li;
3154       d2lnLidlz2 *= inv_Li;
3155
3156       dlnLidlz *= scaler;
3157       d2lnLidlz2 *= scaler;
3158
3159       dlnLdlz  += wrptr[i] * dlnLidlz;
3160       d2lnLdlz2 += wrptr[i] * (d2lnLidlz2 - dlnLidlz * dlnLidlz);
3161     }
3162
3163   *ext_dlnLdlz   = dlnLdlz;
3164   *ext_d2lnLdlz2 = d2lnLdlz2;
3165}
3166
3167
3168static void coreGTRGAMMASECONDARYINVAR(double *gammaRates, double *EIGN, double *sumtable, int upper, int *wrptr,
3169                                       volatile double *ext_dlnLdlz,  volatile double *ext_d2lnLdlz2, double lz, double *frequencies,
3170                                       double propInvar, int *iptr)
3171{
3172  double  *sum, diagptable[256];
3173  int     i, l, j;
3174  double  dlnLdlz = 0;
3175  double d2lnLdlz2 = 0;
3176  double ki, kisqr;
3177  double freqs[16];
3178  double scaler =  0.25 * (1.0 - propInvar);
3179  double tmp;
3180  double inv_Li, dlnLidlz, d2lnLidlz2;
3181
3182  for(i = 0; i < 16; i++)
3183    freqs[i] = frequencies[i] * propInvar;
3184
3185  for(i = 0; i < 4; i++)
3186    {
3187      ki = gammaRates[i];
3188      kisqr = ki * ki;
3189
3190      for(l = 1; l < 16; l++)
3191        {
3192          diagptable[i * 64 + l * 4]     = EXP(EIGN[l-1] * ki * lz);
3193          diagptable[i * 64 + l * 4 + 1] = EIGN[l-1] * ki;
3194          diagptable[i * 64 + l * 4 + 2] = EIGN[l-1] * EIGN[l-1] * kisqr;
3195        }
3196    }
3197
3198
3199   for(i = 0; i < upper; i++)
3200     {
3201       sum = &sumtable[i * 64];
3202       inv_Li   = 0.0;
3203       dlnLidlz = 0.0;
3204       d2lnLidlz2 = 0.0;
3205
3206       for(j = 0; j < 4; j++)
3207        {
3208          inv_Li += sum[j * 16];
3209
3210          for(l = 1; l < 16; l++)
3211            {
3212              inv_Li     += (tmp = diagptable[j * 64 + l * 4] * sum[j * 16 + l]);
3213              dlnLidlz   +=  tmp * diagptable[j * 64 + l * 4 + 1];
3214              d2lnLidlz2 +=  tmp * diagptable[j * 64 + l * 4 + 2];
3215            }
3216        }
3217
3218       inv_Li *= scaler;
3219
3220       if(iptr[i] < 16)
3221         inv_Li += freqs[iptr[i]];
3222
3223       inv_Li = 1.0 / inv_Li;
3224
3225       dlnLidlz   *= inv_Li;
3226       d2lnLidlz2 *= inv_Li;
3227
3228       dlnLidlz *= scaler;
3229       d2lnLidlz2 *= scaler;
3230
3231       dlnLdlz  += wrptr[i] * dlnLidlz;
3232       d2lnLdlz2 += wrptr[i] * (d2lnLidlz2 - dlnLidlz * dlnLidlz);
3233     }
3234
3235   *ext_dlnLdlz   = dlnLdlz;
3236   *ext_d2lnLdlz2 = d2lnLdlz2;
3237}
3238
3239static void coreGTRGAMMASECONDARYINVAR_6(double *gammaRates, double *EIGN, double *sumtable, int upper, int *wrptr,
3240                                         volatile double *ext_dlnLdlz,  volatile double *ext_d2lnLdlz2, double lz, double *frequencies,
3241                                         double propInvar, int *iptr)
3242{
3243  double  *sum, diagptable[96];
3244  int     i, l, j;
3245  double  dlnLdlz = 0;
3246  double d2lnLdlz2 = 0;
3247  double ki, kisqr;
3248  double freqs[6];
3249  double scaler =  0.25 * (1.0 - propInvar);
3250  double tmp;
3251  double inv_Li, dlnLidlz, d2lnLidlz2;
3252
3253  for(i = 0; i < 6; i++)
3254    freqs[i] = frequencies[i] * propInvar;
3255
3256  for(i = 0; i < 4; i++)
3257    {
3258      ki = gammaRates[i];
3259      kisqr = ki * ki;
3260
3261      for(l = 1; l < 6; l++)
3262        {
3263          diagptable[i * 24 + l * 4]     = EXP(EIGN[l-1] * ki * lz);
3264          diagptable[i * 24 + l * 4 + 1] = EIGN[l-1] * ki;
3265          diagptable[i * 24 + l * 4 + 2] = EIGN[l-1] * EIGN[l-1] * kisqr;
3266        }
3267    }
3268
3269
3270   for(i = 0; i < upper; i++)
3271     {
3272       sum = &sumtable[i * 24];
3273       inv_Li   = 0.0;
3274       dlnLidlz = 0.0;
3275       d2lnLidlz2 = 0.0;
3276
3277       for(j = 0; j < 4; j++)
3278        {
3279          inv_Li += sum[j * 6];
3280
3281          for(l = 1; l < 6; l++)
3282            {
3283              inv_Li     += (tmp = diagptable[j * 24 + l * 4] * sum[j * 6 + l]);
3284              dlnLidlz   +=  tmp * diagptable[j * 24 + l * 4 + 1];
3285              d2lnLidlz2 +=  tmp * diagptable[j * 24 + l * 4 + 2];
3286            }
3287        }
3288
3289       inv_Li *= scaler;
3290
3291       if(iptr[i] < 6)
3292         inv_Li += freqs[iptr[i]];
3293
3294       inv_Li = 1.0 / inv_Li;
3295
3296       dlnLidlz   *= inv_Li;
3297       d2lnLidlz2 *= inv_Li;
3298
3299       dlnLidlz *= scaler;
3300       d2lnLidlz2 *= scaler;
3301
3302       dlnLdlz  += wrptr[i] * dlnLidlz;
3303       d2lnLdlz2 += wrptr[i] * (d2lnLidlz2 - dlnLidlz * dlnLidlz);
3304     }
3305
3306   *ext_dlnLdlz   = dlnLdlz;
3307   *ext_d2lnLdlz2 = d2lnLdlz2;
3308}
3309
3310static void coreGTRGAMMASECONDARYINVAR_7(double *gammaRates, double *EIGN, double *sumtable, int upper, int *wrptr,
3311                                         volatile double *ext_dlnLdlz,  volatile double *ext_d2lnLdlz2, double lz, double *frequencies,
3312                                         double propInvar, int *iptr)
3313{
3314  double  *sum, diagptable[112];
3315  int     i, l, j;
3316  double  dlnLdlz = 0;
3317  double d2lnLdlz2 = 0;
3318  double ki, kisqr;
3319  double freqs[7];
3320  double scaler =  0.25 * (1.0 - propInvar);
3321  double tmp;
3322  double inv_Li, dlnLidlz, d2lnLidlz2;
3323
3324  for(i = 0; i < 7; i++)
3325    freqs[i] = frequencies[i] * propInvar;
3326
3327  for(i = 0; i < 4; i++)
3328    {
3329      ki = gammaRates[i];
3330      kisqr = ki * ki;
3331
3332      for(l = 1; l < 7; l++)
3333        {
3334          diagptable[i * 28 + l * 4]     = EXP(EIGN[l-1] * ki * lz);
3335          diagptable[i * 28 + l * 4 + 1] = EIGN[l-1] * ki;
3336          diagptable[i * 28 + l * 4 + 2] = EIGN[l-1] * EIGN[l-1] * kisqr;
3337        }
3338    }
3339
3340
3341   for(i = 0; i < upper; i++)
3342     {
3343       sum = &sumtable[i * 28];
3344       inv_Li   = 0.0;
3345       dlnLidlz = 0.0;
3346       d2lnLidlz2 = 0.0;
3347
3348       for(j = 0; j < 4; j++)
3349        {
3350          inv_Li += sum[j * 7];
3351
3352          for(l = 1; l < 7; l++)
3353            {
3354              inv_Li     += (tmp = diagptable[j * 28 + l * 4] * sum[j * 7 + l]);
3355              dlnLidlz   +=  tmp * diagptable[j * 28 + l * 4 + 1];
3356              d2lnLidlz2 +=  tmp * diagptable[j * 28 + l * 4 + 2];
3357            }
3358        }
3359
3360       inv_Li *= scaler;
3361
3362       if(iptr[i] < 7)
3363         inv_Li += freqs[iptr[i]];
3364
3365       inv_Li = 1.0 / inv_Li;
3366
3367       dlnLidlz   *= inv_Li;
3368       d2lnLidlz2 *= inv_Li;
3369
3370       dlnLidlz *= scaler;
3371       d2lnLidlz2 *= scaler;
3372
3373       dlnLdlz  += wrptr[i] * dlnLidlz;
3374       d2lnLdlz2 += wrptr[i] * (d2lnLidlz2 - dlnLidlz * dlnLidlz);
3375     }
3376
3377   *ext_dlnLdlz   = dlnLdlz;
3378   *ext_d2lnLdlz2 = d2lnLdlz2;
3379}
3380
3381
3382static void getVects(tree *tr, unsigned char **tipX1, unsigned char **tipX2, double **x1_start, double **x2_start, int *tipCase, int model,
3383                     double **x1_gapColumn, double **x2_gapColumn, unsigned int **x1_gap, unsigned int **x2_gap)
3384{
3385  int 
3386    rateHet,
3387    states = tr->partitionData[model].states,   
3388    pNumber = tr->td[0].ti[0].pNumber, 
3389    qNumber = tr->td[0].ti[0].qNumber;
3390
3391  if(tr->rateHetModel == CAT)
3392    rateHet = 1;
3393  else
3394    rateHet = 4;
3395   
3396  *x1_start = (double*)NULL,
3397  *x2_start = (double*)NULL;
3398  *tipX1 = (unsigned char*)NULL,
3399  *tipX2 = (unsigned char*)NULL;
3400
3401  if(isTip(pNumber, tr->mxtips) || isTip(qNumber, tr->mxtips))
3402    {
3403      if(!( isTip(pNumber, tr->mxtips) && isTip(qNumber, tr->mxtips)) )
3404        {
3405          *tipCase = TIP_INNER;
3406          if(isTip(qNumber, tr->mxtips))
3407            {
3408              *tipX1 = tr->partitionData[model].yVector[qNumber];
3409              *x2_start = tr->partitionData[model].xVector[pNumber - tr->mxtips - 1];
3410             
3411              if(tr->saveMemory)
3412                {
3413                  *x2_gap = &(tr->partitionData[model].gapVector[pNumber * tr->partitionData[model].gapVectorLength]);
3414                  *x2_gapColumn   = &tr->partitionData[model].gapColumn[(pNumber - tr->mxtips - 1) * states * rateHet]; 
3415                }
3416            }
3417          else
3418            {
3419              *tipX1 = tr->partitionData[model].yVector[pNumber];
3420              *x2_start = tr->partitionData[model].xVector[qNumber - tr->mxtips - 1];
3421             
3422              if(tr->saveMemory)
3423                {
3424                  *x2_gap = &(tr->partitionData[model].gapVector[qNumber * tr->partitionData[model].gapVectorLength]);
3425                  *x2_gapColumn   = &tr->partitionData[model].gapColumn[(qNumber - tr->mxtips - 1) * states * rateHet];
3426                }
3427            }
3428        }
3429      else
3430        {
3431          *tipCase = TIP_TIP;
3432          *tipX1 = tr->partitionData[model].yVector[pNumber];
3433          *tipX2 = tr->partitionData[model].yVector[qNumber];
3434        }
3435    }
3436  else
3437    {
3438      *tipCase = INNER_INNER;
3439
3440      *x1_start = tr->partitionData[model].xVector[pNumber - tr->mxtips - 1];
3441      *x2_start = tr->partitionData[model].xVector[qNumber - tr->mxtips - 1];
3442     
3443      if(tr->saveMemory)
3444        {
3445          *x1_gap = &(tr->partitionData[model].gapVector[pNumber * tr->partitionData[model].gapVectorLength]);
3446          *x1_gapColumn   = &tr->partitionData[model].gapColumn[(pNumber - tr->mxtips - 1) * states * rateHet]; 
3447     
3448          *x2_gap = &(tr->partitionData[model].gapVector[qNumber * tr->partitionData[model].gapVectorLength]);
3449          *x2_gapColumn   = &tr->partitionData[model].gapColumn[(qNumber - tr->mxtips - 1) * states * rateHet]; 
3450        }
3451    }
3452
3453}
3454
3455
3456
3457
3458void makenewzIterative(tree *tr)
3459{
3460  int 
3461    model, 
3462    tipCase;
3463
3464  double
3465    *x1_start = (double*)NULL,
3466    *x2_start = (double*)NULL,
3467    *x1_gapColumn = (double*)NULL,
3468    *x2_gapColumn = (double*)NULL;
3469
3470  unsigned char
3471    *tipX1,
3472    *tipX2; 
3473 
3474  unsigned int
3475    *x1_gap = (unsigned int*)NULL,
3476    *x2_gap = (unsigned int*)NULL;                           
3477 
3478  newviewIterative(tr);
3479
3480  for(model = 0; model < tr->NumberOfModels; model++)
3481    {
3482      /*printf("MNZIterative %d %d\n", model, tr->executeModel[model]);*/
3483     
3484      if(tr->executeModel[model])
3485        {
3486          int 
3487            width = tr->partitionData[model].width,
3488            states = tr->partitionData[model].states;
3489
3490         
3491          getVects(tr, &tipX1, &tipX2, &x1_start, &x2_start, &tipCase, model, &x1_gapColumn, &x2_gapColumn, &x1_gap, &x2_gap);
3492         
3493
3494          switch(tr->partitionData[model].dataType)
3495            {
3496            case BINARY_DATA:
3497              switch(tr->rateHetModel)
3498                {
3499                case CAT:
3500                  sumCAT_BINARY(tipCase, tr->partitionData[model].sumBuffer, x1_start, x2_start, tr->partitionData[model].tipVector, tipX1, tipX2,
3501                                width);
3502                  break;
3503                case GAMMA:
3504                case GAMMA_I:
3505                  sumGAMMA_BINARY(tipCase, tr->partitionData[model].sumBuffer, x1_start, x2_start, tr->partitionData[model].tipVector, tipX1, tipX2,
3506                                  width);
3507                  break;
3508                default:
3509                  assert(0);
3510                }
3511              break;
3512            case DNA_DATA:
3513              switch(tr->rateHetModel)
3514                {
3515                case CAT:               
3516#ifdef __SIM_SSE3
3517                  if(tr->saveMemory)
3518                    {
3519                     
3520                      sumCAT_SAVE(tipCase, tr->partitionData[model].sumBuffer, x1_start, x2_start, tr->partitionData[model].tipVector, tipX1, tipX2,
3521                                  width, x1_gapColumn, x2_gapColumn, x1_gap, x2_gap);
3522                    }
3523                  else
3524#endif
3525                    sumCAT(tipCase, tr->partitionData[model].sumBuffer, x1_start, x2_start, tr->partitionData[model].tipVector, tipX1, tipX2,
3526                           width);
3527                  break;
3528                case GAMMA:
3529                case GAMMA_I:             
3530                      {
3531                        if(tr->saveMemory)
3532                          sumGAMMA_GAPPED_SAVE(tipCase, tr->partitionData[model].sumBuffer, x1_start, x2_start, tr->partitionData[model].tipVector, tipX1, tipX2,
3533                                               width, x1_gapColumn, x2_gapColumn, x1_gap, x2_gap);
3534                        else                     
3535                          sumGAMMA(tipCase, tr->partitionData[model].sumBuffer, x1_start, x2_start, tr->partitionData[model].tipVector, tipX1, tipX2,
3536                                   width);                       
3537                      }
3538                  break;
3539                default:
3540                  assert(0);
3541                }
3542              break;
3543            case AA_DATA:
3544              switch(tr->rateHetModel)
3545                {
3546                case CAT:       
3547#ifdef __SIM_SSE3
3548                  if(tr->saveMemory)     
3549                    sumGTRCATPROT_SAVE(tipCase, tr->partitionData[model].sumBuffer, x1_start, x2_start, tr->partitionData[model].tipVector,
3550                                       tipX1, tipX2, width, x1_gapColumn, x2_gapColumn, x1_gap, x2_gap);
3551                  else
3552#endif
3553                    sumGTRCATPROT(tipCase, tr->partitionData[model].sumBuffer, x1_start, x2_start, tr->partitionData[model].tipVector,
3554                                  tipX1, tipX2, width);
3555                  break;
3556                case GAMMA:
3557                case GAMMA_I:             
3558#ifdef __SIM_SSE3
3559                  if(tr->saveMemory)
3560                    sumGAMMAPROT_GAPPED_SAVE(tipCase, tr->partitionData[model].sumBuffer, x1_start, x2_start, tr->partitionData[model].tipVector, tipX1, tipX2,
3561                                             width, x1_gapColumn, x2_gapColumn, x1_gap, x2_gap);
3562                  else
3563#endif           
3564                    if(tr->partitionData[model].protModels == LG4 || tr->partitionData[model].protModels == LG4X)                     
3565                      {                 
3566                        sumGAMMAPROT_LG4(tipCase,  tr->partitionData[model].sumBuffer, x1_start, x2_start, tr->partitionData[model].tipVector_LG4,
3567                                         tipX1, tipX2, width);
3568                      }
3569                    else
3570                      sumGAMMAPROT(tipCase,  tr->partitionData[model].sumBuffer, x1_start, x2_start, tr->partitionData[model].tipVector,
3571                                   tipX1, tipX2, width);                       
3572                  break;
3573                default:
3574                  assert(0);
3575                }
3576              break;
3577            case SECONDARY_DATA:
3578              switch(tr->rateHetModel)
3579                {
3580                case CAT:
3581                  sumGTRCATSECONDARY(tipCase, tr->partitionData[model].sumBuffer, x1_start, x2_start, tr->partitionData[model].tipVector,
3582                                     tipX1, tipX2, width);
3583                  break;
3584                case GAMMA:
3585                case GAMMA_I:
3586                  sumGAMMASECONDARY(tipCase,  tr->partitionData[model].sumBuffer, x1_start, x2_start, tr->partitionData[model].tipVector,
3587                                    tipX1, tipX2, width);
3588                  break;
3589                default:
3590                  assert(0);
3591                }
3592              break;
3593            case SECONDARY_DATA_6:
3594              switch(tr->rateHetModel)
3595                {
3596                case CAT:
3597                  sumGTRCATSECONDARY_6(tipCase, tr->partitionData[model].sumBuffer, x1_start, x2_start, tr->partitionData[model].tipVector,
3598                                       tipX1, tipX2, width);
3599                  break;
3600                case GAMMA:
3601                case GAMMA_I:
3602                  sumGAMMASECONDARY_6(tipCase,  tr->partitionData[model].sumBuffer, x1_start, x2_start, tr->partitionData[model].tipVector,
3603                                      tipX1, tipX2, width);
3604                  break;
3605                default:
3606                  assert(0);
3607                }
3608              break;
3609            case SECONDARY_DATA_7:
3610              switch(tr->rateHetModel)
3611                {
3612                case CAT:
3613                  sumGTRCATSECONDARY_7(tipCase, tr->partitionData[model].sumBuffer, x1_start, x2_start, tr->partitionData[model].tipVector,
3614                                       tipX1, tipX2, width);
3615                  break;
3616                case GAMMA:
3617                case GAMMA_I:
3618                  sumGAMMASECONDARY_7(tipCase,  tr->partitionData[model].sumBuffer, x1_start, x2_start, tr->partitionData[model].tipVector,
3619                                      tipX1, tipX2, width);
3620                  break;
3621                default:
3622                  assert(0);
3623                }
3624              break;
3625            case GENERIC_32:
3626               switch(tr->rateHetModel)
3627                {
3628                case CAT:
3629                  sumCatFlex(tipCase, tr->partitionData[model].sumBuffer, x1_start, x2_start, tr->partitionData[model].tipVector,
3630                             tipX1, tipX2, width, states);
3631                  break;
3632                case GAMMA:
3633                case GAMMA_I:
3634                  sumGammaFlex(tipCase,  tr->partitionData[model].sumBuffer, x1_start, x2_start, tr->partitionData[model].tipVector,
3635                               tipX1, tipX2, width, states);
3636                  break;
3637                default:
3638                  assert(0);
3639                }
3640              break;
3641            default:
3642              assert(0);
3643            }
3644        }
3645    }
3646}
3647
3648
3649
3650
3651
3652void execCore(tree *tr, volatile double *_dlnLdlz, volatile double *_d2lnLdlz2)
3653{
3654  int 
3655    model, 
3656    branchIndex;
3657
3658  double 
3659    lz;
3660
3661  _dlnLdlz[0]   = 0.0;
3662  _d2lnLdlz2[0] = 0.0;
3663
3664#ifdef _DEBUG_MULTI_EPA 
3665  printf("MNZC: ");
3666#endif 
3667
3668  for(model = 0; model < tr->NumberOfModels; model++)
3669    {
3670     
3671#ifdef _DEBUG_MULTI_EPA 
3672          printf("%d ", tr->executeModel[model]);
3673#endif
3674      if(tr->executeModel[model])
3675        {
3676          int 
3677            width = tr->partitionData[model].width,
3678            states = tr->partitionData[model].states;
3679         
3680          double 
3681            *weights         = tr->partitionData[model].weights,
3682            *sumBuffer       = (double*)NULL;
3683         
3684
3685          volatile double
3686            dlnLdlz   = 0.0,
3687            d2lnLdlz2 = 0.0;
3688         
3689
3690            sumBuffer = tr->partitionData[model].sumBuffer;
3691
3692
3693          if(tr->multiBranch)
3694            {
3695              branchIndex = model;
3696              lz = tr->coreLZ[model];
3697              _dlnLdlz[model]   = 0.0;
3698              _d2lnLdlz2[model] = 0.0;
3699            }
3700          else
3701            {
3702              branchIndex = 0;
3703              lz = tr->coreLZ[0];
3704            }
3705
3706         
3707
3708          switch(tr->partitionData[model].dataType)
3709            {
3710            case BINARY_DATA:
3711              switch(tr->rateHetModel)
3712                {
3713                case CAT:
3714                  coreGTRCAT_BINARY(width, tr->partitionData[model].numberOfCategories, sumBuffer,
3715                                    &dlnLdlz, &d2lnLdlz2, 
3716                                    tr->partitionData[model].perSiteRates, tr->partitionData[model].EIGN,  tr->partitionData[model].rateCategory, lz, tr->partitionData[model].wgt);
3717                  break;
3718                case GAMMA:
3719                  coreGTRGAMMA_BINARY(width, sumBuffer,
3720                                      &dlnLdlz, &d2lnLdlz2, tr->partitionData[model].EIGN, tr->partitionData[model].gammaRates, lz,
3721                                      tr->partitionData[model].wgt);
3722                  break;
3723                case GAMMA_I:
3724                  coreGTRGAMMAINVAR_BINARY(tr->partitionData[model].propInvariant, tr->partitionData[model].frequencies,
3725                                           tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN,
3726                                           sumBuffer, &dlnLdlz, &d2lnLdlz2,
3727                                           tr->partitionData[model].invariant, tr->partitionData[model].wgt,
3728                                           width, lz);
3729                  break;
3730                default:
3731                  assert(0);
3732                }
3733              break;
3734            case DNA_DATA:           
3735                {
3736                  switch(tr->rateHetModel)
3737                    {
3738                    case CAT:
3739                      coreGTRCAT(width, tr->partitionData[model].numberOfCategories, sumBuffer,
3740                                 &dlnLdlz, &d2lnLdlz2,
3741                                 tr->partitionData[model].perSiteRates, tr->partitionData[model].EIGN,  tr->partitionData[model].rateCategory, lz, tr->partitionData[model].wgt);
3742                      break;
3743                    case GAMMA:             
3744                      coreGTRGAMMA(width, sumBuffer,
3745                                   &dlnLdlz, &d2lnLdlz2, tr->partitionData[model].EIGN, tr->partitionData[model].gammaRates, lz,
3746                                   tr->partitionData[model].wgt);
3747                      break;
3748                    case GAMMA_I:
3749                      coreGTRGAMMAINVAR(tr->partitionData[model].propInvariant, tr->partitionData[model].frequencies,
3750                                        tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN,
3751                                        sumBuffer, &dlnLdlz, &d2lnLdlz2,
3752                                        tr->partitionData[model].invariant, tr->partitionData[model].wgt,
3753                                        width, lz);
3754                      break;
3755                    default:
3756                      assert(0);
3757                    }
3758                }
3759              break;
3760            case AA_DATA:
3761             
3762                {
3763                  switch(tr->rateHetModel)
3764                    {
3765                    case CAT:
3766                      coreGTRCATPROT(tr->partitionData[model].EIGN, lz, tr->partitionData[model].numberOfCategories,  tr->partitionData[model].perSiteRates,
3767                                     tr->partitionData[model].rateCategory, width,                                 
3768                                     &dlnLdlz, &d2lnLdlz2,
3769                                     sumBuffer, tr->partitionData[model].wgt);
3770                      break;
3771                    case GAMMA: 
3772                      if(tr->partitionData[model].protModels == LG4 || tr->partitionData[model].protModels == LG4X)
3773                        {                                                 
3774                          coreGTRGAMMAPROT_LG4(tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN_LG4,
3775                                              sumBuffer, width, tr->partitionData[model].wgt,
3776                                               &dlnLdlz, &d2lnLdlz2, lz, weights);
3777
3778                        }
3779                      else
3780                        coreGTRGAMMAPROT(tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN,
3781                                         sumBuffer, width, tr->partitionData[model].wgt,
3782                                         &dlnLdlz, &d2lnLdlz2, lz);
3783                      break;
3784                    case GAMMA_I:
3785                      coreGTRGAMMAPROTINVAR(tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN,
3786                                            sumBuffer, width, tr->partitionData[model].wgt,
3787                                            &dlnLdlz, &d2lnLdlz2, lz, tr->partitionData[model].frequencies,
3788                                            tr->partitionData[model].propInvariant, tr->partitionData[model].invariant);
3789                      break;
3790                    default:
3791                      assert(0);
3792                    }
3793                }
3794              break;
3795            case SECONDARY_DATA:             
3796              switch(tr->rateHetModel)
3797                {
3798                case CAT:
3799                  coreGTRCATSECONDARY(tr->partitionData[model].EIGN, lz, tr->partitionData[model].numberOfCategories,  tr->partitionData[model].perSiteRates,
3800                                      tr->partitionData[model].rateCategory, width,                                   
3801                                      &dlnLdlz, &d2lnLdlz2,
3802                                      sumBuffer, tr->partitionData[model].wgt);
3803                  break;
3804                case GAMMA:
3805                  coreGTRGAMMASECONDARY(tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN,
3806                                        sumBuffer, width, tr->partitionData[model].wgt,
3807                                        &dlnLdlz, &d2lnLdlz2, lz);
3808                  break;
3809                case GAMMA_I:           
3810                  coreGTRGAMMASECONDARYINVAR(tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN,
3811                                             sumBuffer, width, tr->partitionData[model].wgt,
3812                                             &dlnLdlz, &d2lnLdlz2, lz, tr->partitionData[model].frequencies,
3813                                             tr->partitionData[model].propInvariant, tr->partitionData[model].invariant);
3814                  break;
3815                default:
3816                  assert(0);
3817                }
3818              break;
3819            case SECONDARY_DATA_6:           
3820              switch(tr->rateHetModel)
3821                {
3822                case CAT:
3823                  coreGTRCATSECONDARY_6(tr->partitionData[model].EIGN, lz, tr->partitionData[model].numberOfCategories,  tr->partitionData[model].perSiteRates,
3824                                        tr->partitionData[model].rateCategory, width,                                   
3825                                        &dlnLdlz, &d2lnLdlz2,
3826                                        sumBuffer, tr->partitionData[model].wgt);
3827                  break;
3828                case GAMMA:
3829                  coreGTRGAMMASECONDARY_6(tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN,
3830                                          sumBuffer, width, tr->partitionData[model].wgt,
3831                                          &dlnLdlz, &d2lnLdlz2, lz);
3832                  break;
3833                case GAMMA_I:
3834                  coreGTRGAMMASECONDARYINVAR_6(tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN,
3835                                               sumBuffer, width, tr->partitionData[model].wgt,
3836                                               &dlnLdlz, &d2lnLdlz2, lz, tr->partitionData[model].frequencies,
3837                                               tr->partitionData[model].propInvariant, tr->partitionData[model].invariant);
3838                  break;
3839                default:
3840                  assert(0);
3841                }
3842              break;
3843            case SECONDARY_DATA_7:         
3844              switch(tr->rateHetModel)
3845                {
3846                case CAT:
3847                  coreGTRCATSECONDARY_7(tr->partitionData[model].EIGN, lz, tr->partitionData[model].numberOfCategories,  tr->partitionData[model].perSiteRates,
3848                                        tr->partitionData[model].rateCategory, width,                                 
3849                                        &dlnLdlz, &d2lnLdlz2,
3850                                        sumBuffer, tr->partitionData[model].wgt);
3851                  break;
3852                case GAMMA:
3853                  coreGTRGAMMASECONDARY_7(tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN,
3854                                          sumBuffer, width, tr->partitionData[model].wgt,
3855                                          &dlnLdlz, &d2lnLdlz2, lz);
3856                  break;
3857                case GAMMA_I:
3858                  coreGTRGAMMASECONDARYINVAR_7(tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN,
3859                                               sumBuffer, width, tr->partitionData[model].wgt,
3860                                               &dlnLdlz, &d2lnLdlz2, lz, tr->partitionData[model].frequencies,
3861                                               tr->partitionData[model].propInvariant, tr->partitionData[model].invariant);
3862                  break;
3863                default:
3864                  assert(0);
3865                }
3866              break;
3867            case GENERIC_32:
3868              switch(tr->rateHetModel)
3869                {
3870                case CAT:
3871                  coreCatFlex(tr->partitionData[model].EIGN, lz, tr->partitionData[model].numberOfCategories,  tr->partitionData[model].perSiteRates,
3872                              tr->partitionData[model].rateCategory, width,                           
3873                              &dlnLdlz, &d2lnLdlz2,
3874                              sumBuffer, states, tr->partitionData[model].wgt);
3875                  break;
3876                case GAMMA:
3877                  coreGammaFlex(tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN,
3878                                sumBuffer, width, tr->partitionData[model].wgt,
3879                                &dlnLdlz, &d2lnLdlz2, lz, states);
3880                  break;
3881                case GAMMA_I:
3882                  coreGammaInvarFlex(tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN,
3883                                     sumBuffer, width, tr->partitionData[model].wgt,
3884                                     &dlnLdlz, &d2lnLdlz2, lz, tr->partitionData[model].frequencies,
3885                                     tr->partitionData[model].propInvariant, tr->partitionData[model].invariant, states);
3886                  break;
3887                default:
3888                  assert(0);
3889                }
3890              break;
3891            default:
3892              assert(0);
3893            }
3894
3895          _dlnLdlz[branchIndex]   = _dlnLdlz[branchIndex]   + dlnLdlz;
3896          _d2lnLdlz2[branchIndex] = _d2lnLdlz2[branchIndex] + d2lnLdlz2;
3897        }
3898    }
3899
3900#ifdef _DEBUG_MULTI_EPA
3901  printf("\n");
3902#endif   
3903
3904}
3905
3906
3907
3908
3909static void topLevelMakenewz(tree *tr, double *z0, int _maxiter, double *result)
3910{
3911  double   z[NUM_BRANCHES], zprev[NUM_BRANCHES], zstep[NUM_BRANCHES];
3912  volatile double  dlnLdlz[NUM_BRANCHES], d2lnLdlz2[NUM_BRANCHES];
3913  int i, maxiter[NUM_BRANCHES], numBranches, model;
3914  boolean
3915    firstIteration = TRUE,
3916    outerConverged[NUM_BRANCHES],
3917    loopConverged;
3918
3919  if(tr->multiBranch)
3920    numBranches = tr->NumberOfModels;
3921  else
3922    numBranches = 1;
3923
3924  for(i = 0; i < numBranches; i++)
3925    {
3926      z[i] = z0[i];       
3927
3928      maxiter[i] = _maxiter;
3929     
3930      outerConverged[i] = FALSE;
3931      tr->curvatOK[i]       = TRUE;     
3932    }
3933
3934  if(tr->perPartitionEPA)   
3935    {
3936      if(tr->multiBranch)
3937        for(i = 0; i < numBranches; i++)
3938          if(!(tr->executeModel[i]))
3939            outerConverged[i] = TRUE;
3940    }
3941
3942  do
3943    {
3944      for(i = 0; i < numBranches; i++)
3945        {
3946          if(outerConverged[i] == FALSE && tr->curvatOK[i] == TRUE)
3947            {
3948              tr->curvatOK[i] = FALSE;
3949
3950              zprev[i] = z[i];
3951
3952              zstep[i] = (1.0 - zmax) * z[i] + zmin;
3953            }
3954        }
3955
3956      for(i = 0; i < numBranches; i++)
3957        {
3958          if(outerConverged[i] == FALSE && tr->curvatOK[i] == FALSE)
3959            {
3960              double lz;
3961
3962              if (z[i] < zmin)                   
3963                z[i] = zmin;           
3964              else
3965                {
3966                  if(z[i] > zmax)                   
3967                    z[i] = zmax;                                   
3968                }
3969
3970              lz    = log(z[i]);
3971
3972              tr->coreLZ[i] = lz;
3973            }
3974        }
3975
3976      if(tr->multiBranch)
3977        {
3978          for(model = 0; model < tr->NumberOfModels; model++)
3979            {
3980              if(tr->executeModel[model])
3981                tr->executeModel[model] = !tr->curvatOK[model];
3982
3983            }
3984        }
3985      else
3986        {
3987          for(model = 0; model < tr->NumberOfModels; model++)
3988            {
3989              if(tr->perPartitionEPA)
3990                {
3991                  if(tr->executeModel[model])
3992                    tr->executeModel[model] = !tr->curvatOK[0];
3993                }
3994              else
3995                tr->executeModel[model] = !tr->curvatOK[0];
3996            }
3997        }
3998
3999
4000#ifdef _USE_PTHREADS
4001      if(firstIteration)
4002        {
4003          masterBarrier(THREAD_MAKENEWZ_FIRST, tr);
4004          firstIteration = FALSE;
4005        }
4006      else
4007        masterBarrier(THREAD_MAKENEWZ, tr);
4008
4009      if(!tr->multiBranch)
4010        {
4011          dlnLdlz[0] = 0.0;
4012          d2lnLdlz2[0] = 0.0;
4013          for(i = 0; i < NumberOfThreads; i++)
4014            {
4015              dlnLdlz[0]   += reductionBuffer[i];
4016              d2lnLdlz2[0] += reductionBufferTwo[i];
4017            }
4018        }
4019      else
4020        {
4021          int j;
4022          for(j = 0; j < tr->NumberOfModels; j++)
4023            {
4024              dlnLdlz[j] = 0.0;
4025              d2lnLdlz2[j] = 0.0;
4026              for(i = 0; i < NumberOfThreads; i++)
4027                {
4028                  dlnLdlz[j]   += reductionBuffer[i * tr->NumberOfModels + j];
4029                  d2lnLdlz2[j] += reductionBufferTwo[i * tr->NumberOfModels + j];
4030                }
4031            }
4032        }
4033#else
4034      if(firstIteration)
4035        {
4036          makenewzIterative(tr);
4037          firstIteration = FALSE;
4038        }
4039      execCore(tr, dlnLdlz, d2lnLdlz2);
4040#endif
4041
4042      for(i = 0; i < numBranches; i++)
4043        {
4044          if(outerConverged[i] == FALSE && tr->curvatOK[i] == FALSE)
4045            {
4046              if ((d2lnLdlz2[i] >= 0.0) && (z[i] < zmax))
4047                zprev[i] = z[i] = 0.37 * z[i] + 0.63;  /*  Bad curvature, shorten branch */
4048              else
4049                tr->curvatOK[i] = TRUE;
4050            }
4051        }
4052
4053      for(i = 0; i < numBranches; i++)
4054        {
4055          if(tr->curvatOK[i] == TRUE && outerConverged[i] == FALSE)
4056            {
4057              if (d2lnLdlz2[i] < 0.0)
4058                {
4059                  double 
4060                    tantmp = -dlnLdlz[i] / d2lnLdlz2[i];
4061                 
4062                  if(tantmp < 100)
4063                    {
4064                      z[i] *= EXP(tantmp);
4065                      if (z[i] < zmin)
4066                        z[i] = zmin;
4067
4068                      if (z[i] > 0.25 * zprev[i] + 0.75)
4069                        z[i] = 0.25 * zprev[i] + 0.75;
4070                    }
4071                  else
4072                    z[i] = 0.25 * zprev[i] + 0.75;
4073                }
4074             
4075              if(z[i] > zmax) 
4076                z[i] = zmax;
4077
4078              maxiter[i] = maxiter[i] - 1;
4079
4080              /* the previous termination condition did not guarantee an improvement in LnL in
4081                 cases where the initial branch was very long ! */
4082              //if(maxiter[i] > 0 && (ABS(z[i] - zprev[i]) > zstep[i]))
4083              if((ABS(z[i] - zprev[i]) > zstep[i]))
4084                {
4085                  /* We should make a more informed decision here,
4086                     based on the log like improvement */
4087
4088                  if(maxiter[i] < -20)
4089                    {
4090                      z[i] = z0[i];
4091                      outerConverged[i] = TRUE;
4092                    }
4093                  else
4094                    outerConverged[i] = FALSE;
4095                }
4096              else                             
4097                outerConverged[i] = TRUE;                                               
4098            }
4099        }
4100
4101      loopConverged = TRUE;
4102      for(i = 0; i < numBranches; i++)
4103        loopConverged = loopConverged && outerConverged[i];
4104    }
4105  while (!loopConverged);
4106
4107  for(model = 0; model < tr->NumberOfModels; model++)
4108    tr->executeModel[model] = TRUE;
4109
4110  for(i = 0; i < numBranches; i++)
4111    result[i] = z[i];
4112}
4113
4114#ifdef _USE_PTHREADS
4115
4116static void sumClassify(tree *tr, int tipCase, double *_x1, double *_x2, unsigned char *_tipX1, unsigned char *_tipX2, boolean *executeModel)
4117{
4118  int 
4119    model = 0,
4120    columnCounter = 0, 
4121    offsetCounter = 0;
4122 
4123  double
4124    *x1_start = (double *)NULL,
4125    *x2_start = (double *)NULL;
4126
4127  unsigned char 
4128    *tipX1 = (unsigned char*)NULL,
4129    *tipX2 = (unsigned char*)NULL;
4130
4131#ifdef _DEBUG_MULTI_EPA
4132  if(tr->threadID == THREAD_TO_DEBUG)
4133    printf("MNZS: ");
4134#endif
4135
4136  for(model = 0; model < tr->NumberOfModels; model++)
4137    {
4138      int 
4139        width = tr->partitionData[model].upper - tr->partitionData[model].lower;       
4140     
4141#ifdef _DEBUG_MULTI_EPA
4142      if(tr->threadID == THREAD_TO_DEBUG)
4143        printf("%d", executeModel[model]);
4144#endif       
4145
4146      if(executeModel[model])
4147        {
4148          double *sumBuffer = &tr->temporarySumBuffer[offsetCounter];
4149
4150          switch(tipCase)
4151            {
4152            case TIP_TIP:
4153              tipX1 = &_tipX1[columnCounter];
4154              tipX2 = &_tipX2[columnCounter];
4155              break;
4156            case TIP_INNER:
4157              tipX1 = &_tipX1[columnCounter];
4158              x2_start    = &_x2[offsetCounter];
4159              break;
4160            case INNER_INNER:
4161              x1_start = &_x1[offsetCounter];
4162              x2_start = &_x2[offsetCounter];
4163              break;
4164            default:
4165              assert(0);             
4166            }
4167
4168          switch(tr->partitionData[model].dataType)
4169            {
4170            case BINARY_DATA:
4171              switch(tr->rateHetModel)
4172                {
4173                case CAT:               
4174                  sumCAT_BINARY(tipCase, sumBuffer, x1_start, x2_start, tr->partitionData[model].tipVector, tipX1, tipX2,
4175                                width);
4176                  break;
4177                case GAMMA:
4178                case GAMMA_I:
4179                  sumGAMMA_BINARY(tipCase, sumBuffer, x1_start, x2_start, tr->partitionData[model].tipVector, tipX1, tipX2,
4180                                  width);
4181                  break;
4182                default:
4183                  assert(0);
4184                }
4185              break;
4186            case DNA_DATA:
4187              switch(tr->rateHetModel)
4188                {
4189                case CAT:                 
4190                  sumCAT(tipCase, sumBuffer, x1_start, x2_start, tr->partitionData[model].tipVector, tipX1, tipX2,
4191                         width);
4192                  break;
4193                case GAMMA:
4194                case GAMMA_I:             
4195                  sumGAMMA(tipCase, sumBuffer, x1_start, x2_start, tr->partitionData[model].tipVector, tipX1, tipX2,
4196                           width);
4197                  break;
4198                default:
4199                  assert(0);
4200                }
4201              break;
4202            case AA_DATA:
4203              switch(tr->rateHetModel)
4204                {
4205                case CAT:               
4206                  sumGTRCATPROT(tipCase, sumBuffer, x1_start, x2_start, tr->partitionData[model].tipVector,
4207                                tipX1, tipX2, width);
4208                  break;
4209                case GAMMA:
4210                case GAMMA_I:           
4211                  sumGAMMAPROT(tipCase,  sumBuffer, x1_start, x2_start, tr->partitionData[model].tipVector,
4212                               tipX1, tipX2, width);
4213                  break;
4214                default:
4215                  assert(0);
4216                }
4217              break;
4218            case SECONDARY_DATA:
4219              switch(tr->rateHetModel)
4220                {
4221                case CAT:
4222                  sumGTRCATSECONDARY(tipCase, sumBuffer, x1_start, x2_start, tr->partitionData[model].tipVector,
4223                                     tipX1, tipX2, width);
4224                  break;
4225                case GAMMA:
4226                case GAMMA_I:
4227                  sumGAMMASECONDARY(tipCase,  sumBuffer, x1_start, x2_start, tr->partitionData[model].tipVector,
4228                                    tipX1, tipX2, width);
4229                  break;
4230                default:
4231                  assert(0);
4232                }
4233              break;
4234            case SECONDARY_DATA_6:
4235              switch(tr->rateHetModel)
4236                {
4237                case CAT:
4238                  sumGTRCATSECONDARY_6(tipCase, sumBuffer, x1_start, x2_start, tr->partitionData[model].tipVector,
4239                                       tipX1, tipX2, width);
4240                  break;
4241                case GAMMA:
4242                case GAMMA_I:
4243                  sumGAMMASECONDARY_6(tipCase,  sumBuffer, x1_start, x2_start, tr->partitionData[model].tipVector,
4244                                      tipX1, tipX2, width);
4245                  break;
4246                default:
4247                  assert(0);
4248                }
4249              break;
4250            case SECONDARY_DATA_7:
4251              switch(tr->rateHetModel)
4252                {
4253                case CAT:
4254                  sumGTRCATSECONDARY_7(tipCase, sumBuffer, x1_start, x2_start, tr->partitionData[model].tipVector,
4255                                       tipX1, tipX2, width);
4256                  break;
4257                case GAMMA:
4258                case GAMMA_I:
4259                  sumGAMMASECONDARY_7(tipCase,  sumBuffer, x1_start, x2_start, tr->partitionData[model].tipVector,
4260                                      tipX1, tipX2, width);
4261                  break;
4262                default:
4263                  assert(0);
4264                }
4265              break;
4266            default:
4267              assert(0);
4268            }
4269        }
4270     
4271      columnCounter += width;
4272      offsetCounter += width * tr->partitionData[model].states * tr->discreteRateCategories;
4273    } 
4274#ifdef _DEBUG_MULTI_EPA
4275  if(tr->threadID == THREAD_TO_DEBUG)
4276    printf("\n");
4277#endif
4278}
4279
4280static void coreClassify(tree *tr, volatile double *_dlnLdlz, volatile double *_d2lnLdlz2, double *coreLZ, boolean *executeModel)
4281{
4282  int 
4283    model, 
4284    branchIndex,
4285    offsetCounter = 0,
4286    columnCounter = 0;
4287  double lz;
4288
4289  _dlnLdlz[0]   = 0.0;
4290  _d2lnLdlz2[0] = 0.0;
4291 
4292#ifdef _DEBUG_MULTI_EPA
4293  if(tr->threadID == THREAD_TO_DEBUG)
4294    printf("MNZC: ");
4295#endif
4296 
4297  for(model = 0; model < tr->NumberOfModels; model++)
4298    {
4299      int 
4300        width = tr->partitionData[model].upper - tr->partitionData[model].lower;
4301
4302      if(tr->multiBranch)
4303        branchIndex = model;
4304      else
4305        branchIndex = 0;       
4306
4307     
4308#ifdef _DEBUG_MULTI_EPA
4309  if(tr->threadID == THREAD_TO_DEBUG)
4310    printf("%d", executeModel[model]);
4311#endif   
4312      if(executeModel[model])
4313        {         
4314          double 
4315            *sumBuffer       = &tr->temporarySumBuffer[offsetCounter],
4316            *patrat          = tr->partitionData[model].perSiteRates;
4317         
4318          int 
4319            *rateCategory  = &tr->contiguousRateCategory[columnCounter],
4320            *wgt           = &tr->contiguousWgt[columnCounter],
4321            *invariant     = &tr->contiguousInvariant[columnCounter];
4322         
4323
4324          volatile double
4325            dlnLdlz   = 0.0,
4326            d2lnLdlz2 = 0.0;             
4327
4328          if(tr->multiBranch)
4329            {
4330              branchIndex = model;
4331              lz = coreLZ[model];
4332              _dlnLdlz[model]   = 0.0;
4333              _d2lnLdlz2[model] = 0.0;
4334            }
4335          else
4336            {
4337              branchIndex = 0;
4338              lz = coreLZ[0];
4339            }
4340
4341          switch(tr->partitionData[model].dataType)
4342            {
4343            case BINARY_DATA:
4344              switch(tr->rateHetModel)
4345                {
4346                case CAT:
4347                  coreGTRCAT_BINARY(width, tr->partitionData[model].numberOfCategories, sumBuffer,
4348                                    &dlnLdlz, &d2lnLdlz2,
4349                                    patrat, tr->partitionData[model].EIGN,  rateCategory, lz, wgt);
4350                  break;
4351                case GAMMA:
4352                  coreGTRGAMMA_BINARY(width, sumBuffer,
4353                                      &dlnLdlz, &d2lnLdlz2, tr->partitionData[model].EIGN, tr->partitionData[model].gammaRates, lz,
4354                                      wgt);
4355                  break;
4356                case GAMMA_I:
4357                  coreGTRGAMMAINVAR_BINARY(tr->partitionData[model].propInvariant, tr->partitionData[model].frequencies,
4358                                           tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN,
4359                                           sumBuffer, &dlnLdlz, &d2lnLdlz2,
4360                                           invariant, wgt,
4361                                           width, lz);
4362                  break;
4363                default:
4364                  assert(0);
4365                }
4366              break;
4367            case DNA_DATA:           
4368              switch(tr->rateHetModel)
4369                {
4370                case CAT:
4371                  coreGTRCAT(width, tr->partitionData[model].numberOfCategories, sumBuffer,
4372                             &dlnLdlz, &d2lnLdlz2,
4373                             patrat, tr->partitionData[model].EIGN,  rateCategory, lz, tr->partitionData[model].wgt);
4374                  break;
4375                case GAMMA:
4376                  coreGTRGAMMA(width, sumBuffer,
4377                               &dlnLdlz, &d2lnLdlz2, tr->partitionData[model].EIGN, tr->partitionData[model].gammaRates, lz,
4378                               wgt);
4379                  break;
4380                case GAMMA_I:
4381                  coreGTRGAMMAINVAR(tr->partitionData[model].propInvariant, tr->partitionData[model].frequencies,
4382                                    tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN,
4383                                    sumBuffer, &dlnLdlz, &d2lnLdlz2,
4384                                    invariant, wgt,
4385                                        width, lz);
4386                  break;
4387                default:
4388                  assert(0);
4389                }               
4390              break;
4391            case AA_DATA:           
4392              switch(tr->rateHetModel)
4393                {
4394                case CAT:
4395                  coreGTRCATPROT(tr->partitionData[model].EIGN, lz, tr->partitionData[model].numberOfCategories,  patrat,
4396                                 rateCategory, width,                           
4397                                 &dlnLdlz, &d2lnLdlz2,
4398                                 sumBuffer, wgt);
4399                  break;
4400                case GAMMA:
4401                  coreGTRGAMMAPROT(tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN,
4402                                   sumBuffer, width, wgt,
4403                                   &dlnLdlz, &d2lnLdlz2, lz);
4404                  break;
4405                case GAMMA_I:
4406                  coreGTRGAMMAPROTINVAR(tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN,
4407                                        sumBuffer, width, wgt,
4408                                        &dlnLdlz, &d2lnLdlz2, lz, tr->partitionData[model].frequencies,
4409                                        tr->partitionData[model].propInvariant, invariant);
4410                  break;
4411                default:
4412                  assert(0);
4413                }           
4414              break;
4415            case SECONDARY_DATA:
4416              switch(tr->rateHetModel)
4417                {
4418                case CAT:
4419                  coreGTRCATSECONDARY(tr->partitionData[model].EIGN, lz, tr->partitionData[model].numberOfCategories,  patrat,
4420                                      rateCategory, width,                                   
4421                                      &dlnLdlz, &d2lnLdlz2,
4422                                      sumBuffer, wgt);
4423                  break;
4424                case GAMMA:
4425                  coreGTRGAMMASECONDARY(tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN,
4426                                        sumBuffer, width, wgt,
4427                                        &dlnLdlz, &d2lnLdlz2, lz);
4428                  break;
4429                case GAMMA_I:
4430                  coreGTRGAMMASECONDARYINVAR(tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN,
4431                                             sumBuffer, width, wgt,
4432                                             &dlnLdlz, &d2lnLdlz2, lz, tr->partitionData[model].frequencies,
4433                                             tr->partitionData[model].propInvariant, invariant);
4434                  break;
4435                default:
4436                  assert(0);
4437                }
4438              break;
4439            case SECONDARY_DATA_6:
4440              switch(tr->rateHetModel)
4441                {
4442                case CAT:
4443                  coreGTRCATSECONDARY_6(tr->partitionData[model].EIGN, lz, tr->partitionData[model].numberOfCategories, patrat,
4444                                        rateCategory, width,                                   
4445                                        &dlnLdlz, &d2lnLdlz2,
4446                                        sumBuffer, wgt);
4447                  break;
4448                case GAMMA:
4449                  coreGTRGAMMASECONDARY_6(tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN,
4450                                          sumBuffer, width, wgt,
4451                                          &dlnLdlz, &d2lnLdlz2, lz);
4452                  break;
4453                case GAMMA_I:
4454                  coreGTRGAMMASECONDARYINVAR_6(tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN,
4455                                               sumBuffer, width, wgt,
4456                                               &dlnLdlz, &d2lnLdlz2, lz, tr->partitionData[model].frequencies,
4457                                               tr->partitionData[model].propInvariant, invariant);
4458                  break;
4459                default:
4460                  assert(0);
4461                }
4462              break;
4463            case SECONDARY_DATA_7:
4464              switch(tr->rateHetModel)
4465                {
4466                case CAT:
4467                  coreGTRCATSECONDARY_7(tr->partitionData[model].EIGN, lz, tr->partitionData[model].numberOfCategories,  patrat,
4468                                        rateCategory, width,                                   
4469                                        &dlnLdlz, &d2lnLdlz2,
4470                                        sumBuffer, wgt);
4471                  break;
4472                case GAMMA:
4473                  coreGTRGAMMASECONDARY_7(tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN,
4474                                          sumBuffer, width, wgt,
4475                                          &dlnLdlz, &d2lnLdlz2, lz);
4476                  break;
4477                case GAMMA_I:
4478                  coreGTRGAMMASECONDARYINVAR_7(tr->partitionData[model].gammaRates, tr->partitionData[model].EIGN,
4479                                               sumBuffer, width, wgt,
4480                                               &dlnLdlz, &d2lnLdlz2, lz, tr->partitionData[model].frequencies,
4481                                               tr->partitionData[model].propInvariant, invariant);
4482                  break;
4483                default:
4484                  assert(0);
4485                }
4486              break;
4487            default:
4488              assert(0);
4489            }
4490
4491          _dlnLdlz[branchIndex]   = _dlnLdlz[branchIndex]   + dlnLdlz;
4492          _d2lnLdlz2[branchIndex] = _d2lnLdlz2[branchIndex] + d2lnLdlz2;
4493        }
4494
4495      columnCounter += width;     
4496      offsetCounter += width * tr->partitionData[model].states * tr->discreteRateCategories;
4497    }
4498
4499#ifdef _DEBUG_MULTI_EPA
4500  if(tr->threadID == THREAD_TO_DEBUG)
4501    printf("\n");
4502#endif
4503
4504}
4505
4506void makenewzClassify(tree *tr, int _maxiter, double *result, double *z0, double *x1_start, double *x2_start, unsigned char *tipX1, 
4507                      unsigned char *tipX2, int tipCase, boolean *partitionConverged, int insertion)
4508{
4509  double   
4510    z[NUM_BRANCHES], 
4511    zprev[NUM_BRANCHES], 
4512    zstep[NUM_BRANCHES],
4513    coreLZ[NUM_BRANCHES];
4514  volatile double 
4515    dlnLdlz[NUM_BRANCHES], 
4516    d2lnLdlz2[NUM_BRANCHES];
4517 
4518  int i, maxiter[NUM_BRANCHES], numBranches, model;
4519 
4520  boolean firstIteration = TRUE;
4521  boolean outerConverged[NUM_BRANCHES];
4522  boolean loopConverged,
4523    curvatOK[NUM_BRANCHES],
4524    executeModel[NUM_BRANCHES];
4525
4526 
4527
4528  if(tr->multiBranch)
4529    numBranches = tr->NumberOfModels;
4530  else
4531    numBranches = 1;
4532
4533 
4534
4535  for(i = 0; i < numBranches; i++)
4536    {
4537      z[i] = z0[i];
4538      maxiter[i] = _maxiter;
4539      outerConverged[i] = FALSE;
4540      curvatOK[i]       = TRUE;
4541      if(partitionConverged[i])
4542        executeModel[i] = FALSE;
4543      else
4544        executeModel[i] = TRUE;
4545    }
4546 
4547  if(tr->perPartitionEPA)   
4548    {
4549      setPartitionMask(tr, insertion, executeModel);
4550
4551      if(tr->multiBranch)
4552        for(i = 0; i < numBranches; i++)
4553          if(!executeModel[i])
4554            outerConverged[i] = TRUE;
4555    }
4556
4557     
4558 
4559  do
4560    {
4561      for(i = 0; i < numBranches; i++)
4562        {
4563          if(outerConverged[i] == FALSE && curvatOK[i] == TRUE)
4564            {
4565              curvatOK[i] = FALSE;
4566
4567              zprev[i] = z[i];
4568
4569              zstep[i] = (1.0 - zmax) * z[i] + zmin;
4570            }
4571        }
4572
4573      for(i = 0; i < numBranches; i++)
4574        {
4575          if(outerConverged[i] == FALSE && curvatOK[i] == FALSE)
4576            {
4577              double lz;
4578
4579              if (z[i] < zmin) z[i] = zmin;
4580              else if (z[i] > zmax) z[i] = zmax;
4581              lz    = log(z[i]);
4582
4583              coreLZ[i] = lz;
4584            }
4585        }
4586
4587     
4588      if(tr->multiBranch)
4589        {
4590          for(model = 0; model < tr->NumberOfModels; model++)
4591            {
4592              if(executeModel[model])
4593                executeModel[model] = !curvatOK[model];       
4594            }
4595        }
4596      else
4597        { 
4598         
4599          for(model = 0; model < tr->NumberOfModels; model++)
4600            {
4601              if(tr->perPartitionEPA)
4602                {               
4603                  if(executeModel[model])
4604                    {
4605                      executeModel[model] = !curvatOK[0];                                                   
4606                    }
4607                  /*else
4608                    executeModel[model] = FALSE;*/
4609                }
4610              else
4611                executeModel[model] = !curvatOK[0];
4612            }
4613         
4614        }     
4615
4616     
4617      if(firstIteration)
4618        {         
4619          sumClassify(tr, tipCase, x1_start, x2_start, tipX1, tipX2, executeModel);       
4620          firstIteration = FALSE;
4621        }       
4622       
4623      coreClassify(tr, dlnLdlz, d2lnLdlz2, coreLZ, executeModel);       
4624
4625      for(i = 0; i < numBranches; i++)
4626        {
4627          if(outerConverged[i] == FALSE && curvatOK[i] == FALSE)
4628            {       
4629              if ((d2lnLdlz2[i] >= 0.0) && (z[i] < zmax))
4630                zprev[i] = z[i] = 0.37 * z[i] + 0.63;  /*  Bad curvature, shorten branch */
4631              else
4632                curvatOK[i] = TRUE;
4633            }
4634        }
4635
4636      for(i = 0; i < numBranches; i++)
4637        {
4638          if(curvatOK[i] == TRUE && outerConverged[i] == FALSE)
4639            {
4640              if (d2lnLdlz2[i] < 0.0)
4641                {
4642                  double tantmp = -dlnLdlz[i] / d2lnLdlz2[i];
4643                  if (tantmp < 100)
4644                    {
4645                      z[i] *= EXP(tantmp);
4646                      if (z[i] < zmin)
4647                        z[i] = zmin;
4648
4649                      if (z[i] > 0.25 * zprev[i] + 0.75)
4650                        z[i] = 0.25 * zprev[i] + 0.75;
4651                    }
4652                  else
4653                    z[i] = 0.25 * zprev[i] + 0.75;
4654                }
4655              if (z[i] > zmax) z[i] = zmax;
4656
4657              maxiter[i] = maxiter[i] - 1;
4658              if(maxiter[i] > 0 && (ABS(z[i] - zprev[i]) > zstep[i]))
4659                outerConverged[i] = FALSE;
4660              else
4661                outerConverged[i] = TRUE;
4662            }
4663        }
4664
4665      loopConverged = TRUE;
4666      for(i = 0; i < numBranches; i++)
4667        loopConverged = loopConverged && outerConverged[i];
4668    }
4669  while (!loopConverged);   
4670 
4671  for(i = 0; i < numBranches; i++)   
4672    result[i] = z[i]; 
4673}
4674
4675#endif
4676
4677void makenewzGeneric(tree *tr, nodeptr p, nodeptr q, double *z0, int maxiter, double *result, boolean mask)
4678{
4679  int 
4680    i; 
4681
4682  tr->td[0].ti[0].pNumber = p->number;
4683  tr->td[0].ti[0].qNumber = q->number;
4684
4685  for(i = 0; i < tr->numBranches; i++)
4686    {     
4687      tr->td[0].ti[0].qz[i] =  z0[i];
4688
4689      if(mask)
4690        {           
4691          if(tr->partitionConverged[i])
4692            tr->executeModel[i] = FALSE;
4693          else
4694            tr->executeModel[i] = TRUE;
4695        }
4696    }
4697   
4698  tr->td[0].count = 1;
4699   
4700  if(!p->x)
4701    computeTraversalInfo(p, &(tr->td[0].ti[0]), &(tr->td[0].count), tr->mxtips, tr->numBranches);
4702
4703  if(!q->x)
4704    computeTraversalInfo(q, &(tr->td[0].ti[0]), &(tr->td[0].count), tr->mxtips, tr->numBranches);
4705     
4706  topLevelMakenewz(tr, z0, maxiter, result);
4707 
4708  for(i = 0; i < tr->numBranches; i++)
4709      tr->executeModel[i] = TRUE;
4710}
4711
4712
4713
4714void makenewzGenericDistance(tree *tr, int maxiter, double *z0, double *result, int taxon1, int taxon2)
4715{
4716  int 
4717    i;
4718
4719  assert(taxon1 != taxon2);
4720  assert(0 < taxon1 && taxon1 <= tr->mxtips);
4721  assert(0 < taxon2 && taxon2 <= tr->mxtips);
4722
4723  tr->td[0].ti[0].pNumber = taxon1;
4724  tr->td[0].ti[0].qNumber = taxon2;
4725  tr->td[0].ti[0].tipCase = TIP_TIP;
4726
4727  for(i = 0; i < tr->numBranches; i++)
4728    tr->td[0].ti[0].qz[i] =  defaultz;      /* TODO why defaultz ? */
4729
4730  tr->td[0].count = 1;
4731
4732  topLevelMakenewz(tr, z0, maxiter, result);
4733}
4734
4735
4736
Note: See TracBrowser for help on using the repository browser.