4

The sliding diagonal vector contains 16 elements, each one an 8-bit unsigned integer.

Without SSE and a bit simplified it would have looked like this in C:

int width=1000000; // a big number
uint8_t matrix[width][16];
fill_matrix_with_interesting_values(&matrix);

for (int i=0; i < width - 16; ++i) {
  uint8_t diagonal_vector[16];
  for (int j=0; j<16; ++j) {
    diagonal_vector[j] = matrix[i+j][j];
  }
  do_something(&diagonal_vector);
}

but in my case I can only load column-wise (vertically) from the matrix with the _mm_load_si128 intrinsics function. The sliding diagonal vector is moving horizontally so I need to load 16 column vectors in advance and use one element from each of those column vectors to create the diagonal vector.

Is it possible to make a fast low-memory implementation for this with SSE?

Update Nov 14 2016: Providing some more details. In my case I read single-letter codes from a text file in FASTA format. Each letter represents a certain amino acid. Each amino acid has a specific column vector associated with it. That column vector is looked up from a constant table (a BLOSUM matrix). In C code it would look like this

while (uint8_t c = read_next_letter_from_file()) {
   column_vector = lookup_from_const_table(c)
   uint8_t diagonal_vector[16];
   ... rearrange the values from the latest column
       vectors into the diagonal_vector ...

   do_something(&diagonal_vector)
}
2
  • Do you have to store the data that way? Can you instead store it with the diagonals contiguous? Since one diagonal doesn't overlap with the next, this storage format can still compactly store all the data (with padding at the ends, or something). To get an arbitrary matrix[i][j] element from the diagonal storage, I guess you'd do something like diags[i-j][j]. So you don't need two copies of the data, but you can choose whether diagonal or column-wise access is SIMD-friendly. Nov 13, 2016 at 21:25
  • 1
    The data is actually not stored that way. I simplified the example too much. The column_vector is looked up from a constant table. I've now updated the question with more information. Nov 14, 2016 at 9:10

1 Answer 1

3

The implementation I will present only needs one column load per iteration. First we initialize some variables

const __m128i mask1=_mm_set_epi8(0,0,0,0,0,0,0,0,255,255,255,255,255,255,255,255);
const __m128i mask2=_mm_set_epi8(0,0,0,0,255,255,255,255,0,0,0,0,255,255,255,255);
const __m128i mask3=_mm_set_epi8(0,0,255,255,0,0,255,255,0,0,255,255,0,0,255,255);
const __m128i mask4=_mm_set_epi8(0,255,0,255,0,255,0,255,0,255,0,255,0,255,0,255);
__m128i v0, v1, v2, v3, v4, v5, v6, v7, v8, v9, v10, v11, v12, v13, v14, v15;

Then for each step the variable v_column_load is loaded with the next column.

v15 = v_column_load;
v7 = _mm_blendv_epi8(v7,v15,mask1);
v3 = _mm_blendv_epi8(v3,v7,mask2);
v1 = _mm_blendv_epi8(v1,v3,mask3);
v0 = _mm_blendv_epi8(v0,v1,mask4);
v_diagonal = v0;

In the next step the variable name numbers in v0, v1, v3, v7, v15 are incremented by 1 and adjusted to be in the range 0 to 15. In other words: newnumber = ( oldnumber + 1 ) modulo 16.

v0 = v_column_load;
v8 = _mm_blendv_epi8(v8,v0,mask1);
v4 = _mm_blendv_epi8(v4,v8,mask2);
v2 = _mm_blendv_epi8(v2,v4,mask3);
v1 = _mm_blendv_epi8(v1,v2,mask4);
v_diagonal = v1;

After 16 iterations the v_diagonal will start to contain the correct diagonal values.

Looking at mask1,mask2, mask3, mask4, we see a pattern that can be used to generalize this algorithm for other vector lengths (2^n).

For instance, for vector length 8, we would only need 3 masks and the iteration steps would look like this:

v7 = a a a a a a a a
v6 =
v5 =
v4 =
v3 =         a a a a
v2 =
v1 =             a a
v0 =               a

v0 = b b b b b b b b
v7 = a a a a a a a a
v6 =
v5 =
v4 =         b b b b
v3 =         a a a a
v2 =             b b
v1 =             a b

