13
\$\begingroup\$

Arrays of double can be zero-initialized with the { 0 } initializer, which works even on systems that have a binary representation of 0.0 that is not all bits zero. While the IEEE 754 representation of floating point zero is all bits zero, the C Standard makes no guarantees about the representation of floating point numbers.

Large dynamically allocated arrays of integer types can be zero-initialized by using calloc(), but this method is not portable for floating point types, since calloc() initializes all bits to zero. Assigning a value of 0.0 to each individual double in the segment using a loop seems grossly inefficient.

I wrote the code for the function calloc_d() to solve this problem. The function takes the size_t argument nmemb, and returns a zero-initialized segment of memory large enough to hold nmemb doubles. If there is an allocation error, the function returns NULL. The caller is of course responsible for deallocation.

After using malloc() to allocate the required amount of memory, the first bytes are assigned the value of 0.0. Then memmove() is used to copy the bytes of the first zero, then the bytes of the first two zeros, and so on until at least half of the memory has been initialized. Finally, the remaining memory is initialized to 0.0 by copying bytes from the first half.

The code works, and has no obvious issues, as far as I can tell. I am interested in any comments about shortcomings of this approach, or shortcomings of my code, and suggestions for improvement. I would also be interested in suggestions for alternative, possibly more efficient methods.

I have included the function in a working program below. MAXCOUNT doubles are allocated for, initialized, and displayed. The value of MAXCOUNT is currently set to 100, but I have tested it for values up to 1000000.

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define MAXCOUNT 100      // number of doubles to allocate for

double * calloc_d(size_t nmemb);

int main(void)
{
    double *dbl_store = calloc_d(MAXCOUNT);
    size_t i;

    for(i = 0; i < MAXCOUNT; i++) {
        if (i % 10 == 0) putchar('\n');
        printf("%-8.4f", dbl_store[i]);
    }
    putchar('\n');

    free(dbl_store);

    return 0;
}

double * calloc_d(size_t nmemb)
{
    double *ret;
    double *next;         // pointer to beginning of uninitialized segment
    size_t alloc_bytes = sizeof(*ret) * nmemb;
    size_t init_sz;       // size of initialized segment

    ret = malloc(alloc_bytes);
    next = ret;

    if (ret) {
        ret[0] = 0.0;
        init_sz = sizeof(*ret);
        ++next;

        while (init_sz < (alloc_bytes + sizeof(*ret)) / 2) {
            memmove(next, ret, init_sz);
            init_sz *= 2;
            next = ret + init_sz / sizeof(*ret);    
        }
        memmove(next, ret, alloc_bytes - init_sz);
    }

    return ret;
}
\$\endgroup\$
10
  • 4
    \$\begingroup\$ How long does a naive solution take? Is it actually making a significant difference compared to the runtime of the rest of your code? These questions should be answered before engaging in a serious attempt to optimise this. \$\endgroup\$ Commented Dec 19, 2016 at 9:45
  • 4
    \$\begingroup\$ I strongly suspect that the compiler will optimise the naive loop to an assembly rep stosq (or its equivalent) on most platforms. "Repeat store string <size>" is a highly optimised in the hardware, and I doubt it can be beaten by anything else for this purpose. It will also never be slower than rep movsq(which I believe is usually how memcpy is implemented) - in fact stosq is almost always faster than movsq because it eliminates the need to read memory while writing to it. \$\endgroup\$ Commented Dec 19, 2016 at 11:06
  • 3
    \$\begingroup\$ Don't use memmove() when you can be confident that memcpy() will suffice. That is indeed the case here. \$\endgroup\$
    – PellMel
    Commented Dec 19, 2016 at 14:47
  • \$\begingroup\$ @PellMel: Why not use memmove()? The speed difference compared to memcpy() is likely to be very small, and it adds some robustness to the code. \$\endgroup\$ Commented Dec 19, 2016 at 15:24
  • \$\begingroup\$ @ThomasPadron-McCarthy, I said to use memmove() "when you can be confident that memcpy() will suffice". Robustness is not then a factor. Why in that case would one choose the function likely to be slower -- by a factor that will vary, and which you cannot be certain will be small? \$\endgroup\$
    – PellMel
    Commented Dec 19, 2016 at 16:33

