/**

 * test/demos/gauss.c

 * Copyright (C) 2008 Markus Broeker

 */

#define VERSION 0.8

#ifndef VERBOSE

#define VERBOSE 1

#endif

#ifndef ZEROASSUMPTION

#define ZEROASSUMPTION 1E-3     /* 0.001 = 0.0 */

#endif

#include <stdio.h>

#include <stdlib.h>

void ausgabe (float **A, int MAXX, int MAXY)

{

    int n, k;

    for (k = 0; k < MAXY; k++) {

        for (n = 0; n < MAXX - 1; n++) {

            if ((A[k][n] < ZEROASSUMPTION) && (A[k][n] > -ZEROASSUMPTION))

                A[k][n] = 0.0;

            printf ("%8.2f ", A[k][n]);

        }

        printf ("=%8.2f\n", A[k][n]);

    }

}

void gauss (float **A, int MAXX, int MAXY)

{

    int i, j, k, n;

    float h;

    i = 0;

    for (j = 0; j < MAXY; j++) {

        for (k = j + 1; k < MAXY; k++) {

            if (A[k][i] != 0)

                h = A[j][i] / A[k][i];

            else

                continue;

            for (n = 0; n < MAXX; n++) {

                A[k][n] *= -h;

                A[k][n] += A[j][n];

            }

        }

        i++;

    }

#ifdef PUPILS

    printf ("\n");

    ausgabe (A, MAXX, MAXY);

    printf ("\n");

#endif

    i--;                        /* i = (MAXX-1)-1; */

    for (j = MAXY - 1; j > -1; j--) {

        for (k = j - 1; k > -1; k--) {

            h = A[j][i] / A[k][i];

            for (n = MAXX - 1; n > k - 1; n--) {

                A[k][n] *= -h;

                A[k][n] += (A[j][n]);

                if (A[j][i] == 0.0) {

                    printf ("Divison by Zero\n");

                    exit (-1);

                }

                A[j][n] /= (A[j][i]);

            }

        }

        i--;

    }

    if ((h = A[0][0]) == 0.0) {

        printf ("Division by Zero\n");

        exit (-1);

    }

    /*

     * h=A[0][0];

     */

    A[0][0] /= h;

    A[0][MAXX - 1] /= h;

}

int main (int argc, char **argv)

{

    int i, k;

    int MAXX, MAXY;

    float value;

    float **A;

    if (argc != 2) {

        printf ("Benutzung: %s UNBEKANNTE\n", argv[0]);

        printf ("Berechnet ein lineares Gleichungssystem mit n UNBEKANNTEN und n Zeilen ( nxn ).\n");

        printf ("Praktische Beispiele stehen in der README Datei.\n");

        return EXIT_FAILURE;

    }

    MAXY = atoi (argv[1]);

    MAXX = MAXY + 1;

    if ((A = calloc (MAXY, sizeof (float *))) == NULL) {

        printf ("Nicht genug Speicher verfuegbar\n");

        return EXIT_FAILURE;

    }

    for (i = 0; i < MAXY; i++) {

        if ((A[i] = calloc (MAXX, sizeof (float))) == NULL) {

            printf ("Nicht genug Speicher verfuegbar\n");

            for (k = 0; k < i; k++)

                free (A[k]);

            free (A);

            return EXIT_FAILURE;

        }

    }

    for (k = 0; k < MAXY; k++) {

        printf ("Spalte %d: ", k);

        for (i = 0; i < MAXX; i++) {

            if (scanf ("%f", &value) < 1) {

                printf ("READ ERROR\n");

                return EXIT_FAILURE;

            }

            A[k][i] = value;

        }

    }

    if (VERBOSE) {

        printf ("\n");

        ausgabe (A, MAXX, MAXY);

    }

    gauss (A, MAXX, MAXY);

    printf ("\nLoesung des GL-Systems( %d x %d )\n\n", MAXX - 1, MAXY);

    for (i = 0; i < MAXY; i++)

        printf ("a[%2d] = %8.2f\n", i, A[i][MAXX - 1]);

    for (i = 0; i < MAXY; i++)

        if (A[i] != NULL)

            free (A[i]);

    if (A != NULL)

        free (A);

    return EXIT_SUCCESS;

}