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