Initialize C array fast (to write)

Question

Is there any better solution to initialize C arrays (created on the heap) fast? Like we do with curly brackets

double** matrix_multiply(const double **l_matrix, const double **r_matrix);

foo() {
    double DCT_matrix[8][8] = {
        { 0.3536,  0.3536,  0.3536,  0.3536,  0.3536,  0.3536,  0.3536,  0.3536 },
        { 0.4904,  0.4157,  0.2778,  0.0975, -0.0975, -0.2778, -0.4157, -0.4904 },
        { 0.4619,  0.1913, -0.1913, -0.4619, -0.4619, -0.1913,  0.1913,  0.4619 },
        { 0.4157, -0.0975, -0.4904, -0.2778,  0.2778,  0.4904,  0.0975, -0.4157 },
        { 0.3536, -0.3536, -0.3536,  0.3536,  0.3536, -0.3536, -0.3536,  0.3536 },
        { 0.2778, -0.4904,  0.0975,  0.4157, -0.4157, -0.0975,  0.4904, -0.2778 },
        { 0.1913, -0.4619,  0.4619, -0.1913, -0.1913,  0.4619, -0.4619,  0.1913 },
        { 0.0975, -0.2778,  0.4157, -0.4904,  0.4904, -0.4157,  0.2778, -0.0975 }
    }; 

    const double other_matrix[8][8] = {
        {  26,  -5,  -5, -5, -5, -5, -5,   8 },
        {  64,  52,   8, 26, 26, 26,  8, -18 },
        { 126,  70,  26, 26, 52, 26, -5,  -5 },
        { 111,  52,   8, 52, 52, 38, -5,  -5 },
        {  52,  26,   8, 39, 38, 21,  8,   8 },
        {   0,   8,  -5,  8, 26, 52, 70,  26 },
        {  -5, -23, -18, 21,  8,  8, 52,  38 },
        { -18,   8,  -5, -5, -5,  8, 26,   8 }
    };

    matrix_multiply(DCT_matrix, other_matrix); // Segfault
}