5 Answers 5

28
\$\begingroup\$

So my first hunch when I read your question is:

Okay this guy is making a complicated solution to a simple problem based on a hunch that the trivial solution is "grossly inefficient" without any measurements to back up said claim. I suspect OPs solution is slower than the trivial solution.

So I set out test this, I pit your fancy implementation against my trivial implementation:

double * calloc_d2(size_t nmemb){
    double* ret = malloc(sizeof(double)*nmemb);
    double* end = ret + nmemb;
    
    if(ret){
        double* i = ret;
        while(i != end){
            *i = 0.0;
            ++i;
        }
    }
    return ret;
}

I tested this here on ideone.com (I think it is using clang) which shows that this simple bit of code is consistently around 3-5% faster than your fancy code.

So I have two pieces of advice

  1. Express your intent as clearly as possible, and let the compiler worry about generating high performance code. Compilers are pretty smart...
  2. Always, always, always measure.

Addendum/Edit

Seeing as this answer has gotten quite some attention I've gone ahead and dug deeper.

Using godbolt compiler explorer I compiled (-O3) and disassembled my trivial loop implementation.

Here is the disassembly on GCC 6.2:

calloc_d2(unsigned long):
        push    rbp
        push    rbx
        lea     rbx, [0+rdi*8]
        sub     rsp, 8
        mov     rdi, rbx
        call    malloc
        test    rax, rax
        mov     rbp, rax
        je      .L1
        lea     rax, [rax+rbx]
        cmp     rbp, rax
        je      .L1
        mov     rdx, rbx
        xor     esi, esi
        mov     rdi, rbp
        call    memset
.L1:
        add     rsp, 8
        mov     rax, rbp
        pop     rbx
        pop     rbp
        ret

Notice how the compiler has just gone "Oh, you want to zero that? Cool, I'll just replace your code with a memset". This goes to show that an optimising compiler is quite smart and is able to transform your code in astounding ways.

I also tried different versions of GCC and clang on ARM and MIPS and PowerPC, they all generate similar instructions which essentially boil down to call malloc then memset.

With the sole exception of ICC 17 (-O3 -march=native) which goes above and beyond:

calloc_d2(unsigned long):
        push      rbp                                           #5.1
        mov       rbp, rsp                                      #5.1
        and       rsp, -32                                      #5.1
        push      r15                                           #5.1
        push      rbx                                           #5.1
        sub       rsp, 16                                       #5.1
        mov       rbx, rdi                                      #5.1
        lea       rdi, QWORD PTR [rbx*8]                        #6.50
        call      malloc         
        -- snip --
        call      __intel_avx_rep_memset                        #12.8
        -- snip --

Notice that it calls into __intel_avx_rep_memset instead, this is a 512 bit wide AVX optimized memset that will only function on modern CPUs (it did this because I passed -march=native asking it to throw backward compatibility out the window). There is some branching code at the end which deals with cases where nmemb isn't evenly divisible by 512/64=8.

But here is the key take away, we clearly expressed our intent with the trivial loop. The compiler could then replaced the trivial loop (after deducing it's functionality) with an (512 bit) AVX optimised memset written by Intel's engineers specifically for modern Intel CPUs.

There is this thing called the "as-if rule" in C/C++. Which pretty much says:

As long as the program does the same reads and writes to volatile memory in the same order as-if the program was executed according to the language specification, then the program may do ANY transformation of the code it pleases. Including reordering instructions, reordering stores and loads to non-volatile memory, or removing instructions all together.

Just to show you how awesome some compilers can be look at what clang does here:

auto stupid(){
   auto x = new int[30];
  for(int i = 0; i < 30; i++)
    x[i] = i;
  
  auto sum = 0;
  for(int i = 0; i < 30; i++)
    sum += x[i]*x[i];
  delete [] x;
  return sum;
}

gets compiled into:

stupid():                             # @stupid()
        mov     eax, 8555
        ret

It deduces that the function always returns the same value and it completely removes all code and just returns a constant! It also removes the calls to new and delete as it realises they are pointless and don't affect the output of the function. If that is not impressive, I don't know what is. :)

So I would like to take this opportunity to reiterate my first piece of advice again:

Clearly express what it is you want your code to do. Let the compiler worry about generating good machine code, the compiler is likely much better at it than you are.

