Matrix structure in C

I know of two ways to represent a matrix in plain C. One is to have arrays in arrays. So the matrix is an array, where each entry contains a row of the matrix. You access the entries with something like matrix[row][col].

The other one is to linearize the matrix and calculate the index with row * width + col where width is the width of the matrix.

I wondered which one would be faster, so I tested this. It seems that the linearized version is a bit faster:

Compiler filename time
gcc double-ref.c 1.726 ± 0.003
gcc linear.c 1.517 ± 0.006
gcc linear-inline.c 1.512 ± 0.002
clang double-ref.c 1.729 ± 0.005
clang linear.c 1.635 ± 0.005
clang linear-inline.c 1.64 ± 0.01

I have the linearized version twice. One has the index calculation in a function, the other has this inlined at every use manually. There is no significant difference, the compiler did the inlining itself well. This means that there is no need to abdicate this index function.

The test code

numbers.h

First of all, I need to set the number of iterations to the same values in every program. This is done with a header file:

#pragma once

#define DIM 10000
#define ITER 100000

double-ref.c

This is the variant which uses the double array type, with double pointers. The matrix_alloc function is more complex, since it has to setup up the raw memory given by malloc to have double * in the first dim segments such that they point to the beginning of the correct row.

#include "numbers.h"

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

/**
  Allocates a square matrix.
  */
double **matrix_alloc(int dim) {
    int row;
    /* Allocate space for `dim` pointers (row beginnings) and for `dim^2`
       numbers.
       */
    double **space = malloc(dim * sizeof(double *) + dim*dim * sizeof(double));

    double *first_number = (double *) (space + dim);

    for (row = 0; row < dim; row++) {
        space[row] = first_number + dim * row;
    }

    return space;
}

void matrix_free(double **matrix) {
    free(matrix);
}

/**
  Computes the trace of the matrix.
  */
double matrix_trace(double **matrix, int dim) {
    double sum = 0;
    int row;
    for (row = 0; row < dim; row++) {
        sum += matrix[row][row];
    }
    return sum;
}

int main(int argc, char **argv) {
    int dim = DIM;
    int iter = ITER;
    double **matrix = matrix_alloc(dim);
    double trace_sum;
    while (iter-- != 0) {
        matrix[0][0]++;
        trace_sum = matrix_trace(matrix, dim);
    }

    printf("Trace Sum: %g\n", trace_sum);

    matrix_free(matrix);

    return 0;
}

linear-inline.c

This is with a one-dimensional array and inline index calculation. This makes matrix_alloc way easier.

#include "numbers.h"

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

/**
  Allocates a square matrix.
  */
double *matrix_alloc(int dim) {
    return malloc(dim * dim * sizeof(double));
}

void matrix_free(double *matrix) {
    free(matrix);
}

/**
  Computes the trace of the matrix.
  */
double matrix_trace(double *matrix, int dim) {
    double sum = 0;
    int row;
    for (row = 0; row < dim; row++) {
        sum += matrix[row * dim + row];
    }
    return sum;
}

int main(int argc, char **argv) {
    int dim = DIM;
    int iter = ITER;
    double *matrix = matrix_alloc(dim);
    double trace_sum;
    while (iter-- != 0) {
        matrix[0]++;
        trace_sum = matrix_trace(matrix, dim);
    }

    printf("Trace Sum: %g\n", trace_sum);

    matrix_free(matrix);

    return 0;
}

linear.c

Same as above, but now with an extra index function (idx).

#include "numbers.h"

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

int idx(int row, int col, int dim) {
    return row * dim + col;
}

/**
  Allocates a square matrix.
  */
double *matrix_alloc(int dim) {
    return malloc(dim * dim * sizeof(double));
}

void matrix_free(double *matrix) {
    free(matrix);
}

/**
  Computes the trace of the matrix.
  */
double matrix_trace(double *matrix, int dim) {
    double sum = 0;
    int row;
    for (row = 0; row < dim; row++) {
        sum += matrix[idx(row, row, dim)];
    }
    return sum;
}

int main(int argc, char **argv) {
    int dim = DIM;
    int iter = ITER;
    double *matrix = matrix_alloc(dim);
    double trace_sum = 0;
    while (iter-- != 0) {
        matrix[idx(0, 0, dim)]++;
        trace_sum += matrix_trace(matrix, dim);
    }

    printf("Trace Sum: %g\n", trace_sum);

    matrix_free(matrix);

    return 0;
}