Matrix definition, adding two different data type matrices

Question

I'm writing my own Matrix type version in C programming language. I'm also trying to use information hiding, so I have my code divided into two different files: a header file containing type definitions and function prototypes and a file containing code.

Since I need to work with square matrices, I defined my matrix struct as follow:

struct matrix {
  int dim;
  enum type matrix_type;
  union {
    int ** int_matrix;
    long ** long_matrix;
    float ** float_matrix;
    double ** double_matrix;
  } data;
};

where enum type is defined as follow:

enum type {M_INT = 1, M_LONG = 2, M_FLOAT = 3, M_DOUBLE = 4};

Matrix type is defined in header file as follow:

typedef struct matrix * Matrix;

Build and destroy functions work correctly, a filling function works too, I've found a way to calculate matrix determinant, too. To cope with the choice to have Matrix->data split into different types, I'm using the switch statement intensively, as following code shows:

static int terminate(char * message)
{
  fprintf(stderr, "%s\n", message);
  exit(EXIT_FAILURE);
}

static int set_dim (Matrix m, int dim)
{
  m->dim = dim;
  return 0;
}

static int set_type (Matrix m, int t)
{
  if (t < 1 || t > 4) {
    terminate("Matrix type not valid.");
  }
  m->matrix_type = t;
  return 0;
}


int create_matrix (Matrix * m, int m_dim, int m_type)
{
  if ( ( (*m) = malloc( sizeof(struct matrix) ) ) == NULL )
    terminate("Not possible to create a new matrix type.");

  set_dim(*m, m_dim);
  set_type(*m, m_type);

  switch (m_type) {
    case M_INT:
      if ( !( (*m)->data.int_matrix = (int **) malloc(m_dim * sizeof(*((*m)->data.int_matrix)))) )
        terminate("Not enough space.");
      for(int i = 0; i < m_dim; i++)
        if (!( (*m)->data.int_matrix[i] = (int *) malloc(m_dim * sizeof(*((*m)->data.int_matrix[i])))) )
          terminate("Not enough space.");

      init_int_matrix((*m)->dim, (*m)->data.int_matrix);
      break;

    case M_LONG:
      if ( !( (*m)->data.long_matrix = (long **) malloc(m_dim * sizeof(*((*m)->data.long_matrix)))) )
        terminate("Not enough space.");
      for(int i = 0; i < m_dim; i++)
        if (!( (*m)->data.long_matrix[i] = (long *) malloc(m_dim * sizeof(*((*m)->data.long_matrix[i])))) )
          terminate("Not enough space.");

      init_long_matrix((*m)->dim, (*m)->data.long_matrix);
      break;

    case M_FLOAT:
      if ( !( (*m)->data.float_matrix = (float **) malloc(m_dim * sizeof(*((*m)->data.float_matrix)))) )
        terminate("Not enough space.");
      for(int i = 0; i < m_dim; i++)
        if (!( (*m)->data.float_matrix[i] = (float *) malloc(m_dim * sizeof(*((*m)->data.float_matrix[i])))) )
          terminate("Not enough space.");

      init_float_matrix((*m)->dim, (*m)->data.float_matrix);
      break;

    case M_DOUBLE:
      if ( !( (*m)->data.double_matrix = (double **) malloc(m_dim * sizeof(*((*m)->data.double_matrix)))) )
        terminate("Not enough space.");
      for(int i = 0; i < m_dim; i++)
        if (!( (*m)->data.double_matrix[i] = (double *) malloc(m_dim * sizeof(*((*m)->data.double_matrix[i])))) )
          terminate("Not enough space.");

      init_double_matrix((*m)->dim, (*m)->data.double_matrix);
      break;
  }

  return 0;
}

Now I want to write a function to add two matrices. However, it seems that I'm stuck with my own Matrix definition. The easy choice is to code a function which adds two matrices only when they share the same matrix_type, returning false or -1 otherwise; but what about writing a function to add two matrices despite their inner data definition (int, long, float, double)? Is there any elegant and practical way to do that, instead of writing a function with nested if statements to check any possible type coupling?

