5
\$\begingroup\$

Follow up on this post. I have refactored my code and created a matrix_t struct and implemented various functions to create, delete, modify and perform arithmetic on them.

matrix.h

#ifndef MATRIX_H
#define MATRIX_H

typedef struct matrix_t
{
    size_t rows;
    size_t cols;
    double *data;
} matrix_t;

matrix_t *mx_new(const size_t rows, const size_t cols);
matrix_t *mx_identity(const size_t size);
matrix_t *mx_ones(const size_t rows, const size_t cols);
matrix_t *mx_clone(const matrix_t *matrix);
matrix_t *mx_row_major(const double *data, const size_t rows, const size_t cols);
matrix_t *mx_col_major(const double *data, const size_t rows, const size_t cols);

void mx_free(matrix_t *matrix);

void mx_print(const matrix_t *matrix);
void mx_fprint(const matrix_t *matrix, FILE *stream);

void mx_set(const matrix_t *matrix, const size_t i, const size_t j, const double x);
double mx_get(const matrix_t *matrix, const size_t i, const size_t j);
void mx_set_row(const matrix_t *matrix, const size_t i, const double *elem);
void mx_set_col(const matrix_t *matrix, const size_t j, const double *elem);
double *mx_get_row(const matrix_t *matrix, const size_t i);
double *mx_get_col(const matrix_t *matrix, const size_t j);


matrix_t *mx_add(const matrix_t *matrixA, const matrix_t *matrixB);
matrix_t *mx_sub(const matrix_t *matrixA, const matrix_t *matrixB);
matrix_t *mx_mul(const matrix_t *matrixA, const matrix_t *matrixB);
matrix_t *mx_elem_mul(const matrix_t *matrixA, const matrix_t *matrixB);
matrix_t *mx_elem_div(const matrix_t *matrixA, const matrix_t *matrixB);
matrix_t *mx_scalar_mul(const matrix_t *matrix, const double x);
matrix_t *mx_scalar_div(const matrix_t *matrix, const double x);

matrix_t *mx_transpose(const matrix_t *matrix);

void mx_swap_rows(matrix_t *matrix, const size_t i1, const size_t i2);
void mx_swap_cols(matrix_t *matrix, const size_t j1, const size_t j2);

#endif

matrix.c

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

#ifndef OUT_OF_MEMORY
#define OUT_OF_MEMORY "Insufficient space in memory. Program terminated.\n"
#endif

/* Matrix Memory Management Functions */


matrix_t *mx_new(const size_t rows, const size_t cols)
{
    matrix_t *matrix = malloc(sizeof(matrix_t));
    if (!matrix) {
        fprintf(stderr, OUT_OF_MEMORY);
    }
    matrix->rows = rows;
    matrix->cols = cols;
    matrix->data = malloc(rows * cols * sizeof(double));
    if (!matrix->data) {
        fprintf(stderr, OUT_OF_MEMORY);
    }
    return matrix;
}

matrix_t *mx_identity(const size_t size)
{
    matrix_t *matrix = mx_new(size, size);
    for (size_t i = 0; i < size; i++)
    {
        for (size_t j = 0; j < size; j++)
        {
            mx_set(matrix, i, j, (i == j ? 1 : 0));
        }
    }
    return matrix;
}
matrix_t *mx_ones(const size_t rows, const size_t cols)
{
    matrix_t *matrix = mx_new(rows, cols);
    for (size_t i = 0; i < rows; i++)
    {
        for (size_t j = 0; j < cols; j++)
        {
            mx_set(matrix, i, j, 1);
        }
    }
    return matrix;
}

matrix_t *mx_row_major(const double *data, const size_t rows, const size_t cols) {
    matrix_t *matrix = mx_new(rows, cols);
    for (size_t i = 0; i < rows; i++)
    {
        for (size_t j = 0; j < cols; j++)
        {
            mx_set(matrix, i, j, *(data + i * rows + j));
        }
    }
    return matrix;
}
matrix_t *mx_col_major(const double *data, const size_t rows, const size_t cols) {
    matrix_t *matrix = mx_new(rows, cols);
    for (size_t i = 0; i < rows; i++)
    {
        for (size_t j = 0; j < cols; j++)
        {
            mx_set(matrix, i, j, *(data + j * cols + i));
        }
    }
    return matrix;
}

