# Matrices in C implementation

I recently wrote a C implementation of Matrices with add, subtract, and multiply. I want to expand this out eventually where I can diagonalize, efficiently square, row-reduce, etc...

I was wondering if there are any more efficient ways to do what I'm currently doing (I'm really new to C and pointers/references in general) and without changing too much, getting

*(result.matrix+ i*r  + j)


to work like

result[i][j]


Here's the code below

/**
Matrix Multiplication
matrices.c
Matrix data structure in C.

@author Michael Asper
@version 1.0 3/29/17
*/

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

typedef struct Matrix {
int     rowSize;
int     columnSize;
long int*    matrix;
} Matrix;

/**
Randomizes the elements of a matrix

@param *m pointer to Matrix to randomize;
*/
void randomize(Matrix *m){
int i,j;
for(i = 0; i < m->rowSize ; i++){
for(j = 0; j < m->columnSize; j++){
*(m->matrix + i*m->rowSize  + j)= rand() % 5000;
}
}
}

/**
Returns a r x c Matrix with all 0s.

@param r The row size of the matrix
@param c The column size of the matrix
@return r x c Matrix
*/
Matrix createMatrix(int r, int c){
Matrix temp = {r, c, calloc(r * c, sizeof(long int *))};
return temp;
}

/**
Returns a r x c Matrix with random numbers.

@param r The row size of the matrix
@param c The column size of the matrix
@return r x c Matrix
*/
Matrix createRandMatrix(int r, int c){
Matrix temp = createMatrix(r,c);
randomize(&temp);
return temp;
}

/**
Prints matrix.

@param *m Pointer to Matrix you want to print
*/
void printMatrix(Matrix *m){

int i,j;
for(i = 0; i < m->rowSize ; i++){
for(j = 0; j < m->columnSize; j++){
printf("%li ", *(m->matrix + i*m->rowSize  + j));
}
printf("\n");
}
}

/**

@param *a pointer to first matrix (A);
@param *b pointer to second matrix (B);
@return A+B
*/
//check if matrices are compatible
if(a->rowSize != b->rowSize || a->columnSize != b->columnSize){
fprintf(stderr, "Error: Incompatible sizes");
exit(0);
}
//create result matrix
int r = a->rowSize;
int c = a->columnSize;
Matrix result = createMatrix(r,c);
int i,j;
for(i = 0; i < r ; i++){
for(j = 0; j < c; j++){
//result[i][j] = a[i][j]+b[i][j]
*(result.matrix+ i*r  + j) = *(a->matrix + i*r  + j) + *(b->matrix + i*r  + j);
}
}
return result;
}

/**
Subtracts two matrices together

@param *a pointer to first matrix (A);
@param *b pointer to second matrix (B);
@return A-B
*/
Matrix sub(Matrix *a, Matrix *b){
//check if matrices are compatible
if(a->rowSize != b->rowSize || a->columnSize != b->columnSize){
fprintf(stderr, "Error: Incompatible sizes");
exit(0);
}
//create result matrix
int r = a->rowSize;
int c = a->columnSize;
Matrix result = createMatrix(r,c);
//subtracts matrix
int i,j;
for(i = 0; i < r ; i++){
for(j = 0; j < c; j++){
//result[i][j] = a[i][j]-b[i][j]
*(result.matrix+ i*r  + j) = *(a->matrix + i*r  + j) - *(b->matrix + i*r  + j);
}
}
return result;
}

/**
Multiplies two matrices together

@param *a pointer to first matrix (A);
@param *b pointer to second matrix (B);
@return A*B
*/
Matrix multiply(Matrix *a, Matrix *b){
//check if matrices are compatible
if(a->columnSize != b->rowSize ){
fprintf(stderr, "Error: Incompatible sizes");
exit(0);
}

//initialize return matrix
int r = a->rowSize;
int c = b->columnSize;
Matrix result = createMatrix(r,c);

//multiply matrices
int i,j;
for(i = 0; i < r ; i++){
for(j = 0; j < c; j++){
long int sum = 0;
int k;
for(k = 0; k < a->columnSize; k++){
//sum += a[i][k] * b[k][j]
sum = sum + (*(a->matrix + i*a->rowSize  + k)**(b->matrix + k*b->rowSize  + j));
}
*(result.matrix+ i*r  + j) = sum;
}
}
return result;
}

int main(){
// seed random with time
time_t t;
srand((unsigned) time(&t));

//setup random matrices and multiply
Matrix a = createRandMatrix(3,100);
Matrix b = createRandMatrix(100,3);
Matrix result = multiply(&a,&b);
printMatrix(&result);

return 0;
}

• 3/29/17 don't. Use the international standard which everyone understands: 2017-03-29 – Daniel Jour Mar 30 '17 at 13:31

I was wondering if there are any more efficient ways to do what I'm currently doing (I'm really new to C and pointers/references in general) and without changing too much, getting

*(result.matrix+ i*r + j)

to work like

result[i][j]

One way is to do as @mdfst13 has recommended and allocate an array of arrays. That's a fine way to do it, but it can be a performance issue when accessing the array element by element in a loop. CPUs generally optimize to access the next few bytes past the last access since you're likely to need bytes near the ones you previously accessed. If you have a separate array per row, this can throw off that optimization when you reach the end of each row. You'd want to first profile to verify that's an issue or not. If it is, another option you have is to leave it as a single allocation and simply write an accessor function. Something like:

long int getElement(const Matrix m, const int r, const int c)
{
return *(m.matrix + r * m.rowSize + c);
}


This also allows you the opportunity to do some range checking when in debug mode by doing something like:

long int getElement(const Matrix m, const int r, const int c)
{
#if DEBUG
assert((r >= 0) && (r < m.rowSize));
assert((c >= 0) && (c < m.colSize));
#endif
return *(m.matrix + r * m.rowSize + c);
}