Besides, I consider and think about my choice to define a struct matrix as I did. Is there any real advantage of splitting matrix data by int, long, float and double? What if I define only a double data type and then I cast it back to int or long when I need it? In this second case I might get rid of pointers and unions.

Edit: I'm currently using gcc 9.3.0 compiler, gnu11.

Answer 1

To simplify the handling of all type combinations from quadratic to linear complexity, you could convert one or both arguments to a common appropriate type and only handle operations on identical types.

You would still need to precisely define the semantics in some corner cases, for example

• how do you handle integer overflow?
• do you round values when multiplying a matrix by a scalar or do you change the matrix type?

Unless you have a very good reason to handle value types other than double, I would recommend you first write code for this case.

Answer 2

This problem is approached differently. In your current implementation, you have to explicitly handle all possible combinations of all types. Adding a new type is increasingly hard, as it adds complexity and blows up the number of combinations.

This is a common problem. It is known. So before you start making a big pile of endless, unreadable and unmaintainable switches with many hidden bugs, read up on virtual table, object-oriented programming in C and about generic programming and templates.

Abstract interfaces are done with virtual functions and virtual table. For examples, I can point file_operations from Linux kernel and GNOME glib library design.

There are no good ways to do templates in C - if you want them, move to another programming language, like Rust or C++. In C you can either do big pile of endless macros (that are hard to document and hard to program and have unreadable error messages), or use define/undef tricks (which are a bit better, but you never know to which pass the errors are about), or abstract every operation with function pointers, just like in GNOME glib library or qsort does (which result in runtime cost of calling the function), or use something outside of C for code generation.

Below is the "big pile of endless macros" method that define a template matrix for a custom type with all accompanying functions inside a chosen namespace. There's also a definition for matrix_add function for 3 different types. Finally, a short main shows usage example, with an 2x2 int + 2x2 double matrix addition that result in a float matrix.

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

// ---------------------------------------------------------------

// Put this in a header
#define MATRIX_H(M, TYPE) \
\
/* the matrix type, for namespace access */ \
typedef TYPE M##matrix_type; \
\
/* the structure that represetntes a matrix */ \
typedef struct { \
    size_t dim; \
    TYPE *arr; \
} M##matrix; \
\
int M##matrix_create(M##matrix *m, size_t dim); \
int M##matrix_create_from_array(M##matrix *m, size_t N, TYPE (*arr)[N]); \
void M##matrix_destroy(M##matrix *m); \
TYPE *M##matrix_getp(M##matrix *m, size_t x, size_t y); \
TYPE M##matrix_get(const M##matrix *m, size_t x, size_t y); \
int M##matrix_add(M##matrix *r, const M##matrix *a, const M##matrix *b); \
void M##matrix_print(M##matrix *m); \
\
MATRIX_OP3_H(M, M, M) \
\
/* */


// put this in a source file
#define MATRIX_C(M, TYPE, PRINTF) \
\
int M##matrix_create(M##matrix *m, size_t dim) { \
    m->arr = malloc(sizeof(TYPE) * dim * dim); \
    if (!m->arr) return -ENOMEM; \
    m->dim = dim; \
    return 0; \
} \
\
int M##matrix_create_from_array(M##matrix *m, size_t N, TYPE (*arr)[N]) { \
    int r = M##matrix_create(m, N); \
    if (r) return r; \
    const size_t dim = m->dim; \
    for (size_t x = 0; x < dim; ++x) {  \
        for (size_t y = 0; y < dim; ++y) {  \
            *M##matrix_getp(m, x, y) = arr[x][y]; \
        } \
    } \
    return 0; \
} \
\
void M##matrix_destroy(M##matrix *m) { \
    free(m->arr); \
} \
\
size_t M##matrix_dim(const M##matrix *m) { \
    return m->dim; \
} \
\
TYPE *M##matrix_getp(M##matrix *m, size_t x, size_t y) { \
    return &m->arr[x * m->dim + y]; \
} \
\
TYPE M##matrix_get(const M##matrix *m, size_t x, size_t y) { \
    return *M##matrix_getp((M##matrix *)m, x, y); \
} \
\
void M##matrix_print(M##matrix *m) { \
    const size_t dim = m->dim; \
    for (size_t x = 0; x < dim; ++x) { \
        for (size_t y = 0; y < dim; ++y) { \
            printf("%"PRINTF "%s", \
                M##matrix_get(m, x, y), \
                y + 1 == dim ? "\n" : " "); \
        } \
    } \
} \
\
/* stub implementation */ \
TYPE M##M##M##matrix_elem_add(TYPE a, TYPE b) { \
    return a + b; \
} \
\
MATRIX_OP3_C(M, M, M) \
\
int M##matrix_add(M##matrix *m, const M##matrix *a, const M##matrix *b) { \
    return M##M##M##matrix_add(m, a, b); \
} \
\
/* */

// shortcut for both
#define MATRIX_HC(M, TYPE, PRINTF) \
    MATRIX_H(M, TYPE) \
    MATRIX_C(M, TYPE, PRINTF)

// ---------------------------------------------------------------

#define MATRIX_OP3_H(M, A, B)  \
\
/* callback for adding 2 variables of differnet types */ \
/* must be implemented by user */ \
M##matrix_type M##A##B##matrix_elem_add(A##matrix_type a, B##matrix_type b); \
\
int M##A##B##matrix_add(M##matrix *m, const A##matrix *a, const B##matrix *b); \
\
/* */

#define MATRIX_OP3_C(M, A, B)  \
\
int M##A##B##matrix_add(M##matrix *m, const A##matrix *a, const B##matrix *b) { \
    assert(M##matrix_dim(m) == A##matrix_dim(a)); \
    assert(M##matrix_dim(m) == B##matrix_dim(b)); \
    const size_t dim = M##matrix_dim(m); \
    for (size_t x = 0; x < dim; ++x) { \
        for (size_t y = 0; y < dim; ++y) { \
            *M##matrix_getp(m, x, y) = M##A##B##matrix_elem_add( \
                A##matrix_get(a, x, y), B##matrix_get(b, x, y)); \
        } \
    } \
} \
\
/* */

#define MATRIX_OP3_HC(M, A, B) \
    MATRIX_OP3_H(M, A, B) \
    MATRIX_OP3_C(M, A, B)

// ---------------------------------------------------------------

// iimatrix_* operations are for integer matrix
MATRIX_HC(ii, int, "d")
// ddmatrix_* operations are for double matrix
MATRIX_HC(dd, double, "f")
// asa above
MATRIX_HC(ff, float, "f")

// callback for adding int + double = float
float ffiiddmatrix_elem_add(int a, double b) {
    return a + b;
}
// define integer + double = float matrix operation
MATRIX_OP3_HC(ff, ii, dd)

int main() {
    // create an integer matrix from array
    iimatrix ii;
    iimatrix_create_from_array(&ii, 2, (int[2][2]){{1,2},{3,4}});
    printf("int matrix:\n");
    iimatrix_print(&ii);
    printf("\n");

    // create double matrix from array
    ddmatrix dd;
    ddmatrix_create_from_array(&dd, 2, (double[2][2]){{5,6},{7,8}});
    printf("double matrix:\n");
    ddmatrix_print(&dd);
    printf("\n");

    // create float matrix and add integer to double matrixes
    ffmatrix ff;
    ffmatrix_create(&ff, 2);
    ffiiddmatrix_add(&ff, &ii, &dd);

    printf("resulting float matrix:\n");
    ffmatrix_print(&ff);
    printf("\n");

    iimatrix_destroy(&ii);
    ddmatrix_destroy(&dd);
    ffmatrix_destroy(&ff);
}

The code outputs:

int matrix:
1 2
3 4

double matrix:
5.000000 6.000000
7.000000 8.000000

resulting float matrix:
6.000000 8.000000
10.000000 12.000000