/**
 * 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;
}