Add an inline version of bgzf_read for small reads.

[jkbonfield]

The bgzf_read function is long and not something that's likely to get inlined even if it was in the header.

However most of the time our calls request a small amount of data and they fit within the buffer we've read, so we offer a static inline to do the memcpy when we can, falling back to the long function when we cannot.

In terms of CPU time it's not much difference, but the key thing is that it's often CPU time saved in a main thread given the bulk of the decode is often threaded. An example of test_view -B -@16

develop:

real    0m48.158s
user    6m2.901s
sys     0m28.134s

real    0m48.730s
user    6m3.707s
sys     0m28.473s

real    0m48.653s
user    6m5.215s
sys     0m28.637s

This PR:

real    0m41.731s
user    5m59.780s
sys     0m30.393s

real    0m41.945s
user    6m0.367s
sys     0m30.426s

So we can see it's a consistent win when threading, potentially 10-15% faster throughput depending on work loads.

[jkbonfield]

/usr/bin/time was reporting ~810% CPU utilisation for develop and ~940% for this PR, demonstrating the reduced time in main thread.

Commenting out the sanity checks in bam_read1 for if (c->n_cigar > 0) onwards increases that CPU utilisation to 1060% (1120% if I get rid of the tag2cigar call), leading to 37s real time. This shows that there is room for improvement, as these sort of checks ought to be executed in the multi-threaded decode instead. We thread the bgzf_read and not the bam_read.

This is why multi-threaded decoding of SAM and SAM.gz is sometimes more performant than BAM, despite the string->integer conversion overheads, as SAM parsing is entirely threadable.

[daviesrob]

This shadows the earlier x, and also adds an almost-duplicate block to fill out the BAM structure. A way to avoid both problems might be something like this:

uint8_t tmp[32], *x;
uint32_t new_l_data;

// ...

if (fp->block_length - fp->block_offset > 32) {
    x = (uint8_t *)fp->uncompressed_block + fp->block_offset;
    fp->block_offset += 32;
} else {
    x = tmp;
    if (bgzf_read(fp, x, 32) != 32) return -3;
}
c->tid = le_to_u32(x);
c->pos = le_to_i32(x+4);
// ... etc.

[jkbonfield]

I see what you mean by duplicate now. It's totally different code and not duplicated at all in style nor in implementation, but the end result is still the same thing. I've taken your suggestion and refactored it.

[daviesrob]

The intermediate variable could be avoided by:

c->l_qname = le_to_u8(x + 8);
c->qual = le_to_u8(x + 9);
c->bin = le_to_u16(x + 10);

I'm not sure if it would be any quicker though...

[jkbonfield]

I was mainly mirroring what it did before, but yes it may make things better (or worse). I'll experiment.

[daviesrob]

Similarly, could be:

c->n_cigar = le_to_u16(x + 12);
c->flag = le_to_u16(x + 14);

[daviesrob]

I think it would be less obscure to replace the individual bgzf_read() calls, as there aren't many of them. Also bgzf_read(fp, x, 32) should remain as-is as the amount of data available would be known to be too small in that case.

You could also merge bgzf_read(fp, &block_len, 4) into the one for the rest of the core header. That would eliminate an entire bgzf_read call, and I don't really see a good reason why it wouldn't work as long as we were careful about checking error cases and returning the right values.

[jkbonfield]

I'd previously tested merging the 4 byte and 32 byte read together and concluded it made no difference in speed (once I'd put the inline version in - it did do before that), and it meant more changes which would have made reviewing more work.

I can try it again if you want, but it's more changes to evaluate.

[jkbonfield]

Also, I don't get the comment about bgzf_read(fp, x, 32) remaining as-is. Initially it wasn't because it's small enough to be inlined and a fixed 32-byte read may still cause the memcpy to be optimised away (I don't know), and because we weren't subverting things and copying out the struct directly initially. However since then, it may not matter much either way as that code path will be very rare. I think it's probably an irrelevance at best.

[jkbonfield]

You also have to think about what's actually happening here. Why is it that bgzf_read(fp, x, 4) for example is faster with the bgzf_read_fast inline? Some of it is the removal of a function, but I think a lot of it is because when we inline a memcpy of constant 4 bytes it becomes a straight 32-bit load to register. When we call the function, the length is a variable so it has to call memcpy direct and that in turn has a lot of code in it for dealing with arbitrary length and is probably optimised for copying long things more than short things.

So anywhere we have a bgzf_read with a small-ish constant size, it's generally good practice to call the inline variant so the compiler can optimise away memcpy.

[daviesrob]

The optimisation can never work in the specific case of bgzf_read_small(fp, x, 32) due to the earlier line:

if (fp->block_length - fp->block_offset > 32) {

i.e. the same optimisation (plus a zero-copy one) happens elsewhere.

[jkbonfield]

Ah yes I see what you mean. It makes no difference to CPU from what I can see (maybe the compiler detects it as an impossibility and removes the code?), but yes it's pointless to go via the static inline.

Now we have the same code for both paths in parsing the x[] array/pointer there's no as much gain to be had on the fp->block_length - fp->block_offset > 32 check anyway (if we kept the "small" variant), however it does still seem to be about 5% faster on small small-read datasets so we may as well keep it.

[daviesrob]

Added a couple of extra tests...