matrix_t *mx_clone(const matrix_t *matrix)
{
    size_t rows, cols;
    matrix_t *copy;
    rows = matrix->rows;
    cols = matrix->cols;
    copy = mx_new(rows, cols);
    for (size_t i = 0; i < rows; i++)
    {
        for (size_t j = 0; j < cols; j++)
        {
            mx_set(copy, i, j, mx_get(matrix, i, j));
        }
    }
    return copy;
}

void mx_free(matrix_t *matrix)
{
    free(matrix->data);
    free(matrix);
}


/* Matrix Output Functions */


void mx_fprint(const matrix_t *matrix, FILE *stream)
{
    for (size_t i = 0; i < matrix->rows; i++)
    {
        for (size_t j = 0; j < matrix->cols; j++)
        {
            fprintf(stream, "%5.2f ", mx_get(matrix, i, j));
        }
        fprintf(stream, "\n");
    }
}

void mx_print(const matrix_t *matrix)
{
    mx_fprint(matrix, stdout);
}

/* Setting and Getting individual elements as well as rows and columns*/

void mx_set(const matrix_t *matrix, const size_t i, const size_t j, const double x)
{
    *(matrix->data + i * matrix->rows + j) = x;
}

double mx_get(const matrix_t *matrix, const size_t i, const size_t j)
{
    return *(matrix->data + i * matrix->rows + j);
}

void mx_set_row(const matrix_t *matrix, const size_t i, const double *elem)
{
    for (size_t j = 0; j < matrix->cols; j++)
    {
        mx_set(matrix, i, j, *(elem + j));
    }
}
void mx_set_col(const matrix_t *matrix, const size_t j, const double *elem)
{
    for (size_t i = 0; i < matrix->rows; i++)
    {
        mx_set(matrix, i, j, *(elem + i));
    }
}
double *mx_get_row(const matrix_t *matrix, const size_t i)
{
    double *elem = malloc(matrix->cols  *sizeof(double));
    if (!elem)
    {
        fprintf(stderr, OUT_OF_MEMORY);
        exit(EXIT_FAILURE);
    }
    for (size_t j = 0; j < matrix->cols; j++)
    {
        *(elem + j) = mx_get(matrix, i, j);
    }
    return elem;
}
double *mx_get_col(const matrix_t *matrix, const size_t j)
{
    double *elem = malloc(matrix->rows  *sizeof(double));
    if (!elem)
    {
        fprintf(stderr, OUT_OF_MEMORY);
        exit(EXIT_FAILURE);
    }
    for (size_t i = 0; i < matrix->rows; i++)
    {
        *(elem + i) = mx_get(matrix, i, j);
    }
    return elem;
}

/* Arithmetic on Matrices */

static void check_same_size(const matrix_t *matrixA, const matrix_t *matrixB)
{
    if (matrixA->rows == matrixB->rows && matrixA->cols == matrixB->cols)
    {
        return;
    }
    fprintf(stderr, "Matrices are not compatible. They have different dimensions.\n");
    exit(EXIT_FAILURE);
}

matrix_t *mx_add(const matrix_t *matrixA, const matrix_t *matrixB)
{
    check_same_size(matrixA, matrixB);
    size_t rows = matrixA->rows;
    size_t cols = matrixA->cols;
    matrix_t *sum = mx_new(rows, cols);
    for (size_t i = 0; i < rows; i++)
    {
        for (size_t j = 0; j < cols; j++)
        {
            double a = mx_get(matrixA, i, j);
            double b = mx_get(matrixB, i, j);
            mx_set(sum, i, j, a + b);
        }
    }
    return sum;
}

