Zero-initializing large dynamically allocated arrays of double

Question

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;
}

Answer 1

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.

Answer 2

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.

Answer 3

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.

Answer 4

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;
}

Answer 5

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.