# Gauss elimination

The gauss function takes 3 parameters:

• a pointer to the matrix of coefficients
• the number of variables
• a pointer to matrix of solution

#include <stdio.h>
#include <stdlib.h>
int gauss(double* matrix , int x, double* sol);
int main(){
int x;
scanf("%i",&x);
double *mat = malloc((x+1)*x*sizeof(double)), *sol =(malloc(x*sizeof(double)));
for(int i=0;i<(x+1)*x; i++) scanf("%lf",&mat[i]);
int state = gauss(mat, x, sol);
if (state ==0) for(int i=0;i<x;i++) printf("X%i = %lf\n",i,sol[i]);
else printf("The system has %s%s", state==-1 ? "infinte number of":"no","solution(s)");
free(mat);
free(sol);
return 0;
}
int gauss(double* matrix , int x, double* sol){
for(int i=0;i<x-1;i++)for(int j=i+1; j<x; j++) for(int k =x;k>=0;k--){
matrix[j*(x+1)+k] -= matrix[i*(x+1)+k]*matrix[j*(x+1)+i]/matrix[i*(x+2)];
}
if(matrix[x*(x+1)-2 ] ==0){
if(matrix[x*(x+1)-1]==0) return -1;
return 1;
}
sol[x-1]= matrix[(x-1)*(x+1)+x]/matrix[(x-1)*(x+1)+x-1];
for(int i =x-2; i>=0; i--){
double counter=0;
for(int j=x-1;j>i;j--) counter+= matrix[i*(x+1)+j] *sol[j];
sol[i] = (matrix[i*(x+1)+x]-counter) / matrix[i*(x+2)];
}
return 0;
}


I was wondering what I can do to make the code shorter (keeping the good style), with better performance and with better memory consumption.

I used linear matrix because I simply don't know how to declare multidimensional array with malloc and I don't know how to pass it as a parameter to function.

I used indent -kr -i8 -l1000 -lc1000 for styling and I used preprocessor for simplifying the matrix math but some pieces make no sense as as matrix(x, -2):

#include <stdio.h>
#include <stdlib.h>
#define matrix(i,j) matrix[i*(x+1)+j]
int gauss(double *matrix, int x, double *sol);
int main()
{
int x;
scanf("%i", &x);
double *mat = malloc((x + 1) * x * sizeof(double)), *sol = (malloc(x * sizeof(double)));
for (int i = 0; i < (x + 1) * x; i++)
scanf("%lf", &mat[i]);
int state = gauss(mat, x, sol);
if (state == 0)
for (int i = 0; i < x; i++)
printf("X%i = %lf\n", i, sol[i]);
else
printf("The system has %s solutions", state == -1 ? "infinte number of" : "no");
free(mat);
free(sol);
return 0;
}
int gauss(double *matrix, int x, double *sol)
{
for (int i = 0; i < x - 1; i++)
for (int j = i + 1; j < x; j++)
for (int k = x; k >= 0; k--)
matrix(j, k) -= matrix(i, k) * matrix(j, i) / matrix(i, i);
if (matrix(x, -2) == 0) {
if (matrix(x, -1) == 0)
return -1;
return 1;
}
sol[x - 1] = matrix((x - 1), x) / matrix((x - 1), (x - 1));
for (int i = x - 2; i >= 0; i--) {
double counter = 0;
for (int j = x - 1; j > i; j--)
counter += matrix(i, j) * sol[j];
sol[i] = (matrix(i, x) - counter) / matrix(i, i);
}
return 0;
}


Is there anything I can do to make this code better?

• I meant by shortest code "less instructions" thats it. and I it is my habit to use less space in math equations so that I don't see the equation in 2 lines in my text editor – oddcoder Feb 3 '15 at 12:46
• In what sense are you reinventing-the-wheel? Is there an existing library function that you want to reimplement? – 200_success Feb 3 '15 at 20:47
• not a specific library but of coarse Gauss is implemented in many languages and I am trying to do this in c – oddcoder Feb 4 '15 at 11:39
• reinventing-the-wheel tag removed. – 200_success Feb 4 '15 at 15:39

• The code has numerical problems.

• The elimination loop blindly divides by matrix((x - 1), (x - 1)) which very well could be 0. Note that a mere presence of 0 at the main diagonal doesn't mean that the system is unsolvable.

• Along the same line, division by 0 is not the only way to end up with Inf. Same will happen when dividing a pretty large number by a pretty small one.

• Even if such division doesn't yield Inf, dividing large by small amplifies errors.

To summarize, a correct implementation must have pivot selection.

• The code has logical problems.

The solvability is decided by looking at just one (or two) values. It is absolutely not enough. Later on, while calculating solutions a division by 0 still may occur. You need to test matrix(i, i) and (matrix(i, x) - counter) at line 38 at each iteration.

• Structure

Don't be shy of functions. Identifying functions not just increase readability, but makes the code reusable and unit-testable.

• gauss is logically divided into 2 algorithms: first, calculate the upper triangular form of a matrix, second, solve the triangular form. That is better to be reflected in the code:

