Gauss elimination

Question

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?

Answer 1

• 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;
    }
    

Answer 2

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.