matrix determinant in c using Gauss elimination Core dumped error

Question

I'm trying to make a simple console application in C which will calculate the determinant of a Matrix using the Gauss elimination. after a lot of tests I found out that my program is not working because of the core dumped error.After 2 days of editing and undoing, i could not find the problem. Any help is more than welcomed.

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

int recherche_pivot(int k, int n, float *A)
{
    int i, j;
    if (A[((k - 1) * n + k) - 1] != 0)
    {
        return k;
    }
    else
    { //parcours du reste de la colonne
        for (i = k + 1; i <= n; i++)
        {
            if (A[((k - 1) * n + i) - 1] != 0)
            {
                return i;
            }
        }
        return -1;
    }
}

void fois(int n, float p, int i, float * A, float *b, float * x)
{
    int a;
    for (a = 1; a <= n; a++)
    {
        x[a - 1] = A[((i - 1) * n + a) - 1] * p;
    }
    x[n] = b[i - 1] * p;
}
void afficher_system(int n, float * X, float *b)
{
    int i, j;
    for (i = 1; i <= n; i++)
    {
        for (j = 1; j <= n; j++)
            printf("%f  ", X[((i - 1) * n + j) - 1]);
        printf(" | %f", b[i - 1]);
        printf("nn");
    }
    printf("nnnn");
}

void saisirmatrice(int n, float *A)
{
    int i, j;
    for (i = 1; i <= n; i++)
        for (j = 1; j <= n; j++)
            scanf("%f", &A[((i - 1) * n + j) - 1]);
}

void affichermatrice(int n, float *A)
{
    int i, j;
    for (i = 1; i <= n; i++)
        for (j = 1; j <= n; j++)
            printf("A[%d][%d] = %fn", i, j, A[((i - 1) * n + j) - 1]);
}

void elemination(int n, int k, float *b, float *A)
{
    int i, l, j;
    float * L, piv;
    L = (float *) malloc((n) * sizeof(float));
    for (i = k + 1; i <= n; i++)
    {
        piv = -1 * (A[((i - 1) * n + k) - 1] / A[((k - 1) * n + k) - 1]);
        fois(n, piv, k, A, b, L);
        //afficher_vecteur(n,L);
        for (j = 1; j <= n; j++)
        {
            A[((i - 1) * n + j) - 1] = A[((i - 1) * n + j) - 1] + L[j - 1];
        }
        b[i - 1] = b[i - 1] + L[n];
        afficher_system(n, A, b);
    }
}

void permutter(int n, float * A, int i, int j, float * b)
{
    int a;
    float t[n + 1];
    for (a = 1; a <= n; a++)
    {
        t[a - 1] = A[((i - 1) * n + a) - 1];
        A[((i - 1) * n + a) - 1] = A[((j - 1) * n + a) - 1];
        A[((j - 1) * n + a) - 1] = t[a - 1];
    }
    t[n] = b[i - 1];
    b[i - 1] = b[j - 1];
    b[j - 1] = t[n];
}

void main()
{
    float * A, det, *L, *R, *b, s;
    int i, j, i0, n, k, stop = 0;

    printf("Veuillez donner la taille de la matrice");
    scanf("%d", &n);
    A = (float *) malloc(sizeof(float) * (n * n));
    L = (float*) malloc(n * sizeof(float));
    R = (float*) malloc(n * sizeof(float));
    b = (float*) malloc(n * sizeof(float));
    printf("Veuillez remplir la matrice");
    saisirmatrice(n, A);
    det = 1;
    stop = 0;
    k = 1;
    do
    {
        do
        {
            i0 = recherche_pivot(k, n, A);
            if (i0 == k)
            {
                //Elémination
                elemination(n, k, b, A);
                k++;
            }
            else if (i0 == -1)
            {
                stop = 1;
            }
            else
            { //cas ou ligne pivot=i0 != k
              //permutation
                det = -det;
                permutter(n, A, k, i0, b);
                //elemination
                elemination(n, k, b, A);
                //afficher_matrice(n,A);
                k++;

            }
        } while ((k <= n) && (stop == 0));
    } while (stop == 1 || k == n);

    for (i = 1; i < n; i++)
    {
        det = det * A[((i - 1) * n + i) - 1];
    }

    printf("Le determinant est :%f", det);
    free(A);
    free(L);
    free(R);
    free(b);
}