int gauss(double * matrix, int size, double * solution)
{
triangularize(matrix, size);
return solve_triangular_form(matrix, size, solution);
}

• No raw loops. Each loop is there for a purpose. You must understand what that purpose is, give it a name and make it into a function. In fact, I highly recommend the reversed process. For example, I would write the following snippet before thinking of how I'd implement the pieces:

error_t triangularize(double * matrix, int size)
{
for (int row = 0; row < size; ++row) {
int pivot = select_pivot_column(matrix, row, size);
swap_columns(matrix, row, pivot);
int error_t = eliminate_column(matrix, row, size);
if (rc != ERROR_NONE) {
return rc;
}
}
return ERROR_NONE;
}


I think your style is okay, and the algorithm seems nice and short. Good work!

I find short code (that has not been deliberately obfuscated) is almost always easier to understand than long-winded code. Sometimes when I read code, I will have a quick look at the comments, then strip them out so I can see the code more clearly!

Here's some suggestions to improve it. I am not considering the algorithm.

You seem to be using C99 features such as "declare anywhere" (which is arguably not such a great idea)... but anyway, here is how to do dynamic multi-dimensional arrays in C99 with malloc: https://stackoverflow.com/a/12805980/218294 This is the biggest improvement you can make to your code.

You can get rid of malloc and put your arrays on the stack, like this:

#include <stdio.h>

void print_matrix(int rows, int cols, double M[rows][cols])
{
for (int i=0; i<rows; ++i) {
for (int j=0; j<cols; ++j)
printf("%lf ", M[i][j]);
printf("\n");
}
}

int main(void)
{
int rows = 5;
int cols = rows+1;
double M[rows][cols];
for (int i=0; i<rows; ++i)
for (int j=0; j<cols; ++j)
M[i][j] = i+j;
print_matrix(rows, cols, M);
return 0;
}


If C99 did not include these features for dynamic multi-dimensional arrays, it might be a good idea to write macros to clarify the array access, e.g. #define Matrix(i,j) matrix[i*(x+1)+j] then instead of matrix[j*(x+1)+k] you could write Matrix(j,k).

I suggest adding some spaces at the top level of expressions, this would be easier to read, e.g.:

matrix[i*(x+1)+k] * matrix[j*(x+1)+i] / matrix[i*(x+2)]


sol =(malloc(x*sizeof(double))) - no need for outer parens here. You have inconsistent spacing around = here (and similar elsewhere), just a minor thing.

The second %s in printf("The system has %s%s", ... is unnecessary, as the argument is a string literal - just put that string into the format; also put a space before it.

Since your style includes long lines, you might free both matrices on the same line: free(mat); free(sol);

You can exclude the braces in the triple for loop, as there is only one statement inside the loop. You are doing this elsewhere, so do it consistently.

Try to follow K&R (Kernighan and Ritchie) code style. I used this shell command:

indent -kr -i8 -l1000 -lc1000


here's what this program (GNU Indent) did to your code (no copyright infringement intended!) The code is more readable. I think it put too much space inside the array indexing [brackets], you could reduce that again (but better use the multi-dimensional arrays).

#include <stdio.h>
#include <stdlib.h>
int gauss(double *matrix, int x, double *sol);
int main()
{
int x;
scanf("%i", &x);
double *mat = malloc((x + 1) * x * sizeof(double)), *sol = (malloc(x * sizeof(double)));
for (int i = 0; i < (x + 1) * x; i++)
scanf("%lf", &mat[i]);
int state = gauss(mat, x, sol);
if (state == 0)
for (int i = 0; i < x; i++)
printf("X%i = %lf\n", i, sol[i]);
else
printf("The system has %s%s", state == -1 ? "infinte number of" : "no", "solution(s)");
free(mat);
free(sol);
return 0;
}

int gauss(double *matrix, int x, double *sol)
{
for (int i = 0; i < x - 1; i++)
for (int j = i + 1; j < x; j++)
for (int k = x; k >= 0; k--) {
matrix[j * (x + 1) + k] -= matrix[i * (x + 1) + k] * matrix[j * (x + 1) + i] / matrix[i * (x + 2)];
}
if (matrix[x * (x + 1) - 2] == 0) {
if (matrix[x * (x + 1) - 1] == 0)
return -1;
return 1;
}
sol[x - 1] = matrix[(x - 1) * (x + 1) + x] / matrix[(x - 1) * (x + 1) + x - 1];
for (int i = x - 2; i >= 0; i--) {
double counter = 0;
for (int j = x - 1; j > i; j--)
counter += matrix[i * (x + 1) + j] * sol[j];
sol[i] = (matrix[i * (x + 1) + x] - counter) / matrix[i * (x + 2)];
}
return 0;
}


If you revise your code, please post it again when you have finished.

• "You can exclude the braces in the triple for loop, as there is only one statement inside the loop." – It seems to me that the opposite advice is given more frequently :) – Martin R Feb 3 '15 at 13:49
• putting the array in stack will cause stack overflow incase of solving large number of equations so I used heaps to get rid of this problem. – oddcoder Feb 3 '15 at 13:55