Matrices in C implementation

Question

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

/**
    Adds two matrices together

    @param *a pointer to first matrix (A);
    @param *b pointer to second matrix (B);
    @return A+B
*/
Matrix add(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);
    //add matrices
    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;
}

Answer 1

I just wanted to comment on your question about accessing the array efficiently. You ask:

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.

Answer 2

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;

to a more readable

            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.

Answer 3

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.

Answer 4

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

Answer 5

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.