matrix_t *mx_sub(const matrix_t *matrixA, const matrix_t *matrixB)
{
    check_same_size(matrixA, matrixB);
    size_t rows = matrixA->rows;
    size_t cols = matrixA->cols;
    matrix_t *diff = mx_new(rows, cols);
    for (size_t i = 0; i < rows; i++)
    {
        for (size_t j = 0; j < cols; j++)
        {
            double a = mx_get(matrixA, i, j);
            double b = mx_get(matrixB, i, j);
            mx_set(diff, i, j, a - b);
        }
    }
    return diff;
}
matrix_t *mx_mul(const matrix_t *matrixA, const matrix_t *matrixB)
{
    // check the size
    if (matrixA->cols != matrixB->rows)
    {
        fprintf(stderr, "Matrices are not compatible for multiplication.\n");
        exit(EXIT_FAILURE);
    }
    size_t m = matrixA->rows;
    size_t n = matrixA->cols;
    size_t p = matrixB->cols;
    matrix_t *product = mx_new(m, p);
    // standard general multiplication
    for (size_t i = 0; i < m; i++)
    {
        for (size_t j = 0; j < p; j++)
        {
            double sum = 0;
            for (size_t k = 0; k < n; k++)
            {
                sum += mx_get(matrixA, i, k) * mx_get(matrixB, k, j);
            }
            mx_set(product, i, j, sum);
        }
    }
    return product;
}
matrix_t *mx_elem_mul(const matrix_t *matrixA, const matrix_t *matrixB)
{
    check_same_size(matrixA, matrixB);
    size_t rows = matrixA->rows;
    size_t cols = matrixA->cols;
    matrix_t *result = mx_new(rows, cols);
    for (size_t i = 0; i < rows; i++)
    {
        for (size_t j = 0; j < cols; j++)
        {
            double a = mx_get(matrixA, i, j);
            double b = mx_get(matrixB, i, j);
            mx_set(result, i, j, a * b);
        }
    }
    return result;
}
matrix_t *mx_elem_div(const matrix_t *matrixA, const matrix_t *matrixB)
{
    check_same_size(matrixA, matrixB);
    size_t rows = matrixA->rows;
    size_t cols = matrixA->cols;
    matrix_t *result = mx_new(rows, cols);
    for (size_t i = 0; i < rows; i++)
    {
        for (size_t j = 0; j < cols; j++)
        {
            double a = mx_get(matrixA, i, j);
            double b = mx_get(matrixB, i, j);
            mx_set(result, i, j, a / b);
        }
    }
    return result;
}
matrix_t *mx_scalar_mul(const matrix_t *matrix, const double x)
{
    size_t rows = matrix->rows;
    size_t cols = matrix->cols;
    matrix_t *result = mx_new(rows, cols);
    for (size_t i = 0; i < rows; i++)
    {
        for (size_t j = 0; j < cols; j++)
        {
            mx_set(result, i, j, mx_get(matrix, i, j) * x);
        }
    }
    return result;
}
matrix_t *mx_scalar_div(const matrix_t *matrix, const double x)
{
    size_t rows = matrix->rows;
    size_t cols = matrix->cols;
    matrix_t *result = mx_new(rows, cols);
    for (size_t i = 0; i < rows; i++)
    {
        for (size_t j = 0; j < cols; j++)
        {
            mx_set(result, i, j, mx_get(matrix, i, j) / x);
        }
    }
    return result;
}

/* Special Matrix Functions */

matrix_t *mx_transpose(const matrix_t *matrix)
{
    size_t rows = matrix->cols;
    size_t cols = matrix->rows;
    matrix_t *trans = mx_new(rows, cols);
    for (size_t i = 0; i < rows; i++)
    {
        for (size_t j = 0; j < cols; j++)
        {
            mx_set(trans, i, j, mx_get(matrix, i, j));
        }
    }
}

/* Matrix Mutating Functions */

void mx_swap_rows(matrix_t *matrix, const size_t i1, const size_t i2)
{
    double *row = mx_get_row(matrix, i1);
    mx_set_row(matrix, i1, mx_get_row(matrix, i2));
    mx_set_row(matrix, i2, row);
}
void mx_swap_cols(matrix_t *matrix, const size_t j1, const size_t j2)
{
    double *col = mx_get_col(matrix, j1);
    mx_set_col(matrix, j1, mx_get_col(matrix, j2));
    mx_set_col(matrix, j2, col);
}

I'd like comments and critique on general code style, performance issues, and anything that seems out of place.

\$\endgroup\$

3 Answers 3

Reset to default
4
\$\begingroup\$

Error handling

I do not expect math library functions to decide when program must be terminated, I'd then avoid to call exit() (or even printing output to stderr) instead of the caller.

If you do not really need the error code then you may simply return NULL instead:

double *mx_get_col(const matrix_t *matrix, const size_t j)
{
    double *column = malloc(matrix->rows  *sizeof(double));
    if (!column)
        return NULL;

    // ...
}

The other way is to use a return code:

error_t mx_get_col(const matrix_t *matrix, const size_t j, double **column)
{
    *column = malloc(matrix->rows  *sizeof(double));
    if (!*column)
        return ENOMEM;

    // ...
}

