inverse.c

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

typedef struct {
    size_t width, height;
    double data[];
} matrix;

matrix *matrix_create(size_t width, size_t height);
matrix *matrix_create_identity(size_t n);
matrix *matrix_clone(matrix *mat);
void matrix_delete(matrix *mat);
double *matrix_at(matrix *mat,size_t x, size_t y);
matrix *matrix_invert(matrix *mat);
void matrix_fprint(matrix *mat, FILE *stream);

int main() {
    uint dimension;
    printf("Enter matrix size\n");
    scanf("%i", &dimension);

    // Create our matrix and populate it.
    matrix *mat = matrix_create(dimension, dimension);
    for(size_t y = 0; y < dimension; y++) for(size_t x = 0; x < dimension; x++) {
        printf("Enter matrix coefficient a[%i][%i]\n", (int) y + 1, (int) x + 1);
        scanf("%lf", matrix_at(mat, x, y));
    }

    // Invert matrix and print it.
    matrix *inv = matrix_invert(mat);
    if(inv == NULL) {
        printf("Matrix uninvertable.\n");
    } else {
        printf("Output matrix: \n");
        matrix_fprint(inv, stdout);
        matrix_delete(inv);
    }

    matrix_delete(mat);
}

// Create a matrix width by height, contents are uninitialised.
matrix *matrix_create(size_t width, size_t height) {
    matrix *mat = malloc(sizeof(matrix) + sizeof(double) * width * height);
    mat -> width = width;
    mat -> height = height;
    return mat;
}

// Create an n by n identity matrix.
matrix *matrix_create_identity(size_t n) {
    matrix *mat = matrix_create(n, n);
    memset(mat -> data, 0, sizeof(double) * n * n);
    for(size_t i = 0; i < n; i++) *matrix_at(mat, i, i) = 1;
    return mat;
}

// Clone the matrix in mat into a new matrix and return this matrix.
matrix *matrix_clone(matrix *mat) {
    matrix *clone = matrix_create(mat -> width, mat -> height);
    memcpy(clone -> data, mat -> data, sizeof(double) * mat -> width * mat -> height);
    return clone;
}

// Delete a matrix, freeing the memory.
void matrix_delete(matrix *mat) {
    free(mat);
}

// Return the memory location of the element (x, y) for the matrix mat.
double *matrix_at(matrix *mat, size_t x, size_t y) {
    return &(mat -> data[y * mat -> width + x]);
}

// Invert matrix input using Gaussian Elimination, creating a new matrix. The original matrix is unaffected.
matrix *matrix_invert(matrix *input) {
    if(input -> width != input -> height) {
        fprintf(stderr, "Matrix to be inverted has invalid dimensions (width != height).\n");
        return NULL;
    }
    matrix *mat = matrix_clone(input), *inv = matrix_create_identity(input -> width), *tmp = matrix_create(input -> width, 1);
    for(size_t i = 0; i < input -> width; i++) {
        if(*matrix_at(mat, i, i) == 0) {
            // If our matrix has a zero in (i, i) we must chose a new y such that (i, y) is non-zero and swap all x in (x, y) and (x, i)
            size_t swap = i;
            for(;*matrix_at(mat, i, swap) == 0; swap++);
            if(swap == input -> width) {
                // No inverse
                matrix_delete(mat);
                matrix_delete(tmp);
                matrix_delete(inv);
                return NULL;
            }
            // Swap using the temporary matrix row tmp.
            memcpy(matrix_at(tmp, 0, 0),    matrix_at(mat, 0, i),    sizeof(double) * input -> width);
            memcpy(matrix_at(mat, 0, i),    matrix_at(mat, 0, swap), sizeof(double) * input -> width);
            memcpy(matrix_at(mat, 0, swap), matrix_at(tmp, 0, 0),    sizeof(double) * input -> width);
            memcpy(matrix_at(tmp, 0, 0),    matrix_at(inv, 0, i),    sizeof(double) * input -> width);
            memcpy(matrix_at(inv, 0, i),    matrix_at(inv, 0, swap), sizeof(double) * input -> width);
            memcpy(matrix_at(inv, 0, swap), matrix_at(tmp, 0, 0),    sizeof(double) * input -> width);
        }
        // Divide row i by the number in (i, i), normalising it
        for(size_t x = 0; x < input -> width; x++) {
            if(x > i) *matrix_at(mat, x, i) /= *matrix_at(mat, i, i);
            *matrix_at(inv, x, i) /= *matrix_at(mat, i, i);
        }
        *matrix_at(mat, i, i) = 1;
        // Set all numbers in column i to be zero, excluding (i, i) by adding row i times - (i, y)
        for(size_t y = 0; y < input -> width; y++) {
            if(y == i) continue;
            if(*matrix_at(mat, i, y) == 0) continue;
            double mul = - *matrix_at(mat, i, y);
            for(size_t x = 0; x < input -> width; x++) {
                if(x >= i) *matrix_at(mat, x, y) += mul * *matrix_at(mat, x, i);
                *matrix_at(inv, x, y) += mul * *matrix_at(inv, x, i);
            }
        }
    }
    matrix_delete(tmp);
    matrix_delete(mat);
    return inv;
}

// Print the contents of a matrix mat onto the stream provided.
void matrix_fprint(matrix *mat, FILE *stream) {
    for(size_t y = 0; y < mat -> height; y++) {
        for(size_t x = 0; x < mat -> width; x++) {
            fprintf(stream, "%+6lf\t", *matrix_at(mat, x, y));
        }
        fprintf(stream, "\n");
    }
}