Answer 1

There are many problems in the above code. Since arrays are zero-indexed in C, you should count the rows and columns of your matrices starting from zero, instead of counting from 1 and then attempting to convert when array-indexing. There is no need to cast the result of malloc(), and it is better to use an identifier rather than an explicit type as the argument for the sizeof operator:

A = malloc(sizeof(*A) * n * n));

You allocate space for L and R in main(), and then never use these pointers until the end of the program when they are freed. Then you allocate for L within the elemination() function; but you never free this memory, so you have a memory leak. You also allocate space for b in main(), but you don't store any values in b before passing it to the elemination() function. This is bound to cause problems.

There is no need for dynamic allocation here in the first place; I suggest using a variable length array to store the elements of the matrix. These have been available since C99, and will allow you to avoid all of the allocation issues.

There is a problem in the recherche_pivot() function, where you compare:

 if(A[((k - 1) * n + i) - 1] != 0) {}

This is a problem because the array element is a floating point value which is the result of arithmetic operations; this value should not be directly compared with 0. I suggest selecting an appropriate DELTA value to represent a zero range, and instead comparing:

#define DELTA  0.000001
...
if (fabs(A[((k - 1) * n + i) - 1]) < DELTA) {}

In the permutter() function you use an array, float t[n];, to hold temporary values. But an array is unnecessary here since you don't need to save these temporary values after the swap; instead just use float t;. Further, when you interchange the values in b[], you use t[n] to store the temporary value, but this is out of bounds.

The elemination() function should probably iterate over all of the rows (excepting the kth row), rather that starting from the kth row, or it should start at the k+1th row. As it is, the kth row is used to eliminate itself. Finally, the actual algorithm that you use to perform the Gaussian elimination in main() is broken. Among other things, the call permutter(n, A, k, i0, b); swaps the kth row with the i0th row, but i0 is the pivot column of the kth row. This makes no sense.

It actually looks like you want to do more than just calculate determinants with this code, since you have b, which is the constant vector of a linear system. This is not needed for the task alluded to in the title of your question. Also, it appears that your code gives a result of 1 for any 1X1 determinant. This is incorrect; it should be the value of the single number in this case.

The Gaussian elimination method for calculating the determinant requires that you keep track of how many row-interchanges are performed, and that you keep a running product of any factors by which individual rows are multiplied. Adding a multiple of one row to another row to replace that row does not change the value of the determinant, and this is the operation used in the reduce() function below. The final result is the product of the diagonal entries in the reduced matrix, multiplied by -1 once for every row-interchange operation, divided by the product of all of the factors used to scale individual rows. In this case, there are no such factors, so the result is simply the product of the diagonal elements of the reduced matrix, with the sign correction. This is the method used by the code posted in the original question.

There were so many issues here that I just wrote a fresh program that implements this algorithm. I think that it is close, at least in spirit, to what you were trying to accomplish. I did add some input validation for the size of the matrix, checking to be sure that the user inputs a positive number, and prompting for re-entry if the input is bad. The input loop that fills the matrix would benefit from similar input validation. Also note that the input size is stored in a signed int, to allow checks for negative input, and a successful input is cast and stored in a variable of type size_t, which is an unsigned integer type guaranteed to hold any array index. This is the correct type to use when indexing arrays, and you will note that size_t is used throughout the program.

#include <stdio.h>
#include <math.h>
#include <stdbool.h>

#define DELTA 0.000001