Hmmm. If you have all the initialization values available at compile time, why do you need to allocate on heap? — Gerhardh, Jun 29 '22 at 08:59
@Gerhardh function takes `double**` it won't work with stack allocated arrays — Kanony, Jun 29 '22 at 09:01
What function? How does your `int` array match to `double**`. Without more context we cannot provide best solution. Also, I did not mention stack allocated memory. You could use a `static int arr[] = {...};` which would not place it on the stack, — Gerhardh, Jun 29 '22 at 09:04
@Kanony In other words, you have given us an [XY-problem](https://meta.stackexchange.com/questions/66377/what-is-the-xy-problem) ;) — klutt, Jun 29 '22 at 09:07
no it's the same analogy, the part I need is how to initialize it faster. Ok I will edit it — Kanony, Jun 29 '22 at 09:09
A 2D array is something very different than a pointer to a pointer. Allocating a 1D array on the heap as your initial code did, does not solve that. — Gerhardh, Jun 29 '22 at 09:18
@Gerhardh I could have created it using `malloc` and pass and it would work but that would take too much time to do — Kanony, Jun 29 '22 at 09:19
What about `const static double a0[]={... }; ... const static double a7[]={... }; const static double*DCT_matrix[] = {a0, ... a7};` ? — Gerhardh, Jun 29 '22 at 09:21
Your code is probably already as good as it can be. Look at this compiler output: https://www.godbolt.org/z/T15qxeK5j — Jabberwocky, Jun 29 '22 at 09:23
It's unclear if you mean "fast to write" or if you have a performance issue. — klutt, Jun 29 '22 at 09:26
The question does not really make sense in it's current state. You're talking about heap allocation, but there is no heap allocation in your code. And have not posted the actual error (which is not an error, but a warning) — klutt, Jun 29 '22 at 09:33
I get the feeling that the main problem is that the function matrix_multiply is poorly designed — klutt, Jun 29 '22 at 09:36
@klutt well you get it quite wrong actually if the parameter list was `double[8][8]` instead of `double**` it would work perfectly. Your warning is actually a segfault — Kanony, Jun 29 '22 at 09:40
Yes, that is true. I only pointed out that it's not a compiler error. It's a compiler warning. Your post gives the impression that this does not compile. — klutt, Jun 29 '22 at 09:41
But I still wonder why matrix_multiply takes `double**` as arguments. That seem like poor design. `double*` or `double (*)[8]` would make more sense — klutt, Jun 29 '22 at 09:43
@klutt some matrices will be created on the heap so it's better to write `double**` but I have a matrix which is 8*8 and it can't be changed and to write its code for the heap allocation would take so much time (DCT_matrix) — Kanony, Jun 29 '22 at 09:46
If the compiler supports variable length arrays, you could have `void matrix_multiply(size_t n, double p_matrix[][n], const double l_matrix[][n], const double r_matrix[][n]);`. It still seems like an XY problem – the product of the two constant matrices is known before compiling. Why don't you just provide that? — Weather Vane, Jun 29 '22 at 09:47
Why do you assume that `double**` is more suited for heap allocation? You can do that with the two other types too. — klutt, Jun 29 '22 at 09:47
And it also smells bad design that matrix_multiply does not have three arguments, two for input and one for output. — klutt, Jun 29 '22 at 09:48
Yes, but the common way of doing that is still to provide an output argument, and then you can return a pointer to that if you want just to make it possible to use chaining. — klutt, Jun 29 '22 at 09:50
@klutt is it common a convention to design function like that? — Kanony, Jun 29 '22 at 09:52
Yes it is. You can see a lot of examples of that in the standard library. `memcpy` for instance — klutt, Jun 29 '22 at 09:53
@klutt yes they are 8*8 but some are in heap and some are in stack — Kanony, Jun 29 '22 at 09:55
A pointer array with `double**` will by definition be needlessly slow. See [Correctly allocating multi-dimensional arrays](https://stackoverflow.com/questions/42094465/correctly-allocating-multi-dimensional-arrays). As for initialization, you can't do that more efficiently than with memcpy etc since dynamic allocation is done in run-time. It's overall slower than static allocation. — Lundin, Jun 29 '22 at 10:02
@Kanony I updated the answer with a workaround for your current code — klutt, Jun 29 '22 at 10:46
@klutt thanks man for the effort. But I better redesign all as you suggested to do — Kanony, Jun 29 '22 at 10:50
@Kanony I think the most valuable knowledge you gained from this question was about the XY-problem ;) — klutt, Jun 29 '22 at 10:51
@Kanony, I suggest using a pointer to VLA, `typeof(double[])* matrix_multiply(size_t rows, size_t cols, double l_matrix[static rows][cols], double r_matrix[static rows][cols]);`. It works for both stack allocated matrices like `double A[8][10];` and for dynamically allocated `double (*A)[cols] = calloc(rows, sizeof *A);` — tstanisl, Jun 29 '22 at 11:22

klutt · Accepted Answer · 2022-06-29T11:52:54.433

In the comment section, it has become clear that the main flaw is the design of matrix_multiply. There's no reason to pass double** to it.

Firstly, the layout of double** indicates that the matrix is created like this:

double **mat = malloc(8 * sizeof *mat);
for(int i=0; i<8; i++) 
    mat[i] = malloc(8 * sizeof *mat[0]);

This is how this is taught in introductory classes when learning how to deal with pointers. But this causes many unnecessary malloc calls, which also will make it clunky to free all the memory. On top of that, it is slow. Both while allocating/freeing and while working with the memory since the matrix will not be cache friendly. Here is an answer I wrote that is about the cache

I have a few others here

So, in your function, I assume you will have a lot of stuff like l_matrix[x][y]. Replace this with l_matrix[x + 8*y] when going from double** to double*.

On top of that, since you have no output argument, I assume that you allocate the output inside the function. Don't do that. There are many situations when you want to reuse memory allocated for matrices. Instead, provide an output argument and allocate the memory outside the function. This gives us the signature:

void matrix_multiply(const double *A, const double *B, double *output)

Also, since in your case the matrix size is fixed, I would recommend renaming the function in some way that reflects that. Maybe DCT_matrix_multiply but that's up to you.

If you want, you can change the return type to double* and end the function with return output;. That's not mandatory, but it gives the option of doing this:

matrix_multiply(matrix_multiply(A, B, C), 
                matrix_multiply(A, D, E),
                F);

However, be VERY careful if you do that, because the evaluation order of arguments is not specified.

If you really want a function that does the allocation, then write a wrapper:

double *wrapper(const double *A, const double *B) {
    double *output = malloc(8*8 * sizeof *output);
    if(!output) return NULL;
    return matrix_multiply(A, B, output);
}

If you want to work with matrices with non fixed size, it's a bit trickier, but in general, I would still recommend passing double* with added information about sizes. Something like:

double *mul(const double *A, const double *B, double *C, size_t x, size_t y);

If you want to support multiplication with matrices and vectors, the signature will become even more complicated, and the complete solution for this is outside the scope of this answer. But I would write a separate function in such a case.

Another option you might want to investigate is having this signature.

typedef dm double[8][8]; // Short for dct matrix
dm *matrix_multiply(const dm *A, const *dm B, dm *output);

It makes some sense here since you're working with fixed sizes, but in general, this approach is not very common. The reason is simply that the benefits of it doesn't really outweigh the drawbacks of having to deal with pointers to 2d-arrays.

Workaround to your current code

If you want to postpone the rewriting but still be able to initialize the way you want, write a converter.

double **convert (double mat[8][8])
{
    double **ret = malloc (8 * sizeof *ret);
  
    if (!ret) return NULL;

    for (int i = 0; i < 8; i++)
        ret[i] = &mat[i][0];

    return ret;
}

Then you can do

double **A = convert(DCT_matrix);
double **B = convert(other_matrix);

// Use A and B

free(A);
free(B);

As written, `matrix_multiply` has no way to know the sizes of the matrices, so they must be hardcoded. Therefore the better signature is something like `void matrix_mulyiply(const double[8][8], const double[8][8], double[8][8])`. No need to muck with handcrafted 2d indexing inside. — n. m. could be an AI, Jun 29 '22 at 10:26
@n.1.8e9-where's-my-sharem. Yes I know. In OP's case, the size is fixed. Edited the answer a bit to clarify this. — klutt, Jun 29 '22 at 10:31

Initialize C array fast (to write)

1 Answers1

Workaround to your current code