\$\endgroup\$
4
  • \$\begingroup\$ I quote: "Assigning a value of 0.0 to each individual double in the segment using a loop seems grossly inefficient." The values I used are in the link, but I'll post them here for reference I used 1 million doubles and I ran 10 trials. \$\endgroup\$
    – Emily L.
    Commented Dec 19, 2016 at 1:33
  • \$\begingroup\$ @DavidBowling Yes, memcpy/memset are often highly optimised for whatever architecture you're using. But note that they are library calls, you are making subroutine calls and returns where as my code will be inlined, loop unrolled and the compiler will work its magic. Also know that filling a memory region with a constant pattern is a very common thing to do and compilers can detect this and generate optimal code for it. \$\endgroup\$
    – Emily L.
    Commented Dec 19, 2016 at 1:43
  • \$\begingroup\$ @DavidBowling Think of it this way the goal is to set all bytes to some pattern. You do some work in addition to setting bytes (calculating the sizes of the memmove blocks, looping etc) so you have overhead in addition to setting all bytes. The compiler sees my code and just generates the byte setting without the overhead of your code. The compiler is good enough to generate byte setting code that is on par with what memcpy and memset do. \$\endgroup\$
    – Emily L.
    Commented Dec 19, 2016 at 1:56
  • 2
    \$\begingroup\$ @DavidBowling It doesn't make sense that memcpy() would be fast to initialize memory because memcpy() requires both reads and writes, and the reads are just wasted cycles. What would be fastest would be an optimized "memsetDouble" function that took a double sized bit pattern and set a block of memory to that pattern. The code that Emily L wrote probably is close to optimal if you use an optimizing compiler. But the fastest way would involve the same kinds of target specific optimizations that are used for memset() (which is possibly hand coded). \$\endgroup\$
    – JS1
    Commented Dec 19, 2016 at 3:14
12
\$\begingroup\$

In the 99.9% of cases, your architecture will use IEEE floats/doubles, and using calloc() is ideal for both simplicity and speed. calloc() is far more efficent when allocating large zeroed memory since the operating system will often simply allocate a blank zero page of memory from the page table and not need to actually set each word's value.

So you can test when this is safe to do by simply noticing if you're using IEEE floats (std::numeric_limits::is_iec55 will be TRUE) and if so, use calloc(). Otherwise fall back to Emily's method.

\$\endgroup\$
1
  • \$\begingroup\$ You can even check dynamically whether 0.0 has all bits zero, which may be true even for implementations that don't use IEEE-754 floating point representation. \$\endgroup\$
    – PellMel
    Commented Dec 19, 2016 at 14:46
3
\$\begingroup\$

This exercise does appear to be trying to optimize a niche part of code and suspect focus on speed improvements are more useful elsewhere. Yet given that, here are some ideas:

  1. Using memcpy()/memove() versus a simple loop is certainty only faster when nmemb exceeds some platform dependent threshold n. Consider steering code based on that value. Use a single simple loop to replicate n elements and then after that, double with memcpy()/memove(). This echoes the idea of @Emily L. to measure performance. With this approach, OP may have discovered the value of n is very large rendering the exercise moot.

  2. Expand functionality: I'd would expect calloc_d(size_t nmemb) to perform that same, with a ret[0] = 0.0;, ret[0] = 42.0; or ret[0] = DBL_MAX;. So using calloc_d(size_t nmemb, double initial_value) { ... ret[0] = initial_value; makes this function far more useful.

  3. memcpy() is usually a bit faster than memmove(), but that difference is more significant with smaller n and #1 above suggest only using memcpy()/memmove() when n is large. Do not expect significant difference in using either of these functions when combined with step #1.


In the end, I would go with usingcalloc(), which zeros the bit data and a while() loop for rare machines that do not encode 0.0 as all zero bits.

\$\endgroup\$
0
3
\$\begingroup\$

Footnote

After reading all the answers, I decided to see what would be the optimal way to initialize the array and how close the various answers came to the optimal time. Unfortunately, I'm working on a 32-bit target so I used float and uint32_t to test instead of double and uint64_t. But my code can easily be modified if someone wants to test it on a 64-bit target.

