5
\$\begingroup\$

As one of my first C programs I want to read in a matrix and compute its determinant. I don't pose limits on the size of the matrix and this makes things more complicated.

Version 0

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>

#define ROW_LENGTH 8
#define CHUNK 32

typedef struct {
    size_t size;
    double *elements;
} Matrix;

double determinant(Matrix *);
Matrix *parse_input();
Matrix *create_matrix(size_t size);
void free_matrix(Matrix *);
double at(Matrix *, int, int);
char *readline();


int main(void) {
    Matrix *M = parse_input();
    printf("%f\n", determinant(M));
    free_matrix(M);
    return 0;
}


double determinant(Matrix *M) {
    if (M->size == 1) {
        return M->elements[0];
    } else if (M->size == 2) {
        return at(M, 0, 0) * at(M, 1, 1) - at(M, 0, 1) * at(M, 1, 0);
    }

    // Make the matrix triangular
    size_t i, j, t;
    double r = 1;
    for (j = 0; j < M->size; j++) {
        if (!at(M, j, j)) return 0;
        for (i = j + 1; i < M->size; i++) {
            double ratio = at(M, i, j) / at(M, j, j);
            for (t = 0; t < M->size; t++) {
                M->elements[i * M->size + t] -= ratio * at(M, j, t);
            }
        }
    }
    for (i = 0; i < M->size; i++) {
        r *= at(M, i, i);
    }
    return r;
}


Matrix *parse_input() {
    char *row = readline();
    size_t t;
    size_t N = 0, P = 0;
    size_t i = 1, j = 0;

    double *first_row;
    if (!(first_row = malloc(ROW_LENGTH * sizeof first_row))) {
        puts("Could not allocate memory.");
        exit(EXIT_FAILURE);
    }
    char *number = strtok(row, " ");
    while (number) {
        if (N == ROW_LENGTH) {
            if (!(first_row = realloc(first_row, 2 * N * sizeof first_row))) {
                puts("Could not allocate memory.");
                free(first_row);
                exit(EXIT_FAILURE);
            }
        }
        first_row[N++] = atof(number);
        number = strtok(NULL, " ");
    }
    Matrix *M = create_matrix(N);
    for (t = 0; t < N; t++) {
        M->elements[t] = first_row[t];
    }
    free(row);
    free(first_row);
    while (++P < N) {
        j = 0;
        row = readline();
        char *number = strtok(row, " ");
        while (number) {
            M->elements[i * M->size + j++] = atof(number);
            number = strtok(NULL, " ");
        }
        i++;
        free(row);
    }
    return M;
}


Matrix *create_matrix(size_t size) {
    Matrix *M = malloc(sizeof(Matrix));
    M->size = size;
    M->elements = calloc(size * size, sizeof(double));
    return M;
}


void free_matrix(Matrix *matrix) {
    free(matrix->elements);
    free(matrix);
}


double at(Matrix *M, int i, int j) {
    return M->elements[i * M->size + j];
}


char *readline() {
    char *input = NULL;
    char tmpbuf[CHUNK];
    size_t inputlen = 0, tmplen = 0;

    do {
        fgets(tmpbuf, CHUNK, stdin);
        tmplen = strlen(tmpbuf);
        inputlen += tmplen;
        input = realloc(input, inputlen + 1);
        if (!input) {
            puts("Could not allocate memory.");
            exit(EXIT_FAILURE);
        }
        strcat(input, tmpbuf);
    } while (tmplen == CHUNK - 1 && tmpbuf[CHUNK - 2] != '\n');

    return input;
}

Version 1

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>

#define ROW_LENGTH 8
#define CHUNK 32

typedef struct {
    size_t size;
    double *elements;
} Matrix;

double determinant(Matrix *);
Matrix *parse_input();
Matrix *create_matrix(size_t);
void free_matrix(Matrix *);
double at(Matrix *, int, int);
char *readline();


int main(void) {
    Matrix *M = parse_input();
    printf("%f\n", determinant(M));
    free_matrix(M);
    return 0;
}