(Or if you want to pass a pointer in instead of passing by value, as suggested in mdfst13's answer, you could change the prototype to take a pointer to a Matrix and dereference the fields via pointer.)

Now when you want to access an element of the array, you would write:

long int x = getElement(result, i, j);


It's not as concise as just result[i][j], but it's better than writing out the math every time and potentially getting it wrong in some subtle way.

• I'll eventually being using OpenMP/MPI to parallelize this code for my university's cluster, would either accessing method affect parallelization? – c0rruptbytes Mar 30 '17 at 5:16
• It would depend entirely on the access patterns. Will different threads be reading and writing data to/from the same matrix at the same time? (Sounds like a bad idea, but I don't know your requirements.) Do you tend to work on columns or more on rows? The best thing to do is put it in a function that you can change later. That way you can come up with a (deterministic - not random!) test and run it under both conditions to see which is faster. – user1118321 Mar 30 '17 at 5:21
• There's no pass by reference in C. I know this may be nit picking, but IMO one should use correct terminology throughout to avoid confusion. (+1 for the rest of the answer though) – Daniel Jour Mar 30 '17 at 13:37
• @DanielJour Fair point! I guess I should have said pass a pointer to the data in. I'll update it. – user1118321 Mar 30 '17 at 16:23
• Would "throw off that optimization" slow down the performance? – fu DL Sep 25 '19 at 8:26

You can read about the general handling of multi-dimension dynamic sized arrays at Stack Overflow. But specific to your code

    long int*    matrix;


Replace this with

    long int **  matrix;


And then change

Matrix createMatrix(int r, int c){
Matrix temp = {r, c, calloc(r * c, sizeof(long int *))};
return temp;
}


to something like

Matrix createMatrix(int r, int c){
Matrix temp = {r, c, calloc(r, sizeof(long int *))};

if (temp.matrix == NULL) {
/* panic */
}

for (int i = 0; i < r; i++) {
temp.matrix[i] = calloc(c, sizeof temp.matrix[i][0]);

if (temp.matrix[i] == NULL) {
/* panic */
}
}

return temp;
}


You should check the result of calloc before using. Replace /* panic */ with any number of alternatives. E.g. exit, log an error, print an error message, or all or some of these. The failure is rare, but it helps to respond to it early.

Your original code used the wrong type in the sizeof command. Presumably this didn't matter because the pointer type was at least as big as a long int.

To avoid that kind of bug, consider doing it the way that I did in the second calloc. Give it an example variable of the right type and let it figure out the type. That way if you update the type in the future, you don't have to update every calloc to match. When using sizeof this way, you don't need the parentheses.

Then you can change lines like

            *(m->matrix + i*m->rowSize  + j)= rand() % 5000;


            m->matrix[i][j] = rand() % 5000;


This also gives you the ability to have different sized rows, but you don't seem to need that.

It looks odd to me to pass struct objects by value. I would normally expect functions to return pointers to struct objects. Or pass a pointer to the struct into the function, which populates it.

Most of the time people just write long and not long int. There's no harm to writing it out; that's just not what I normally see.

And you should free the memory that you calloc when you're done with it. Note that you have to do this row by row, just as you created it. Then the overall variable.

• Why do you think that this would be better? (Apart from fixing the bug, of course.) – Carsten S Mar 30 '17 at 8:53
• I would suggest allocating the array as calloc(rows, sizeof(long[cols]));, so you have a contiguous block of memory. Then when you need to access it (in a readable way) just write long (*matrix)[m->c] = m->matrix;, so you can access it like matrix[i][j]; /* i < m->rowSize && j < m->colSize */. – ammut Mar 30 '17 at 10:54

A few notes:

1. Implement multiplication of a matrix by a scalar (i.e., a real number). Then your subtraction routine becomes:

1. Multiply matrix B by scalar -1. Let this resultant matrix be C
2. Add A and C

More useful and less duplicated code. (Why did you not think of implementing scalar multiplication in the first place? Closure under addition and scalar multiplication (and the presence of a null element) define vector spaces in general, and vector spaces of matrices (real or complex) in particular.)

2. Some might disagree with me, but I think all indexed accesses and size declarations should use parameters of type size_t, as that is the type deemed fit for this purpose by the standard. This means that all those rs and cs should be of type size_t, instead of int. Also, so should rowSize and columnSize. I would prefer calling them either rows and columns or numRows and numColumns respectively, as there is less ambiguity about their purpose that way.

3. Matrices are defined on either real or complex, not integer, scalar fields. I guess that the type of an element of a matrix should be double, if not, at least float. This will come in useful when calculating matrix inverses, as they involve division, and the field of integers is not closed under division.

4. Everything that can be const should be const. That includes all function parameters, as you do not modify them or their contents.

• 1) is definitely a good idea, and 2) you're completely corrupt about rowSize and columnSize being ambiguous as I got really confused writing this partial pivot row reduce function (confused rowLength as rowSize when rowSize is numRows). 3) Just converted to double as well :) – c0rruptbytes Mar 30 '17 at 7:18
• @kllrshrk I guess you meant "correct", not "corrupt" ;) – Tamoghna Chowdhury Mar 30 '17 at 10:53
• As an aside, I've recently been trying out functional programming, so I've been working on a FP-compatible matrix implementation in Scala. If you've got any interest that way, you might want to take a look at github.com/tamchow/ScalaStuff/blob/master/src/in/tamchow/… – Tamoghna Chowdhury Mar 30 '17 at 10:55

If your code runs on a x86/ x64 platform, you could also check if sse is available (which is true on all modern cpu's, but can be verified by cpuid) and use sse to further boost your program's performance. See http://download.intel.com/design/PentiumIII/sml/24504501.pdf

• While this link may answer the question, it is better to include the essential parts of the answer here and provide the link for reference. Link-only answers can become invalid if the linked page changes. - From Review – Graipher Mar 30 '17 at 15:26
• Right, that's reasonable, I will remember this next time. – Lambda Mar 30 '17 at 18:34
• It's not too late for this answer, yet ;-) – Graipher Mar 30 '17 at 19:15

### Macro suggestion

I wanted to quickly suggest something I haven't seen in the other answers. You could use a macro to make your code look nicer:

#define M(m, i, j)    *((m)->matrix + (i)*((m)->rowSize) + (j))


So code like this:

*(result.matrix+ i*r  + j) = *(a->matrix + i*r  + j) + *(b->matrix + i*r  + j);


becomes:

M(&result, i, j) = M(a, i, j) + M(b, i, j);


The advantage of this macro over something like a getElement() function is that you can use M() as an lvalue, whereas getElement() can only be an rvalue.