Tutaj masz zarys algorytmu. Aha, mój kod działa in-place, więc nie da się zastosować bezpośrednio do Twojego przypadku. Relatywnie niewielka modyfikacja powinna Ci też umożliwić rozwiązywanie układu równań.
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
/* Row-major order */
#define IDX(i, j, n) ((i) * (n) + (j))
#define ERR_MATRIX_SINGULAR 1
#ifdef WITH_PIVOTING
static void
swap_rows(double *matrix, ssize_t n, ssize_t p, ssize_t q)
{
/*
* Swaps row p-th and q-th starting on first non-zero entry.
*/
ssize_t ii = p < q ? p : q;
for (; ii < n; ii++) {
double tmp;
tmp = matrix[IDX(p, ii, n)];
matrix[IDX(p, ii, n)] = matrix[IDX(q, ii, n)];
matrix[IDX(q, ii, n)] = tmp;
}
}
#endif
int
gauss(double *matrix, ssize_t n)
{
ssize_t ii, jj, kk;
double a_inv;
for (ii = 0; ii < n; ii++) {
#ifdef WITH_PIVOTING
ssize_t pivot_row;
double pivot;
pivot_row = ii;
pivot = matrix[IDX(ii, ii, n)];
for (jj = ii + 1; jj < n; jj++) {
if (fabs(matrix[IDX(jj, ii, n)]) > fabs(pivot)) {
pivot_row = jj;
pivot = matrix[IDX(pivot_row, ii, n)];
}
}
/*
* Swap rows pivot_row and ii-th.
*/
if (ii != pivot_row)
swap_rows(matrix, n, ii, pivot_row);
#endif
if (matrix[IDX(ii, ii, n)] == 0.0)
return (ERR_MATRIX_SINGULAR);
a_inv = 1.0 / matrix[IDX(ii, ii, n)];
/* Scale rows below i-th */
for (kk = ii + 1; kk < n; kk++) {
double scale;
scale = a_inv * matrix[IDX(kk, ii, n)];
for (jj = ii; jj < n; jj++) {
matrix[IDX(kk, jj, n)] -= scale *
matrix[IDX(ii, jj, n)];
}
}
/* Scale i-th row. */
for (jj = ii; jj < n; jj++)
matrix[IDX(ii, jj, n)] *= a_inv;
}
return (0);
}
int
main()
{
size_t ii, jj;
size_t test_n = 4;
double test_matrix[] = {
4.0, -2.0, 4.0, -2.0,
3.0, 1.0, 4.0, 2.0,
2.0, 4.0, 2.0, 1.0,
2.0, -2.0, 4.0, 2.0
};
if (gauss(test_matrix, test_n) == ERR_MATRIX_SINGULAR) {
printf("Invalid matrix.\n");
}
for (ii = 0; ii < test_n; ii++) {
for (jj = 0; jj < test_n; jj++) {
printf("%02.5lf ", test_matrix[IDX(ii, jj, test_n)]);
}
printf("\n");
}
return (0);
}