double determinant(Matrix *M) {
    if (M->size == 1) {
        return M->elements[0];
    } else if (M->size == 2) {
        return at(M, 0, 0) * at(M, 1, 1) - at(M, 0, 1) * at(M, 1, 0);
    }

    // Make the matrix triangular
    size_t i, j, t;
    double r = 1;
    for (j = 0; j < M->size; j++) {
        if (!at(M, j, j)) return 0;
        for (i = j + 1; i < M->size; i++) {
            double ratio = at(M, i, j) / at(M, j, j);
            for (t = 0; t < M->size; t++) {
                M->elements[i * M->size + t] -= ratio * at(M, j, t);
            }
        }
    }
    for (i = 0; i < M->size; i++) {
        r *= at(M, i, i);
    }
    return r;
}


Matrix *parse_input() {
    char *row = readline();
    size_t t;
    size_t N = 0;
    size_t i = 1, j = 0;
    size_t row_length = ROW_LENGTH;

    double *first_row;
    if (!(first_row = malloc(ROW_LENGTH * sizeof *first_row))) {
        puts("Could not allocate memory.");
        exit(EXIT_FAILURE);
    }
    char *number = strtok(row, " ");
    while (number) {
        if (N == row_length) {
            row_length *= 2;
            if (!(first_row = realloc(first_row, row_length * sizeof *first_row))) {
                puts("Could not allocate memory.");
                free(first_row);
                exit(EXIT_FAILURE);
            }
        }
        first_row[N++] = atof(number);
        number = strtok(NULL, " ");
    }
    Matrix *M = create_matrix(N);
    for (t = 0; t < N; t++) {
        M->elements[t] = first_row[t];
    }
    free(row);
    free(first_row);
    while (i < N) {
        j = 0;
        row = readline();
        char *number = strtok(row, " ");
        while (number && j < N) {
            M->elements[i * M->size + j++] = atof(number);
            number = strtok(NULL, " ");
        }
        i++;
        free(row);
    }
    return M;
}


Matrix *create_matrix(size_t size) {
    Matrix *M;
    if (!(M = malloc(sizeof *M))) {
        puts("Could not allocate memory.");
        exit(EXIT_FAILURE);
    }
    M->size = size;
    M->elements = calloc(size * size, sizeof(double));
    return M;
}


void free_matrix(Matrix *matrix) {
    free(matrix->elements);
    free(matrix);
}


double at(Matrix *M, int i, int j) {
    return M->elements[i * M->size + j];
}


char *readline() {
    char *input = calloc(CHUNK, 1);
    char tmpbuf[CHUNK];
    size_t inputlen = 0, tmplen = 0;

    do {
        fgets(tmpbuf, CHUNK, stdin);
        tmplen = strlen(tmpbuf);
        inputlen += tmplen;
        input = realloc(input, inputlen + 1);
        if (!input) {
            puts("Could not allocate memory.");
            exit(EXIT_FAILURE);
        }
        strcat(input, tmpbuf);
    } while (tmplen == CHUNK - 1 && tmpbuf[CHUNK - 2] != '\n');

    return input;
}

Version 2

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>

#define _M(i, j) (M->elements[(i) * M->size + (j)])
#define ROW_LENGTH 8
#define CHUNK 32

typedef struct {
    size_t size;
    double *elements;
} Matrix;

double determinant(Matrix *);
signed char find_pivot(Matrix *, int);
Matrix *parse_input();
Matrix *create_matrix(size_t);
void free_matrix(Matrix *);
int readline(char **, size_t *, FILE *);


int main(void) {
    Matrix *M = parse_input();
    printf("%f\n", determinant(M));
    free_matrix(M);
    return 0;
}


double determinant(Matrix *M) {
    if (M->size == 1) {
        return M->elements[0];
    } else if (M->size == 2) {
        return _M(0, 0) * _M(1, 1) - _M(0, 1) * _M(1, 0);
    }

    // Make the matrix triangular
    size_t i, j, t;
    signed char sign = 1;
    double ratio, r = 1;
    for (j = 0; j < M->size; j++) {
        if (!_M(j, j)) {
            if (!find_pivot(M, j)) {
                return 0;
            }
            sign *= -1;
        }
        for (i = j + 1; i < M->size; i++) {
            ratio = _M(i, j) / _M(j, j);
            for (t = 0; t < M->size; t++) {
                _M(i, t) -= ratio * _M(j, t);
            }
        }
    }
    for (i = 0; i < M->size; i++) {
        r *= _M(i, i);
    }
    return sign * r;
}


