Created
February 11, 2016 12:25
-
-
Save meiamsome/0919e0b04562a6ae668b to your computer and use it in GitHub Desktop.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#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"); | |
} | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment