Function to solve a system using Gauss-Seidel

Question

I'm trying to write a function in order to solve a system of equations using the Gauss-Seidel method. The function I have so far doesn't do the job. The code stops giving me output after "Enter an initial guess..." and doesn't solve the function I'm inputting. Is this an issue with my function or some other part of my code?

/* code to solve a nxn system using the Gauss-Seidel method */
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define MAX_DIM 100
#define MAX_ITER 500
#define TOLERANCE 1.e-6

void gauss_seidel(double a[][MAX_DIM], double b[], double x[], int n);
void main()
{       
    int i, j, n;
    int violation_counter, answer;
    int violation_rows[MAX_DIM];
    double sum;
    double a[MAX_DIM][MAX_DIM];
    double b[MAX_DIM], x[MAX_DIM];

    /* read in data */
    n = MAX_DIM + 1;
    while (n > MAX_DIM) {
        printf("Enter the dimension of the system to be solved: ");
            scanf_s("%d", &n);
    }
    printf("\n");
    for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++) {
            printf("Enter a[%d][%d] element of the system matrix: ",i,j);
                scanf_s("%lf", &a[i][j]);
        }
        printf("Enter b[%d] of the right-hand-side vector: ", i);
        scanf_s("%lf", &b[i]);
        printf("\n");
    }

    /* test the convergence criterion */
    violation_counter = 0;
    for (i = 0; i < n; i++) {
        sum = 0.0;
        for (j = 0; j < n; j++)
            if (i != j)
                sum = sum + fabs(a[i][j]);
        if (fabs(a[i][i]) < sum) {
            violation_rows[violation_counter] = i;
            violation_counter = violation_counter + 1;
        }
        if (a[i][i] == 0.0) {
            printf("Found diagonal element equal to zero; rearrange equations; exiting ...\n");
                exit(0);
        }
    }
    if (violation_counter > 0) {
        printf("The Gauss-Seidel convergence criterion is violated in %d rows out of %d\n", violation_counter, n);
            printf("Specifically, it was violated in rows:\n");
        for (i = 0; i < violation_counter; i++)
            printf("%d ", violation_rows[i]);
        printf("\n");
        printf("Enter 1 if you want to continue; any other number to abort : ");
            scanf_s("%d", &answer);
        if (answer != 1)
            exit(1);
        printf("Check results carefully\n\n");
    }

    /* initialize the solution vector -- initial guesses */
    for (i = 0; i < n; i++) {
        printf("Enter an initial guess for x[%d] of the solution vector : ",i);
            scanf_s("%lf", &x[i]);
    }

    /* solve the system */
    gauss_seidel(a, b, x, n);

    /* output solution */
    for (i = 0; i < n; i++)
        printf("x[%d]=%f\n", i, x[i]);
    printf("\n");
}

/* function to solve a system using Gauss-Seidel */
void gauss_seidel(double a[][MAX_DIM], double b[], double x[], int n)
{
    double maxerror = 1.0;
    double iteration_error;
    double e, sum, temp;
    int i, j;

    while (maxerror > 1.e-6) {
        iteration_error = 0.0;

        for (i = 0; i < n; i++) {
            sum = 0.0;
            for (j = 0; j < n; j++) {
                if (i != j)
                    sum = sum + (a[i][j] * x[j]);
            }
        }

        temp = (a[i][n] - sum) / a[i][i];
        e = fabs((temp - x[i]) / x[i]);
        x[i] = temp;

        if (e > iteration_error)
            iteration_error = e;
    }
    maxerror = iteration_error;
}

Please improve the indentation of your code. In the current state it's hard to read. What do you mean by "It stops at..." ? Did you run it in a debugger? Did you get some error message? — Dan Mašek, Apr 03 '16 at 23:54
In gauss_seidel(), you're only changing maxerror after the while loop, so the while loop will never finish. — Georg, Apr 04 '16 at 04:56
@Georg is that the reason why the output stops at "Enter an initial guess..."?. I'm assuming that's a fault in the function then? — DanielRossi, Apr 04 '16 at 05:22
You should check the return value from `scanf_s()` each time you use it. Yes, if the loop with 'Enter an initial guess' completes, the next step jumps into `gauss_seidel()` and, as Georg says, the function doesn't return because 1.0E0 > 1.0E-6 permanently. — Jonathan Leffler, Apr 04 '16 at 06:01

Function to solve a system using Gauss-Seidel

0 Answers0