v1 = c c c c c c c c
v0 = b b b b b b b b
v7 = a a a a a a a a
v6 =
v5 =         c c c c
v4 =         b b b b
v3 =         a a c c
v2 =           a b c

v2 = d d d d d d d d
v1 = c c c c c c c c
v0 = b b b b b b b b
v7 = a a a a a a a a
v6 =         d d d d
v5 =         c c c c
v4 =         b b d d
v3 =         a a c d

v3 = e e e e e e e e
v2 = d d d d d d d d
v1 = c c c c c c c c
v0 = b b b b b b b b
v7 = a a a a e e e e
v6 =         d d d d
v5 =     a a c c e e
v4 =       a b b d a

v4 = f f f f f f f f
v3 = e e e e e e e e
v2 = d d d d d d d d
v1 = c c c c c c c c
v0 = b b b b f f f f
v7 = a a a a e e e e
v6 =     b b d d f f
v5 =     a b c d e f

v5 = g g g g g g g g
v4 = f f f f f f f f
v3 = e e e e e e e e
v2 = d d d d d d d d
v1 = c c c c g g g g
v0 = b b b b f f f f
v7 = a a c c e e g g
v6 =   a b c d e f g

v6 = h h h h h h h h
v5 = g g g g g g g g
v4 = f f f f f f f f
v3 = e e e e e e e e
v2 = d d d d h h h h
v1 = c c c c g g g g
v0 = b b d d f f h h
v7 = a b c d e f g h  <-- this vector now contains the diagonal

v7 = i i i i i i i i
v6 = h h h h h h h h
v5 = g g g g g g g g
v4 = f f f f f f f f
v3 = e e e e i i i i
v2 = d d d d h h h h
v1 = c c e e g g i i
v0 = b c d e f g h i  <-- this vector now contains the diagonal

v0 = j j j j j j j j
v7 = i i i i i i i i
v6 = h h h h h h h h
v5 = g g g g g g g g
v4 = f f f f j j j j
v3 = e e e e i i i i
v2 = d d f f h h j j
v1 = c d e f g h i j  <-- this vector now contains the diagonal

Sidenote: I discovered this way of loading a diagonal vector when I was working on an implementation of the Smith-Waterman algorithm. Some more information can be found on the old SourceForge project web page.

16
  • I've always been curious, does anyone try to store sequence data in memory with multiple elements packed into a byte? An ACGT nucleotide only requires two bits, so you can pack either 4 per byte, or a codon of 3 base pairs into a byte with padding. Unfortunately there are more than 16 amino-acids, so an AA model can't pack more than one into a byte without insane bit-string crap. You could pack three 5-bit AAs into 16 bits, but that's probably not worth the overhead for random-access. Nov 13, 2016 at 21:33
  • (I used to work with the phylogenetics group at Dalhousie, mostly on finding maximum-likelihood trees for evolutionary models, so the real work each iteration was on matrices of floats; evolutionary distances were just inputs. e.g. software like phyml, IIRC, not on finding matches in giant sequence databases.) Nov 13, 2016 at 21:35
  • Since your masks are compile-time constants, your code would run faster if you used constant blends when possible (i.e. for the steps where granularity >= 16 bits). PBLENDW runs faster than PBLENDVB (see agner.org/optimize). The first step could maybe also be done with MOVHPS or something, to merge a new upper half of a vector. Nov 13, 2016 at 21:39
  • If you need multiple consecutive diagonals, isn't there some redundancy you can take advantage of here, maybe with PALIGNR and PSHUFB? It seems like you're doing 16 * 4 PBLENDVB operations for every 16 diagonal vectors you produce. Nov 13, 2016 at 21:43
  • I would imagine that PUNPCKLBW to interleave bytes from (the low half of) two vectors might be useful, but it still requires one shuffle-port uop to combine a pair of registers, and may not lead to any fewer than 4 ops per result vector. Hmm. If you have AVX2, VPBLENDD can do immediate blends (with 32b granularity) with better throughput than PBLENDW on Intel Haswell, but it's still the same operation. Nov 13, 2016 at 21:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.