Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

avxLikelihood.c

Visit:

Last change on this file was 13845, checked in by westram, 11 years ago
remove executable flag
File size: 104.3 KB

Line
1	#include <unistd.h>
2
3	#include <math.h>
4	#include <time.h>
5	#include <stdlib.h>
6	#include <stdio.h>
7	#include <ctype.h>
8	#include <string.h>
9	#include <stdint.h>
10	#include <limits.h>
11	#include "axml.h"
12	#include <stdint.h>
13	#include <xmmintrin.h>
14	#include <pmmintrin.h>
15	#include <immintrin.h>
16
17	#ifdef _FMA
18	#include <x86intrin.h>
19	#define FMAMACC(a,b,c) _mm256_macc_pd(b,c,a)
20	#endif
21
22	extern const unsigned int mask32[32];
23
24	const union __attribute__ ((aligned (BYTE_ALIGNMENT)))
25	{
26	uint64_t i[4];
27	__m256d m;
28
29	} absMask_AVX = {{0x7fffffffffffffffULL, 0x7fffffffffffffffULL, 0x7fffffffffffffffULL, 0x7fffffffffffffffULL}};
30
31
32
33	static inline __m256d hadd4(__m256d v, __m256d u)
34	{
35	__m256d
36	a, b;
37
38	v = _mm256_hadd_pd(v, v);
39	a = _mm256_permute2f128_pd(v, v, 1);
40	v = _mm256_add_pd(a, v);
41
42	u = _mm256_hadd_pd(u, u);
43	b = _mm256_permute2f128_pd(u, u, 1);
44	u = _mm256_add_pd(b, u);
45
46	v = _mm256_mul_pd(v, u);
47
48	return v;
49	}
50
51	static inline __m256d hadd3(__m256d v)
52	{
53	__m256d
54	a;
55
56	v = _mm256_hadd_pd(v, v);
57	a = _mm256_permute2f128_pd(v, v, 1);
58	v = _mm256_add_pd(a, v);
59
60	return v;
61	}
62
63
64	void newviewGTRGAMMA_AVX(int tipCase,
65	double x1, double x2, double *x3,
66	double extEV, double tipVector,
67	int ex3, unsigned char tipX1, unsigned char *tipX2,
68	const int n, double left, double right, int wgt, int scalerIncrement, const boolean useFastScaling
69	)
70	{
71
72	int
73	i,
74	k,
75	scale,
76	addScale = 0;
77
78	__m256d
79	minlikelihood_avx = _mm256_set1_pd( minlikelihood ),
80	twoto = _mm256_set1_pd(twotothe256);
81
82
83	switch(tipCase)
84	{
85	case TIP_TIP:
86	{
87	double
88	*uX1,
89	umpX1[1024] __attribute__ ((aligned (BYTE_ALIGNMENT))),
90	*uX2,
91	umpX2[1024] __attribute__ ((aligned (BYTE_ALIGNMENT)));
92
93	for (i = 1; i < 16; i++)
94	{
95	__m256d
96	tv = _mm256_load_pd(&(tipVector[i * 4]));
97
98	int
99	j;
100
101	for (j = 0; j < 4; j++)
102	for (k = 0; k < 4; k++)
103	{
104	__m256d
105	left1 = _mm256_load_pd(&left[j * 16 + k * 4]);
106
107	left1 = _mm256_mul_pd(left1, tv);
108	left1 = hadd3(left1);
109
110	_mm256_store_pd(&umpX1[i * 64 + j * 16 + k * 4], left1);
111	}
112
113	for (j = 0; j < 4; j++)
114	for (k = 0; k < 4; k++)
115	{
116	__m256d
117	left1 = _mm256_load_pd(&right[j * 16 + k * 4]);
118
119	left1 = _mm256_mul_pd(left1, tv);
120	left1 = hadd3(left1);
121
122	_mm256_store_pd(&umpX2[i * 64 + j * 16 + k * 4], left1);
123	}
124	}
125
126
127	for(i = 0; i < n; i++)
128	{
129	uX1 = &umpX1[64 * tipX1[i]];
130	uX2 = &umpX2[64 * tipX2[i]];
131
132	for(k = 0; k < 4; k++)
133	{
134	__m256d
135	xv = _mm256_setzero_pd();
136
137	int
138	l;
139
140	for(l = 0; l < 4; l++)
141	{
142	__m256d
143	x1v = _mm256_mul_pd(_mm256_load_pd(&uX1[k * 16 + l * 4]), _mm256_load_pd(&uX2[k * 16 + l * 4]));
144
145	__m256d
146	evv = _mm256_load_pd(&extEV[l * 4]);
147	#ifdef _FMA
148	xv = FMAMACC(xv,x1v,evv);
149	#else
150	xv = _mm256_add_pd(xv, _mm256_mul_pd(x1v, evv));
151	#endif
152	}
153
154	_mm256_store_pd(&x3[16 * i + 4 * k], xv);
155	}
156	}
157	}
158	break;
159	case TIP_INNER:
160	{
161	double
162	*uX1,
163	umpX1[1024] __attribute__ ((aligned (BYTE_ALIGNMENT)));
164
165	for (i = 1; i < 16; i++)
166	{
167	__m256d
168	tv = _mm256_load_pd(&(tipVector[i*4]));
169
170	int
171	j;
172
173	for (j = 0; j < 4; j++)
174	for (k = 0; k < 4; k++)
175	{
176	__m256d
177	left1 = _mm256_load_pd(&left[j * 16 + k * 4]);
178
179	left1 = _mm256_mul_pd(left1, tv);
180	left1 = hadd3(left1);
181
182	_mm256_store_pd(&umpX1[i * 64 + j * 16 + k * 4], left1);
183	}
184	}
185
186	for(i = 0; i < n; i++)
187	{
188	__m256d
189	xv[4];
190
191	scale = 1;
192	uX1 = &umpX1[64 * tipX1[i]];
193
194	for(k = 0; k < 4; k++)
195	{
196	__m256d
197	xvr = _mm256_load_pd(&(x2[i * 16 + k * 4]));
198
199	int
200	l;
201
202	xv[k] = _mm256_setzero_pd();
203
204	for(l = 0; l < 4; l++)
205	{
206	__m256d
207	x1v = _mm256_load_pd(&uX1[k * 16 + l * 4]),
208	x2v = _mm256_mul_pd(xvr, _mm256_load_pd(&right[k * 16 + l * 4]));
209
210	x2v = hadd3(x2v);
211	x1v = _mm256_mul_pd(x1v, x2v);
212
213	__m256d
214	evv = _mm256_load_pd(&extEV[l * 4]);
215
216	#ifdef _FMA
217	xv[k] = FMAMACC(xv[k],x1v,evv);
218	#else
219	xv[k] = _mm256_add_pd(xv[k], _mm256_mul_pd(x1v, evv));
220	#endif
221	}
222
223	if(scale)
224	{
225	__m256d
226	v1 = _mm256_and_pd(xv[k], absMask_AVX.m);
227
228	v1 = _mm256_cmp_pd(v1, minlikelihood_avx, _CMP_LT_OS);
229
230	if(_mm256_movemask_pd( v1 ) != 15)
231	scale = 0;
232	}
233	}
234
235	if(scale)
236	{
237	xv[0] = _mm256_mul_pd(xv[0], twoto);
238	xv[1] = _mm256_mul_pd(xv[1], twoto);
239	xv[2] = _mm256_mul_pd(xv[2], twoto);
240	xv[3] = _mm256_mul_pd(xv[3], twoto);
241
242	if(useFastScaling)
243	addScale += wgt[i];
244	else
245	ex3[i] += 1;
246	}
247
248	_mm256_store_pd(&x3[16 * i], xv[0]);
249	_mm256_store_pd(&x3[16 * i + 4], xv[1]);
250	_mm256_store_pd(&x3[16 * i + 8], xv[2]);
251	_mm256_store_pd(&x3[16 * i + 12], xv[3]);
252	}
253	}
254	break;
255	case INNER_INNER:
256	{
257	for(i = 0; i < n; i++)
258	{
259	__m256d
260	xv[4];
261
262	scale = 1;
263
264	for(k = 0; k < 4; k++)
265	{
266	__m256d
267
268	xvl = _mm256_load_pd(&(x1[i * 16 + k * 4])),
269	xvr = _mm256_load_pd(&(x2[i * 16 + k * 4]));
270
271	int
272	l;
273
274	xv[k] = _mm256_setzero_pd();
275
276	for(l = 0; l < 4; l++)
277	{
278	__m256d
279	x1v = _mm256_mul_pd(xvl, _mm256_load_pd(&left[k * 16 + l * 4])),
280	x2v = _mm256_mul_pd(xvr, _mm256_load_pd(&right[k * 16 + l * 4]));
281
282	x1v = hadd4(x1v, x2v);
283
284	__m256d
285	evv = _mm256_load_pd(&extEV[l * 4]);
286
287	xv[k] = _mm256_add_pd(xv[k], _mm256_mul_pd(x1v, evv));
288	}
289
290	if(scale)
291	{
292	__m256d
293	v1 = _mm256_and_pd(xv[k], absMask_AVX.m);
294
295	v1 = _mm256_cmp_pd(v1, minlikelihood_avx, _CMP_LT_OS);
296
297	if(_mm256_movemask_pd( v1 ) != 15)
298	scale = 0;
299	}
300	}
301
302	if(scale)
303	{
304	xv[0] = _mm256_mul_pd(xv[0], twoto);
305	xv[1] = _mm256_mul_pd(xv[1], twoto);
306	xv[2] = _mm256_mul_pd(xv[2], twoto);
307	xv[3] = _mm256_mul_pd(xv[3], twoto);
308
309	if(useFastScaling)
310	addScale += wgt[i];
311	else
312	ex3[i] += 1;
313	}
314
315	_mm256_store_pd(&x3[16 * i], xv[0]);
316	_mm256_store_pd(&x3[16 * i + 4], xv[1]);
317	_mm256_store_pd(&x3[16 * i + 8], xv[2]);
318	_mm256_store_pd(&x3[16 * i + 12], xv[3]);
319	}
320	}
321	break;
322	default:
323	assert(0);
324	}
325
326	if(useFastScaling)
327	*scalerIncrement = addScale;
328
329	}
330
331	void newviewGTRGAMMA_AVX_GAPPED_SAVE(int tipCase,
332	double x1_start, double x2_start, double *x3_start,
333	double extEV, double tipVector,
334	int ex3, unsigned char tipX1, unsigned char *tipX2,
335	const int n, double left, double right, int wgt, int scalerIncrement, const boolean useFastScaling,
336	unsigned int x1_gap, unsigned int x2_gap, unsigned int *x3_gap,
337	double x1_gapColumn, double x2_gapColumn, double *x3_gapColumn
338	)
339	{
340
341	int
342	i,
343	k,
344	scale,
345	scaleGap,
346	addScale = 0;
347
348	__m256d
349	minlikelihood_avx = _mm256_set1_pd( minlikelihood ),
350	twoto = _mm256_set1_pd(twotothe256);
351
352	double
353	*x1,
354	*x2,
355	*x3,
356	*x1_ptr = x1_start,
357	*x2_ptr = x2_start;
358
359	switch(tipCase)
360	{
361	case TIP_TIP:
362	{
363	double
364	*uX1,
365	umpX1[1024] __attribute__ ((aligned (BYTE_ALIGNMENT))),
366	*uX2,
367	umpX2[1024] __attribute__ ((aligned (BYTE_ALIGNMENT)));
368
369	for (i = 1; i < 16; i++)
370	{
371	__m256d
372	tv = _mm256_load_pd(&(tipVector[i * 4]));
373
374	int
375	j;
376
377	for (j = 0; j < 4; j++)
378	for (k = 0; k < 4; k++)
379	{
380	__m256d
381	left1 = _mm256_load_pd(&left[j * 16 + k * 4]);
382
383	left1 = _mm256_mul_pd(left1, tv);
384	left1 = hadd3(left1);
385
386	_mm256_store_pd(&umpX1[i * 64 + j * 16 + k * 4], left1);
387	}
388
389	for (j = 0; j < 4; j++)
390	for (k = 0; k < 4; k++)
391	{
392	__m256d
393	left1 = _mm256_load_pd(&right[j * 16 + k * 4]);
394
395	left1 = _mm256_mul_pd(left1, tv);
396	left1 = hadd3(left1);
397
398	_mm256_store_pd(&umpX2[i * 64 + j * 16 + k * 4], left1);
399	}
400	}
401
402	x3 = x3_gapColumn;
403
404	{
405	uX1 = &umpX1[960];
406	uX2 = &umpX2[960];
407
408	for(k = 0; k < 4; k++)
409	{
410	__m256d
411	xv = _mm256_setzero_pd();
412
413	int
414	l;
415
416	for(l = 0; l < 4; l++)
417	{
418	__m256d
419	x1v = _mm256_mul_pd(_mm256_load_pd(&uX1[k * 16 + l * 4]), _mm256_load_pd(&uX2[k * 16 + l * 4]));
420
421	__m256d
422	evv = _mm256_load_pd(&extEV[l * 4]);
423	#ifdef _FMA
424	xv = FMAMACC(xv,x1v,evv);
425	#else
426	xv = _mm256_add_pd(xv, _mm256_mul_pd(x1v, evv));
427	#endif
428	}
429
430	_mm256_store_pd(&x3[4 * k], xv);
431	}
432	}
433
434	x3 = x3_start;
435
436	for(i = 0; i < n; i++)
437	{
438	if(!(x3_gap[i / 32] & mask32[i % 32]))
439	{
440	uX1 = &umpX1[64 * tipX1[i]];
441	uX2 = &umpX2[64 * tipX2[i]];
442
443	for(k = 0; k < 4; k++)
444	{
445	__m256d
446	xv = _mm256_setzero_pd();
447
448	int
449	l;
450
451	for(l = 0; l < 4; l++)
452	{
453	__m256d
454	x1v = _mm256_mul_pd(_mm256_load_pd(&uX1[k * 16 + l * 4]), _mm256_load_pd(&uX2[k * 16 + l * 4]));
455
456	__m256d
457	evv = _mm256_load_pd(&extEV[l * 4]);
458	#ifdef _FMA
459	xv = FMAMACC(xv,x1v,evv);
460	#else
461	xv = _mm256_add_pd(xv, _mm256_mul_pd(x1v, evv));
462	#endif
463	}
464
465	_mm256_store_pd(&x3[4 * k], xv);
466	}
467
468	x3 += 16;
469	}
470	}
471	}
472	break;
473	case TIP_INNER:
474	{
475	double
476	*uX1,
477	umpX1[1024] __attribute__ ((aligned (BYTE_ALIGNMENT)));
478
479	for (i = 1; i < 16; i++)
480	{
481	__m256d
482	tv = _mm256_load_pd(&(tipVector[i*4]));
483
484	int
485	j;
486
487	for (j = 0; j < 4; j++)
488	for (k = 0; k < 4; k++)
489	{
490	__m256d
491	left1 = _mm256_load_pd(&left[j * 16 + k * 4]);
492
493	left1 = _mm256_mul_pd(left1, tv);
494	left1 = hadd3(left1);
495
496	_mm256_store_pd(&umpX1[i * 64 + j * 16 + k * 4], left1);
497	}
498	}
499
500	{
501	__m256d
502	xv[4];
503
504	scaleGap = 1;
505	uX1 = &umpX1[960];
506
507	x2 = x2_gapColumn;
508	x3 = x3_gapColumn;
509
510	for(k = 0; k < 4; k++)
511	{
512	__m256d
513	xvr = _mm256_load_pd(&(x2[k * 4]));
514
515	int
516	l;
517
518	xv[k] = _mm256_setzero_pd();
519
520	for(l = 0; l < 4; l++)
521	{
522	__m256d
523	x1v = _mm256_load_pd(&uX1[k * 16 + l * 4]),
524	x2v = _mm256_mul_pd(xvr, _mm256_load_pd(&right[k * 16 + l * 4]));
525
526	x2v = hadd3(x2v);
527	x1v = _mm256_mul_pd(x1v, x2v);
528
529	__m256d
530	evv = _mm256_load_pd(&extEV[l * 4]);
531
532	#ifdef _FMA
533	xv[k] = FMAMACC(xv[k],x1v,evv);
534	#else
535	xv[k] = _mm256_add_pd(xv[k], _mm256_mul_pd(x1v, evv));
536	#endif
537	}
538
539	if(scaleGap)
540	{
541	__m256d
542	v1 = _mm256_and_pd(xv[k], absMask_AVX.m);
543
544	v1 = _mm256_cmp_pd(v1, minlikelihood_avx, _CMP_LT_OS);
545
546	if(_mm256_movemask_pd( v1 ) != 15)
547	scaleGap = 0;
548	}
549	}
550
551	if(scaleGap)
552	{
553	xv[0] = _mm256_mul_pd(xv[0], twoto);
554	xv[1] = _mm256_mul_pd(xv[1], twoto);
555	xv[2] = _mm256_mul_pd(xv[2], twoto);
556	xv[3] = _mm256_mul_pd(xv[3], twoto);
557	}
558
559	_mm256_store_pd(&x3[0], xv[0]);
560	_mm256_store_pd(&x3[4], xv[1]);
561	_mm256_store_pd(&x3[8], xv[2]);
562	_mm256_store_pd(&x3[12], xv[3]);
563	}
564
565	x3 = x3_start;
566
567	for(i = 0; i < n; i++)
568	{
569	if((x3_gap[i / 32] & mask32[i % 32]))
570	{
571	if(scaleGap)
572	{
573	if(useFastScaling)
574	addScale += wgt[i];
575	else
576	ex3[i] += 1;
577	}
578	}
579	else
580	{
581	if(x2_gap[i / 32] & mask32[i % 32])
582	x2 = x2_gapColumn;
583	else
584	{
585	x2 = x2_ptr;
586	x2_ptr += 16;
587	}
588
589	__m256d
590	xv[4];
591
592	scale = 1;
593	uX1 = &umpX1[64 * tipX1[i]];
594
595	for(k = 0; k < 4; k++)
596	{
597	__m256d
598	xvr = _mm256_load_pd(&(x2[k * 4]));
599
600	int
601	l;
602
603	xv[k] = _mm256_setzero_pd();
604
605	for(l = 0; l < 4; l++)
606	{
607	__m256d
608	x1v = _mm256_load_pd(&uX1[k * 16 + l * 4]),
609	x2v = _mm256_mul_pd(xvr, _mm256_load_pd(&right[k * 16 + l * 4]));
610
611	x2v = hadd3(x2v);
612	x1v = _mm256_mul_pd(x1v, x2v);
613
614	__m256d
615	evv = _mm256_load_pd(&extEV[l * 4]);
616
617	#ifdef _FMA
618	xv[k] = FMAMACC(xv[k],x1v,evv);
619	#else
620	xv[k] = _mm256_add_pd(xv[k], _mm256_mul_pd(x1v, evv));
621	#endif
622	}
623
624	if(scale)
625	{
626	__m256d
627	v1 = _mm256_and_pd(xv[k], absMask_AVX.m);
628
629	v1 = _mm256_cmp_pd(v1, minlikelihood_avx, _CMP_LT_OS);
630
631	if(_mm256_movemask_pd( v1 ) != 15)
632	scale = 0;
633	}
634	}
635
636	if(scale)
637	{
638	xv[0] = _mm256_mul_pd(xv[0], twoto);
639	xv[1] = _mm256_mul_pd(xv[1], twoto);
640	xv[2] = _mm256_mul_pd(xv[2], twoto);
641	xv[3] = _mm256_mul_pd(xv[3], twoto);
642
643	if(useFastScaling)
644	addScale += wgt[i];
645	else
646	ex3[i] += 1;
647	}
648
649	_mm256_store_pd(&x3[0], xv[0]);
650	_mm256_store_pd(&x3[4], xv[1]);
651	_mm256_store_pd(&x3[8], xv[2]);
652	_mm256_store_pd(&x3[12], xv[3]);
653
654	x3 += 16;
655	}
656	}
657	}
658	break;
659	case INNER_INNER:
660	{
661	{
662	x1 = x1_gapColumn;
663	x2 = x2_gapColumn;
664	x3 = x3_gapColumn;
665
666	__m256d
667	xv[4];
668
669	scaleGap = 1;
670
671	for(k = 0; k < 4; k++)
672	{
673	__m256d
674
675	xvl = _mm256_load_pd(&(x1[k * 4])),
676	xvr = _mm256_load_pd(&(x2[k * 4]));
677
678	int
679	l;
680
681	xv[k] = _mm256_setzero_pd();
682
683	for(l = 0; l < 4; l++)
684	{
685	__m256d
686	x1v = _mm256_mul_pd(xvl, _mm256_load_pd(&left[k * 16 + l * 4])),
687	x2v = _mm256_mul_pd(xvr, _mm256_load_pd(&right[k * 16 + l * 4]));
688
689	x1v = hadd4(x1v, x2v);
690
691	__m256d
692	evv = _mm256_load_pd(&extEV[l * 4]);
693
694	xv[k] = _mm256_add_pd(xv[k], _mm256_mul_pd(x1v, evv));
695	}
696
697	if(scaleGap)
698	{
699	__m256d
700	v1 = _mm256_and_pd(xv[k], absMask_AVX.m);
701
702	v1 = _mm256_cmp_pd(v1, minlikelihood_avx, _CMP_LT_OS);
703
704	if(_mm256_movemask_pd( v1 ) != 15)
705	scaleGap = 0;
706	}
707	}
708
709	if(scaleGap)
710	{
711	xv[0] = _mm256_mul_pd(xv[0], twoto);
712	xv[1] = _mm256_mul_pd(xv[1], twoto);
713	xv[2] = _mm256_mul_pd(xv[2], twoto);
714	xv[3] = _mm256_mul_pd(xv[3], twoto);
715	}
716
717	_mm256_store_pd(&x3[0], xv[0]);
718	_mm256_store_pd(&x3[4], xv[1]);
719	_mm256_store_pd(&x3[8], xv[2]);
720	_mm256_store_pd(&x3[12], xv[3]);
721	}
722
723	x3 = x3_start;
724
725	for(i = 0; i < n; i++)
726	{
727	if(x3_gap[i / 32] & mask32[i % 32])
728	{
729	if(scaleGap)
730	{
731	if(useFastScaling)
732	addScale += wgt[i];
733	else
734	ex3[i] += 1;
735	}
736	}
737	else
738	{
739	if(x1_gap[i / 32] & mask32[i % 32])
740	x1 = x1_gapColumn;
741	else
742	{
743	x1 = x1_ptr;
744	x1_ptr += 16;
745	}
746
747	if(x2_gap[i / 32] & mask32[i % 32])
748	x2 = x2_gapColumn;
749	else
750	{
751	x2 = x2_ptr;
752	x2_ptr += 16;
753	}
754
755	__m256d
756	xv[4];
757
758	scale = 1;
759
760	for(k = 0; k < 4; k++)
761	{
762	__m256d
763
764	xvl = _mm256_load_pd(&(x1[k * 4])),
765	xvr = _mm256_load_pd(&(x2[k * 4]));
766
767	int
768	l;
769
770	xv[k] = _mm256_setzero_pd();
771
772	for(l = 0; l < 4; l++)
773	{
774	__m256d
775	x1v = _mm256_mul_pd(xvl, _mm256_load_pd(&left[k * 16 + l * 4])),
776	x2v = _mm256_mul_pd(xvr, _mm256_load_pd(&right[k * 16 + l * 4]));
777
778	x1v = hadd4(x1v, x2v);
779
780	__m256d
781	evv = _mm256_load_pd(&extEV[l * 4]);
782
783	xv[k] = _mm256_add_pd(xv[k], _mm256_mul_pd(x1v, evv));
784	}
785
786	if(scale)
787	{
788	__m256d
789	v1 = _mm256_and_pd(xv[k], absMask_AVX.m);
790
791	v1 = _mm256_cmp_pd(v1, minlikelihood_avx, _CMP_LT_OS);
792
793	if(_mm256_movemask_pd( v1 ) != 15)
794	scale = 0;
795	}
796	}
797
798	if(scale)
799	{
800	xv[0] = _mm256_mul_pd(xv[0], twoto);
801	xv[1] = _mm256_mul_pd(xv[1], twoto);
802	xv[2] = _mm256_mul_pd(xv[2], twoto);
803	xv[3] = _mm256_mul_pd(xv[3], twoto);
804
805	if(useFastScaling)
806	addScale += wgt[i];
807	else
808	ex3[i] += 1;
809	}
810
811	_mm256_store_pd(&x3[0], xv[0]);
812	_mm256_store_pd(&x3[4], xv[1]);
813	_mm256_store_pd(&x3[8], xv[2]);
814	_mm256_store_pd(&x3[12], xv[3]);
815
816	x3 += 16;
817	}
818	}
819	}
820	break;
821	default:
822	assert(0);
823	}
824
825	if(useFastScaling)
826	*scalerIncrement = addScale;
827
828	}
829
830
831
832
833	void newviewGTRCAT_AVX(int tipCase, double EV, int cptr,
834	double x1_start, double x2_start, double x3_start, double tipVector,
835	int ex3, unsigned char tipX1, unsigned char *tipX2,
836	int n, double left, double right, int wgt, int scalerIncrement, const boolean useFastScaling)
837	{
838	double
839	*le,
840	*ri,
841	*x1,
842	*x2;
843
844	int
845	i,
846	addScale = 0;
847
848	__m256d
849	minlikelihood_avx = _mm256_set1_pd( minlikelihood ),
850	twoto = _mm256_set1_pd(twotothe256);
851
852	switch(tipCase)
853	{
854	case TIP_TIP:
855	for (i = 0; i < n; i++)
856	{
857	int
858	l;
859
860	le = &left[cptr[i] * 16];
861	ri = &right[cptr[i] * 16];
862
863	x1 = &(tipVector[4 * tipX1[i]]);
864	x2 = &(tipVector[4 * tipX2[i]]);
865
866	__m256d
867	vv = _mm256_setzero_pd();
868
869	for(l = 0; l < 4; l++)
870	{
871	__m256d
872	x1v = _mm256_mul_pd(_mm256_load_pd(x1), _mm256_load_pd(&le[l * 4])),
873	x2v = _mm256_mul_pd(_mm256_load_pd(x2), _mm256_load_pd(&ri[l * 4]));
874
875	x1v = hadd4(x1v, x2v);
876
877	__m256d
878	evv = _mm256_load_pd(&EV[l * 4]);
879	#ifdef _FMA
880	vv = FMAMACC(vv,x1v,evv);
881	#else
882	vv = _mm256_add_pd(vv, _mm256_mul_pd(x1v, evv));
883	#endif
884	}
885
886	_mm256_store_pd(&x3_start[4 * i], vv);
887	}
888	break;
889	case TIP_INNER:
890	for (i = 0; i < n; i++)
891	{
892	int
893	l;
894
895	x1 = &(tipVector[4 * tipX1[i]]);
896	x2 = &x2_start[4 * i];
897
898	le = &left[cptr[i] * 16];
899	ri = &right[cptr[i] * 16];
900
901	__m256d
902	vv = _mm256_setzero_pd();
903
904	for(l = 0; l < 4; l++)
905	{
906	__m256d
907	x1v = _mm256_mul_pd(_mm256_load_pd(x1), _mm256_load_pd(&le[l * 4])),
908	x2v = _mm256_mul_pd(_mm256_load_pd(x2), _mm256_load_pd(&ri[l * 4]));
909
910	x1v = hadd4(x1v, x2v);
911
912	__m256d
913	evv = _mm256_load_pd(&EV[l * 4]);
914
915	#ifdef _FMA
916	vv = FMAMACC(vv,x1v,evv);
917	#else
918	vv = _mm256_add_pd(vv, _mm256_mul_pd(x1v, evv));
919	#endif
920	}
921
922
923	__m256d
924	v1 = _mm256_and_pd(vv, absMask_AVX.m);
925
926	v1 = _mm256_cmp_pd(v1, minlikelihood_avx, _CMP_LT_OS);
927
928	if(_mm256_movemask_pd( v1 ) == 15)
929	{
930	vv = _mm256_mul_pd(vv, twoto);
931
932	if(useFastScaling)
933	addScale += wgt[i];
934	else
935	ex3[i] += 1;
936	}
937
938	_mm256_store_pd(&x3_start[4 * i], vv);
939	}
940	break;
941	case INNER_INNER:
942	for (i = 0; i < n; i++)
943	{
944	int
945	l;
946
947	x1 = &x1_start[4 * i];
948	x2 = &x2_start[4 * i];
949
950
951	le = &left[cptr[i] * 16];
952	ri = &right[cptr[i] * 16];
953
954	__m256d
955	vv = _mm256_setzero_pd();
956
957	for(l = 0; l < 4; l++)
958	{
959	__m256d
960	x1v = _mm256_mul_pd(_mm256_load_pd(x1), _mm256_load_pd(&le[l * 4])),
961	x2v = _mm256_mul_pd(_mm256_load_pd(x2), _mm256_load_pd(&ri[l * 4]));
962
963	x1v = hadd4(x1v, x2v);
964
965	__m256d
966	evv = _mm256_load_pd(&EV[l * 4]);
967	#ifdef _FMA
968	vv = FMAMACC(vv,x1v,evv);
969	#else
970	vv = _mm256_add_pd(vv, _mm256_mul_pd(x1v, evv));
971	#endif
972	}
973
974
975	__m256d
976	v1 = _mm256_and_pd(vv, absMask_AVX.m);
977
978	v1 = _mm256_cmp_pd(v1, minlikelihood_avx, _CMP_LT_OS);
979
980	if(_mm256_movemask_pd( v1 ) == 15)
981	{
982	vv = _mm256_mul_pd(vv, twoto);
983
984	if(useFastScaling)
985	addScale += wgt[i];
986	else
987	ex3[i] += 1;
988	}
989
990	_mm256_store_pd(&x3_start[4 * i], vv);
991
992	}
993	break;
994	default:
995	assert(0);
996	}
997
998	if(useFastScaling)
999	*scalerIncrement = addScale;
1000	}
1001
1002
1003	void newviewGTRCAT_AVX_GAPPED_SAVE(int tipCase, double EV, int cptr,
1004	double x1_start, double x2_start, double x3_start, double tipVector,
1005	int ex3, unsigned char tipX1, unsigned char *tipX2,
1006	int n, double left, double right, int wgt, int scalerIncrement, const boolean useFastScaling,
1007	unsigned int x1_gap, unsigned int x2_gap, unsigned int *x3_gap,
1008	double x1_gapColumn, double x2_gapColumn, double *x3_gapColumn, const int maxCats)
1009	{
1010	double
1011	*le,
1012	*ri,
1013	*x1,
1014	*x2,
1015	*x3,
1016	*x1_ptr = x1_start,
1017	*x2_ptr = x2_start,
1018	*x3_ptr = x3_start;
1019
1020	int
1021	i,
1022	scaleGap = 0,
1023	addScale = 0;
1024
1025	__m256d
1026	minlikelihood_avx = _mm256_set1_pd( minlikelihood ),
1027	twoto = _mm256_set1_pd(twotothe256);
1028
1029
1030	{
1031	int
1032	l;
1033
1034	x1 = x1_gapColumn;
1035	x2 = x2_gapColumn;
1036	x3 = x3_gapColumn;
1037
1038	le = &left[maxCats * 16];
1039	ri = &right[maxCats * 16];
1040
1041	__m256d
1042	vv = _mm256_setzero_pd();
1043
1044	for(l = 0; l < 4; l++)
1045	{
1046	__m256d
1047	x1v = _mm256_mul_pd(_mm256_load_pd(x1), _mm256_load_pd(&le[l * 4])),
1048	x2v = _mm256_mul_pd(_mm256_load_pd(x2), _mm256_load_pd(&ri[l * 4]));
1049
1050	x1v = hadd4(x1v, x2v);
1051
1052	__m256d
1053	evv = _mm256_load_pd(&EV[l * 4]);
1054	#ifdef _FMA
1055	vv = FMAMACC(vv,x1v,evv);
1056	#else
1057	vv = _mm256_add_pd(vv, _mm256_mul_pd(x1v, evv));
1058	#endif
1059	}
1060
1061	if(tipCase != TIP_TIP)
1062	{
1063	__m256d
1064	v1 = _mm256_and_pd(vv, absMask_AVX.m);
1065
1066	v1 = _mm256_cmp_pd(v1, minlikelihood_avx, _CMP_LT_OS);
1067
1068	if(_mm256_movemask_pd( v1 ) == 15)
1069	{
1070	vv = _mm256_mul_pd(vv, twoto);
1071	scaleGap = 1;
1072	}
1073	}
1074
1075	_mm256_store_pd(x3, vv);
1076	}
1077
1078	switch(tipCase)
1079	{
1080	case TIP_TIP:
1081	for (i = 0; i < n; i++)
1082	{
1083	if(noGap(x3_gap, i))
1084	{
1085	int
1086	l;
1087
1088	x1 = &(tipVector[4 * tipX1[i]]);
1089	x2 = &(tipVector[4 * tipX2[i]]);
1090
1091	x3 = x3_ptr;
1092
1093	if(isGap(x1_gap, i))
1094	le = &left[maxCats * 16];
1095	else
1096	le = &left[cptr[i] * 16];
1097
1098	if(isGap(x2_gap, i))
1099	ri = &right[maxCats * 16];
1100	else
1101	ri = &right[cptr[i] * 16];
1102
1103	__m256d
1104	vv = _mm256_setzero_pd();
1105
1106	for(l = 0; l < 4; l++)
1107	{
1108	__m256d
1109	x1v = _mm256_mul_pd(_mm256_load_pd(x1), _mm256_load_pd(&le[l * 4])),
1110	x2v = _mm256_mul_pd(_mm256_load_pd(x2), _mm256_load_pd(&ri[l * 4]));
1111
1112	x1v = hadd4(x1v, x2v);
1113
1114	__m256d
1115	evv = _mm256_load_pd(&EV[l * 4]);
1116	#ifdef _FMA
1117	vv = FMAMACC(vv,x1v,evv);
1118	#else
1119	vv = _mm256_add_pd(vv, _mm256_mul_pd(x1v, evv));
1120	#endif
1121	}
1122
1123	_mm256_store_pd(x3, vv);
1124
1125	x3_ptr += 4;
1126	}
1127	}
1128	break;
1129	case TIP_INNER:
1130	for (i = 0; i < n; i++)
1131	{
1132	if(isGap(x3_gap, i))
1133	{
1134	if(scaleGap)
1135	{
1136	if(useFastScaling)
1137	addScale += wgt[i];
1138	else
1139	ex3[i] += 1;
1140	}
1141	}
1142	else
1143	{
1144	int
1145	l;
1146
1147	x1 = &(tipVector[4 * tipX1[i]]);
1148	x3 = x3_ptr;
1149
1150	if(isGap(x1_gap, i))
1151	le = &left[maxCats * 16];
1152	else
1153	le = &left[cptr[i] * 16];
1154
1155	if(isGap(x2_gap, i))
1156	{
1157	ri = &right[maxCats * 16];
1158	x2 = x2_gapColumn;
1159	}
1160	else
1161	{
1162	ri = &right[cptr[i] * 16];
1163	x2 = x2_ptr;
1164	x2_ptr += 4;
1165	}
1166
1167	__m256d
1168	vv = _mm256_setzero_pd();
1169
1170	for(l = 0; l < 4; l++)
1171	{
1172	__m256d
1173	x1v = _mm256_mul_pd(_mm256_load_pd(x1), _mm256_load_pd(&le[l * 4])),
1174	x2v = _mm256_mul_pd(_mm256_load_pd(x2), _mm256_load_pd(&ri[l * 4]));
1175
1176	x1v = hadd4(x1v, x2v);
1177
1178	__m256d
1179	evv = _mm256_load_pd(&EV[l * 4]);
1180
1181	#ifdef _FMA
1182	vv = FMAMACC(vv,x1v,evv);
1183	#else
1184	vv = _mm256_add_pd(vv, _mm256_mul_pd(x1v, evv));
1185	#endif
1186	}
1187
1188
1189	__m256d
1190	v1 = _mm256_and_pd(vv, absMask_AVX.m);
1191
1192	v1 = _mm256_cmp_pd(v1, minlikelihood_avx, _CMP_LT_OS);
1193
1194	if(_mm256_movemask_pd( v1 ) == 15)
1195	{
1196	vv = _mm256_mul_pd(vv, twoto);
1197
1198	if(useFastScaling)
1199	addScale += wgt[i];
1200	else
1201	ex3[i] += 1;
1202	}
1203
1204	_mm256_store_pd(x3, vv);
1205
1206	x3_ptr += 4;
1207	}
1208	}
1209	break;
1210	case INNER_INNER:
1211	for (i = 0; i < n; i++)
1212	{
1213	if(isGap(x3_gap, i))
1214	{
1215	if(scaleGap)
1216	{
1217	if(useFastScaling)
1218	addScale += wgt[i];
1219	else
1220	ex3[i] += 1;
1221	}
1222	}
1223	else
1224	{
1225	int
1226	l;
1227
1228	x3 = x3_ptr;
1229
1230	if(isGap(x1_gap, i))
1231	{
1232	x1 = x1_gapColumn;
1233	le = &left[maxCats * 16];
1234	}
1235	else
1236	{
1237	le = &left[cptr[i] * 16];
1238	x1 = x1_ptr;
1239	x1_ptr += 4;
1240	}
1241
1242	if(isGap(x2_gap, i))
1243	{
1244	x2 = x2_gapColumn;
1245	ri = &right[maxCats * 16];
1246	}
1247	else
1248	{
1249	ri = &right[cptr[i] * 16];
1250	x2 = x2_ptr;
1251	x2_ptr += 4;
1252	}
1253
1254	__m256d
1255	vv = _mm256_setzero_pd();
1256
1257	for(l = 0; l < 4; l++)
1258	{
1259	__m256d
1260	x1v = _mm256_mul_pd(_mm256_load_pd(x1), _mm256_load_pd(&le[l * 4])),
1261	x2v = _mm256_mul_pd(_mm256_load_pd(x2), _mm256_load_pd(&ri[l * 4]));
1262
1263	x1v = hadd4(x1v, x2v);
1264
1265	__m256d
1266	evv = _mm256_load_pd(&EV[l * 4]);
1267	#ifdef _FMA
1268	vv = FMAMACC(vv,x1v,evv);
1269	#else
1270	vv = _mm256_add_pd(vv, _mm256_mul_pd(x1v, evv));
1271	#endif
1272	}
1273
1274
1275	__m256d
1276	v1 = _mm256_and_pd(vv, absMask_AVX.m);
1277
1278	v1 = _mm256_cmp_pd(v1, minlikelihood_avx, _CMP_LT_OS);
1279
1280	if(_mm256_movemask_pd( v1 ) == 15)
1281	{
1282	vv = _mm256_mul_pd(vv, twoto);
1283
1284	if(useFastScaling)
1285	addScale += wgt[i];
1286	else
1287	ex3[i] += 1;
1288	}
1289
1290	_mm256_store_pd(x3, vv);
1291
1292	x3_ptr += 4;
1293	}
1294	}
1295	break;
1296	default:
1297	assert(0);
1298	}
1299
1300	if(useFastScaling)
1301	*scalerIncrement = addScale;
1302	}
1303
1304	void newviewGTRCATPROT_AVX(int tipCase, double *extEV,
1305	int *cptr,
1306	double x1, double x2, double x3, double tipVector,
1307	int ex3, unsigned char tipX1, unsigned char *tipX2,
1308	int n, double left, double right, int wgt, int scalerIncrement, const boolean useFastScaling)
1309	{
1310	double
1311	le, ri, v, vl, *vr;
1312
1313	int i, l, scale, addScale = 0;
1314
1315	#ifdef _FMA
1316	int k;
1317	#endif
1318
1319	switch(tipCase)
1320	{
1321	case TIP_TIP:
1322	{
1323	for (i = 0; i < n; i++)
1324	{
1325	le = &left[cptr[i] * 400];
1326	ri = &right[cptr[i] * 400];
1327
1328	vl = &(tipVector[20 * tipX1[i]]);
1329	vr = &(tipVector[20 * tipX2[i]]);
1330	v = &x3[20 * i];
1331
1332	__m256d vv[5];
1333
1334	vv[0] = _mm256_setzero_pd();
1335	vv[1] = _mm256_setzero_pd();
1336	vv[2] = _mm256_setzero_pd();
1337	vv[3] = _mm256_setzero_pd();
1338	vv[4] = _mm256_setzero_pd();
1339
1340	for(l = 0; l < 20; l++)
1341	{
1342	__m256d
1343	x1v = _mm256_setzero_pd(),
1344	x2v = _mm256_setzero_pd();
1345
1346	double
1347	ev = &extEV[l 20],
1348	lv = &le[l 20],
1349	rv = &ri[l 20];
1350
1351	#ifdef _FMA
1352	for(k = 0; k < 20; k += 4)
1353	{
1354	__m256d vlv = _mm256_load_pd(&vl[k]);
1355	__m256d lvv = _mm256_load_pd(&lv[k]);
1356	x1v = FMAMACC(x1v,vlv,lvv);
1357	__m256d vrv = _mm256_load_pd(&vr[k]);
1358	__m256d rvv = _mm256_load_pd(&rv[k]);
1359	x2v = FMAMACC(x2v,vrv,rvv);
1360	}
1361	#else
1362	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[0]), _mm256_load_pd(&lv[0])));
1363	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[4]), _mm256_load_pd(&lv[4])));
1364	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[8]), _mm256_load_pd(&lv[8])));
1365	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[12]), _mm256_load_pd(&lv[12])));
1366	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[16]), _mm256_load_pd(&lv[16])));
1367
1368	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[0]), _mm256_load_pd(&rv[0])));
1369	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[4]), _mm256_load_pd(&rv[4])));
1370	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[8]), _mm256_load_pd(&rv[8])));
1371	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[12]), _mm256_load_pd(&rv[12])));
1372	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[16]), _mm256_load_pd(&rv[16])));
1373	#endif
1374
1375	x1v = hadd4(x1v, x2v);
1376	#ifdef _FMA
1377	for(k = 0; k < 5; k++)
1378	{
1379	__m256d evv = _mm256_load_pd(&ev[k*4]);
1380	vv[k] = FMAMACC(vv[k],x1v,evv);
1381	}
1382	#else
1383	__m256d
1384	evv[5];
1385
1386	evv[0] = _mm256_load_pd(&ev[0]);
1387	evv[1] = _mm256_load_pd(&ev[4]);
1388	evv[2] = _mm256_load_pd(&ev[8]);
1389	evv[3] = _mm256_load_pd(&ev[12]);
1390	evv[4] = _mm256_load_pd(&ev[16]);
1391
1392	vv[0] = _mm256_add_pd(vv[0], _mm256_mul_pd(x1v, evv[0]));
1393	vv[1] = _mm256_add_pd(vv[1], _mm256_mul_pd(x1v, evv[1]));
1394	vv[2] = _mm256_add_pd(vv[2], _mm256_mul_pd(x1v, evv[2]));
1395	vv[3] = _mm256_add_pd(vv[3], _mm256_mul_pd(x1v, evv[3]));
1396	vv[4] = _mm256_add_pd(vv[4], _mm256_mul_pd(x1v, evv[4]));
1397	#endif
1398	}
1399	_mm256_store_pd(&v[0], vv[0]);
1400	_mm256_store_pd(&v[4], vv[1]);
1401	_mm256_store_pd(&v[8], vv[2]);
1402	_mm256_store_pd(&v[12], vv[3]);
1403	_mm256_store_pd(&v[16], vv[4]);
1404	}
1405	}
1406	break;
1407	case TIP_INNER:
1408	for (i = 0; i < n; i++)
1409	{
1410	le = &left[cptr[i] * 400];
1411	ri = &right[cptr[i] * 400];
1412
1413	vl = &(tipVector[20 * tipX1[i]]);
1414	vr = &x2[20 * i];
1415	v = &x3[20 * i];
1416
1417	__m256d vv[5];
1418
1419	vv[0] = _mm256_setzero_pd();
1420	vv[1] = _mm256_setzero_pd();
1421	vv[2] = _mm256_setzero_pd();
1422	vv[3] = _mm256_setzero_pd();
1423	vv[4] = _mm256_setzero_pd();
1424
1425
1426
1427	for(l = 0; l < 20; l++)
1428	{
1429	__m256d
1430	x1v = _mm256_setzero_pd(),
1431	x2v = _mm256_setzero_pd();
1432
1433	double
1434	ev = &extEV[l 20],
1435	lv = &le[l 20],
1436	rv = &ri[l 20];
1437	#ifdef _FMA
1438	for(k = 0; k < 20; k += 4)
1439	{
1440	__m256d vlv = _mm256_load_pd(&vl[k]);
1441	__m256d lvv = _mm256_load_pd(&lv[k]);
1442	x1v = FMAMACC(x1v,vlv,lvv);
1443	__m256d vrv = _mm256_load_pd(&vr[k]);
1444	__m256d rvv = _mm256_load_pd(&rv[k]);
1445	x2v = FMAMACC(x2v,vrv,rvv);
1446	}
1447	#else
1448	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[0]), _mm256_load_pd(&lv[0])));
1449	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[4]), _mm256_load_pd(&lv[4])));
1450	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[8]), _mm256_load_pd(&lv[8])));
1451	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[12]), _mm256_load_pd(&lv[12])));
1452	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[16]), _mm256_load_pd(&lv[16])));
1453
1454	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[0]), _mm256_load_pd(&rv[0])));
1455	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[4]), _mm256_load_pd(&rv[4])));
1456	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[8]), _mm256_load_pd(&rv[8])));
1457	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[12]), _mm256_load_pd(&rv[12])));
1458	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[16]), _mm256_load_pd(&rv[16])));
1459	#endif
1460
1461	x1v = hadd4(x1v, x2v);
1462
1463	__m256d
1464	evv[5];
1465
1466	evv[0] = _mm256_load_pd(&ev[0]);
1467	evv[1] = _mm256_load_pd(&ev[4]);
1468	evv[2] = _mm256_load_pd(&ev[8]);
1469	evv[3] = _mm256_load_pd(&ev[12]);
1470	evv[4] = _mm256_load_pd(&ev[16]);
1471
1472	#ifdef _FMA
1473	for(k = 0; k < 5; k++)
1474	vv[k] = FMAMACC(vv[k],x1v,evv[k]);
1475	#else
1476	vv[0] = _mm256_add_pd(vv[0], _mm256_mul_pd(x1v, evv[0]));
1477	vv[1] = _mm256_add_pd(vv[1], _mm256_mul_pd(x1v, evv[1]));
1478	vv[2] = _mm256_add_pd(vv[2], _mm256_mul_pd(x1v, evv[2]));
1479	vv[3] = _mm256_add_pd(vv[3], _mm256_mul_pd(x1v, evv[3]));
1480	vv[4] = _mm256_add_pd(vv[4], _mm256_mul_pd(x1v, evv[4]));
1481	#endif
1482	}
1483
1484
1485	__m256d minlikelihood_avx = _mm256_set1_pd( minlikelihood );
1486
1487	scale = 1;
1488
1489	for(l = 0; scale && (l < 20); l += 4)
1490	{
1491	__m256d
1492	v1 = _mm256_and_pd(vv[l / 4], absMask_AVX.m);
1493	v1 = _mm256_cmp_pd(v1, minlikelihood_avx, _CMP_LT_OS);
1494
1495	if(_mm256_movemask_pd( v1 ) != 15)
1496	scale = 0;
1497	}
1498
1499
1500	if(scale)
1501	{
1502	__m256d
1503	twoto = _mm256_set1_pd(twotothe256);
1504
1505	for(l = 0; l < 20; l += 4)
1506	vv[l / 4] = _mm256_mul_pd(vv[l / 4] , twoto);
1507
1508	if(useFastScaling)
1509	addScale += wgt[i];
1510	else
1511	ex3[i] += 1;
1512	}
1513
1514	_mm256_store_pd(&v[0], vv[0]);
1515	_mm256_store_pd(&v[4], vv[1]);
1516	_mm256_store_pd(&v[8], vv[2]);
1517	_mm256_store_pd(&v[12], vv[3]);
1518	_mm256_store_pd(&v[16], vv[4]);
1519	}
1520	break;
1521	case INNER_INNER:
1522	for(i = 0; i < n; i++)
1523	{
1524	le = &left[cptr[i] * 400];
1525	ri = &right[cptr[i] * 400];
1526
1527	vl = &x1[20 * i];
1528	vr = &x2[20 * i];
1529	v = &x3[20 * i];
1530
1531	__m256d vv[5];
1532
1533	vv[0] = _mm256_setzero_pd();
1534	vv[1] = _mm256_setzero_pd();
1535	vv[2] = _mm256_setzero_pd();
1536	vv[3] = _mm256_setzero_pd();
1537	vv[4] = _mm256_setzero_pd();
1538
1539	for(l = 0; l < 20; l++)
1540	{
1541	__m256d
1542	x1v = _mm256_setzero_pd(),
1543	x2v = _mm256_setzero_pd();
1544
1545	double
1546	ev = &extEV[l 20],
1547	lv = &le[l 20],
1548	rv = &ri[l 20];
1549
1550	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[0]), _mm256_load_pd(&lv[0])));
1551	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[4]), _mm256_load_pd(&lv[4])));
1552	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[8]), _mm256_load_pd(&lv[8])));
1553	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[12]), _mm256_load_pd(&lv[12])));
1554	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[16]), _mm256_load_pd(&lv[16])));
1555
1556	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[0]), _mm256_load_pd(&rv[0])));
1557	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[4]), _mm256_load_pd(&rv[4])));
1558	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[8]), _mm256_load_pd(&rv[8])));
1559	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[12]), _mm256_load_pd(&rv[12])));
1560	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[16]), _mm256_load_pd(&rv[16])));
1561
1562	x1v = hadd4(x1v, x2v);
1563	#ifdef _FMA
1564	for(k = 0; k < 5; k++)
1565	{
1566	__m256d evv = _mm256_load_pd(&ev[k*4]);
1567	vv[k] = FMAMACC(vv[k],x1v,evv);
1568	}
1569	#else
1570	__m256d
1571	evv[5];
1572
1573	evv[0] = _mm256_load_pd(&ev[0]);
1574	evv[1] = _mm256_load_pd(&ev[4]);
1575	evv[2] = _mm256_load_pd(&ev[8]);
1576	evv[3] = _mm256_load_pd(&ev[12]);
1577	evv[4] = _mm256_load_pd(&ev[16]);
1578
1579	vv[0] = _mm256_add_pd(vv[0], _mm256_mul_pd(x1v, evv[0]));
1580	vv[1] = _mm256_add_pd(vv[1], _mm256_mul_pd(x1v, evv[1]));
1581	vv[2] = _mm256_add_pd(vv[2], _mm256_mul_pd(x1v, evv[2]));
1582	vv[3] = _mm256_add_pd(vv[3], _mm256_mul_pd(x1v, evv[3]));
1583	vv[4] = _mm256_add_pd(vv[4], _mm256_mul_pd(x1v, evv[4]));
1584	#endif
1585	}
1586
1587
1588	__m256d minlikelihood_avx = _mm256_set1_pd( minlikelihood );
1589
1590	scale = 1;
1591
1592	for(l = 0; scale && (l < 20); l += 4)
1593	{
1594	__m256d
1595	v1 = _mm256_and_pd(vv[l / 4], absMask_AVX.m);
1596	v1 = _mm256_cmp_pd(v1, minlikelihood_avx, _CMP_LT_OS);
1597
1598	if(_mm256_movemask_pd( v1 ) != 15)
1599	scale = 0;
1600	}
1601
1602	if(scale)
1603	{
1604	__m256d
1605	twoto = _mm256_set1_pd(twotothe256);
1606
1607	for(l = 0; l < 20; l += 4)
1608	vv[l / 4] = _mm256_mul_pd(vv[l / 4] , twoto);
1609
1610	if(useFastScaling)
1611	addScale += wgt[i];
1612	else
1613	ex3[i] += 1;
1614	}
1615
1616	_mm256_store_pd(&v[0], vv[0]);
1617	_mm256_store_pd(&v[4], vv[1]);
1618	_mm256_store_pd(&v[8], vv[2]);
1619	_mm256_store_pd(&v[12], vv[3]);
1620	_mm256_store_pd(&v[16], vv[4]);
1621
1622	}
1623	break;
1624	default:
1625	assert(0);
1626	}
1627
1628	if(useFastScaling)
1629	*scalerIncrement = addScale;
1630	}
1631
1632	void newviewGTRCATPROT_AVX_GAPPED_SAVE(int tipCase, double *extEV,
1633	int *cptr,
1634	double x1, double x2, double x3, double tipVector,
1635	int ex3, unsigned char tipX1, unsigned char *tipX2,
1636	int n, double left, double right, int wgt, int scalerIncrement, const boolean useFastScaling,
1637	unsigned int x1_gap, unsigned int x2_gap, unsigned int *x3_gap,
1638	double x1_gapColumn, double x2_gapColumn, double *x3_gapColumn, const int maxCats)
1639	{
1640	double
1641	*le,
1642	*ri,
1643	*v,
1644	*vl,
1645	*vr,
1646	*x1_ptr = x1,
1647	*x2_ptr = x2,
1648	*x3_ptr = x3;
1649
1650	int
1651	i,
1652	l,
1653	scale,
1654	addScale = 0,
1655	scaleGap = 0;
1656
1657	#ifdef _FMA
1658	int k;
1659	#endif
1660
1661	{
1662	le = &left[maxCats * 400];
1663	ri = &right[maxCats * 400];
1664
1665	vl = x1_gapColumn;
1666	vr = x2_gapColumn;
1667	v = x3_gapColumn;
1668
1669	__m256d vv[5];
1670
1671	vv[0] = _mm256_setzero_pd();
1672	vv[1] = _mm256_setzero_pd();
1673	vv[2] = _mm256_setzero_pd();
1674	vv[3] = _mm256_setzero_pd();
1675	vv[4] = _mm256_setzero_pd();
1676
1677	for(l = 0; l < 20; l++)
1678	{
1679	__m256d
1680	x1v = _mm256_setzero_pd(),
1681	x2v = _mm256_setzero_pd();
1682
1683	double
1684	ev = &extEV[l 20],
1685	lv = &le[l 20],
1686	rv = &ri[l 20];
1687
1688	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[0]), _mm256_load_pd(&lv[0])));
1689	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[4]), _mm256_load_pd(&lv[4])));
1690	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[8]), _mm256_load_pd(&lv[8])));
1691	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[12]), _mm256_load_pd(&lv[12])));
1692	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[16]), _mm256_load_pd(&lv[16])));
1693
1694	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[0]), _mm256_load_pd(&rv[0])));
1695	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[4]), _mm256_load_pd(&rv[4])));
1696	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[8]), _mm256_load_pd(&rv[8])));
1697	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[12]), _mm256_load_pd(&rv[12])));
1698	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[16]), _mm256_load_pd(&rv[16])));
1699
1700	x1v = hadd4(x1v, x2v);
1701	#ifdef _FMA
1702	for(k = 0; k < 5; k++)
1703	{
1704	__m256d evv = _mm256_load_pd(&ev[k*4]);
1705	vv[k] = FMAMACC(vv[k],x1v,evv);
1706	}
1707	#else
1708	__m256d
1709	evv[5];
1710
1711	evv[0] = _mm256_load_pd(&ev[0]);
1712	evv[1] = _mm256_load_pd(&ev[4]);
1713	evv[2] = _mm256_load_pd(&ev[8]);
1714	evv[3] = _mm256_load_pd(&ev[12]);
1715	evv[4] = _mm256_load_pd(&ev[16]);
1716
1717	vv[0] = _mm256_add_pd(vv[0], _mm256_mul_pd(x1v, evv[0]));
1718	vv[1] = _mm256_add_pd(vv[1], _mm256_mul_pd(x1v, evv[1]));
1719	vv[2] = _mm256_add_pd(vv[2], _mm256_mul_pd(x1v, evv[2]));
1720	vv[3] = _mm256_add_pd(vv[3], _mm256_mul_pd(x1v, evv[3]));
1721	vv[4] = _mm256_add_pd(vv[4], _mm256_mul_pd(x1v, evv[4]));
1722	#endif
1723	}
1724
1725
1726	if(tipCase != TIP_TIP)
1727	{
1728	__m256d minlikelihood_avx = _mm256_set1_pd( minlikelihood );
1729
1730	scale = 1;
1731
1732	for(l = 0; scale && (l < 20); l += 4)
1733	{
1734	__m256d
1735	v1 = _mm256_and_pd(vv[l / 4], absMask_AVX.m);
1736	v1 = _mm256_cmp_pd(v1, minlikelihood_avx, _CMP_LT_OS);
1737
1738	if(_mm256_movemask_pd( v1 ) != 15)
1739	scale = 0;
1740	}
1741
1742	if(scale)
1743	{
1744	__m256d
1745	twoto = _mm256_set1_pd(twotothe256);
1746
1747	for(l = 0; l < 20; l += 4)
1748	vv[l / 4] = _mm256_mul_pd(vv[l / 4] , twoto);
1749
1750	scaleGap = 1;
1751	}
1752	}
1753
1754	_mm256_store_pd(&v[0], vv[0]);
1755	_mm256_store_pd(&v[4], vv[1]);
1756	_mm256_store_pd(&v[8], vv[2]);
1757	_mm256_store_pd(&v[12], vv[3]);
1758	_mm256_store_pd(&v[16], vv[4]);
1759	}
1760
1761
1762
1763	switch(tipCase)
1764	{
1765	case TIP_TIP:
1766	{
1767	for (i = 0; i < n; i++)
1768	{
1769	if(noGap(x3_gap, i))
1770	{
1771	vl = &(tipVector[20 * tipX1[i]]);
1772	vr = &(tipVector[20 * tipX2[i]]);
1773	v = x3_ptr;
1774
1775	if(isGap(x1_gap, i))
1776	le = &left[maxCats * 400];
1777	else
1778	le = &left[cptr[i] * 400];
1779
1780	if(isGap(x2_gap, i))
1781	ri = &right[maxCats * 400];
1782	else
1783	ri = &right[cptr[i] * 400];
1784
1785	__m256d vv[5];
1786
1787	vv[0] = _mm256_setzero_pd();
1788	vv[1] = _mm256_setzero_pd();
1789	vv[2] = _mm256_setzero_pd();
1790	vv[3] = _mm256_setzero_pd();
1791	vv[4] = _mm256_setzero_pd();
1792
1793	for(l = 0; l < 20; l++)
1794	{
1795	__m256d
1796	x1v = _mm256_setzero_pd(),
1797	x2v = _mm256_setzero_pd();
1798
1799	double
1800	ev = &extEV[l 20],
1801	lv = &le[l 20],
1802	rv = &ri[l 20];
1803
1804	#ifdef _FMA
1805	for(k = 0; k < 20; k += 4)
1806	{
1807	__m256d vlv = _mm256_load_pd(&vl[k]);
1808	__m256d lvv = _mm256_load_pd(&lv[k]);
1809	x1v = FMAMACC(x1v,vlv,lvv);
1810	__m256d vrv = _mm256_load_pd(&vr[k]);
1811	__m256d rvv = _mm256_load_pd(&rv[k]);
1812	x2v = FMAMACC(x2v,vrv,rvv);
1813	}
1814	#else
1815	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[0]), _mm256_load_pd(&lv[0])));
1816	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[4]), _mm256_load_pd(&lv[4])));
1817	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[8]), _mm256_load_pd(&lv[8])));
1818	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[12]), _mm256_load_pd(&lv[12])));
1819	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[16]), _mm256_load_pd(&lv[16])));
1820
1821	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[0]), _mm256_load_pd(&rv[0])));
1822	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[4]), _mm256_load_pd(&rv[4])));
1823	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[8]), _mm256_load_pd(&rv[8])));
1824	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[12]), _mm256_load_pd(&rv[12])));
1825	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[16]), _mm256_load_pd(&rv[16])));
1826	#endif
1827
1828	x1v = hadd4(x1v, x2v);
1829	#ifdef _FMA
1830	for(k = 0; k < 5; k++)
1831	{
1832	__m256d evv = _mm256_load_pd(&ev[k*4]);
1833	vv[k] = FMAMACC(vv[k],x1v,evv);
1834	}
1835	#else
1836	__m256d
1837	evv[5];
1838
1839	evv[0] = _mm256_load_pd(&ev[0]);
1840	evv[1] = _mm256_load_pd(&ev[4]);
1841	evv[2] = _mm256_load_pd(&ev[8]);
1842	evv[3] = _mm256_load_pd(&ev[12]);
1843	evv[4] = _mm256_load_pd(&ev[16]);
1844
1845	vv[0] = _mm256_add_pd(vv[0], _mm256_mul_pd(x1v, evv[0]));
1846	vv[1] = _mm256_add_pd(vv[1], _mm256_mul_pd(x1v, evv[1]));
1847	vv[2] = _mm256_add_pd(vv[2], _mm256_mul_pd(x1v, evv[2]));
1848	vv[3] = _mm256_add_pd(vv[3], _mm256_mul_pd(x1v, evv[3]));
1849	vv[4] = _mm256_add_pd(vv[4], _mm256_mul_pd(x1v, evv[4]));
1850	#endif
1851	}
1852
1853	_mm256_store_pd(&v[0], vv[0]);
1854	_mm256_store_pd(&v[4], vv[1]);
1855	_mm256_store_pd(&v[8], vv[2]);
1856	_mm256_store_pd(&v[12], vv[3]);
1857	_mm256_store_pd(&v[16], vv[4]);
1858
1859	x3_ptr += 20;
1860	}
1861	}
1862	}
1863	break;
1864	case TIP_INNER:
1865	for (i = 0; i < n; i++)
1866	{
1867	if(isGap(x3_gap, i))
1868	{
1869	if(scaleGap)
1870	{
1871	if(useFastScaling)
1872	addScale += wgt[i];
1873	else
1874	ex3[i] += 1;
1875	}
1876	}
1877	else
1878	{
1879	vl = &(tipVector[20 * tipX1[i]]);
1880
1881	vr = x2_ptr;
1882	v = x3_ptr;
1883
1884	if(isGap(x1_gap, i))
1885	le = &left[maxCats * 400];
1886	else
1887	le = &left[cptr[i] * 400];
1888
1889	if(isGap(x2_gap, i))
1890	{
1891	ri = &right[maxCats * 400];
1892	vr = x2_gapColumn;
1893	}
1894	else
1895	{
1896	ri = &right[cptr[i] * 400];
1897	vr = x2_ptr;
1898	x2_ptr += 20;
1899	}
1900
1901	__m256d vv[5];
1902
1903	vv[0] = _mm256_setzero_pd();
1904	vv[1] = _mm256_setzero_pd();
1905	vv[2] = _mm256_setzero_pd();
1906	vv[3] = _mm256_setzero_pd();
1907	vv[4] = _mm256_setzero_pd();
1908
1909	for(l = 0; l < 20; l++)
1910	{
1911	__m256d
1912	x1v = _mm256_setzero_pd(),
1913	x2v = _mm256_setzero_pd();
1914
1915	double
1916	ev = &extEV[l 20],
1917	lv = &le[l 20],
1918	rv = &ri[l 20];
1919	#ifdef _FMA
1920	for(k = 0; k < 20; k += 4)
1921	{
1922	__m256d vlv = _mm256_load_pd(&vl[k]);
1923	__m256d lvv = _mm256_load_pd(&lv[k]);
1924	x1v = FMAMACC(x1v,vlv,lvv);
1925	__m256d vrv = _mm256_load_pd(&vr[k]);
1926	__m256d rvv = _mm256_load_pd(&rv[k]);
1927	x2v = FMAMACC(x2v,vrv,rvv);
1928	}
1929	#else
1930	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[0]), _mm256_load_pd(&lv[0])));
1931	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[4]), _mm256_load_pd(&lv[4])));
1932	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[8]), _mm256_load_pd(&lv[8])));
1933	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[12]), _mm256_load_pd(&lv[12])));
1934	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[16]), _mm256_load_pd(&lv[16])));
1935
1936	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[0]), _mm256_load_pd(&rv[0])));
1937	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[4]), _mm256_load_pd(&rv[4])));
1938	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[8]), _mm256_load_pd(&rv[8])));
1939	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[12]), _mm256_load_pd(&rv[12])));
1940	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[16]), _mm256_load_pd(&rv[16])));
1941	#endif
1942
1943	x1v = hadd4(x1v, x2v);
1944
1945	__m256d
1946	evv[5];
1947
1948	evv[0] = _mm256_load_pd(&ev[0]);
1949	evv[1] = _mm256_load_pd(&ev[4]);
1950	evv[2] = _mm256_load_pd(&ev[8]);
1951	evv[3] = _mm256_load_pd(&ev[12]);
1952	evv[4] = _mm256_load_pd(&ev[16]);
1953
1954	#ifdef _FMA
1955	for(k = 0; k < 5; k++)
1956	vv[k] = FMAMACC(vv[k],x1v,evv[k]);
1957	#else
1958	vv[0] = _mm256_add_pd(vv[0], _mm256_mul_pd(x1v, evv[0]));
1959	vv[1] = _mm256_add_pd(vv[1], _mm256_mul_pd(x1v, evv[1]));
1960	vv[2] = _mm256_add_pd(vv[2], _mm256_mul_pd(x1v, evv[2]));
1961	vv[3] = _mm256_add_pd(vv[3], _mm256_mul_pd(x1v, evv[3]));
1962	vv[4] = _mm256_add_pd(vv[4], _mm256_mul_pd(x1v, evv[4]));
1963	#endif
1964	}
1965
1966
1967	__m256d minlikelihood_avx = _mm256_set1_pd( minlikelihood );
1968
1969	scale = 1;
1970
1971	for(l = 0; scale && (l < 20); l += 4)
1972	{
1973	__m256d
1974	v1 = _mm256_and_pd(vv[l / 4], absMask_AVX.m);
1975	v1 = _mm256_cmp_pd(v1, minlikelihood_avx, _CMP_LT_OS);
1976
1977	if(_mm256_movemask_pd( v1 ) != 15)
1978	scale = 0;
1979	}
1980
1981	if(scale)
1982	{
1983	__m256d
1984	twoto = _mm256_set1_pd(twotothe256);
1985
1986	for(l = 0; l < 20; l += 4)
1987	vv[l / 4] = _mm256_mul_pd(vv[l / 4] , twoto);
1988
1989	if(useFastScaling)
1990	addScale += wgt[i];
1991	else
1992	ex3[i] += 1;
1993	}
1994
1995	_mm256_store_pd(&v[0], vv[0]);
1996	_mm256_store_pd(&v[4], vv[1]);
1997	_mm256_store_pd(&v[8], vv[2]);
1998	_mm256_store_pd(&v[12], vv[3]);
1999	_mm256_store_pd(&v[16], vv[4]);
2000
2001	x3_ptr += 20;
2002	}
2003	}
2004	break;
2005	case INNER_INNER:
2006	for(i = 0; i < n; i++)
2007	{
2008	if(isGap(x3_gap, i))
2009	{
2010	if(scaleGap)
2011	{
2012	if(useFastScaling)
2013	addScale += wgt[i];
2014	else
2015	ex3[i] += 1;
2016	}
2017	}
2018	else
2019	{
2020
2021	v = x3_ptr;
2022
2023	if(isGap(x1_gap, i))
2024	{
2025	vl = x1_gapColumn;
2026	le = &left[maxCats * 400];
2027	}
2028	else
2029	{
2030	le = &left[cptr[i] * 400];
2031	vl = x1_ptr;
2032	x1_ptr += 20;
2033	}
2034
2035	if(isGap(x2_gap, i))
2036	{
2037	vr = x2_gapColumn;
2038	ri = &right[maxCats * 400];
2039	}
2040	else
2041	{
2042	ri = &right[cptr[i] * 400];
2043	vr = x2_ptr;
2044	x2_ptr += 20;
2045	}
2046
2047	__m256d vv[5];
2048
2049	vv[0] = _mm256_setzero_pd();
2050	vv[1] = _mm256_setzero_pd();
2051	vv[2] = _mm256_setzero_pd();
2052	vv[3] = _mm256_setzero_pd();
2053	vv[4] = _mm256_setzero_pd();
2054
2055	for(l = 0; l < 20; l++)
2056	{
2057	__m256d
2058	x1v = _mm256_setzero_pd(),
2059	x2v = _mm256_setzero_pd();
2060
2061	double
2062	ev = &extEV[l 20],
2063	lv = &le[l 20],
2064	rv = &ri[l 20];
2065
2066	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[0]), _mm256_load_pd(&lv[0])));
2067	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[4]), _mm256_load_pd(&lv[4])));
2068	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[8]), _mm256_load_pd(&lv[8])));
2069	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[12]), _mm256_load_pd(&lv[12])));
2070	x1v = _mm256_add_pd(x1v, _mm256_mul_pd(_mm256_load_pd(&vl[16]), _mm256_load_pd(&lv[16])));
2071
2072	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[0]), _mm256_load_pd(&rv[0])));
2073	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[4]), _mm256_load_pd(&rv[4])));
2074	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[8]), _mm256_load_pd(&rv[8])));
2075	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[12]), _mm256_load_pd(&rv[12])));
2076	x2v = _mm256_add_pd(x2v, _mm256_mul_pd(_mm256_load_pd(&vr[16]), _mm256_load_pd(&rv[16])));
2077
2078	x1v = hadd4(x1v, x2v);
2079	#ifdef _FMA
2080	for(k = 0; k < 5; k++)
2081	{
2082	__m256d evv = _mm256_load_pd(&ev[k*4]);
2083	vv[k] = FMAMACC(vv[k],x1v,evv);
2084	}
2085	#else
2086	__m256d
2087	evv[5];
2088
2089	evv[0] = _mm256_load_pd(&ev[0]);
2090	evv[1] = _mm256_load_pd(&ev[4]);
2091	evv[2] = _mm256_load_pd(&ev[8]);
2092	evv[3] = _mm256_load_pd(&ev[12]);
2093	evv[4] = _mm256_load_pd(&ev[16]);
2094
2095	vv[0] = _mm256_add_pd(vv[0], _mm256_mul_pd(x1v, evv[0]));
2096	vv[1] = _mm256_add_pd(vv[1], _mm256_mul_pd(x1v, evv[1]));
2097	vv[2] = _mm256_add_pd(vv[2], _mm256_mul_pd(x1v, evv[2]));
2098	vv[3] = _mm256_add_pd(vv[3], _mm256_mul_pd(x1v, evv[3]));
2099	vv[4] = _mm256_add_pd(vv[4], _mm256_mul_pd(x1v, evv[4]));
2100	#endif
2101	}
2102
2103
2104	__m256d minlikelihood_avx = _mm256_set1_pd( minlikelihood );
2105
2106	scale = 1;
2107
2108	for(l = 0; scale && (l < 20); l += 4)
2109	{
2110	__m256d
2111	v1 = _mm256_and_pd(vv[l / 4], absMask_AVX.m);
2112	v1 = _mm256_cmp_pd(v1, minlikelihood_avx, _CMP_LT_OS);
2113
2114	if(_mm256_movemask_pd( v1 ) != 15)
2115	scale = 0;
2116	}
2117
2118	if(scale)
2119	{
2120	__m256d
2121	twoto = _mm256_set1_pd(twotothe256);
2122
2123	for(l = 0; l < 20; l += 4)
2124	vv[l / 4] = _mm256_mul_pd(vv[l / 4] , twoto);
2125
2126	if(useFastScaling)
2127	addScale += wgt[i];
2128	else
2129	ex3[i] += 1;
2130	}
2131
2132	_mm256_store_pd(&v[0], vv[0]);
2133	_mm256_store_pd(&v[4], vv[1]);
2134	_mm256_store_pd(&v[8], vv[2]);
2135	_mm256_store_pd(&v[12], vv[3]);
2136	_mm256_store_pd(&v[16], vv[4]);
2137
2138	x3_ptr += 20;
2139	}
2140	}
2141	break;
2142	default:
2143	assert(0);
2144	}
2145
2146	if(useFastScaling)
2147	*scalerIncrement = addScale;
2148	}
2149
2150
2151
2152
2153
2154
2155	void newviewGTRGAMMAPROT_AVX(int tipCase,
2156	double x1, double x2, double x3, double extEV, double *tipVector,
2157	int ex3, unsigned char tipX1, unsigned char *tipX2, int n,
2158	double left, double right, int wgt, int scalerIncrement, const boolean useFastScaling)
2159	{
2160	double
2161	*uX1,
2162	*uX2,
2163	*v,
2164	x1px2,
2165	*vl,
2166	*vr;
2167
2168	int
2169	i,
2170	j,
2171	l,
2172	k,
2173	scale,
2174	addScale = 0;
2175
2176
2177	#ifndef GCC_VERSION
2178	#define GCC_VERSION (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__)
2179	#endif
2180
2181
2182	#if GCC_VERSION < 40500
2183	__m256d
2184	bitmask = _mm256_set_pd(0,0,0,-1);
2185	#else
2186	__m256i
2187	bitmask = _mm256_set_epi32(0, 0, 0, 0, 0, 0, -1, -1);
2188	#endif
2189
2190	switch(tipCase)
2191	{
2192	case TIP_TIP:
2193	{
2194
2195	double
2196	umpX1[1840] __attribute__ ((aligned (BYTE_ALIGNMENT))),
2197	umpX2[1840] __attribute__ ((aligned (BYTE_ALIGNMENT)));
2198
2199	for(i = 0; i < 23; i++)
2200	{
2201	v = &(tipVector[20 * i]);
2202
2203	for(k = 0; k < 80; k++)
2204	{
2205	double
2206	ll = &left[k 20],
2207	rr = &right[k 20];
2208
2209	__m256d
2210	umpX1v = _mm256_setzero_pd(),
2211	umpX2v = _mm256_setzero_pd();
2212
2213	for(l = 0; l < 20; l+=4)
2214	{
2215	__m256d vv = _mm256_load_pd(&v[l]);
2216	#ifdef _FMA
2217	__m256d llv = _mm256_load_pd(&ll[l]);
2218	umpX1v = FMAMACC(umpX1v,vv,llv);
2219	__m256d rrv = _mm256_load_pd(&rr[l]);
2220	umpX2v = FMAMACC(umpX2v,vv,rrv);
2221	#else
2222	umpX1v = _mm256_add_pd(umpX1v,_mm256_mul_pd(vv,_mm256_load_pd(&ll[l])));
2223	umpX2v = _mm256_add_pd(umpX2v,_mm256_mul_pd(vv,_mm256_load_pd(&rr[l])));
2224	#endif
2225	}
2226
2227	umpX1v = hadd3(umpX1v);
2228	umpX2v = hadd3(umpX2v);
2229	_mm256_maskstore_pd(&umpX1[80 * i + k], bitmask, umpX1v);
2230	_mm256_maskstore_pd(&umpX2[80 * i + k], bitmask, umpX2v);
2231	}
2232	}
2233
2234	for(i = 0; i < n; i++)
2235	{
2236	uX1 = &umpX1[80 * tipX1[i]];
2237	uX2 = &umpX2[80 * tipX2[i]];
2238
2239	for(j = 0; j < 4; j++)
2240	{
2241	__m256d vv[5];
2242
2243	v = &x3[i * 80 + j * 20];
2244
2245	vv[0] = _mm256_setzero_pd();
2246	vv[1] = _mm256_setzero_pd();
2247	vv[2] = _mm256_setzero_pd();
2248	vv[3] = _mm256_setzero_pd();
2249	vv[4] = _mm256_setzero_pd();
2250
2251	for(k = 0; k < 20; k++)
2252	{
2253	x1px2 = uX1[j * 20 + k] * uX2[j * 20 + k];
2254
2255	__m256d x1px2v = _mm256_set1_pd(x1px2);
2256
2257	__m256d extEvv = _mm256_load_pd(&extEV[20 * k]);
2258	#ifdef _FMA
2259	vv[0] = FMAMACC(vv[0],x1px2v,extEvv);
2260	#else
2261	vv[0] = _mm256_add_pd(vv[0],_mm256_mul_pd(x1px2v,extEvv));
2262	#endif
2263	_mm256_store_pd(&v[0],vv[0]);
2264
2265	extEvv = _mm256_load_pd(&extEV[20 * k + 4]);
2266	#ifdef _FMA
2267	vv[1] = FMAMACC(vv[1],x1px2v,extEvv);
2268	#else
2269	vv[1] = _mm256_add_pd(vv[1],_mm256_mul_pd(x1px2v,extEvv));
2270	#endif
2271	_mm256_store_pd(&v[4],vv[1]);
2272
2273	extEvv = _mm256_load_pd(&extEV[20 * k + 8]);
2274	#ifdef _FMA
2275	vv[2] = FMAMACC(vv[2],x1px2v,extEvv);
2276	#else
2277	vv[2] = _mm256_add_pd(vv[2],_mm256_mul_pd(x1px2v,extEvv));
2278	#endif
2279	_mm256_store_pd(&v[8],vv[2]);
2280
2281	extEvv = _mm256_load_pd(&extEV[20 * k + 12]);
2282	#ifdef _FMA
2283	vv[3] = FMAMACC(vv[3],x1px2v,extEvv);
2284	#else
2285	vv[3] = _mm256_add_pd(vv[3],_mm256_mul_pd(x1px2v,extEvv));
2286	#endif
2287	_mm256_store_pd(&v[12],vv[3]);
2288
2289	extEvv = _mm256_load_pd(&extEV[20 * k + 16]);
2290	#ifdef _FMA
2291	vv[4] = FMAMACC(vv[4],x1px2v,extEvv);
2292	#else
2293	vv[4] = _mm256_add_pd(vv[4],_mm256_mul_pd(x1px2v,extEvv));
2294	#endif
2295	_mm256_store_pd(&v[16],vv[4]);
2296	}
2297	}
2298	}
2299	}
2300	break;
2301	case TIP_INNER:
2302	{
2303
2304	double
2305	umpX1[1840] __attribute__ ((aligned (BYTE_ALIGNMENT))),
2306	ump_x2[20] __attribute__ ((aligned (BYTE_ALIGNMENT)));
2307
2308	for(i = 0; i < 23; i++)
2309	{
2310	v = &(tipVector[20 * i]);
2311
2312	for(k = 0; k < 80; k++)
2313	{
2314	__m256d umpX1v = _mm256_setzero_pd();
2315	for(l = 0; l < 20; l+=4)
2316	{
2317	__m256d vv = _mm256_load_pd(&v[l]);
2318	__m256d leftv = _mm256_load_pd(&left[k * 20 + l]);
2319	#ifdef _FMA
2320
2321	umpX1v = FMAMACC(umpX1v, vv, leftv);
2322	#else
2323	umpX1v = _mm256_add_pd(umpX1v, _mm256_mul_pd(vv, leftv));
2324	#endif
2325	}
2326	umpX1v = hadd3(umpX1v);
2327	_mm256_maskstore_pd(&umpX1[80 * i + k], bitmask, umpX1v);
2328	}
2329	}
2330
2331	for (i = 0; i < n; i++)
2332	{
2333	uX1 = &umpX1[80 * tipX1[i]];
2334
2335	for(k = 0; k < 4; k++)
2336	{
2337	v = &(x2[80 * i + k * 20]);
2338
2339	for(l = 0; l < 20; l++)
2340	{
2341	__m256d ump_x2v = _mm256_setzero_pd();
2342
2343	__m256d vv = _mm256_load_pd(&v[0]);
2344	__m256d rightv = _mm256_load_pd(&right[k400+l20+0]);
2345	#ifdef _FMA
2346	ump_x2v = FMAMACC(ump_x2v,vv,rightv);
2347	#else
2348	ump_x2v = _mm256_add_pd(ump_x2v, _mm256_mul_pd(vv, rightv));
2349	#endif
2350
2351	vv = _mm256_load_pd(&v[4]);
2352	rightv = _mm256_load_pd(&right[k400+l20+4]);
2353	#ifdef _FMA
2354	ump_x2v = FMAMACC(ump_x2v,vv,rightv);
2355	#else
2356	ump_x2v = _mm256_add_pd(ump_x2v, _mm256_mul_pd(vv, rightv));
2357	#endif
2358
2359	vv = _mm256_load_pd(&v[8]);
2360	rightv = _mm256_load_pd(&right[k400+l20+8]);
2361	#ifdef _FMA
2362	ump_x2v = FMAMACC(ump_x2v,vv,rightv);
2363	#else
2364	ump_x2v = _mm256_add_pd(ump_x2v, _mm256_mul_pd(vv, rightv));
2365	#endif
2366
2367	vv = _mm256_load_pd(&v[12]);
2368	rightv = _mm256_load_pd(&right[k400+l20+12]);
2369	#ifdef _FMA
2370	ump_x2v = FMAMACC(ump_x2v,vv,rightv);
2371	#else
2372	ump_x2v = _mm256_add_pd(ump_x2v, _mm256_mul_pd(vv, rightv));
2373	#endif
2374
2375	vv = _mm256_load_pd(&v[16]);
2376	rightv = _mm256_load_pd(&right[k400+l20+16]);
2377	#ifdef _FMA
2378	ump_x2v = FMAMACC(ump_x2v,vv,rightv);
2379	#else
2380	ump_x2v = _mm256_add_pd(ump_x2v, _mm256_mul_pd(vv, rightv));
2381	#endif
2382
2383	ump_x2v = hadd3(ump_x2v);
2384	_mm256_maskstore_pd(&ump_x2[l], bitmask, ump_x2v);
2385	}
2386
2387	v = &(x3[80 * i + 20 * k]);
2388
2389
2390	__m256d vv[5];
2391
2392	vv[0] = _mm256_setzero_pd();
2393	vv[1] = _mm256_setzero_pd();
2394	vv[2] = _mm256_setzero_pd();
2395	vv[3] = _mm256_setzero_pd();
2396	vv[4] = _mm256_setzero_pd();
2397
2398	for(l = 0; l < 20; l++)
2399	{
2400	x1px2 = uX1[k * 20 + l] * ump_x2[l];
2401	__m256d x1px2v = _mm256_set1_pd(x1px2);
2402
2403	#ifdef _FMA
2404	__m256d ev = _mm256_load_pd(&extEV[l * 20 + 0]);
2405	vv[0] = FMAMACC(vv[0],x1px2v, ev);
2406	#else
2407	vv[0] = _mm256_add_pd(vv[0],_mm256_mul_pd(x1px2v, _mm256_load_pd(&extEV[l * 20 + 0])));
2408	#endif
2409	_mm256_store_pd(&v[0],vv[0]);
2410
2411	#ifdef _FMA
2412	ev = _mm256_load_pd(&extEV[l * 20 + 4]);
2413	vv[1] = FMAMACC(vv[1],x1px2v, ev);
2414	#else
2415	vv[1] = _mm256_add_pd(vv[1],_mm256_mul_pd(x1px2v, _mm256_load_pd(&extEV[l * 20 + 4])));
2416	#endif
2417	_mm256_store_pd(&v[4],vv[1]);
2418
2419	#ifdef _FMA
2420	ev = _mm256_load_pd(&extEV[l * 20 + 8]);
2421	vv[2] = FMAMACC(vv[2],x1px2v, ev);
2422	#else
2423	vv[2] = _mm256_add_pd(vv[2],_mm256_mul_pd(x1px2v, _mm256_load_pd(&extEV[l * 20 + 8])));
2424	#endif
2425	_mm256_store_pd(&v[8],vv[2]);
2426
2427	#ifdef _FMA
2428	ev = _mm256_load_pd(&extEV[l * 20 + 12]);
2429	vv[3] = FMAMACC(vv[3],x1px2v, ev);
2430	#else
2431	vv[3] = _mm256_add_pd(vv[3],_mm256_mul_pd(x1px2v, _mm256_load_pd(&extEV[l * 20 + 12])));
2432	#endif
2433	_mm256_store_pd(&v[12],vv[3]);
2434
2435
2436	#ifdef _FMA
2437	ev = _mm256_load_pd(&extEV[l * 20 + 16]);
2438	vv[4] = FMAMACC(vv[4],x1px2v, ev);
2439	#else
2440	vv[4] = _mm256_add_pd(vv[4],_mm256_mul_pd(x1px2v, _mm256_load_pd(&extEV[l * 20 + 16])));
2441	#endif
2442	_mm256_store_pd(&v[16],vv[4]);
2443
2444	}
2445	}
2446
2447	v = &x3[80 * i];
2448	__m256d minlikelihood_avx = _mm256_set1_pd(minlikelihood);
2449	scale = 1;
2450	for(l = 0; scale && (l < 80); l += 4)
2451	{
2452	__m256d vv = _mm256_load_pd(&v[l]);
2453	__m256d vv_abs = _mm256_and_pd(vv,absMask_AVX.m);
2454	vv_abs = _mm256_cmp_pd(vv_abs,minlikelihood_avx,_CMP_LT_OS);
2455	if(_mm256_movemask_pd(vv_abs) != 15)
2456	scale = 0;
2457	}
2458
2459	if(scale)
2460	{
2461	__m256d twotothe256v = _mm256_set_pd(twotothe256,twotothe256,twotothe256,twotothe256);
2462	for(l = 0; l < 80; l += 4)
2463	{
2464	__m256d vv = _mm256_load_pd(&v[l]);
2465	_mm256_store_pd(&v[l],_mm256_mul_pd(vv,twotothe256v));
2466	}
2467	if(useFastScaling)
2468	addScale += wgt[i];
2469	else
2470	ex3[i] += 1;
2471	}
2472	}
2473	}
2474	break;
2475	case INNER_INNER:
2476	for(i = 0; i < n; i++)
2477	{
2478	scale = 1;
2479
2480	for(k = 0; k < 4; k++)
2481	{
2482	vl = &(x1[80 * i + 20 * k]);
2483	vr = &(x2[80 * i + 20 * k]);
2484	v = &(x3[80 * i + 20 * k]);
2485
2486	__m256d vv[5];
2487
2488	vv[0] = _mm256_setzero_pd();
2489	vv[1] = _mm256_setzero_pd();
2490	vv[2] = _mm256_setzero_pd();
2491	vv[3] = _mm256_setzero_pd();
2492	vv[4] = _mm256_setzero_pd();
2493
2494	for(l = 0; l < 20; l++)
2495	{
2496	__m256d al = _mm256_setzero_pd();
2497	__m256d ar = _mm256_setzero_pd();
2498
2499	__m256d leftv = _mm256_load_pd(&left[k * 400 + l * 20 + 0]);
2500	__m256d rightv = _mm256_load_pd(&right[k * 400 + l * 20 + 0]);
2501	__m256d vlv = _mm256_load_pd(&vl[0]);
2502	__m256d vrv = _mm256_load_pd(&vr[0]);
2503
2504	#ifdef _FMA
2505
2506	al = FMAMACC(al, vlv, leftv);
2507	ar = FMAMACC(ar, vrv, rightv);
2508	#else
2509	al = _mm256_add_pd(al,_mm256_mul_pd(vlv,leftv));
2510	ar = _mm256_add_pd(ar,_mm256_mul_pd(vrv,rightv));
2511	#endif
2512
2513	leftv = _mm256_load_pd(&left[k * 400 + l * 20 + 4]);
2514	rightv = _mm256_load_pd(&right[k * 400 + l * 20 + 4]);
2515	vlv = _mm256_load_pd(&vl[4]);
2516	vrv = _mm256_load_pd(&vr[4]);
2517	#ifdef _FMA
2518
2519	al = FMAMACC(al, vlv, leftv);
2520	ar = FMAMACC(ar, vrv, rightv);
2521	#else
2522	al = _mm256_add_pd(al,_mm256_mul_pd(vlv,leftv));
2523	ar = _mm256_add_pd(ar,_mm256_mul_pd(vrv,rightv));
2524	#endif
2525
2526	leftv = _mm256_load_pd(&left[k * 400 + l * 20 + 8]);
2527	rightv = _mm256_load_pd(&right[k * 400 + l * 20 + 8]);
2528	vlv = _mm256_load_pd(&vl[8]);
2529	vrv = _mm256_load_pd(&vr[8]);
2530	#ifdef _FMA
2531
2532	al = FMAMACC(al, vlv, leftv);
2533	ar = FMAMACC(ar, vrv, rightv);
2534	#else
2535	al = _mm256_add_pd(al,_mm256_mul_pd(vlv,leftv));
2536	ar = _mm256_add_pd(ar,_mm256_mul_pd(vrv,rightv));
2537	#endif
2538
2539	leftv = _mm256_load_pd(&left[k * 400 + l * 20 + 12]);
2540	rightv = _mm256_load_pd(&right[k * 400 + l * 20 + 12]);
2541	vlv = _mm256_load_pd(&vl[12]);
2542	vrv = _mm256_load_pd(&vr[12]);
2543	#ifdef _FMA
2544
2545	al = FMAMACC(al, vlv, leftv);
2546	ar = FMAMACC(ar, vrv, rightv);
2547	#else
2548	al = _mm256_add_pd(al,_mm256_mul_pd(vlv,leftv));
2549	ar = _mm256_add_pd(ar,_mm256_mul_pd(vrv,rightv));
2550	#endif
2551
2552	leftv = _mm256_load_pd(&left[k * 400 + l * 20 + 16]);
2553	rightv = _mm256_load_pd(&right[k * 400 + l * 20 + 16]);
2554	vlv = _mm256_load_pd(&vl[16]);
2555	vrv = _mm256_load_pd(&vr[16]);
2556
2557	#ifdef _FMA
2558	al = FMAMACC(al, vlv, leftv);
2559	ar = FMAMACC(ar, vrv, rightv);
2560	#else
2561	al = _mm256_add_pd(al,_mm256_mul_pd(vlv,leftv));
2562	ar = _mm256_add_pd(ar,_mm256_mul_pd(vrv,rightv));
2563	#endif
2564
2565	/**************************************************************************************************************/
2566
2567	al = hadd3(al);
2568	ar = hadd3(ar);
2569	al = _mm256_mul_pd(ar,al);
2570
2571	/************************************************************************************************************/
2572	#ifdef _FMA
2573	__m256d ev = _mm256_load_pd(&extEV[20 * l + 0]);
2574	vv[0] = FMAMACC(vv[0], al, ev);
2575	#else
2576	vv[0] = _mm256_add_pd(vv[0],_mm256_mul_pd(al, _mm256_load_pd(&extEV[20 * l + 0])));
2577	#endif
2578	_mm256_store_pd(&v[0],vv[0]);
2579
2580	#ifdef _FMA
2581	ev = _mm256_load_pd(&extEV[20 * l + 4]);
2582	vv[1] = FMAMACC(vv[1], al, ev);
2583	#else
2584	vv[1] = _mm256_add_pd(vv[1],_mm256_mul_pd(al, _mm256_load_pd(&extEV[20 * l + 4])));
2585	#endif
2586	_mm256_store_pd(&v[4],vv[1]);
2587
2588	#ifdef _FMA
2589	ev = _mm256_load_pd(&extEV[20 * l + 8]);
2590	vv[2] = FMAMACC(vv[2], al, ev);
2591	#else
2592	vv[2] = _mm256_add_pd(vv[2],_mm256_mul_pd(al, _mm256_load_pd(&extEV[20 * l + 8])));
2593	#endif
2594	_mm256_store_pd(&v[8],vv[2]);
2595
2596	#ifdef _FMA
2597	ev = _mm256_load_pd(&extEV[20 * l + 12]);
2598	vv[3] = FMAMACC(vv[3], al, ev);
2599	#else
2600	vv[3] = _mm256_add_pd(vv[3],_mm256_mul_pd(al, _mm256_load_pd(&extEV[20 * l + 12])));
2601	#endif
2602	_mm256_store_pd(&v[12],vv[3]);
2603
2604	#ifdef _FMA
2605	ev = _mm256_load_pd(&extEV[20 * l + 16]);
2606	vv[4] = FMAMACC(vv[4], al, ev);
2607	#else
2608	vv[4] = _mm256_add_pd(vv[4],_mm256_mul_pd(al, _mm256_load_pd(&extEV[20 * l + 16])));
2609	#endif
2610	_mm256_store_pd(&v[16],vv[4]);
2611	}
2612	}
2613	v = &(x3[80 * i]);
2614	scale = 1;
2615	__m256d minlikelihood_avx = _mm256_set1_pd(minlikelihood);
2616
2617	for(l = 0; scale && (l < 80); l += 4)
2618	{
2619	__m256d vv = _mm256_load_pd(&v[l]);
2620	__m256d vv_abs = _mm256_and_pd(vv,absMask_AVX.m);
2621	vv_abs = _mm256_cmp_pd(vv_abs,minlikelihood_avx,_CMP_LT_OS);
2622	if(_mm256_movemask_pd(vv_abs) != 15)
2623	scale = 0;
2624	}
2625
2626	if(scale)
2627	{
2628	__m256d twotothe256v = _mm256_set_pd(twotothe256,twotothe256,twotothe256,twotothe256);
2629	for(l = 0; l < 80; l += 4)
2630	{
2631	__m256d vv = _mm256_load_pd(&v[l]);
2632	_mm256_store_pd(&v[l],_mm256_mul_pd(vv,twotothe256v));
2633	}
2634	if(useFastScaling)
2635	addScale += wgt[i];
2636	else
2637	ex3[i] += 1;
2638	}
2639	}
2640	break;
2641	default:
2642	assert(0);
2643	}
2644
2645	if(useFastScaling)
2646	*scalerIncrement = addScale;
2647	}
2648
2649	void newviewGTRGAMMAPROT_AVX_LG4(int tipCase,
2650	double x1, double x2, double x3, double extEV[4], double *tipVector[4],
2651	int ex3, unsigned char tipX1, unsigned char *tipX2, int n,
2652	double left, double right, int wgt, int scalerIncrement, const boolean useFastScaling)
2653	{
2654	double
2655	*uX1,
2656	*uX2,
2657	*v,
2658	x1px2,
2659	*vl,
2660	*vr;
2661
2662	int
2663	i,
2664	j,
2665	l,
2666	k,
2667	scale,
2668	addScale = 0;
2669
2670
2671	#ifndef GCC_VERSION
2672	#define GCC_VERSION (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__)
2673	#endif
2674
2675
2676	#if GCC_VERSION < 40500
2677	__m256d
2678	bitmask = _mm256_set_pd(0,0,0,-1);
2679	#else
2680	__m256i
2681	bitmask = _mm256_set_epi32(0, 0, 0, 0, 0, 0, -1, -1);
2682	#endif
2683
2684	switch(tipCase)
2685	{
2686	case TIP_TIP:
2687	{
2688
2689	double
2690	umpX1[1840] __attribute__ ((aligned (BYTE_ALIGNMENT))),
2691	umpX2[1840] __attribute__ ((aligned (BYTE_ALIGNMENT)));
2692
2693
2694	for(i = 0; i < 23; i++)
2695	{
2696	for(k = 0; k < 80; k++)
2697	{
2698	double
2699	ll = &left[k 20],
2700	rr = &right[k 20];
2701
2702	__m256d
2703	umpX1v = _mm256_setzero_pd(),
2704	umpX2v = _mm256_setzero_pd();
2705
2706	v = &(tipVector[k / 20][20 * i]);
2707
2708	for(l = 0; l < 20; l+=4)
2709	{
2710	__m256d vv = _mm256_load_pd(&v[l]);
2711	#ifdef _FMA
2712	__m256d llv = _mm256_load_pd(&ll[l]);
2713	umpX1v = FMAMACC(umpX1v,vv,llv);
2714	__m256d rrv = _mm256_load_pd(&rr[l]);
2715	umpX2v = FMAMACC(umpX2v,vv,rrv);
2716	#else
2717	umpX1v = _mm256_add_pd(umpX1v,_mm256_mul_pd(vv,_mm256_load_pd(&ll[l])));
2718	umpX2v = _mm256_add_pd(umpX2v,_mm256_mul_pd(vv,_mm256_load_pd(&rr[l])));
2719	#endif
2720	}
2721
2722	umpX1v = hadd3(umpX1v);
2723	umpX2v = hadd3(umpX2v);
2724	_mm256_maskstore_pd(&umpX1[80 * i + k], bitmask, umpX1v);
2725	_mm256_maskstore_pd(&umpX2[80 * i + k], bitmask, umpX2v);
2726	}
2727	}
2728
2729	for(i = 0; i < n; i++)
2730	{
2731	uX1 = &umpX1[80 * tipX1[i]];
2732	uX2 = &umpX2[80 * tipX2[i]];
2733
2734	for(j = 0; j < 4; j++)
2735	{
2736	__m256d vv[5];
2737
2738	v = &x3[i * 80 + j * 20];
2739
2740	vv[0] = _mm256_setzero_pd();
2741	vv[1] = _mm256_setzero_pd();
2742	vv[2] = _mm256_setzero_pd();
2743	vv[3] = _mm256_setzero_pd();
2744	vv[4] = _mm256_setzero_pd();
2745
2746	for(k = 0; k < 20; k++)
2747	{
2748	x1px2 = uX1[j * 20 + k] * uX2[j * 20 + k];
2749
2750	__m256d x1px2v = _mm256_set1_pd(x1px2);
2751
2752	__m256d extEvv = _mm256_load_pd(&extEV[j][20 * k]);
2753	#ifdef _FMA
2754	vv[0] = FMAMACC(vv[0],x1px2v,extEvv);
2755	#else
2756	vv[0] = _mm256_add_pd(vv[0],_mm256_mul_pd(x1px2v,extEvv));
2757	#endif
2758	_mm256_store_pd(&v[0],vv[0]);
2759
2760	extEvv = _mm256_load_pd(&extEV[j][20 * k + 4]);
2761	#ifdef _FMA
2762	vv[1] = FMAMACC(vv[1],x1px2v,extEvv);
2763	#else
2764	vv[1] = _mm256_add_pd(vv[1],_mm256_mul_pd(x1px2v,extEvv));
2765	#endif
2766	_mm256_store_pd(&v[4],vv[1]);
2767
2768	extEvv = _mm256_load_pd(&extEV[j][20 * k + 8]);
2769	#ifdef _FMA
2770	vv[2] = FMAMACC(vv[2],x1px2v,extEvv);
2771	#else
2772	vv[2] = _mm256_add_pd(vv[2],_mm256_mul_pd(x1px2v,extEvv));
2773	#endif
2774	_mm256_store_pd(&v[8],vv[2]);
2775
2776	extEvv = _mm256_load_pd(&extEV[j][20 * k + 12]);
2777	#ifdef _FMA
2778	vv[3] = FMAMACC(vv[3],x1px2v,extEvv);
2779	#else
2780	vv[3] = _mm256_add_pd(vv[3],_mm256_mul_pd(x1px2v,extEvv));
2781	#endif
2782	_mm256_store_pd(&v[12],vv[3]);
2783
2784	extEvv = _mm256_load_pd(&extEV[j][20 * k + 16]);
2785	#ifdef _FMA
2786	vv[4] = FMAMACC(vv[4],x1px2v,extEvv);
2787	#else
2788	vv[4] = _mm256_add_pd(vv[4],_mm256_mul_pd(x1px2v,extEvv));
2789	#endif
2790	_mm256_store_pd(&v[16],vv[4]);
2791	}
2792	}
2793	}
2794	}
2795	break;
2796	case TIP_INNER:
2797	{
2798
2799	double
2800	umpX1[1840] __attribute__ ((aligned (BYTE_ALIGNMENT))),
2801	ump_x2[20] __attribute__ ((aligned (BYTE_ALIGNMENT)));
2802
2803	for(i = 0; i < 23; i++)
2804	{
2805	for(k = 0; k < 80; k++)
2806	{
2807	__m256d umpX1v = _mm256_setzero_pd();
2808
2809	v = &(tipVector[k / 20][20 * i]);
2810
2811	for(l = 0; l < 20; l+=4)
2812	{
2813	__m256d vv = _mm256_load_pd(&v[l]);
2814	__m256d leftv = _mm256_load_pd(&left[k * 20 + l]);
2815	#ifdef _FMA
2816
2817	umpX1v = FMAMACC(umpX1v, vv, leftv);
2818	#else
2819	umpX1v = _mm256_add_pd(umpX1v, _mm256_mul_pd(vv, leftv));
2820	#endif
2821	}
2822	umpX1v = hadd3(umpX1v);
2823	_mm256_maskstore_pd(&umpX1[80 * i + k], bitmask, umpX1v);
2824	}
2825	}
2826
2827	for (i = 0; i < n; i++)
2828	{
2829	uX1 = &umpX1[80 * tipX1[i]];
2830
2831	for(k = 0; k < 4; k++)
2832	{
2833	v = &(x2[80 * i + k * 20]);
2834
2835	for(l = 0; l < 20; l++)
2836	{
2837	__m256d ump_x2v = _mm256_setzero_pd();
2838
2839	__m256d vv = _mm256_load_pd(&v[0]);
2840	__m256d rightv = _mm256_load_pd(&right[k400+l20+0]);
2841	#ifdef _FMA
2842	ump_x2v = FMAMACC(ump_x2v,vv,rightv);
2843	#else
2844	ump_x2v = _mm256_add_pd(ump_x2v, _mm256_mul_pd(vv, rightv));
2845	#endif
2846
2847	vv = _mm256_load_pd(&v[4]);
2848	rightv = _mm256_load_pd(&right[k400+l20+4]);
2849	#ifdef _FMA
2850	ump_x2v = FMAMACC(ump_x2v,vv,rightv);
2851	#else
2852	ump_x2v = _mm256_add_pd(ump_x2v, _mm256_mul_pd(vv, rightv));
2853	#endif
2854
2855	vv = _mm256_load_pd(&v[8]);
2856	rightv = _mm256_load_pd(&right[k400+l20+8]);
2857	#ifdef _FMA
2858	ump_x2v = FMAMACC(ump_x2v,vv,rightv);
2859	#else
2860	ump_x2v = _mm256_add_pd(ump_x2v, _mm256_mul_pd(vv, rightv));
2861	#endif
2862
2863	vv = _mm256_load_pd(&v[12]);
2864	rightv = _mm256_load_pd(&right[k400+l20+12]);
2865	#ifdef _FMA
2866	ump_x2v = FMAMACC(ump_x2v,vv,rightv);
2867	#else
2868	ump_x2v = _mm256_add_pd(ump_x2v, _mm256_mul_pd(vv, rightv));
2869	#endif
2870
2871	vv = _mm256_load_pd(&v[16]);
2872	rightv = _mm256_load_pd(&right[k400+l20+16]);
2873	#ifdef _FMA
2874	ump_x2v = FMAMACC(ump_x2v,vv,rightv);
2875	#else
2876	ump_x2v = _mm256_add_pd(ump_x2v, _mm256_mul_pd(vv, rightv));
2877	#endif
2878
2879	ump_x2v = hadd3(ump_x2v);
2880	_mm256_maskstore_pd(&ump_x2[l], bitmask, ump_x2v);
2881	}
2882
2883	v = &(x3[80 * i + 20 * k]);
2884
2885
2886	__m256d vv[5];
2887
2888	vv[0] = _mm256_setzero_pd();
2889	vv[1] = _mm256_setzero_pd();
2890	vv[2] = _mm256_setzero_pd();
2891	vv[3] = _mm256_setzero_pd();
2892	vv[4] = _mm256_setzero_pd();
2893
2894	for(l = 0; l < 20; l++)
2895	{
2896	x1px2 = uX1[k * 20 + l] * ump_x2[l];
2897	__m256d x1px2v = _mm256_set1_pd(x1px2);
2898
2899	#ifdef _FMA
2900	__m256d ev = _mm256_load_pd(&extEV[l * 20 + 0]);
2901	vv[0] = FMAMACC(vv[0],x1px2v, ev);
2902	#else
2903	vv[0] = _mm256_add_pd(vv[0],_mm256_mul_pd(x1px2v, _mm256_load_pd(&extEV[k][l * 20 + 0])));
2904	#endif
2905	_mm256_store_pd(&v[0],vv[0]);
2906
2907	#ifdef _FMA
2908	ev = _mm256_load_pd(&extEV[l * 20 + 4]);
2909	vv[1] = FMAMACC(vv[1],x1px2v, ev);
2910	#else
2911	vv[1] = _mm256_add_pd(vv[1],_mm256_mul_pd(x1px2v, _mm256_load_pd(&extEV[k][l * 20 + 4])));
2912	#endif
2913	_mm256_store_pd(&v[4],vv[1]);
2914
2915	#ifdef _FMA
2916	ev = _mm256_load_pd(&extEV[l * 20 + 8]);
2917	vv[2] = FMAMACC(vv[2],x1px2v, ev);
2918	#else
2919	vv[2] = _mm256_add_pd(vv[2],_mm256_mul_pd(x1px2v, _mm256_load_pd(&extEV[k][l * 20 + 8])));
2920	#endif
2921	_mm256_store_pd(&v[8],vv[2]);
2922
2923	#ifdef _FMA
2924	ev = _mm256_load_pd(&extEV[l * 20 + 12]);
2925	vv[3] = FMAMACC(vv[3],x1px2v, ev);
2926	#else
2927	vv[3] = _mm256_add_pd(vv[3],_mm256_mul_pd(x1px2v, _mm256_load_pd(&extEV[k][l * 20 + 12])));
2928	#endif
2929	_mm256_store_pd(&v[12],vv[3]);
2930
2931
2932	#ifdef _FMA
2933	ev = _mm256_load_pd(&extEV[l * 20 + 16]);
2934	vv[4] = FMAMACC(vv[4],x1px2v, ev);
2935	#else
2936	vv[4] = _mm256_add_pd(vv[4],_mm256_mul_pd(x1px2v, _mm256_load_pd(&extEV[k][l * 20 + 16])));
2937	#endif
2938	_mm256_store_pd(&v[16],vv[4]);
2939
2940	}
2941	}
2942
2943	v = &x3[80 * i];
2944	__m256d minlikelihood_avx = _mm256_set1_pd(minlikelihood);
2945	scale = 1;
2946	for(l = 0; scale && (l < 80); l += 4)
2947	{
2948	__m256d vv = _mm256_load_pd(&v[l]);
2949	__m256d vv_abs = _mm256_and_pd(vv,absMask_AVX.m);
2950	vv_abs = _mm256_cmp_pd(vv_abs,minlikelihood_avx,_CMP_LT_OS);
2951	if(_mm256_movemask_pd(vv_abs) != 15)
2952	scale = 0;
2953	}
2954
2955	if(scale)
2956	{
2957	__m256d twotothe256v = _mm256_set_pd(twotothe256,twotothe256,twotothe256,twotothe256);
2958	for(l = 0; l < 80; l += 4)
2959	{
2960	__m256d vv = _mm256_load_pd(&v[l]);
2961	_mm256_store_pd(&v[l],_mm256_mul_pd(vv,twotothe256v));
2962	}
2963	if(useFastScaling)
2964	addScale += wgt[i];
2965	else
2966	ex3[i] += 1;
2967	}
2968	}
2969	}
2970	break;
2971	case INNER_INNER:
2972	for(i = 0; i < n; i++)
2973	{
2974	scale = 1;
2975
2976	for(k = 0; k < 4; k++)
2977	{
2978	vl = &(x1[80 * i + 20 * k]);
2979	vr = &(x2[80 * i + 20 * k]);
2980	v = &(x3[80 * i + 20 * k]);
2981
2982	__m256d vv[5];
2983
2984	vv[0] = _mm256_setzero_pd();
2985	vv[1] = _mm256_setzero_pd();
2986	vv[2] = _mm256_setzero_pd();
2987	vv[3] = _mm256_setzero_pd();
2988	vv[4] = _mm256_setzero_pd();
2989
2990	for(l = 0; l < 20; l++)
2991	{
2992	__m256d al = _mm256_setzero_pd();
2993	__m256d ar = _mm256_setzero_pd();
2994
2995	__m256d leftv = _mm256_load_pd(&left[k * 400 + l * 20 + 0]);
2996	__m256d rightv = _mm256_load_pd(&right[k * 400 + l * 20 + 0]);
2997	__m256d vlv = _mm256_load_pd(&vl[0]);
2998	__m256d vrv = _mm256_load_pd(&vr[0]);
2999
3000	#ifdef _FMA
3001
3002	al = FMAMACC(al, vlv, leftv);
3003	ar = FMAMACC(ar, vrv, rightv);
3004	#else
3005	al = _mm256_add_pd(al,_mm256_mul_pd(vlv,leftv));
3006	ar = _mm256_add_pd(ar,_mm256_mul_pd(vrv,rightv));
3007	#endif
3008
3009	leftv = _mm256_load_pd(&left[k * 400 + l * 20 + 4]);
3010	rightv = _mm256_load_pd(&right[k * 400 + l * 20 + 4]);
3011	vlv = _mm256_load_pd(&vl[4]);
3012	vrv = _mm256_load_pd(&vr[4]);
3013	#ifdef _FMA
3014
3015	al = FMAMACC(al, vlv, leftv);
3016	ar = FMAMACC(ar, vrv, rightv);
3017	#else
3018	al = _mm256_add_pd(al,_mm256_mul_pd(vlv,leftv));
3019	ar = _mm256_add_pd(ar,_mm256_mul_pd(vrv,rightv));
3020	#endif
3021
3022	leftv = _mm256_load_pd(&left[k * 400 + l * 20 + 8]);
3023	rightv = _mm256_load_pd(&right[k * 400 + l * 20 + 8]);
3024	vlv = _mm256_load_pd(&vl[8]);
3025	vrv = _mm256_load_pd(&vr[8]);
3026	#ifdef _FMA
3027
3028	al = FMAMACC(al, vlv, leftv);
3029	ar = FMAMACC(ar, vrv, rightv);
3030	#else
3031	al = _mm256_add_pd(al,_mm256_mul_pd(vlv,leftv));
3032	ar = _mm256_add_pd(ar,_mm256_mul_pd(vrv,rightv));
3033	#endif
3034
3035	leftv = _mm256_load_pd(&left[k * 400 + l * 20 + 12]);
3036	rightv = _mm256_load_pd(&right[k * 400 + l * 20 + 12]);
3037	vlv = _mm256_load_pd(&vl[12]);
3038	vrv = _mm256_load_pd(&vr[12]);
3039	#ifdef _FMA
3040
3041	al = FMAMACC(al, vlv, leftv);
3042	ar = FMAMACC(ar, vrv, rightv);
3043	#else
3044	al = _mm256_add_pd(al,_mm256_mul_pd(vlv,leftv));
3045	ar = _mm256_add_pd(ar,_mm256_mul_pd(vrv,rightv));
3046	#endif
3047
3048	leftv = _mm256_load_pd(&left[k * 400 + l * 20 + 16]);
3049	rightv = _mm256_load_pd(&right[k * 400 + l * 20 + 16]);
3050	vlv = _mm256_load_pd(&vl[16]);
3051	vrv = _mm256_load_pd(&vr[16]);
3052
3053	#ifdef _FMA
3054	al = FMAMACC(al, vlv, leftv);
3055	ar = FMAMACC(ar, vrv, rightv);
3056	#else
3057	al = _mm256_add_pd(al,_mm256_mul_pd(vlv,leftv));
3058	ar = _mm256_add_pd(ar,_mm256_mul_pd(vrv,rightv));
3059	#endif
3060
3061	/**************************************************************************************************************/
3062
3063	al = hadd3(al);
3064	ar = hadd3(ar);
3065	al = _mm256_mul_pd(ar,al);
3066
3067	/************************************************************************************************************/
3068	#ifdef _FMA
3069	__m256d ev = _mm256_load_pd(&extEV[20 * l + 0]);
3070	vv[0] = FMAMACC(vv[0], al, ev);
3071	#else
3072	vv[0] = _mm256_add_pd(vv[0],_mm256_mul_pd(al, _mm256_load_pd(&extEV[k][20 * l + 0])));
3073	#endif
3074	_mm256_store_pd(&v[0],vv[0]);
3075
3076	#ifdef _FMA
3077	ev = _mm256_load_pd(&extEV[20 * l + 4]);
3078	vv[1] = FMAMACC(vv[1], al, ev);
3079	#else
3080	vv[1] = _mm256_add_pd(vv[1],_mm256_mul_pd(al, _mm256_load_pd(&extEV[k][20 * l + 4])));
3081	#endif
3082	_mm256_store_pd(&v[4],vv[1]);
3083
3084	#ifdef _FMA
3085	ev = _mm256_load_pd(&extEV[20 * l + 8]);
3086	vv[2] = FMAMACC(vv[2], al, ev);
3087	#else
3088	vv[2] = _mm256_add_pd(vv[2],_mm256_mul_pd(al, _mm256_load_pd(&extEV[k][20 * l + 8])));
3089	#endif
3090	_mm256_store_pd(&v[8],vv[2]);
3091
3092	#ifdef _FMA
3093	ev = _mm256_load_pd(&extEV[20 * l + 12]);
3094	vv[3] = FMAMACC(vv[3], al, ev);
3095	#else
3096	vv[3] = _mm256_add_pd(vv[3],_mm256_mul_pd(al, _mm256_load_pd(&extEV[k][20 * l + 12])));
3097	#endif
3098	_mm256_store_pd(&v[12],vv[3]);
3099
3100	#ifdef _FMA
3101	ev = _mm256_load_pd(&extEV[20 * l + 16]);
3102	vv[4] = FMAMACC(vv[4], al, ev);
3103	#else
3104	vv[4] = _mm256_add_pd(vv[4],_mm256_mul_pd(al, _mm256_load_pd(&extEV[k][20 * l + 16])));
3105	#endif
3106	_mm256_store_pd(&v[16],vv[4]);
3107	}
3108	}
3109	v = &(x3[80 * i]);
3110	scale = 1;
3111	__m256d minlikelihood_avx = _mm256_set1_pd(minlikelihood);
3112
3113	for(l = 0; scale && (l < 80); l += 4)
3114	{
3115	__m256d vv = _mm256_load_pd(&v[l]);
3116	__m256d vv_abs = _mm256_and_pd(vv,absMask_AVX.m);
3117	vv_abs = _mm256_cmp_pd(vv_abs,minlikelihood_avx,_CMP_LT_OS);
3118	if(_mm256_movemask_pd(vv_abs) != 15)
3119	scale = 0;
3120	}
3121
3122	if(scale)
3123	{
3124	__m256d twotothe256v = _mm256_set_pd(twotothe256,twotothe256,twotothe256,twotothe256);
3125	for(l = 0; l < 80; l += 4)
3126	{
3127	__m256d vv = _mm256_load_pd(&v[l]);
3128	_mm256_store_pd(&v[l],_mm256_mul_pd(vv,twotothe256v));
3129	}
3130	if(useFastScaling)
3131	addScale += wgt[i];
3132	else
3133	ex3[i] += 1;
3134	}
3135	}
3136	break;
3137	default:
3138	assert(0);
3139	}
3140
3141	if(useFastScaling)
3142	*scalerIncrement = addScale;
3143	}
3144
3145
3146	void newviewGTRGAMMAPROT_AVX_GAPPED_SAVE(int tipCase,
3147	double x1_start, double x2_start, double x3_start, double extEV, double *tipVector,
3148	int ex3, unsigned char tipX1, unsigned char *tipX2, int n,
3149	double left, double right, int wgt, int scalerIncrement, const boolean useFastScaling,
3150	unsigned int x1_gap, unsigned int x2_gap, unsigned int *x3_gap,
3151	double x1_gapColumn, double x2_gapColumn, double *x3_gapColumn)
3152	{
3153	double
3154	*x1 = x1_start,
3155	*x2 = x2_start,
3156	*x3_ptr = x3_start,
3157	*x2_ptr = x2_start,
3158	*x1_ptr = x1_start,
3159	*uX1,
3160	*uX2,
3161	*v,
3162	x1px2,
3163	*vl,
3164	*vr;
3165
3166	int
3167	i,
3168	j,
3169	l,
3170	k,
3171	gapScaling = 0,
3172	scale,
3173	addScale = 0;
3174
3175
3176	#ifndef GCC_VERSION
3177	#define GCC_VERSION (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__)
3178	#endif
3179
3180
3181	#if GCC_VERSION < 40500
3182	__m256d
3183	bitmask = _mm256_set_pd(0,0,0,-1);
3184	#else
3185	__m256i
3186	bitmask = _mm256_set_epi32(0, 0, 0, 0, 0, 0, -1, -1);
3187	#endif
3188
3189	switch(tipCase)
3190	{
3191	case TIP_TIP:
3192	{
3193	double
3194	umpX1[1840] __attribute__ ((aligned (BYTE_ALIGNMENT))),
3195	umpX2[1840] __attribute__ ((aligned (BYTE_ALIGNMENT)));
3196
3197
3198
3199	for(i = 0; i < 23; i++)
3200	{
3201	v = &(tipVector[20 * i]);
3202
3203	for(k = 0; k < 80; k++)
3204	{
3205	double
3206	ll = &left[k 20],
3207	rr = &right[k 20];
3208
3209	__m256d
3210	umpX1v = _mm256_setzero_pd(),
3211	umpX2v = _mm256_setzero_pd();
3212
3213	for(l = 0; l < 20; l+=4)
3214	{
3215	__m256d vv = _mm256_load_pd(&v[l]);
3216	#ifdef _FMA
3217	__m256d llv = _mm256_load_pd(&ll[l]);
3218	umpX1v = FMAMACC(umpX1v,vv,llv);
3219	__m256d rrv = _mm256_load_pd(&rr[l]);
3220	umpX2v = FMAMACC(umpX2v,vv,rrv);
3221	#else
3222	umpX1v = _mm256_add_pd(umpX1v,_mm256_mul_pd(vv,_mm256_load_pd(&ll[l])));
3223	umpX2v = _mm256_add_pd(umpX2v,_mm256_mul_pd(vv,_mm256_load_pd(&rr[l])));
3224	#endif
3225	}
3226
3227	umpX1v = hadd3(umpX1v);
3228	umpX2v = hadd3(umpX2v);
3229	_mm256_maskstore_pd(&umpX1[80 * i + k], bitmask, umpX1v);
3230	_mm256_maskstore_pd(&umpX2[80 * i + k], bitmask, umpX2v);
3231	}
3232	}
3233
3234
3235	{
3236	uX1 = &umpX1[1760];
3237	uX2 = &umpX2[1760];
3238
3239	for(j = 0; j < 4; j++)
3240	{
3241	__m256d vv[5];
3242
3243	v = &x3_gapColumn[j * 20];
3244
3245	vv[0] = _mm256_setzero_pd();
3246	vv[1] = _mm256_setzero_pd();
3247	vv[2] = _mm256_setzero_pd();
3248	vv[3] = _mm256_setzero_pd();
3249	vv[4] = _mm256_setzero_pd();
3250
3251	for(k = 0; k < 20; k++)
3252	{
3253	x1px2 = uX1[j * 20 + k] * uX2[j * 20 + k];
3254
3255	__m256d x1px2v = _mm256_set1_pd(x1px2);
3256
3257	__m256d extEvv = _mm256_load_pd(&extEV[20 * k]);
3258	#ifdef _FMA
3259	vv[0] = FMAMACC(vv[0],x1px2v,extEvv);
3260	#else
3261	vv[0] = _mm256_add_pd(vv[0],_mm256_mul_pd(x1px2v,extEvv));
3262	#endif
3263	_mm256_store_pd(&v[0],vv[0]);
3264
3265	extEvv = _mm256_load_pd(&extEV[20 * k + 4]);
3266	#ifdef _FMA
3267	vv[1] = FMAMACC(vv[1],x1px2v,extEvv);
3268	#else
3269	vv[1] = _mm256_add_pd(vv[1],_mm256_mul_pd(x1px2v,extEvv));
3270	#endif
3271	_mm256_store_pd(&v[4],vv[1]);
3272
3273	extEvv = _mm256_load_pd(&extEV[20 * k + 8]);
3274	#ifdef _FMA
3275	vv[2] = FMAMACC(vv[2],x1px2v,extEvv);
3276	#else
3277	vv[2] = _mm256_add_pd(vv[2],_mm256_mul_pd(x1px2v,extEvv));
3278	#endif
3279	_mm256_store_pd(&v[8],vv[2]);
3280
3281	extEvv = _mm256_load_pd(&extEV[20 * k + 12]);
3282	#ifdef _FMA
3283	vv[3] = FMAMACC(vv[3],x1px2v,extEvv);
3284	#else
3285	vv[3] = _mm256_add_pd(vv[3],_mm256_mul_pd(x1px2v,extEvv));
3286	#endif
3287	_mm256_store_pd(&v[12],vv[3]);
3288
3289	extEvv = _mm256_load_pd(&extEV[20 * k + 16]);
3290	#ifdef _FMA
3291	vv[4] = FMAMACC(vv[4],x1px2v,extEvv);
3292	#else
3293	vv[4] = _mm256_add_pd(vv[4],_mm256_mul_pd(x1px2v,extEvv));
3294	#endif
3295	_mm256_store_pd(&v[16],vv[4]);
3296	}
3297	}
3298	}
3299
3300
3301	for(i = 0; i < n; i++)
3302	{
3303	if(!(x3_gap[i / 32] & mask32[i % 32]))
3304	{
3305	uX1 = &umpX1[80 * tipX1[i]];
3306	uX2 = &umpX2[80 * tipX2[i]];
3307
3308	for(j = 0; j < 4; j++)
3309	{
3310	__m256d vv[5];
3311
3312	v = &x3_ptr[j * 20];
3313
3314	vv[0] = _mm256_setzero_pd();
3315	vv[1] = _mm256_setzero_pd();
3316	vv[2] = _mm256_setzero_pd();
3317	vv[3] = _mm256_setzero_pd();
3318	vv[4] = _mm256_setzero_pd();
3319
3320	for(k = 0; k < 20; k++)
3321	{
3322	x1px2 = uX1[j * 20 + k] * uX2[j * 20 + k];
3323
3324	__m256d x1px2v = _mm256_set1_pd(x1px2);
3325
3326	__m256d extEvv = _mm256_load_pd(&extEV[20 * k]);
3327	#ifdef _FMA
3328	vv[0] = FMAMACC(vv[0],x1px2v,extEvv);
3329	#else
3330	vv[0] = _mm256_add_pd(vv[0],_mm256_mul_pd(x1px2v,extEvv));
3331	#endif
3332	_mm256_store_pd(&v[0],vv[0]);
3333
3334	extEvv = _mm256_load_pd(&extEV[20 * k + 4]);
3335	#ifdef _FMA
3336	vv[1] = FMAMACC(vv[1],x1px2v,extEvv);
3337	#else
3338	vv[1] = _mm256_add_pd(vv[1],_mm256_mul_pd(x1px2v,extEvv));
3339	#endif
3340	_mm256_store_pd(&v[4],vv[1]);
3341
3342	extEvv = _mm256_load_pd(&extEV[20 * k + 8]);
3343	#ifdef _FMA
3344	vv[2] = FMAMACC(vv[2],x1px2v,extEvv);
3345	#else
3346	vv[2] = _mm256_add_pd(vv[2],_mm256_mul_pd(x1px2v,extEvv));
3347	#endif
3348	_mm256_store_pd(&v[8],vv[2]);
3349
3350	extEvv = _mm256_load_pd(&extEV[20 * k + 12]);
3351	#ifdef _FMA
3352	vv[3] = FMAMACC(vv[3],x1px2v,extEvv);
3353	#else
3354	vv[3] = _mm256_add_pd(vv[3],_mm256_mul_pd(x1px2v,extEvv));
3355	#endif
3356	_mm256_store_pd(&v[12],vv[3]);
3357
3358	extEvv = _mm256_load_pd(&extEV[20 * k + 16]);
3359	#ifdef _FMA
3360	vv[4] = FMAMACC(vv[4],x1px2v,extEvv);
3361	#else
3362	vv[4] = _mm256_add_pd(vv[4],_mm256_mul_pd(x1px2v,extEvv));
3363	#endif
3364	_mm256_store_pd(&v[16],vv[4]);
3365	}
3366	}
3367	x3_ptr += 80;
3368	}
3369	}
3370	}
3371	break;
3372	case TIP_INNER:
3373	{
3374	double
3375	umpX1[1840] __attribute__ ((aligned (BYTE_ALIGNMENT))),
3376	ump_x2[20] __attribute__ ((aligned (BYTE_ALIGNMENT)));
3377
3378
3379
3380	for(i = 0; i < 23; i++)
3381	{
3382	v = &(tipVector[20 * i]);
3383
3384	for(k = 0; k < 80; k++)
3385	{
3386	__m256d umpX1v = _mm256_setzero_pd();
3387	for(l = 0; l < 20; l+=4)
3388	{
3389	__m256d vv = _mm256_load_pd(&v[l]);
3390	__m256d leftv = _mm256_load_pd(&left[k * 20 + l]);
3391	#ifdef _FMA
3392
3393	umpX1v = FMAMACC(umpX1v, vv, leftv);
3394	#else
3395	umpX1v = _mm256_add_pd(umpX1v, _mm256_mul_pd(vv, leftv));
3396	#endif
3397	}
3398	umpX1v = hadd3(umpX1v);
3399	_mm256_maskstore_pd(&umpX1[80 * i + k], bitmask, umpX1v);
3400	}
3401	}
3402
3403	{
3404	uX1 = &umpX1[1760];
3405
3406	for(k = 0; k < 4; k++)
3407	{
3408	v = &(x2_gapColumn[k * 20]);
3409
3410	for(l = 0; l < 20; l++)
3411	{
3412	__m256d ump_x2v = _mm256_setzero_pd();
3413
3414	__m256d vv = _mm256_load_pd(&v[0]);
3415	__m256d rightv = _mm256_load_pd(&right[k400+l20+0]);
3416	#ifdef _FMA
3417	ump_x2v = FMAMACC(ump_x2v,vv,rightv);
3418	#else
3419	ump_x2v = _mm256_add_pd(ump_x2v, _mm256_mul_pd(vv, rightv));
3420	#endif
3421
3422	vv = _mm256_load_pd(&v[4]);
3423	rightv = _mm256_load_pd(&right[k400+l20+4]);
3424	#ifdef _FMA
3425	ump_x2v = FMAMACC(ump_x2v,vv,rightv);
3426	#else
3427	ump_x2v = _mm256_add_pd(ump_x2v, _mm256_mul_pd(vv, rightv));
3428	#endif
3429
3430	vv = _mm256_load_pd(&v[8]);
3431	rightv = _mm256_load_pd(&right[k400+l20+8]);
3432	#ifdef _FMA
3433	ump_x2v = FMAMACC(ump_x2v,vv,rightv);
3434	#else
3435	ump_x2v = _mm256_add_pd(ump_x2v, _mm256_mul_pd(vv, rightv));
3436	#endif
3437
3438	vv = _mm256_load_pd(&v[12]);
3439	rightv = _mm256_load_pd(&right[k400+l20+12]);
3440	#ifdef _FMA
3441	ump_x2v = FMAMACC(ump_x2v,vv,rightv);
3442	#else
3443	ump_x2v = _mm256_add_pd(ump_x2v, _mm256_mul_pd(vv, rightv));
3444	#endif
3445
3446	vv = _mm256_load_pd(&v[16]);
3447	rightv = _mm256_load_pd(&right[k400+l20+16]);
3448	#ifdef _FMA
3449	ump_x2v = FMAMACC(ump_x2v,vv,rightv);
3450	#else
3451	ump_x2v = _mm256_add_pd(ump_x2v, _mm256_mul_pd(vv, rightv));
3452	#endif
3453
3454	ump_x2v = hadd3(ump_x2v);
3455	_mm256_maskstore_pd(&ump_x2[l], bitmask, ump_x2v);
3456	}
3457
3458	v = &x3_gapColumn[20 * k];
3459
3460	__m256d vv[5];
3461
3462	vv[0] = _mm256_setzero_pd();
3463	vv[1] = _mm256_setzero_pd();
3464	vv[2] = _mm256_setzero_pd();
3465	vv[3] = _mm256_setzero_pd();
3466	vv[4] = _mm256_setzero_pd();
3467
3468	for(l = 0; l < 20; l++)
3469	{
3470	x1px2 = uX1[k * 20 + l] * ump_x2[l];
3471	__m256d x1px2v = _mm256_set1_pd(x1px2);
3472
3473	#ifdef _FMA
3474	__m256d ev = _mm256_load_pd(&extEV[l * 20 + 0]);
3475	vv[0] = FMAMACC(vv[0],x1px2v, ev);
3476	#else
3477	vv[0] = _mm256_add_pd(vv[0],_mm256_mul_pd(x1px2v, _mm256_load_pd(&extEV[l * 20 + 0])));
3478	#endif
3479	_mm256_store_pd(&v[0],vv[0]);
3480
3481	#ifdef _FMA
3482	ev = _mm256_load_pd(&extEV[l * 20 + 4]);
3483	vv[1] = FMAMACC(vv[1],x1px2v, ev);
3484	#else
3485	vv[1] = _mm256_add_pd(vv[1],_mm256_mul_pd(x1px2v, _mm256_load_pd(&extEV[l * 20 + 4])));
3486	#endif
3487	_mm256_store_pd(&v[4],vv[1]);
3488
3489	#ifdef _FMA
3490	ev = _mm256_load_pd(&extEV[l * 20 + 8]);
3491	vv[2] = FMAMACC(vv[2],x1px2v, ev);
3492	#else
3493	vv[2] = _mm256_add_pd(vv[2],_mm256_mul_pd(x1px2v, _mm256_load_pd(&extEV[l * 20 + 8])));
3494	#endif
3495	_mm256_store_pd(&v[8],vv[2]);
3496
3497	#ifdef _FMA
3498	ev = _mm256_load_pd(&extEV[l * 20 + 12]);
3499	vv[3] = FMAMACC(vv[3],x1px2v, ev);
3500	#else
3501	vv[3] = _mm256_add_pd(vv[3],_mm256_mul_pd(x1px2v, _mm256_load_pd(&extEV[l * 20 + 12])));
3502	#endif
3503	_mm256_store_pd(&v[12],vv[3]);
3504
3505
3506	#ifdef _FMA
3507	ev = _mm256_load_pd(&extEV[l * 20 + 16]);
3508	vv[4] = FMAMACC(vv[4],x1px2v, ev);
3509	#else
3510	vv[4] = _mm256_add_pd(vv[4],_mm256_mul_pd(x1px2v, _mm256_load_pd(&extEV[l * 20 + 16])));
3511	#endif
3512	_mm256_store_pd(&v[16],vv[4]);
3513
3514	}
3515	}
3516
3517	v = x3_gapColumn;
3518	__m256d minlikelihood_avx = _mm256_set1_pd(minlikelihood);
3519	scale = 1;
3520	for(l = 0; scale && (l < 80); l += 4)
3521	{
3522	__m256d vv = _mm256_load_pd(&v[l]);
3523	__m256d vv_abs = _mm256_and_pd(vv,absMask_AVX.m);
3524	vv_abs = _mm256_cmp_pd(vv_abs,minlikelihood_avx,_CMP_LT_OS);
3525	if(_mm256_movemask_pd(vv_abs) != 15)
3526	scale = 0;
3527	}
3528
3529	if(scale)
3530	{
3531	__m256d twotothe256v = _mm256_set_pd(twotothe256,twotothe256,twotothe256,twotothe256);
3532	gapScaling = 1;
3533
3534	for(l = 0; l < 80; l += 4)
3535	{
3536	__m256d vv = _mm256_load_pd(&v[l]);
3537	_mm256_store_pd(&v[l],_mm256_mul_pd(vv,twotothe256v));
3538	}
3539	}
3540	}
3541
3542	for (i = 0; i < n; i++)
3543	{
3544	if((x3_gap[i / 32] & mask32[i % 32]))
3545	{
3546	if(gapScaling)
3547	{
3548	if(useFastScaling)
3549	addScale += wgt[i];
3550	else
3551	ex3[i] += 1;
3552	}
3553	}
3554	else
3555	{
3556	uX1 = &umpX1[80 * tipX1[i]];
3557
3558	if(x2_gap[i / 32] & mask32[i % 32])
3559	x2 = x2_gapColumn;
3560	else
3561	{
3562	x2 = x2_ptr;
3563	x2_ptr += 80;
3564	}
3565
3566	for(k = 0; k < 4; k++)
3567	{
3568	v = &(x2[k * 20]);
3569
3570	for(l = 0; l < 20; l++)
3571	{
3572	__m256d ump_x2v = _mm256_setzero_pd();
3573
3574	__m256d vv = _mm256_load_pd(&v[0]);
3575	__m256d rightv = _mm256_load_pd(&right[k400+l20+0]);
3576	#ifdef _FMA
3577	ump_x2v = FMAMACC(ump_x2v,vv,rightv);
3578	#else
3579	ump_x2v = _mm256_add_pd(ump_x2v, _mm256_mul_pd(vv, rightv));
3580	#endif
3581
3582	vv = _mm256_load_pd(&v[4]);
3583	rightv = _mm256_load_pd(&right[k400+l20+4]);
3584	#ifdef _FMA
3585	ump_x2v = FMAMACC(ump_x2v,vv,rightv);
3586	#else
3587	ump_x2v = _mm256_add_pd(ump_x2v, _mm256_mul_pd(vv, rightv));
3588	#endif
3589
3590	vv = _mm256_load_pd(&v[8]);
3591	rightv = _mm256_load_pd(&right[k400+l20+8]);
3592	#ifdef _FMA
3593	ump_x2v = FMAMACC(ump_x2v,vv,rightv);
3594	#else
3595	ump_x2v = _mm256_add_pd(ump_x2v, _mm256_mul_pd(vv, rightv));
3596	#endif
3597
3598	vv = _mm256_load_pd(&v[12]);
3599	rightv = _mm256_load_pd(&right[k400+l20+12]);
3600	#ifdef _FMA
3601	ump_x2v = FMAMACC(ump_x2v,vv,rightv);
3602	#else
3603	ump_x2v = _mm256_add_pd(ump_x2v, _mm256_mul_pd(vv, rightv));
3604	#endif
3605
3606	vv = _mm256_load_pd(&v[16]);
3607	rightv = _mm256_load_pd(&right[k400+l20+16]);
3608	#ifdef _FMA
3609	ump_x2v = FMAMACC(ump_x2v,vv,rightv);
3610	#else
3611	ump_x2v = _mm256_add_pd(ump_x2v, _mm256_mul_pd(vv, rightv));
3612	#endif
3613
3614	ump_x2v = hadd3(ump_x2v);
3615	_mm256_maskstore_pd(&ump_x2[l], bitmask, ump_x2v);
3616	}
3617
3618
3619	v = &x3_ptr[k * 20];
3620
3621	__m256d vv[5];
3622
3623	vv[0] = _mm256_setzero_pd();
3624	vv[1] = _mm256_setzero_pd();
3625	vv[2] = _mm256_setzero_pd();
3626	vv[3] = _mm256_setzero_pd();
3627	vv[4] = _mm256_setzero_pd();
3628
3629	for(l = 0; l < 20; l++)
3630	{
3631	x1px2 = uX1[k * 20 + l] * ump_x2[l];
3632	__m256d x1px2v = _mm256_set1_pd(x1px2);
3633
3634	#ifdef _FMA
3635	__m256d ev = _mm256_load_pd(&extEV[l * 20 + 0]);
3636	vv[0] = FMAMACC(vv[0],x1px2v, ev);
3637	#else
3638	vv[0] = _mm256_add_pd(vv[0],_mm256_mul_pd(x1px2v, _mm256_load_pd(&extEV[l * 20 + 0])));
3639	#endif
3640	_mm256_store_pd(&v[0],vv[0]);
3641
3642	#ifdef _FMA
3643	ev = _mm256_load_pd(&extEV[l * 20 + 4]);
3644	vv[1] = FMAMACC(vv[1],x1px2v, ev);
3645	#else
3646	vv[1] = _mm256_add_pd(vv[1],_mm256_mul_pd(x1px2v, _mm256_load_pd(&extEV[l * 20 + 4])));
3647	#endif
3648	_mm256_store_pd(&v[4],vv[1]);
3649
3650	#ifdef _FMA
3651	ev = _mm256_load_pd(&extEV[l * 20 + 8]);
3652	vv[2] = FMAMACC(vv[2],x1px2v, ev);
3653	#else
3654	vv[2] = _mm256_add_pd(vv[2],_mm256_mul_pd(x1px2v, _mm256_load_pd(&extEV[l * 20 + 8])));
3655	#endif
3656	_mm256_store_pd(&v[8],vv[2]);
3657
3658	#ifdef _FMA
3659	ev = _mm256_load_pd(&extEV[l * 20 + 12]);
3660	vv[3] = FMAMACC(vv[3],x1px2v, ev);
3661	#else
3662	vv[3] = _mm256_add_pd(vv[3],_mm256_mul_pd(x1px2v, _mm256_load_pd(&extEV[l * 20 + 12])));
3663	#endif
3664	_mm256_store_pd(&v[12],vv[3]);
3665
3666
3667	#ifdef _FMA
3668	ev = _mm256_load_pd(&extEV[l * 20 + 16]);
3669	vv[4] = FMAMACC(vv[4],x1px2v, ev);
3670	#else
3671	vv[4] = _mm256_add_pd(vv[4],_mm256_mul_pd(x1px2v, _mm256_load_pd(&extEV[l * 20 + 16])));
3672	#endif
3673	_mm256_store_pd(&v[16],vv[4]);
3674
3675	}
3676	}
3677
3678	v = x3_ptr;
3679	__m256d minlikelihood_avx = _mm256_set1_pd(minlikelihood);
3680	scale = 1;
3681	for(l = 0; scale && (l < 80); l += 4)
3682	{
3683	__m256d vv = _mm256_load_pd(&v[l]);
3684	__m256d vv_abs = _mm256_and_pd(vv,absMask_AVX.m);
3685	vv_abs = _mm256_cmp_pd(vv_abs,minlikelihood_avx,_CMP_LT_OS);
3686	if(_mm256_movemask_pd(vv_abs) != 15)
3687	scale = 0;
3688	}
3689
3690	if(scale)
3691	{
3692	__m256d twotothe256v = _mm256_set_pd(twotothe256,twotothe256,twotothe256,twotothe256);
3693	for(l = 0; l < 80; l += 4)
3694	{
3695	__m256d vv = _mm256_load_pd(&v[l]);
3696	_mm256_store_pd(&v[l],_mm256_mul_pd(vv,twotothe256v));
3697	}
3698	if(useFastScaling)
3699	addScale += wgt[i];
3700	else
3701	ex3[i] += 1;
3702	}
3703	x3_ptr += 80;
3704	}
3705	}
3706	}
3707	break;
3708	case INNER_INNER:
3709	for(k = 0; k < 4; k++)
3710	{
3711	vl = &(x1_gapColumn[20 * k]);
3712	vr = &(x2_gapColumn[20 * k]);
3713	v = &(x3_gapColumn[20 * k]);
3714
3715	__m256d vv[5];
3716
3717	vv[0] = _mm256_setzero_pd();
3718	vv[1] = _mm256_setzero_pd();
3719	vv[2] = _mm256_setzero_pd();
3720	vv[3] = _mm256_setzero_pd();
3721	vv[4] = _mm256_setzero_pd();
3722
3723	for(l = 0; l < 20; l++)
3724	{
3725	__m256d al = _mm256_setzero_pd();
3726	__m256d ar = _mm256_setzero_pd();
3727
3728	__m256d leftv = _mm256_load_pd(&left[k * 400 + l * 20 + 0]);
3729	__m256d rightv = _mm256_load_pd(&right[k * 400 + l * 20 + 0]);
3730	__m256d vlv = _mm256_load_pd(&vl[0]);
3731	__m256d vrv = _mm256_load_pd(&vr[0]);
3732
3733	#ifdef _FMA
3734
3735	al = FMAMACC(al, vlv, leftv);
3736	ar = FMAMACC(ar, vrv, rightv);
3737	#else
3738	al = _mm256_add_pd(al,_mm256_mul_pd(vlv,leftv));
3739	ar = _mm256_add_pd(ar,_mm256_mul_pd(vrv,rightv));
3740	#endif
3741
3742	leftv = _mm256_load_pd(&left[k * 400 + l * 20 + 4]);
3743	rightv = _mm256_load_pd(&right[k * 400 + l * 20 + 4]);
3744	vlv = _mm256_load_pd(&vl[4]);
3745	vrv = _mm256_load_pd(&vr[4]);
3746	#ifdef _FMA
3747
3748	al = FMAMACC(al, vlv, leftv);
3749	ar = FMAMACC(ar, vrv, rightv);
3750	#else
3751	al = _mm256_add_pd(al,_mm256_mul_pd(vlv,leftv));
3752	ar = _mm256_add_pd(ar,_mm256_mul_pd(vrv,rightv));
3753	#endif
3754
3755	leftv = _mm256_load_pd(&left[k * 400 + l * 20 + 8]);
3756	rightv = _mm256_load_pd(&right[k * 400 + l * 20 + 8]);
3757	vlv = _mm256_load_pd(&vl[8]);
3758	vrv = _mm256_load_pd(&vr[8]);
3759	#ifdef _FMA
3760
3761	al = FMAMACC(al, vlv, leftv);
3762	ar = FMAMACC(ar, vrv, rightv);
3763	#else
3764	al = _mm256_add_pd(al,_mm256_mul_pd(vlv,leftv));
3765	ar = _mm256_add_pd(ar,_mm256_mul_pd(vrv,rightv));
3766	#endif
3767
3768	leftv = _mm256_load_pd(&left[k * 400 + l * 20 + 12]);
3769	rightv = _mm256_load_pd(&right[k * 400 + l * 20 + 12]);
3770	vlv = _mm256_load_pd(&vl[12]);
3771	vrv = _mm256_load_pd(&vr[12]);
3772	#ifdef _FMA
3773
3774	al = FMAMACC(al, vlv, leftv);
3775	ar = FMAMACC(ar, vrv, rightv);
3776	#else
3777	al = _mm256_add_pd(al,_mm256_mul_pd(vlv,leftv));
3778	ar = _mm256_add_pd(ar,_mm256_mul_pd(vrv,rightv));
3779	#endif
3780
3781	leftv = _mm256_load_pd(&left[k * 400 + l * 20 + 16]);
3782	rightv = _mm256_load_pd(&right[k * 400 + l * 20 + 16]);
3783	vlv = _mm256_load_pd(&vl[16]);
3784	vrv = _mm256_load_pd(&vr[16]);
3785
3786	#ifdef _FMA
3787	al = FMAMACC(al, vlv, leftv);
3788	ar = FMAMACC(ar, vrv, rightv);
3789	#else
3790	al = _mm256_add_pd(al,_mm256_mul_pd(vlv,leftv));
3791	ar = _mm256_add_pd(ar,_mm256_mul_pd(vrv,rightv));
3792	#endif
3793
3794	/**************************************************************************************************************/
3795
3796	al = hadd3(al);
3797	ar = hadd3(ar);
3798	al = _mm256_mul_pd(ar,al);
3799
3800	/************************************************************************************************************/
3801	#ifdef _FMA
3802	__m256d ev = _mm256_load_pd(&extEV[20 * l + 0]);
3803	vv[0] = FMAMACC(vv[0], al, ev);
3804	#else
3805	vv[0] = _mm256_add_pd(vv[0],_mm256_mul_pd(al, _mm256_load_pd(&extEV[20 * l + 0])));
3806	#endif
3807	_mm256_store_pd(&v[0],vv[0]);
3808
3809	#ifdef _FMA
3810	ev = _mm256_load_pd(&extEV[20 * l + 4]);
3811	vv[1] = FMAMACC(vv[1], al, ev);
3812	#else
3813	vv[1] = _mm256_add_pd(vv[1],_mm256_mul_pd(al, _mm256_load_pd(&extEV[20 * l + 4])));
3814	#endif
3815	_mm256_store_pd(&v[4],vv[1]);
3816
3817	#ifdef _FMA
3818	ev = _mm256_load_pd(&extEV[20 * l + 8]);
3819	vv[2] = FMAMACC(vv[2], al, ev);
3820	#else
3821	vv[2] = _mm256_add_pd(vv[2],_mm256_mul_pd(al, _mm256_load_pd(&extEV[20 * l + 8])));
3822	#endif
3823	_mm256_store_pd(&v[8],vv[2]);
3824
3825	#ifdef _FMA
3826	ev = _mm256_load_pd(&extEV[20 * l + 12]);
3827	vv[3] = FMAMACC(vv[3], al, ev);
3828	#else
3829	vv[3] = _mm256_add_pd(vv[3],_mm256_mul_pd(al, _mm256_load_pd(&extEV[20 * l + 12])));
3830	#endif
3831	_mm256_store_pd(&v[12],vv[3]);
3832
3833	#ifdef _FMA
3834	ev = _mm256_load_pd(&extEV[20 * l + 16]);
3835	vv[4] = FMAMACC(vv[4], al, ev);
3836	#else
3837	vv[4] = _mm256_add_pd(vv[4],_mm256_mul_pd(al, _mm256_load_pd(&extEV[20 * l + 16])));
3838	#endif
3839	_mm256_store_pd(&v[16],vv[4]);
3840	}
3841	}
3842
3843	v = x3_gapColumn;
3844	scale = 1;
3845	__m256d minlikelihood_avx = _mm256_set1_pd(minlikelihood);
3846
3847	for(l = 0; scale && (l < 80); l += 4)
3848	{
3849	__m256d vv = _mm256_load_pd(&v[l]);
3850	__m256d vv_abs = _mm256_and_pd(vv,absMask_AVX.m);
3851	vv_abs = _mm256_cmp_pd(vv_abs,minlikelihood_avx,_CMP_LT_OS);
3852	if(_mm256_movemask_pd(vv_abs) != 15)
3853	scale = 0;
3854	}
3855
3856	if(scale)
3857	{
3858	__m256d twotothe256v = _mm256_set_pd(twotothe256,twotothe256,twotothe256,twotothe256);
3859	gapScaling = 1;
3860
3861	for(l = 0; l < 80; l += 4)
3862	{
3863	__m256d vv = _mm256_load_pd(&v[l]);
3864	_mm256_store_pd(&v[l],_mm256_mul_pd(vv,twotothe256v));
3865	}
3866
3867	}
3868
3869
3870
3871	for(i = 0; i < n; i++)
3872	{
3873
3874	if(x3_gap[i / 32] & mask32[i % 32])
3875	{
3876	if(gapScaling)
3877	{
3878	if(useFastScaling)
3879	addScale += wgt[i];
3880	else
3881	ex3[i] += 1;
3882	}
3883	}
3884	else
3885	{
3886	if(x1_gap[i / 32] & mask32[i % 32])
3887	x1 = x1_gapColumn;
3888	else
3889	{
3890	x1 = x1_ptr;
3891	x1_ptr += 80;
3892	}
3893
3894	if(x2_gap[i / 32] & mask32[i % 32])
3895	x2 = x2_gapColumn;
3896	else
3897	{
3898	x2 = x2_ptr;
3899	x2_ptr += 80;
3900	}
3901
3902	for(k = 0; k < 4; k++)
3903	{
3904	vl = &(x1[20 * k]);
3905	vr = &(x2[20 * k]);
3906	v = &(x3_ptr[20 * k]);
3907
3908	__m256d vv[5];
3909
3910	vv[0] = _mm256_setzero_pd();
3911	vv[1] = _mm256_setzero_pd();
3912	vv[2] = _mm256_setzero_pd();
3913	vv[3] = _mm256_setzero_pd();
3914	vv[4] = _mm256_setzero_pd();
3915
3916	for(l = 0; l < 20; l++)
3917	{
3918	__m256d al = _mm256_setzero_pd();
3919	__m256d ar = _mm256_setzero_pd();
3920
3921	__m256d leftv = _mm256_load_pd(&left[k * 400 + l * 20 + 0]);
3922	__m256d rightv = _mm256_load_pd(&right[k * 400 + l * 20 + 0]);
3923	__m256d vlv = _mm256_load_pd(&vl[0]);
3924	__m256d vrv = _mm256_load_pd(&vr[0]);
3925
3926	#ifdef _FMA
3927
3928	al = FMAMACC(al, vlv, leftv);
3929	ar = FMAMACC(ar, vrv, rightv);
3930	#else
3931	al = _mm256_add_pd(al,_mm256_mul_pd(vlv,leftv));
3932	ar = _mm256_add_pd(ar,_mm256_mul_pd(vrv,rightv));
3933	#endif
3934
3935	leftv = _mm256_load_pd(&left[k * 400 + l * 20 + 4]);
3936	rightv = _mm256_load_pd(&right[k * 400 + l * 20 + 4]);
3937	vlv = _mm256_load_pd(&vl[4]);
3938	vrv = _mm256_load_pd(&vr[4]);
3939	#ifdef _FMA
3940
3941	al = FMAMACC(al, vlv, leftv);
3942	ar = FMAMACC(ar, vrv, rightv);
3943	#else
3944	al = _mm256_add_pd(al,_mm256_mul_pd(vlv,leftv));
3945	ar = _mm256_add_pd(ar,_mm256_mul_pd(vrv,rightv));
3946	#endif
3947
3948	leftv = _mm256_load_pd(&left[k * 400 + l * 20 + 8]);
3949	rightv = _mm256_load_pd(&right[k * 400 + l * 20 + 8]);
3950	vlv = _mm256_load_pd(&vl[8]);
3951	vrv = _mm256_load_pd(&vr[8]);
3952	#ifdef _FMA
3953
3954	al = FMAMACC(al, vlv, leftv);
3955	ar = FMAMACC(ar, vrv, rightv);
3956	#else
3957	al = _mm256_add_pd(al,_mm256_mul_pd(vlv,leftv));
3958	ar = _mm256_add_pd(ar,_mm256_mul_pd(vrv,rightv));
3959	#endif
3960
3961	leftv = _mm256_load_pd(&left[k * 400 + l * 20 + 12]);
3962	rightv = _mm256_load_pd(&right[k * 400 + l * 20 + 12]);
3963	vlv = _mm256_load_pd(&vl[12]);
3964	vrv = _mm256_load_pd(&vr[12]);
3965	#ifdef _FMA
3966
3967	al = FMAMACC(al, vlv, leftv);
3968	ar = FMAMACC(ar, vrv, rightv);
3969	#else
3970	al = _mm256_add_pd(al,_mm256_mul_pd(vlv,leftv));
3971	ar = _mm256_add_pd(ar,_mm256_mul_pd(vrv,rightv));
3972	#endif
3973
3974	leftv = _mm256_load_pd(&left[k * 400 + l * 20 + 16]);
3975	rightv = _mm256_load_pd(&right[k * 400 + l * 20 + 16]);
3976	vlv = _mm256_load_pd(&vl[16]);
3977	vrv = _mm256_load_pd(&vr[16]);
3978
3979	#ifdef _FMA
3980	al = FMAMACC(al, vlv, leftv);
3981	ar = FMAMACC(ar, vrv, rightv);
3982	#else
3983	al = _mm256_add_pd(al,_mm256_mul_pd(vlv,leftv));
3984	ar = _mm256_add_pd(ar,_mm256_mul_pd(vrv,rightv));
3985	#endif
3986
3987	/**************************************************************************************************************/
3988
3989	al = hadd3(al);
3990	ar = hadd3(ar);
3991	al = _mm256_mul_pd(ar,al);
3992
3993	/************************************************************************************************************/
3994	#ifdef _FMA
3995	__m256d ev = _mm256_load_pd(&extEV[20 * l + 0]);
3996	vv[0] = FMAMACC(vv[0], al, ev);
3997	#else
3998	vv[0] = _mm256_add_pd(vv[0],_mm256_mul_pd(al, _mm256_load_pd(&extEV[20 * l + 0])));
3999	#endif
4000	_mm256_store_pd(&v[0],vv[0]);
4001
4002	#ifdef _FMA
4003	ev = _mm256_load_pd(&extEV[20 * l + 4]);
4004	vv[1] = FMAMACC(vv[1], al, ev);
4005	#else
4006	vv[1] = _mm256_add_pd(vv[1],_mm256_mul_pd(al, _mm256_load_pd(&extEV[20 * l + 4])));
4007	#endif
4008	_mm256_store_pd(&v[4],vv[1]);
4009
4010	#ifdef _FMA
4011	ev = _mm256_load_pd(&extEV[20 * l + 8]);
4012	vv[2] = FMAMACC(vv[2], al, ev);
4013	#else
4014	vv[2] = _mm256_add_pd(vv[2],_mm256_mul_pd(al, _mm256_load_pd(&extEV[20 * l + 8])));
4015	#endif
4016	_mm256_store_pd(&v[8],vv[2]);
4017
4018	#ifdef _FMA
4019	ev = _mm256_load_pd(&extEV[20 * l + 12]);
4020	vv[3] = FMAMACC(vv[3], al, ev);
4021	#else
4022	vv[3] = _mm256_add_pd(vv[3],_mm256_mul_pd(al, _mm256_load_pd(&extEV[20 * l + 12])));
4023	#endif
4024	_mm256_store_pd(&v[12],vv[3]);
4025
4026	#ifdef _FMA
4027	ev = _mm256_load_pd(&extEV[20 * l + 16]);
4028	vv[4] = FMAMACC(vv[4], al, ev);
4029	#else
4030	vv[4] = _mm256_add_pd(vv[4],_mm256_mul_pd(al, _mm256_load_pd(&extEV[20 * l + 16])));
4031	#endif
4032	_mm256_store_pd(&v[16],vv[4]);
4033	}
4034	}
4035
4036	v = x3_ptr;
4037	scale = 1;
4038
4039	__m256d minlikelihood_avx = _mm256_set1_pd(minlikelihood);
4040
4041	for(l = 0; scale && (l < 80); l += 4)
4042	{
4043	__m256d vv = _mm256_load_pd(&v[l]);
4044	__m256d vv_abs = _mm256_and_pd(vv,absMask_AVX.m);
4045	vv_abs = _mm256_cmp_pd(vv_abs,minlikelihood_avx,_CMP_LT_OS);
4046	if(_mm256_movemask_pd(vv_abs) != 15)
4047	scale = 0;
4048	}
4049
4050	if(scale)
4051	{
4052	__m256d twotothe256v = _mm256_set_pd(twotothe256,twotothe256,twotothe256,twotothe256);
4053	for(l = 0; l < 80; l += 4)
4054	{
4055	__m256d vv = _mm256_load_pd(&v[l]);
4056	_mm256_store_pd(&v[l],_mm256_mul_pd(vv,twotothe256v));
4057	}
4058	if(useFastScaling)
4059	addScale += wgt[i];
4060	else
4061	ex3[i] += 1;
4062	}
4063	x3_ptr += 80;
4064	}
4065	}
4066	break;
4067	default:
4068	assert(0);
4069	}
4070
4071	if(useFastScaling)
4072	*scalerIncrement = addScale;
4073	}

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: branches/species/GDE/RAxML/avxLikelihood.c

Download in other formats: