« back to all changes in this revision
Viewing changes to src/C/sse41/xalip_sse41.c
- browse files at revision 1
- compare with another revision
- download diff
- download tarball
- view history from revision 1
- Committer: Package Import Robot
- Author(s): Andreas Tille
- Date: 2017-03-31 14:01:39 UTC
- Revision ID: package-import@ubuntu.com-20170331140139-povxg86r196gejfa
Tags: upstream-3+dfsg
Import upstream version 3+dfsg
- files added:
- data
- doc
- man
- src
- src/C
- src/C/generic
- src/C/include
- src/C/prg
- src/C/sse2
- src/C/sse41
- src/C/utils
- src/Fortran
added
removed
src/C/sse41/xalip_sse41.c
1
/*******************************************************
2
PFTOOLS
3
*******************************************************
4
Sep 26, 2011 xalip_sse41.c
5
*******************************************************
6
(C) 2011 Swiss Institute of Bioinformatics
7
Thierry Schuepbach (thierry.schuepbach@isb-sib.ch)
8
*******************************************************/
9
#include <stdlib.h>
10
#ifdef XALIP_DEBUG
11
# include <stdio.h>
12
#endif
13
#include <smmintrin.h>
14
#include "profile.h"
15
16
17
18
19
#define MAX(a,b) (a>b) ? a : b
20
#define MIN(a,b) (a<b) ? a : b
21
22
struct Minima {
23
int _a;
24
int _b;
25
int _c;
26
};
27
28
static void InitR(const int iseq, const size_t N1, const size_t N2, const size_t bseq, const size_t lseq,
29
union lScores * const restrict iop, union Positions * const restrict iom,
30
union Positions * const restrict ioi, const struct Profile * const restrict prf)
31
{
32
int KOPD;
33
#ifdef XALIP_DEBUG
34
fprintf(stdout,"InitR called with index %i bseq %lu\n",iseq, bseq);
35
#endif
36
// Are we treating sequence index below given start?
37
register const ScoreTuple * restrict FirstSequenceProtein = ( iseq < bseq )
38
? &(prf->Scores.Insertion.FirstSequenceProtein[0])
39
: &(prf->Scores.Insertion.Transitions->From[EXTRA]);
40
const size_t FirstSequenceProteinAlignStep = (iseq < bseq) ? 1 : 4;
41
42
register __m128i * restrict pIOP = &(iop[0].xmm);
43
44
// Load FirstSequenceProtein
45
__m128i __FirstSequenceProtein = _mm_loadl_epi64((__m128i*) FirstSequenceProtein);
46
47
// Convert signed WORD into signed DWORD
48
__FirstSequenceProtein = _mm_cvtepi16_epi32(__FirstSequenceProtein);
49
50
// Store into iop
51
*pIOP = __FirstSequenceProtein;
52
53
// Set KOPD ( this is SSE 4.1 )
54
KOPD = _mm_extract_epi32(__FirstSequenceProtein, DELETION);
55
56
register const TransitionScores * const restrict Transitions = prf->Scores.Insertion.Transitions;
57
register const short int * restrict pMatch = &prf->Scores.Match.Alphabet[_D];
58
const size_t AlignStep = prf->Scores.Match.AlignStep;
59
60
// Move to next profile First Sequence
61
FirstSequenceProtein += FirstSequenceProteinAlignStep;
62
pIOP++;
63
64
for (unsigned int iprf=1; iprf<=(unsigned int) prf->Length; ++iprf) {
65
register const int KD = KOPD + (int) *pMatch;
66
pMatch += AlignStep;
67
68
// Transform KD into a vector
69
__m128i __KD = _mm_set1_epi32(KD);
70
// Load Transitions
71
__m128i __Transitions = _mm_loadl_epi64((__m128i*) &(Transitions[iprf].From[DELETION]));
72
// Convert signed WORD into signed DWORD
73
__Transitions = _mm_cvtepi16_epi32(__Transitions);
74
// Add KD to Transitions
75
__Transitions = _mm_add_epi32(__Transitions, __KD);
76
77
// Load FirstSequenceProtein
78
__m128i __FirstSequenceProtein = _mm_loadl_epi64((__m128i*) FirstSequenceProtein);
79
// Convert signed WORD into signed DWORD
80
__FirstSequenceProtein = _mm_cvtepi16_epi32(__FirstSequenceProtein);
81
82
// Move to next profile First Sequence
83
FirstSequenceProtein += FirstSequenceProteinAlignStep;
84
85
// Get maximum ( this is SSE 4.1 )
86
__m128i __max = _mm_max_epi32(__Transitions, __FirstSequenceProtein);
87
88
// Store IOPI and IOPM
89
*pIOP = __max;
90
pIOP++;
91
92
// Set KOPD ( this is SSE 4.1 )
93
KOPD = _mm_extract_epi32(__max, DELETION);
94
}
95
96
97
union Positions TPOS __attribute__((aligned(16)));;
98
TPOS.Element.One = lseq + 1;
99
TPOS.Element.Two = 0;
100
TPOS.Element.B = iseq + 1;
101
TPOS.Element.dummy = 0;
102
103
for (unsigned int iprf=0; iprf<(unsigned int) (N1-1); ++iprf) {
104
_mm_store_si128(&(iom[iprf].xmm), TPOS.xmm);
105
_mm_store_si128(&(ioi[iprf].xmm), TPOS.xmm);
106
}
107
108
if (N1 == 0) { fputs("BUG HERE N1 is NULL\n", stderr); exit(1);}
109
ioi[N1-1].xmm = TPOS.xmm; // Warning N1 > 0 ?
110
TPOS.Element.One = TPOS.Element.B;
111
TPOS.Element.Two = TPOS.Element.B;
112
iom[N1-1].xmm = TPOS.xmm;
113
114
for (unsigned int iprf=N1; iprf<(unsigned int)N2; ++iprf) {
115
_mm_store_si128(&(iom[iprf].xmm), TPOS.xmm);
116
_mm_store_si128(&(ioi[iprf].xmm), TPOS.xmm);
117
}
118
119
TPOS.Element.One = lseq + 1;
120
TPOS.Element.Two = 0;
121
122
for (unsigned int iprf=(unsigned int) N2; iprf<=(unsigned int)prf->Length; ++iprf) {
123
_mm_store_si128(&(iom[iprf].xmm), TPOS.xmm);
124
_mm_store_si128(&(ioi[iprf].xmm), TPOS.xmm);
125
}
126
}
127
128
static void nextR(const struct Profile * const restrict prf, const unsigned char * const restrict Sequence,
129
const int iseq, union lScores * const restrict iop, union Positions * const restrict iom,
130
union Positions * const restrict ioi,const int lseq, struct Alignment * const restrict alignment,
131
struct Minima * const restrict ifer, const _Bool Lock, const size_t N1, const size_t N2)
132
{
133
#ifdef XALIP_DEBUG
134
fprintf(stdout,"NextR called with iseq %i\n",iseq);
135
#endif
136
// Initialization
137
const unsigned int In = iseq + 1;
138
// WARNING: Fortran uses a 1 based index for sequence
139
const unsigned int SequenceIndex = (unsigned int) Sequence[iseq-1];
140
141
if ( iseq >= lseq) {
142
fputs("nextR_last should have been called\n", stderr);
143
exit(1);
144
}
145
146
// Disable match and insert vertices of protected region
147
if (Lock) {
148
iop[N1-1].Element[MATCH] = NLOW;
149
for (size_t iprf=N1; iprf<N2; ++iprf) {
150
iop[iprf].Element[MATCH] = NLOW;
151
iop[iprf].Element[INSERTION] = NLOW;
152
}
153
}
154
////////////////////////////////////////////////////////////////////////////////////////////
155
// Profile position 0
156
////////////////////////////////////////////////////////////////////////////////////////////
157
158
// Save previous match position
159
int Kopm = iop[0].Element[MATCH];
160
union Positions Kpos __attribute__((aligned(16)));
161
Kpos.xmm = iom[0].xmm;
162
163
const union Positions TEMPpos __attribute__((aligned(16))) = { lseq+1, 0, In, 0 };
164
165
const union Positions * restrict PTRpos[4];
166
union Positions Kiod;
167
PTRpos[0] = &TEMPpos;
168
PTRpos[1] = &Kpos;
169
PTRpos[3] = &Kiod;
170
171
// Get pointers to score data
172
const TransitionScores * const restrict Transitions = prf->Scores.Insertion.Transitions;
173
const short int * restrict Insertion = prf->Scores.Insertion.Alphabet;
174
const size_t AlignStep = prf->Scores.Insertion.AlignStep;
175
176
int Ki = iop[0].Element[INSERTION] + (int) Insertion[SequenceIndex];
177
178
// Match position
179
int Ji = Ki + (int) Transitions[0].From[INSERTION].To[MATCH];
180
int itmp = (int) Transitions[0].From[EXTRA].To[MATCH];
181
if ( Ji > itmp) {
182
iop[0].Element[MATCH] = Ji;
183
iom[0].xmm = ioi[0].xmm;
184
} else {
185
iop[0].Element[MATCH] = itmp;
186
iom[0].xmm = TEMPpos.xmm;
187
}
188
189
// Deletion position
190
int Kopd;
191
Ji = Ki + (int) Transitions[0].From[INSERTION].To[DELETION];
192
itmp = (int) Transitions[0].From[EXTRA].To[DELETION];
193
if ( Ji > itmp ) {
194
Kopd = Ji;
195
Kiod.xmm = ioi[0].xmm;
196
} else {
197
Kopd = itmp;
198
Kiod.xmm = TEMPpos.xmm;
199
}
200
201
// Insertion position
202
Ji = Ki + (int) Transitions[0].From[INSERTION].To[INSERTION];
203
itmp = (int) Transitions[0].From[EXTRA].To[INSERTION];
204
if ( Ji > itmp ) {
205
iop[0].Element[INSERTION] = Ji;
206
} else {
207
iop[0].Element[INSERTION] = itmp;
208
ioi[0].xmm = TEMPpos.xmm;
209
}
210
211
// Initialize minima
212
ifer->_a = iseq;
213
ifer->_b = iseq;
214
itmp = MIN(ioi[0].Element.B, iom[0].Element.B);
215
ifer->_c = MIN(itmp, Kiod.Element.B);
216
217
// Initialize alignment
218
union Positions Fpos __attribute__((aligned(16)));;
219
Fpos.Element.One = alignment->JAL1;
220
Fpos.Element.Two = alignment->JAL2;
221
Fpos.Element.B = alignment->JALB;
222
Fpos.Element.dummy = 0;
223
224
////////////////////////////////////////////////////////////////////////////////////////////
225
// Loop through rest of profile
226
////////////////////////////////////////////////////////////////////////////////////////////
227
const short int * restrict Match = prf->Scores.Match.Alphabet;
228
Insertion += AlignStep;
229
230
for (int iprf=1; iprf<=prf->Length; ++iprf) {
231
/////////////////////////////////////////////////////////////////////////////////////////
232
// Match
233
const register int KM = Kopm + (int) Match[SequenceIndex];
234
Kopm = iop[iprf].Element[MATCH];
235
236
__m128i __KM = _mm_set1_epi32(KM);
237
// Load Transitions
238
__m128i __TransitionsM = _mm_loadl_epi64((__m128i*) &(Transitions[iprf].From[MATCH]));
239
// Convert signed WORD into signed DWORD
240
__TransitionsM = _mm_cvtepi16_epi32(__TransitionsM);
241
// Add KM to Transitions
242
__TransitionsM = _mm_add_epi32(__TransitionsM, __KM);
243
244
245
/////////////////////////////////////////////////////////////////////////////////////////
246
// Insertion
247
const register int KI = iop[iprf].Element[INSERTION] + (int) Insertion[SequenceIndex];
248
// one could move on the seq index instead
249
250
__m128i __KI = _mm_set1_epi32(KI);
251
// Load Transitions
252
__m128i __TransitionsI = _mm_loadl_epi64((__m128i*) &(Transitions[iprf].From[INSERTION]));
253
// Convert signed WORD into signed DWORD
254
__TransitionsI = _mm_cvtepi16_epi32(__TransitionsI);
255
// Add KM to Transitions
256
__TransitionsI = _mm_add_epi32(__TransitionsI, __KI);
257
258
/////////////////////////////////////////////////////////////////////////////////////////
259
// Deletion
260
const register int KD = Kopd + (int) Match[_D];
261
262
__m128i __KD = _mm_set1_epi32(KD);
263
// Load Transitions
264
__m128i __TransitionsD = _mm_loadl_epi64((__m128i*) &(Transitions[iprf].From[DELETION]));
265
// Convert signed WORD into signed DWORD
266
__TransitionsD = _mm_cvtepi16_epi32(__TransitionsD);
267
// Add KM to Transitions
268
__TransitionsD = _mm_add_epi32(__TransitionsD, __KD);
269
270
/////////////////////////////////////////////////////////////////////////////////////////
271
// Extensions
272
// Load Transitions
273
__m128i __TransitionsX = _mm_loadl_epi64((__m128i*) &(Transitions[iprf].From[EXTRA]));
274
// Convert signed WORD into signed DWORD
275
__TransitionsX = _mm_cvtepi16_epi32(__TransitionsX);
276
277
// Insert NLOW into Extension vector -> done in io.c
278
//__TransitionsX = _mm_insert_epi32(__TransitionsX, NLOW, DUMMY);
279
280
// Move to next profile index
281
Match += AlignStep;
282
Insertion += AlignStep;
283
284
/////////////////////////////////////////////////////////////////////////////////////////
285
// Compute Maxima (Fortran Nstep function)
286
__m128i __Index = (__m128i) _mm_setzero_si128();
287
__m128i __three = _mm_set_epi32(3,3,3,3);
288
289
__m128i __Mask = _mm_cmpgt_epi32(__TransitionsD, __TransitionsX);
290
__TransitionsX = _mm_blendv_epi8(__TransitionsX, __TransitionsD, __Mask);
291
__Index = _mm_blendv_epi8(__Index, __three, __Mask);
292
293
__m128i __One = _mm_set_epi32(1,1,1,1);
294
__Mask = _mm_cmpgt_epi32(__TransitionsM, __TransitionsX);
295
__TransitionsX = _mm_blendv_epi8(__TransitionsX, __TransitionsM, __Mask);
296
__Index = _mm_blendv_epi8(__Index, __One, __Mask);
297
298
__m128i __Two = _mm_add_epi32(__One, __One);
299
__Mask = _mm_cmpgt_epi32(__TransitionsI, __TransitionsX);
300
__TransitionsX = _mm_blendv_epi8(__TransitionsX, __TransitionsI, __Mask);
301
__Index = _mm_blendv_epi8(__Index, __Two, __Mask);
302
303
// Set new data
304
iop[iprf].xmm = __TransitionsX;
305
//StoreMatchInsertion((__m64*) &iop[iprf],(__m128) __TransitionsX);
306
Kopd = _mm_extract_epi32(__TransitionsX, DELETION);
307
308
/////////////////////////////////////////////////////////////////////////////////////////
309
// Check for new maxima
310
const int KE = _mm_extract_epi32(__TransitionsX, DUMMY);
311
if (KE > alignment->JALS) {
312
alignment->JALS = KE;
313
alignment->JALE = iseq;
314
Fpos.xmm = ioi[iprf].xmm;
315
const unsigned int Id = (unsigned int) _mm_extract_epi32(__Index, DUMMY);
316
if (Id == 1) { // KM is max
317
Fpos.xmm = Kpos.xmm;
318
} else if (Id == 3) { // KD is max
319
Fpos.xmm = Kiod.xmm;
320
}
321
}
322
#ifdef XALIP_DEBUG
323
printf("FPOS %8i %8i %8i\n", Fpos.Element.One, Fpos.Element.Two, Fpos.Element.B);
324
#endif
325
/////////////////////////////////////////////////////////////////////////////////////////
326
// Update alignment positions
327
union Positions Jpos __attribute__((aligned(16)));
328
Jpos.xmm = iom[iprf].xmm;
329
PTRpos[2] = &ioi[iprf];
330
331
const int NewM = _mm_extract_epi32(__Index, MATCH);
332
iom[iprf].xmm = PTRpos[NewM]->xmm;
333
334
const int NewD = _mm_extract_epi32(__Index, DELETION);
335
Kiod.xmm = PTRpos[NewD]->xmm;
336
337
const int NewI = _mm_extract_epi32(__Index, INSERTION);
338
ioi[iprf].xmm = PTRpos[NewI]->xmm;
339
340
Kpos.xmm = Jpos.xmm;
341
342
/////////////////////////////////////////////////////////////////////////////////////////
343
// Update minima
344
345
const int t1 = MIN(ioi[iprf].Element.One, iom[iprf].Element.One);
346
const int t2 = MIN(t1, Kiod.Element.One);
347
ifer->_a = MIN(ifer->_a, t2);
348
349
if (iprf > N1) {
350
ifer->_b = MIN(ifer->_b, t2);
351
}
352
const int t3 = MIN(ioi[iprf].Element.B, iom[iprf].Element.B);
353
const int t4 = MIN(t3, Kiod.Element.B);
354
ifer->_c = MIN(ifer->_c, t4);
355
}
356
357
////////////////////////////////////////////////////////////////////////////////////////////
358
// Epilogue
359
////////////////////////////////////////////////////////////////////////////////////////////
360
361
// Finish updating alignment positions
362
alignment->JAL1 = Fpos.Element.One;
363
alignment->JAL2 = Fpos.Element.Two;
364
alignment->JALB = Fpos.Element.B;
365
366
// Entry and exit from protected regions
367
iom[N1-1].Element.One = MIN(iom[N1-1].Element.One, In);
368
iom[N1-1].Element.Two = In;
369
370
for (int iprf=N1; iprf<N2; ++iprf) {
371
iom[iprf].Element.One = MIN(iom[iprf].Element.One, In);
372
iom[iprf].Element.Two = In;
373
374
ioi[iprf].Element.One = MIN(ioi[iprf].Element.One, In);
375
ioi[iprf].Element.Two = In;
376
}
377
378
#ifdef XALIP_DEBUG
379
for (int iprf=0; iprf<=prf->Length; ++iprf) {
380
printf("NEXTR IOP %4.4i %15i %15i %15i\n", iprf, iop[iprf].Element[MATCH], iop[iprf].Element[INSERTION],
381
iop[iprf].Element[DELETION]);
382
printf("NEXTR IOM %4.4i %15i %15i %15i\n", iprf, iom[iprf].Element.One, iom[iprf].Element.Two,
383
iom[iprf].Element.B);
384
printf("NEXTR IOI %4.4i %15i %15i %15i\n", iprf, ioi[iprf].Element.One, ioi[iprf].Element.Two,
385
ioi[iprf].Element.B);
386
387
}
388
389
printf("NEXTR ALIGN %4i %4i %4i %4i %4i\n",
390
alignment->JAL1, alignment->JAL2, alignment->JALS, alignment->JALB, alignment->JALE);
391
#endif
392
}
393
394
static void nextR_last(const struct Profile * const restrict prf, const unsigned char * const restrict Sequence,
395
const int iseq, union lScores * const restrict iop, union Positions * const restrict iom,
396
union Positions * const restrict ioi,const int lseq, struct Alignment * const restrict alignment,
397
struct Minima * const restrict ifer, const _Bool Lock, const size_t N1, const size_t N2)
398
{
399
#ifdef XALIP_DEBUG
400
fprintf(stdout,"NextR_last called with iseq %i\n",iseq);
401
#endif
402
// Initialization
403
const unsigned int In = iseq + 1;
404
// WARNING: Fortran uses a 1 based index for sequence
405
const unsigned int SequenceIndex = (unsigned int) Sequence[iseq-1];
406
407
if ( iseq < lseq) {
408
fputs("nextR should have been called\n", stderr);
409
exit(1);
410
}
411
412
// Disable match and insert vertices of protected region
413
if (Lock) {
414
iop[N1-1].Element[MATCH] = NLOW;
415
for (int iprf=N1; iprf<N2; ++iprf) {
416
iop[iprf].Element[MATCH] = NLOW;
417
iop[iprf].Element[INSERTION] = NLOW;
418
}
419
}
420
////////////////////////////////////////////////////////////////////////////////////////////
421
// Profile position 0
422
////////////////////////////////////////////////////////////////////////////////////////////
423
424
// Save previous match position
425
int Kopm = iop[0].Element[MATCH];
426
union Positions Kpos __attribute__((aligned(16)));
427
Kpos.xmm = iom[0].xmm;
428
429
const union Positions TEMPpos __attribute__((aligned(16))) = { lseq+1, 0, In, 0 };
430
431
const union Positions * restrict PTRpos[4];
432
union Positions Kiod;
433
PTRpos[0] = &TEMPpos;
434
PTRpos[1] = &Kpos;
435
PTRpos[3] = &Kiod;
436
437
// Get pointers to score data
438
const TransitionScores * const restrict Transitions = prf->Scores.Insertion.Transitions;
439
const short int * restrict Insertion = prf->Scores.Insertion.Alphabet;
440
const size_t AlignStep = prf->Scores.Insertion.AlignStep;
441
442
int Ki = iop[0].Element[INSERTION] + (int) Insertion[SequenceIndex];
443
444
// Match position
445
int Ji = Ki + (int) Transitions[0].From[INSERTION].To[MATCH];
446
int itmp = (int) Transitions[0].From[EXTRA].To[MATCH];
447
if ( Ji > itmp) {
448
iop[0].Element[MATCH] = Ji;
449
iom[0].xmm = ioi[0].xmm;
450
} else {
451
iop[0].Element[MATCH] = itmp;
452
iom[0].xmm = TEMPpos.xmm;
453
}
454
455
// Deletion position
456
int Kopd;
457
Ji = Ki + (int) Transitions[0].From[INSERTION].To[DELETION];
458
itmp = (int) Transitions[0].From[EXTRA].To[DELETION];
459
if ( Ji > itmp ) {
460
Kopd = Ji;
461
Kiod.xmm = ioi[0].xmm;
462
} else {
463
Kopd = itmp;
464
Kiod.xmm = TEMPpos.xmm;
465
}
466
467
// Insertion position
468
Ji = Ki + (int) Transitions[0].From[INSERTION].To[INSERTION];
469
itmp = (int) Transitions[0].From[EXTRA].To[INSERTION];
470
if ( Ji > itmp ) {
471
iop[0].Element[INSERTION] = Ji;
472
} else {
473
iop[0].Element[INSERTION] = itmp;
474
ioi[0].xmm = TEMPpos.xmm;
475
}
476
477
// Initialize minima
478
ifer->_a = iseq;
479
ifer->_b = iseq;
480
itmp = MIN(ioi[0].Element.B, iom[0].Element.B);
481
ifer->_c = MIN(itmp, Kiod.Element.B);
482
483
// Initialize alignment
484
union Positions Fpos __attribute__((aligned(16)));
485
Fpos.Element.One = alignment->JAL1;
486
Fpos.Element.Two = alignment->JAL2;
487
Fpos.Element.B = alignment->JALB;
488
Fpos.Element.dummy = 0;
489
490
////////////////////////////////////////////////////////////////////////////////////////////
491
// Loop through rest of profile
492
////////////////////////////////////////////////////////////////////////////////////////////
493
const short int * restrict Match = prf->Scores.Match.Alphabet;
494
const ScoreTuple * const restrict LastProteinSequence = prf->Scores.Insertion.LastSequenceProtein;
495
Insertion += AlignStep;
496
497
for (int iprf=1; iprf<=prf->Length; ++iprf) {
498
/////////////////////////////////////////////////////////////////////////////////////////
499
// Match
500
const register int KM = Kopm + (int) Match[SequenceIndex];
501
Kopm = iop[iprf].Element[MATCH];
502
503
__m128i __KM = _mm_set1_epi32(KM);
504
// Load Transitions
505
__m128i __TransitionsM = _mm_loadl_epi64((__m128i*) &(Transitions[iprf].From[MATCH]));
506
// Convert signed WORD into signed DWORD
507
__TransitionsM = _mm_cvtepi16_epi32(__TransitionsM);
508
// Insert LastProteinSequence
509
__TransitionsM = _mm_insert_epi32(__TransitionsM, (int) LastProteinSequence[iprf].From[MATCH], DUMMY);
510
// Add KM to Transitions
511
__TransitionsM = _mm_add_epi32(__TransitionsM, __KM);
512
513
514
/////////////////////////////////////////////////////////////////////////////////////////
515
// Insertion
516
const register int KI = iop[iprf].Element[INSERTION] + (int) Insertion[SequenceIndex];
517
// one could move on the seq index instead
518
519
__m128i __KI = _mm_set1_epi32(KI);
520
// Load Transitions
521
__m128i __TransitionsI = _mm_loadl_epi64((__m128i*) &(Transitions[iprf].From[INSERTION]));
522
// Convert signed WORD into signed DWORD
523
__TransitionsI = _mm_cvtepi16_epi32(__TransitionsI);
524
// Insert LastProteinSequence
525
__TransitionsI = _mm_insert_epi32(__TransitionsI, (int) LastProteinSequence[iprf].From[INSERTION], DUMMY);
526
// Add KM to Transitions
527
__TransitionsI = _mm_add_epi32(__TransitionsI, __KI);
528
529
/////////////////////////////////////////////////////////////////////////////////////////
530
// Deletion
531
const register int KD = Kopd + (int) Match[_D];
532
533
__m128i __KD = _mm_set1_epi32(KD);
534
// Load Transitions
535
__m128i __TransitionsD = _mm_loadl_epi64((__m128i*) &(Transitions[iprf].From[DELETION]));
536
// Convert signed WORD into signed DWORD
537
__TransitionsD = _mm_cvtepi16_epi32(__TransitionsD);
538
// Insert LastProteinSequence
539
__TransitionsD = _mm_insert_epi32(__TransitionsD, (int) LastProteinSequence[iprf].From[DELETION], DUMMY);
540
// Add KM to Transitions
541
__TransitionsD = _mm_add_epi32(__TransitionsD, __KD);
542
543
/////////////////////////////////////////////////////////////////////////////////////////
544
// Extensions
545
// Load Transitions
546
__m128i __TransitionsX = _mm_loadl_epi64((__m128i*) &(Transitions[iprf].From[EXTRA]));
547
// Convert signed WORD into signed DWORD
548
__TransitionsX = _mm_cvtepi16_epi32(__TransitionsX);
549
550
// Insert NLOW into Extension vector -> done in io.c
551
//__TransitionsX = _mm_insert_epi32(__TransitionsX, NLOW, DUMMY);
552
553
// Move to next profile index
554
Match += AlignStep;
555
Insertion += AlignStep;
556
557
/////////////////////////////////////////////////////////////////////////////////////////
558
// Compute Maxima (Fortran Nstep function)
559
__m128i __Index = (__m128i) _mm_setzero_ps();
560
__m128i __three = _mm_set_epi32(3,3,3,3);
561
562
__m128i __Mask = _mm_cmpgt_epi32(__TransitionsD, __TransitionsX);
563
__TransitionsX = _mm_blendv_epi8(__TransitionsX, __TransitionsD, __Mask);
564
__Index = _mm_blendv_epi8(__Index, __three, __Mask);
565
566
__m128i __One = _mm_set_epi32(1,1,1,1);
567
__Mask = _mm_cmpgt_epi32(__TransitionsM, __TransitionsX);
568
__TransitionsX = _mm_blendv_epi8(__TransitionsX, __TransitionsM, __Mask);
569
__Index = _mm_blendv_epi8(__Index, __One, __Mask);
570
571
__m128i __Two = _mm_add_epi32(__One, __One);
572
__Mask = _mm_cmpgt_epi32(__TransitionsI, __TransitionsX);
573
__TransitionsX = _mm_blendv_epi8(__TransitionsX, __TransitionsI, __Mask);
574
__Index = _mm_blendv_epi8(__Index, __Two, __Mask);
575
576
// Set new data
577
iop[iprf].xmm = __TransitionsX;
578
//_mm_storel_pi((__m64*) &iop[iprf],(__m128) __TransitionsX);
579
Kopd = _mm_extract_epi32(__TransitionsX, DELETION);
580
581
/////////////////////////////////////////////////////////////////////////////////////////
582
// Check for new maxima
583
const int KE = _mm_extract_epi32(__TransitionsX, DUMMY);
584
if (KE > alignment->JALS) {
585
alignment->JALS = KE;
586
alignment->JALE = iseq;
587
Fpos.xmm = ioi[iprf].xmm;
588
const unsigned int Id = (unsigned int) _mm_extract_epi32(__Index, DUMMY);
589
if (Id == 1) { // KM is max
590
Fpos.xmm = Kpos.xmm;
591
} else if (Id == 3) { // KD is max
592
Fpos.xmm = Kiod.xmm;
593
}
594
}
595
#ifdef XALIP_DEBUG
596
printf("FPOS %8i %8i %8i\n", Fpos.Element.One, Fpos.Element.Two, Fpos.Element.B);
597
#endif
598
/////////////////////////////////////////////////////////////////////////////////////////
599
// Update alignment positions
600
union Positions Jpos __attribute__((aligned(16)));
601
Jpos.xmm = iom[iprf].xmm;
602
PTRpos[2] = &ioi[iprf];
603
604
const int NewM = _mm_extract_epi32(__Index, MATCH);
605
iom[iprf].xmm = PTRpos[NewM]->xmm;
606
607
const int NewD = _mm_extract_epi32(__Index, DELETION);
608
Kiod.xmm = PTRpos[NewD]->xmm;
609
610
const int NewI = _mm_extract_epi32(__Index, INSERTION);
611
ioi[iprf].xmm = PTRpos[NewI]->xmm;
612
613
Kpos.xmm = Jpos.xmm;
614
615
/////////////////////////////////////////////////////////////////////////////////////////
616
// Update minima
617
618
const int t1 = MIN(ioi[iprf].Element.One, iom[iprf].Element.One);
619
const int t2 = MIN(t1, Kiod.Element.One);
620
ifer->_a = MIN(ifer->_a, t2);
621
622
if (iprf > N1) {
623
ifer->_b = MIN(ifer->_b, t2);
624
}
625
const int t3 = MIN(ioi[iprf].Element.B, iom[iprf].Element.B);
626
const int t4 = MIN(t3, Kiod.Element.B);
627
ifer->_c = MIN(ifer->_c, t4);
628
}
629
630
////////////////////////////////////////////////////////////////////////////////////////////
631
// Epilogue
632
////////////////////////////////////////////////////////////////////////////////////////////
633
634
// Finish updating alignment positions
635
alignment->JAL1 = Fpos.Element.One;
636
alignment->JAL2 = Fpos.Element.Two;
637
alignment->JALB = Fpos.Element.B;
638
639
// Entry and exit from protected regions
640
iom[N1-1].Element.One = MIN(iom[N1-1].Element.One, In);
641
iom[N1-1].Element.Two = In;
642
643
for (int iprf=N1; iprf<N2; ++iprf) {
644
iom[iprf].Element.One = MIN(iom[iprf].Element.One, In);
645
iom[iprf].Element.Two = In;
646
647
ioi[iprf].Element.One = MIN(ioi[iprf].Element.One, In);
648
ioi[iprf].Element.Two = In;
649
}
650
651
#ifdef XALIP_DEBUG
652
for (int iprf=0; iprf<=prf->Length; ++iprf) {
653
printf("NEXTR IOP %4.4i %15i %15i %15i\n", iprf, iop[iprf].Element[MATCH], iop[iprf].Element[INSERTION],
654
iop[iprf].Element[DELETION]);
655
printf("NEXTR IOM %4.4i %15i %15i %15i\n", iprf, iom[iprf].Element.One, iom[iprf].Element.Two,
656
iom[iprf].Element.B);
657
printf("NEXTR IOI %4.4i %15i %15i %15i\n", iprf, ioi[iprf].Element.One, ioi[iprf].Element.Two,
658
ioi[iprf].Element.B);
659
660
}
661
662
printf("NEXTR ALIGN %4i %4i %4i %4i %4i\n",
663
alignment->JAL1, alignment->JAL2, alignment->JALS, alignment->JALB, alignment->JALE);
664
#endif
665
}
666
667
int xalip_sse41( const struct Profile * const restrict prf, const unsigned char * const restrict Sequence,
668
union lScores * const restrict iop, union Positions * const restrict iom,
669
union Positions * const restrict ioi, const size_t bseq, const size_t lseq,
670
struct Alignment * const restrict alignment,
671
_Bool * const restrict Lock, const size_t N1, const size_t N2, const _Bool Lopt,
672
const int kcut, const size_t NALI)
673
{
674
int iseq;
675
////////////////////////////////////////////////////////////////////////////////////////////
676
// Prelogue
677
////////////////////////////////////////////////////////////////////////////////////////////
678
679
// Alignment list
680
int iopt = NLOW;
681
size_t nali = 0;
682
683
// Search control fields
684
int ibeg = bseq-1;
685
size_t jlcp = prf->Length / 2;
686
int nsca = 0;
687
688
// Stack Memory
689
struct Minima ifer __attribute__((aligned(16)));
690
691
////////////////////////////////////////////////////////////////////////////////////////////
692
// Two step forward one step backward loop
693
////////////////////////////////////////////////////////////////////////////////////////////
694
MajorLoop:
695
696
iseq = ibeg;
697
struct Alignment lAlignment;
698
lAlignment.JALS = NLOW;
699
lAlignment.JAL1 = 0;
700
lAlignment.JAL2 = 0;
701
lAlignment.JALB = 0;
702
lAlignment.JALE = 0;
703
704
// Initiate work array
705
InitR(iseq, N1, N2, bseq, lseq, iop, iom, ioi, prf);
706
#ifdef XALIP_DEBUG
707
for (size_t i=0; i<=prf->Length; ++i) {
708
printf("IOP %8i %8i %8i %8i\n",
709
iop[i].Element[MATCH], iop[i].Element[INSERTION], iop[i].Element[DELETION],
710
iop[i].Element[DUMMY]);
711
printf("IOM %8i %8i %8i %8i\n",
712
iom[i].Element.One, iom[i].Element.Two, iom[i].Element.B,
713
iom[i].Element.dummy);
714
printf("IOI %8i %8i %8i %8i\n",
715
ioi[i].Element.One, ioi[i].Element.Two, ioi[i].Element.B,
716
ioi[i].Element.dummy);
717
}
718
#endif
719
720
// Initiate search control values
721
int ilcp = iseq;
722
int ifcp = iseq+1;
723
int nlcp = ilcp + jlcp;
724
725
// Move one sequence position forward
726
++iseq;
727
728
////////////////////////////////////////////////////////////////////////////////////////////
729
// Loop over sequence positions
730
////////////////////////////////////////////////////////////////////////////////////////////
731
SeqPosLoop:
732
{
733
++nsca;
734
/////////////////////////////////////////////////////////////////////////////////////////
735
// Update work array
736
if (iseq < lseq ) {
737
#pragma noinline
738
nextR(prf, Sequence, iseq, iop, iom, ioi, lseq, &lAlignment,&ifer, Lock[iseq], N1, N2);
739
} else {
740
#pragma noinline
741
nextR_last(prf, Sequence, iseq, iop, iom, ioi, lseq, &lAlignment,&ifer, Lock[iseq], N1, N2);
742
}
743
/////////////////////////////////////////////////////////////////////////////////////////
744
// If Match found
745
if (lAlignment.JALS >= kcut) {
746
// Determine firdst entry of current row
747
if ( (ifer._a > lAlignment.JAL2) || (iseq == lseq) ) {
748
// Fill in missing alignment data
749
lAlignment.JAL1 = lAlignment.JAL1 == 0 ? lAlignment.JALB : lAlignment.JAL1;
750
lAlignment.JAL2 = lAlignment.JAL2 == 0 ? iseq : lAlignment.JAL2;
751
752
// Check for errors
753
if (lAlignment.JAL2 < lAlignment.JAL1) {
754
fprintf(stderr, "Error: Illegal alignment found in alignment %lu - no list produced.\n"
755
" Alignement should be from %i to %i\n",1+nali, lAlignment.JAL1, lAlignment.JAL2 );
756
return -1;
757
}
758
759
if (++nali > NALI) {
760
fputs("Warning: Too many alignments found - list may be incomplete.\n", stderr);
761
return -2;
762
}
763
764
// Accept alignment
765
struct Alignment * pAlignment = &alignment[nali];
766
pAlignment->JALS = lAlignment.JALS;
767
pAlignment->JALB = lAlignment.JALB;
768
pAlignment->JAL1 = lAlignment.JAL1;
769
pAlignment->JAL2 = lAlignment.JAL2;
770
pAlignment->JALE = lAlignment.JALE;
771
#ifdef XALIP_DEBUG
772
printf("XALIP ALIGN %lu %4.4i %4.4i %4.4i %4.4i %4.4i\n",nali,
773
lAlignment.JAL1, lAlignment.JAL2, lAlignment.JALS, lAlignment.JALB, lAlignment.JALE);
774
#endif
775
// Protect sequence region
776
for (int jseq=lAlignment.JAL1; jseq<=lAlignment.JAL2; ++jseq) {
777
Lock[jseq] = true;
778
}
779
780
// Exit if only searching for optimal alignment
781
if (nali>0 && Lopt)
782
return nali;
783
else
784
goto MajorLoop;
785
} else {
786
if ( ++iseq <= lseq ) goto SeqPosLoop;
787
}
788
} else {
789
// Have we reached next check point ?
790
if (iseq >= nlcp) {
791
// Determine firdst entry of current row
792
if (ifer._b >= ilcp) {
793
ibeg = ifcp - 1;
794
ifcp = ifer._c;
795
ilcp = iseq;
796
}
797
798
// Calculate next check point
799
nlcp += jlcp;
800
}
801
802
// Move one sequence position forward
803
if ( ++iseq <= lseq ) goto SeqPosLoop;
804
}
805
}
806
return nali;
807
}
808
Loggerhead is a web-based interface for Breezy
Version: 2.0.1