This lets me introduce another problem: you're not validating your inputs. Yes, it comes with a performance penalty but you should be ready to pay for it (in very performance critical scenarios at least when compiling debug builds).

error_t mx_get_col(const matrix_t *matrix, const size_t j, double **column)
{
    if (matrix == NULL)
        return EFAULT;

    if (j >= matrix->cols)
        return EINVAL;

    // ...
}

Note that there more preconditions you might want to check, even just in debug:

assert(matrix->data != NULL);

Do not underestimate error handling favoring performance, you will pay a huge price during debugging (if you're lucky). You already have a struct instead of a typedef double** matrix then you're already paying for some overhead.

Features

When working with matrices you do not always want to create a new matrix for each operation (it hurts performance and consumes memory). You may want to consider to add in-place variants:

void mx_scalar_mul_inplace(matrix_t *matrix, const double x);

Note that more often than not you can simply implement mx_scalar_mul() just calling mx_scalar_mul_inplace(). Also note that for some operations (when it makes sense!) it's pretty handy to directly change one of the inputs without allocating new memory (remember that we're talking about a math library).

Sometimes you may also need to reuse an existing structure (again to save multiple allocations):

void mx_scalar_mul_to(const matrix_t *matrix, const double x, matrix_t *result);

You may consider to introduce a vector class.

Naming

It has been said many times but it is worth to repeat again. POSIX reserves _t suffix then you should not use it for your own types.

I, personally, think you should not define a typedef for matrix. What for? It makes code much more clear:

struct matrix *mx_scalar_mul(const struct matrix *matrix, const double x);

Note that naming a struct matrix has very few chances to conflict with anything else, in C identifiers have different name spaces (labels, struct/union/enum, struct members and other identifiers; see 6.2.3).

Misc

As I said I'd avoid to print error messages from library functions but when you have to then you should use fputs() instead of fprintf when printing simple strings without formatting:

fputs(OUT_OF_MEMORY, stderr);

It's slightly faster, safer (because you do not have to worry about formatting characters in your strings) and it better communicates the intent (to print a raw unformatted string). Note that compiler may detect this usage pattern and use fputs() for you.

You may want to declare your parameters as const matrix * const (it's not useful for the caller but as you're already doing with const double it helps to keep your own code correct).

\$\endgroup\$
2
  • \$\begingroup\$ I had enquired about the practice of typedef structs earlier here. \$\endgroup\$ Commented Nov 2, 2017 at 18:14
  • \$\begingroup\$ Yes, it's definitely an opinion. My rule of thumb is that if I have to add suffixes/prefixes because of them then I'm using them wrongly. \$\endgroup\$ Commented Nov 2, 2017 at 18:21
3
\$\begingroup\$

I like it, don't have any criticisms about style or anything like that, but I can't use it personally. I'll explain why and a possible solution.

Performance

I work in performance-critical fields like raytracing and, by far, I work with square matrices that are 4x4. For that your library is a bit overkill in my case considering that it requires two heap allocations and deallocations to create and destroy each matrix and creates a new one with every single operation.