signed char find_pivot(Matrix *M, int j) {
    size_t i;
    for (i = j + 1; i < M->size; i++) {
        if (_M(i, j)) {
            size_t t;
            double tmp;
            for (t = 0; t < M->size; t++) {
                tmp = _M(i, t);
                _M(i, t) = M->elements[j * M->size + t];
                _M(j, t) = tmp;
            }
            return 1;
        }
    }
    return 0;
}


Matrix *parse_input() {
    char *row;
    size_t reading_size = CHUNK;
    if (!(row = malloc(reading_size))) {
        puts("Could not allocate memory.");
        exit(EXIT_FAILURE);
    }
    readline(&row, &reading_size, stdin);
    size_t t;
    size_t N = 0;
    size_t i = 1, j = 0;
    size_t row_length = ROW_LENGTH;

    double *first_row;
    if (!(first_row = malloc(ROW_LENGTH * sizeof *first_row))) {
        puts("Could not allocate memory.");
        exit(EXIT_FAILURE);
    }
    char *number = strtok(row, " ");
    while (number) {
        if (N == row_length) {
            row_length *= 2;
            if (!(first_row = realloc(first_row, row_length * sizeof *first_row))) {
                puts("Could not allocate memory.");
                free(row);
                free(first_row);
                exit(EXIT_FAILURE);
            }
        }
        first_row[N++] = atof(number);
        number = strtok(NULL, " ");
    }
    Matrix *M = create_matrix(N);
    for (t = 0; t < N; t++) {
        M->elements[t] = first_row[t];
    }
    free(first_row);
    while (i < N) {
        j = 0;
        readline(&row, &reading_size, stdin);
        char *number = strtok(row, " ");
        while (number && j < N) {
            M->elements[i * M->size + j++] = atof(number);
            number = strtok(NULL, " ");
        }
        i++;
    }
    free(row);
    return M;
}


Matrix *create_matrix(size_t size) {
    Matrix *M;
    if (!(M = malloc(sizeof *M))) {
        puts("Could not allocate memory.");
        exit(EXIT_FAILURE);
    }
    M->size = size;
    if (!(M->elements = calloc(size * size, sizeof(double)))) {
        puts("Could not allocate memory.");
        exit(EXIT_FAILURE);
    }
    return M;
}


void free_matrix(Matrix *matrix) {
    free(matrix->elements);
    free(matrix);
}


int readline(char **input, size_t *size, FILE *file) {
    char *offset;
    char *p;
    size_t old_size;

    // Already at the end of file
    if (!fgets(*input, *size, file)) {
        return EOF;
    }

    // Check if input already contains a newline
    if (p = strchr(*input, '\n')) {
        *p = 0;
        return 0;
    }

    do {
        old_size = *size;
        *size *= 2;
        if (!(*input = realloc(*input, *size))) {
            puts("Could not allocate memory.");
            free(*input);
            exit(EXIT_FAILURE);
        }
        offset = &((*input)[old_size - 1]);
    } while (fgets(offset, old_size + 1, file) &&
             offset[strlen(offset) - 1] != '\n');
    return 0;
}

How it works
The user inputs the first row of the matrix, the program counts the number of elements and determines how many rows remain (the matrix must be a square one). Then the determinant is computed by first reducing the matrix to triangular form and then computing the product of the elements on the main diagonal.

Example session:

./det
4 23 4
2 -5 2
45 2 40
330.000000

I'm quite proud of how I managed to get it finally correct (after many problems with allocations and memory leaks). However, I was wondering how can I improve readline() and parse_input(). The latter looks quite a mess.