void show_matrix(size_t mx_sz, double mx[mx_sz][mx_sz]);
void interchange(size_t r1, size_t r2, size_t mx_sz, double mx[mx_sz][mx_sz]);
void reduce(double factor, size_t r1, size_t r2,
            size_t mx_sz, double mx[mx_sz][mx_sz]);
size_t get_pivot(size_t row, size_t mx_sz, double mx[mx_sz][mx_sz]);
double find_det(size_t mx_sz, double mx[mx_sz][mx_sz]);

int main(void)
{
    size_t n;
    int read_val, c;

    printf("Enter size of matrix: ");
    while (scanf("%d", &read_val) != 1 || read_val < 1) {
        while ((c = getchar()) != '\n' && c != EOF) {
            continue;               // discard extra characters
        }
        printf("Enter size of matrix: ");
    }
    n = (size_t) read_val;

    double matrix[n][n];

    printf("Enter matrix elements:\n");
    for (size_t i = 0; i < n; i++) {
        for (size_t j = 0; j < n; j++) {
            scanf("%lf", &matrix[i][j]);
        }
    }

    printf("You entered:\n");
    show_matrix(n, matrix);
    putchar('\n');

    double result = find_det(n, matrix);
    show_matrix(n, matrix);
    putchar('\n');
    printf("Determinant: %f\n", result);

    return 0;
}

void show_matrix(size_t n, double mx[n][n])
{
    for (size_t i = 0; i < n; i++) {
        for (size_t j = 0; j < n; j++) {
            printf("%7.2f", mx[i][j]);
        }
        putchar('\n');
    }
}

/* interchange rows r1 and r2 */
void interchange(size_t r1, size_t r2, size_t mx_sz, double mx[mx_sz][mx_sz])
{
    double temp;

    for (size_t j = 0; j < mx_sz; j++) {
        temp = mx[r1][j];
        mx[r1][j] = mx[r2][j];
        mx[r2][j] = temp;
    }
}

/* add factor * row r1 to row r2 to replace row r2 */
void reduce(double factor, size_t r1, size_t r2,
            size_t mx_sz, double mx[mx_sz][mx_sz])
{
    for (size_t j = 0; j < mx_sz; j++) {
        mx[r2][j] += (factor * mx[r1][j]);
    }
}

/* returns pivot column, or mx_sz if there is no pivot */
size_t get_pivot(size_t row, size_t mx_sz, double mx[mx_sz][mx_sz])
{
    size_t j = 0;

    while (j < mx_sz && fabs(mx[row][j]) < DELTA) {
        ++j;
    }

    return j;
}

double find_det(size_t mx_sz, double mx[mx_sz][mx_sz])
{
    size_t pivot1, pivot2;
    size_t row;
    double factor;
    bool finished = false;
    double result = 1.0;

    while (!finished) {
        finished = true;
        row = 1;
        while (row < mx_sz) {
            // determinant is zero if there is a zero row
            if ((pivot1 = get_pivot(row - 1, mx_sz, mx)) == mx_sz ||
                (pivot2 = get_pivot(row, mx_sz, mx)) == mx_sz) {
                return 0.0;
            }
            if (pivot1 == pivot2) {
                factor = -mx[row][pivot1] / mx[row - 1][pivot1];
                reduce(factor, row - 1, row, mx_sz, mx);
                finished = false;
            } else if (pivot2 < pivot1) {
                interchange(row - 1, row, mx_sz, mx);
                result = -result;
                finished = false;
            }
            ++row;
        }
    }

    for (size_t j = 0; j < mx_sz; j++) {
        result *= mx[j][j];
    }

    return result;
}

Sample session:

Enter size of matrix: oops
Enter size of matrix: 0
Enter size of matrix: -1
Enter size of matrix: 3
Enter matrix elements:
0 1 3
1 2 0
0 3 4
You entered:
   0.00   1.00   3.00
   1.00   2.00   0.00
   0.00   3.00   4.00

   1.00   2.00   0.00
  -0.00  -3.00  -9.00
   0.00   0.00  -5.00

Determinant: 5.000000