In my case I probably benefit most from a library specifically designed for 4x4 matrices (and that's what I have), but if there was any possibility of me using your library, I have some performance requests.

An immediate performance improvement would be to use output parameters whenever possible so that I could reuse existing matrices, already hot in the cache, when possible. That is, instead of:

matrix_t *mx_mul(const matrix_t *matrixA, const matrix_t *matrixB);

We have:

void mx_mul(const matrix_t *matrixA, const matrix_t *matrixB, 
            matrix_t *outMatrix);

.. or something to this effect. Syntactically it might seem a little bit more cumbersome, but it might actually result in reduced misuse and human error since it could result in clients requiring fewer matrices to have to free manually at the end of a function, especially in the event of an error. What might have required 16 matrices to manually destroy might now only require 4.

If you really want to make me happy from a performance standpoint, you might do this:

typedef struct matrix_t
{
    size_t rows;
    size_t cols;

    /* Points to `mem` when the matrix fits, otherwise
     * points to dynamic array. */
    double *data;

    /* Used when the matrix data fits. */
    double mem[4*4];
} matrix_t;

void mx_new(const size_t rows, const size_t cols, 
            matrix_t *outMesult);

/* returns 1 if matrix was successfully created. */
int mx_valid(const matrix_t *matrix);

... at which point it would be used like this:

matrix_t m1, m2, m3;
mx_new(4,4, &m1);
mx_new(4,4, &m2);
mx_new(4,4, &m3);

if (mx_valid(&m1) && mx_valid(&m2) && mx_valid(&m3))
{
    /* mx_new above would require no heap allocation since
     * a 4x4 matrix fits into the `mem` field */

    /* do stuff with m1 and m2 and m3 */
    ...

}

/* these calls to destroy the matrices also require no
 * heap deallocation, since 'data' points to 'mem'. These
 * functions do nothing if the matrix is not in a valid
 * state (data points to null). */
mx_free(&m1);
mx_free(&m2);
mx_free(&m3);

mx_free would check if data == mem. If so, it doesn't need to free anything. Otherwise, it only needs to free the data field. The matrix itself is never dynamically allocated by the library itself.

This does introduce a safety concern where matrix_t would no longer be safe to copy by value without invaliding the data pointer if it points to mem, but hey, this is C, and I want performance with my matrix library for small matrices.

Override Memory Allocation

The next one would be to allow overriding the memory allocation functions used by the matrix library, so that instead of malloc and free, I can provide my own with, say, a fixed allocator to be used whenever matrices are encountered of an anticipated size, or a general-purpose allocator which keeps track of memory usage for leak detection (something I lean on a lot in my codebase a lot as part of the unit testing). It should definitely widen the applicability of your library, and most people might be fine just keeping it with the defaults at which point the function pointers can simply point to default functions which use malloc and free.

typedef void* mx_malloc_matrix(void *data, size_t memSize);
typedef void mx_free_matrix(void *data, , size_t memSize, void *mem);

/* Specifies functions to use to allocate and free matrix data. 
 * 'data' gets passed through to these callbacks. */
void mx_memory(void *data, 
               mx_malloc_matrix *mallocFunc, 
               mx_free_matrix *freeFunc);

If you don't like the safety issue mentioned above with cases where data points to the mem field, then you can get rid of the mem field optimization for small matrices and I can get around that by supplying the library with a fixed allocator for those 4x4 cases if your library provides these capabilities.

\$\endgroup\$
1
\$\begingroup\$

I don't have any specific complaints about your code, just some ideas for future directions or alternate ways of doing some things.

Naming

I like the typedefed structs (but avoid the _t as the other answer advises). In fact I often go further and put the pointer in there too. Some people don't like it. Most of their complaints I have found fall down to simple dislike or unfamiliarity, except one. The real problem that can arise is breaking the abstraction.

I think it's defensible to typedef <<whatever gobbledegook>> to make an opaque handle that gets passed around and returned from functions. That's what typedefs are for. Making abstractions.

The problem happens when you use it in an expression that peeks into its type defintion. If you have a simple handler function

mything squelch( mything it ){
    ...
}

It should not use it in some freaky manner requiring you to look up the typedef to figure out if it's correct.

mything squelch( mything it ){
    return  it->glorf( it->florg[ it->grrrf ] ); // whaaa??
}

That would be breaking the abstraction. If it just passed it along, perhaps with additional parameters to a function with a more informative name, then it would okay IMO.

Math functions on matrices

You've probably started to experience already the tedium of writing the same for loops over and over again with a different math operator in each one. I have a solution for that, 2 in fact over on SO. One uses macros, the other less extreme example uses enums. But either way, if you pass-in the operator then you only need one binary_operation() function.

3D matrices?

Further down the rabbit hole, I have some explanation of a multidimensional array structure that can hold 2D, 3D, ..., nD arrays. I'm not saying that you need to do this or that your code needs this capability or anything, just that the path is there. :)

\$\endgroup\$
3
  • \$\begingroup\$ Just to clarify, your suggestion for the naming issue is to simply drop the _t suffix from the name, right? \$\endgroup\$ Commented Nov 4, 2017 at 12:10
  • \$\begingroup\$ @Astrobleme that's correct. \$\endgroup\$ Commented Nov 4, 2017 at 15:08
  • 1
    \$\begingroup\$ @Astrobleme for the rest of my comments on naming, I respect the opinion of the other reviewer. I merely meant to offer an alternative viewpoint. You should feel free to make up your own mind on what constitutes a sensible set of rules for when to typedef and when not to typedef. I hope that's clear. Happy coding! \$\endgroup\$ Commented Nov 22, 2017 at 2:48

You must log in to answer this question.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.