\$\endgroup\$
  • \$\begingroup\$ I think your code will say that the matrix [0,0,1; 0,1,0; 1,0,0] has determinant 0, but it has determinant -1. To fix this you could try swapping rows/columns, or use the QR decomposition. \$\endgroup\$ – dmuir Apr 4 '15 at 11:15
  • \$\begingroup\$ @dmuir Damn you are right! Why doesn't my method work? I thought reducing the matrix to triangular form always worked. \$\endgroup\$ – rubik Apr 4 '15 at 11:26
  • \$\begingroup\$ Well not all matrices M have a LU decomposition in the sense that there are lower and upper triangular matrices L and U such that M = LU. They do however have a decomposition M = PL*U where P is a permutation matrix. This amounts to saying that if you meet a zero on the diagonal, you should swap the row you are on with a lower row that doesn't have a zero in the column. If you can't find such a row then the matrix has zero determinant. You might want to read up about LU decomposition, and an alternative QR decomposition. \$\endgroup\$ – dmuir Apr 4 '15 at 12:41
  • \$\begingroup\$ @dmuir Very interesting. I think the most straightforward way is to do the row swapping. I think now (version 2) the program is correct. \$\endgroup\$ – rubik Apr 4 '15 at 15:47
1
\$\begingroup\$

Bugs

I spotted a few bugs in your program.

  1. You are allocating the wrong size.

     if (!(first_row = malloc(ROW_LENGTH * sizeof first_row))) {
    

    Should be:

     if (!(first_row = malloc(ROW_LENGTH * sizeof *first_row))) {
    

    The same goes for the call to realloc later.

  2. Your program can't handle a row bigger than 16 elements. Right now, you only call realloc under this condition:

     if (N == ROW_LENGTH) {
    

    But that can only happen once. To allow for an infinite size, you need to keep track of the current allocation size and realloc every time N reaches the current allocation max.

  3. The first time you allocate input using realloc(NULL, size), you don't clear it before you use it. Remember that realloc(NULL, size) is equivalent to calling malloc. So input will be uninitialized but you call strcat on it right after.

  4. If the first row has N elements, you will create a NxN matrix. But when you read the second and subsequent rows, you don't limit the row length to N elements when you fill in the matrix. In particular, if the last row contains more than N elements, you will overflow your matrix.

Other things

  1. You check the return value of malloc in most places, but not in create_matrix().

  2. In readline, you are expanding the size of the reallocated buffer by a constant 32 bytes at a time. If the input is really huge, then this will take a long time, as this is an O(N^2) operation. This could be improved by doubling the size of the buffer each time. Also, you can simplify the function by getting rid of tmpbuf. Instead, call fgets directly on the unfilled part of the buffer (the part that realloc just expanded).

  3. Your at() function is an interesting way of accessing matrix elements. I think an even better way would be like this:

    #define _M(i,j) (M->elements[(i) * M->size + (j)])
    

    This assumes you will always name your matrix M, but you could modify it to be more general if you need to. The good thing about this macro is that you can use it as an lvalue, which you couldn't with your function. For example:

    _M(i,j) -= ratio * _M(j,t);
    
  4. You have a variable P in parse_input. I don't know what it stands for but it looks like it has the same exact value as i, which is the current row number. I think you could eliminate P.

\$\endgroup\$
  • \$\begingroup\$ Thanks for the review! I fixed all the issues except 2) and 3) in the second section. I also edited the question with an updated version. I read the documentation for fgets() but I didn't understand how can I get rid of tmpbuf. Can't I add the matrix name in the macro? Sort of like I did with the at() function. \$\endgroup\$ – rubik Apr 4 '15 at 8:57
  • \$\begingroup\$ @rubik Sure you can add the matrix name in the macro if you want. As far as fgets, imagine you just have one buffer. You read 32 bytes into it. Then you resize it to 64 bytes. You can now call fgets(&buf[32], 32, stdin) to read into the "tail" of the buffer. Now you resize to 128 bytes and do fgets(&buf[64], 64, stdin). Of course, I'm off by one each time because of the null terminator. So really it should be fgets(&buf[31], 33, stdin), etc. \$\endgroup\$ – JS1 Apr 4 '15 at 9:49
  • \$\begingroup\$ I rewrote readline() and changed its interface, mixing various inputs I found on the internet. I think now it's faster (I don't have to free every time) and the reading size doubles every time as you were saying. Thanks again! \$\endgroup\$ – rubik Apr 4 '15 at 15:48
  • \$\begingroup\$ @rubik You still have one calloc left to check, the one where you allocate M->elements. Everything else is looking good. \$\endgroup\$ – JS1 Apr 4 '15 at 17:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.