Essentially, all I did was to hand code a version that used rep stosl to do the memset. This turned out to be the fastest time:

Time to initialize array of 200 million floats:

memcpy (OP's function)   : 236 msec
loop (Emily L's function): 201 msec (187 msec with -sse3)
rep stosl                : 178 msec

Update: 64-bit target, sse optimizations

For the above, numbers, I found out that the gcc I was using did not turn on sse optimizations by default. When I turned those on, the loop version dropped to 187 msec.

I then found a 64-bit x86 machine to test on. What I found was that the simple loop actually outperformed the "rep stosq" method. I used gcc -O3 (and sse optimizations were on by default on this machine). Here are the results, using doubles this time:

memcpy (OP's function)   : 845 msec
loop (Emily L's function): 615 msec
rep stosq                : 682 msec

FYI, the assembly code for the simple loop was a 256-byte at a time unrolled loop that looked something like this:

loop_start:
    movupd %xmm1, -240(%rdx)
    movupd %xmm1, -224(%rdx)
    movupd %xmm1, -208(%rdx)
    ... snip ...
    movupd %xmm1, (%rdx)
    addq   $256, %rdx
    addq   $-32, %rcx
    jne    loop_start

So I guess trusting the compiler is probably the best thing to do.

Benchmarking code

Here is the program I used to generate the above benchmarks:

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>
#include <sys/time.h>

void memset_f0(float *array, float val, size_t nmemb)
{
    size_t init_sz;       // size of initialized segment
    float *next = array; // pointer to beginning of uninitialized segment
    size_t numBytes = sizeof(*array) * nmemb;

    array[0] = val;
    init_sz = sizeof(*array);
    ++next;

    while (init_sz < (numBytes + sizeof(*array)) / 2) {
        memmove(next, array, init_sz);
        init_sz *= 2;
        next = array + init_sz / sizeof(*array);
    }
    memmove(next, array, numBytes - init_sz);
}

void memset_f1(float *array, float val, size_t nmemb)
{
    size_t i;

    for (i=0;i<nmemb;i++)
        array[i] = val;
}

void memset_f2(float *array, float val, size_t nmemb)
{
    union {
        uint32_t intVal;
        float    fltVal;
    } u;

    u.fltVal = val;
    __asm__ __volatile__ ("cld\n"    \
                          "rep\n"    \
                          "stosl\n"  \
                          : "+a"(u.intVal), "+D"(array), "+c"(nmemb));
}

int getTimeMsec(void)
{
    struct timeval tv;

    gettimeofday(&tv, NULL);
    return (tv.tv_sec * 1000 + tv.tv_usec / 1000);
}

int main(int argc, char *argv[])
{
    int size = 200000000;
    int method = 0;
    float *array;
    int begin, end;

    if (argc > 1)
        method = atoi(argv[1]);

    if (argc > 2)
        size = atoi(argv[1]);

    array = calloc(size, sizeof(float));

    begin = getTimeMsec();
    switch (method) {
        case 0:
            memset_f0(array, (float) size, size);
            break;
        case 1:
            memset_f1(array, (float) size, size);
            break;
        case 2:
            memset_f2(array, (float) size, size);
            break;
    }
    end = getTimeMsec();
    printf("Elapsed time for size %d method %d: %d msec\n", size, method,
            end - begin);
    return 0;
}
\$\endgroup\$
1
  • \$\begingroup\$ As a note, the code you're testing isn't equivalent to OPs. As val is not a constant 0 the optimisations performed by the compiler will be different. Still an interesting test. An added benefit of the loop is that if the compiler improves so might the performance of the loop while the other two will have the same performance. \$\endgroup\$
    – Emily L.
    Commented Dec 21, 2016 at 19:32
2
\$\begingroup\$

If you are worried that 0.0 might not be all bits zero, then it's absolutely fine to write a function like calloc_d that will return a double* or float*, with all array elements set to zero. If you are also worried about putting in extra work that isn't needed, you can make that function just call calloc - if it breaks, there is a single function to change, and not possibly 100 calls all over the place, so you have already gained a lot.

If you implement filling memory with memcpy, don't copy too many bytes at a time. I'd say copy no more than 20,000 bytes at a time, so that the bytes you are copying are in the fastest cache.

\$\endgroup\$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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