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

#define SPEICHER_F "Kein Platz mehr für Speicher\n"


// Funktion showMatrix
// zeigt Matrix an
void showMatrix(double **mx, int rows, int cols) {

	int i, j;
	printf("\n");
	for(i = 0; i < rows; i++) {
		printf("  |  ");
		for(j = 0; j < cols; j++) {
			printf("%9f ", mx[i][j]);
			}
		printf("|\n");
		}
	printf("\n");

}


// Funktion showVector
// zeigt Vektor an
void showVector(double *ve, int n) {

	int i;
	printf("\n");
	for(i = 0; i < n; i++)
		printf("  |  %9f  |\n", ve[i]);
	printf("\n");

}


// Funktion makeHilbert
// liefert Hilbert-Matrix zurück
double **makeHilbert(int n) {

	double **hmx;
	hmx = (double **) calloc(n, sizeof(double *));
	if(hmx == NULL) printf(SPEICHER_F);

	int i, j;
	for(i = 1; i <= n; i++) {
		hmx[i-1] = (double *) calloc(n, sizeof(double));
		if(hmx[i-1] == NULL) printf(SPEICHER_F);

		for(j = 1; j <= n; j++) {
			hmx[i-1][j-1] = 1 / (double) (i + j - 1);
			}
		}

	return hmx;

}


// Funktion makeB
// liefert Vektor b (für rechte Seite) zurück
double *makeB(double **hmx, int n) {

	double *b;
	b = (double *) calloc(n, sizeof(double));
	if(b == NULL) printf(SPEICHER_F);

	int i, j;
	for(i = 1; i <= n; i++) {
		b[i-1] = 0.0;
		for(j = 1; j <= n; j++) {
			b[i-1] += hmx[i-1][j-1];
			}
		}

	return b;

}


// Funktion matrixProduct
// liefert Matrizenprodukt aus Matrizen ma und mb
double **matrixProduct(double **ma, double **mb, int rows, int cols) {

	// i Zeilenindex, j Spaltenindex, k Laufindex über Produkte von ma und mb
	int i, j, k;
	double **mc;
	mc = (double **) calloc(rows, sizeof(double *));
	// Zeile i von ma durchlaufen
	for(i = 1; i <= rows; i++) {
		mc[i-1] = (double *) calloc(cols, sizeof(double));
		// Spalte j von mb durchlaufen
		for(j = 1; j <= cols; j++) {
			mc[i-1][j-1] = 0.0;
			// Produkte aufsummieren
			for(k = 1; k <= rows; k++)
				mc[i-1][j-1] += ma[i-1][k-1] * mb[k-1][j-1];
			}
		}

	return mc;

}


// Funktion matrixProductVector
// liefert Produkt aus nxn-Matrix und n-Vektor (weil Vektoren vom Typ double * sind)
double *matrixProductVector(double **mx, double *ve, int n) {

	int i, j;
	double *ve2, help;
	ve2 = (double *) calloc(n, sizeof(double));
	for(i = 0; i < n; i++) {
		help = 0.0;
		for(j = 0; j < n; j++)
			help += mx[i][j] * ve[j];
		ve2[i] = help;
		}

	return ve2;

}


// Funktion swapRow
// vertauscht in Matrix mx Zeile i mit Zeile j (beginnend bei 1)
void swapRow(double **mx, int i, int j, int n) {

	if(i > n || j > n) {
		printf("swapRow: Zeile(n) außerhalb des gültigen Bereichs\n");
		return;
		}

	double *help;
	help = (double *) calloc(n, sizeof(double));
	if(help == NULL) printf(SPEICHER_F);

	help = mx[i-1];
	mx[i-1] = mx[j-1];
	mx[j-1] = help;

}


// Funktion columnMax
// gibt die Zeile aus dem Bereich [start, end] mit größtem Element in Spalte col aus
int columnMax(double **mx, int col, int start, int end) {

	int i, max = start;
	for(i = start+1; i <= end; i++)
		if(fabs(mx[max-1][col-1]) < fabs(mx[i-1][col-1]))
			max = i;

	return max;

}


// Funktion matrixInverse
// berechnet Matrixinverse von mx
double **matrixInverse(double **mx, int n) {

	int i, j, k, max;
	double **g, **mx_inverse, help;
	g = (double **) calloc(n, sizeof(double *));
	mx_inverse = (double **) calloc(n, sizeof(double *));
	if(g == NULL) printf(SPEICHER_F);
	if(mx_inverse == NULL) printf(SPEICHER_F);

	// g = mx|Id initialisieren
	for(i = 0; i < n; i++) {
		g[i] = (double *) calloc(2*n, sizeof(double));
		for(j = 0; j < 2*n; j++) {
			if(j < n)
				g[i][j] = mx[i][j];
			else
				g[i][j] = (i+n == j) ? 1.0 : 0.0;
			}
		}

	// g' = a|b mit a untere Dreiecksmatrix konstruieren
	for(j = 1; j <= n; j++) {
		max = columnMax(g, j, j, n);
		if(j != max)
			swapRow(g, j, max, 2*n);
		for(i = j+1; i <= n; i++) {
			help = g[i-1][j-1];
			for(k = 1; k <= 2*n; k++)
				g[i-1][k-1] = g[i-1][k-1] - (help / g[j-1][j-1]) * g[j-1][k-1];
			}
		}

	// g'' = c|d mit c Diagonalmatrix konstruieren
	for(j = 1; j <= n; j++) {
		for(i = j+1; i <= n; i++) {
			help = g[j-1][i-1];
			for(k = 1; k <= 2*n; k++)
				g[j-1][k-1] = g[j-1][k-1] - (help /g[i-1][i-1]) * g[i-1][k-1];
			}
		}

	// g''' = e|mx^(-1) konstruieren bzw. gleich in mx_inverse schreiben
	// (Diagonaleinträge aus c müssen noch skaliert werden)
	for(i = 1; i <= n; i++) {
		mx_inverse[i-1] = (double *) calloc(n, sizeof(double));
		help = g[i-1][i-1];
		for(j = 1; j <= n; j++)
			mx_inverse[i-1][j-1] = g[i-1][j-1+n] / help;
		}

	return mx_inverse;

}


// Funktion infinityNorm
// gibt die Unendlich-Norm einer Matrix aus
double infinityNorm(double **mx, int rows, int cols) {

	int k, j;
	double sum, max = 0.0;
	for(k = 1; k <= rows; k++) {
		sum = 0.0;
		for(j = 1; j <= cols; j++)
			sum += fabs(mx[k-1][j-1]);
		if(max < sum) max = sum;
		}

	return max;

}



int main(int argc, char **argv) {

	// n für die Dimension des Gleichungssystems; i, j, k Laufindexe;
	// max maximales Element in Spalte
	int n, i, j, k, max;

	n = 8;

	// hilbert ist Hilbertmatrix; hilbertInverse ist Inverse der Hilbertmatrix;
	// mx_R und mx_L sind die Matrizen aus der LR-Zerlegung;
	// mx_P ist die Permutationsmatrix für Pivotisierung,
	// ll ist jeweils l_ij für Zeilenumformung, kappa ist die Konditionszahl
	// b ist der Vektor der rechten Seite der Gleichung Hx = b,
	// z ist der Vektor der Gleichung Lz = b
	double **hilbert, **hilbertInverse, **mx_R, **mx_L, **mx_P, ll, kappa, *b, *z, *x, help;



	/*

	Algorithmus für Dimensionen 4, 8, 16, 32

	*/

	for(n = 4; n <= 32; n *= 2) {



	printf("\n\n_______________________\n");
	printf("|\n|  Durchlauf für n = %d\n|\n_______________________\n\n\n", n);



	// Speicherbelegungen und Matrix-Initialisierungen
	hilbert = makeHilbert(n);
	hilbertInverse = matrixInverse(hilbert, n);
	mx_R = makeHilbert(n);
	b = makeB(hilbert, n);

	z = (double *) calloc(n, sizeof(double));
	if(z == NULL) printf(SPEICHER_F);
	x = (double *) calloc(n, sizeof(double));
	if(x == NULL) printf(SPEICHER_F);

	mx_L = (double **) calloc(n, sizeof(double *));
	mx_P = (double **) calloc(n, sizeof(double *));

	for(i = 0; i < n; i++) {
		mx_L[i] = (double *) calloc(n, sizeof(double));
		mx_P[i] = (double *) calloc(n, sizeof(double));
		for(j = 0; j < n; j++) {
			mx_L[i][j] = 0.0;
			mx_P[i][j] = (i == j) ? 1.0 : 0.0;
			}
		}



	/************************
			LR-Zerlegung
	************************/


	for(j = 1; j <= n; j++) {
		// Pivot-Suche: Zeilennummer mit maximalem Element in Spalte j suchen und
		// diese Zeile mit Zeile j vertauschen
		max = columnMax(mx_R, j, j, n);
		if(max != j) {
			//Permutation in Permutationsmatrix eintragen und Zeile in Hilbertmatrix vertauschen
			swapRow(mx_P, j, max, n);
			swapRow(mx_R, j, max, n);
			}
		}

	b = matrixProductVector(matrixInverse(mx_P, n), b, n);

	// alle Spalten durchlaufen
	for(j = 1; j <= n; j++) {
		// ab Zeile i >= j Werte in R verändern und
		// Koeffizienten in Frobeniusmatrix L eintragen
		for(i = j; i <= n; i++) {
			ll = mx_R[i-1][j-1] / mx_R[j-1][j-1];
			mx_L[i-1][j-1] = ll;

			if(i > j) {
				// obere Dreiecksmatrix R konstruieren
				for(k = 1; k <= n; k++) {
					help = ll * mx_R[j-1][k-1];
					mx_R[i-1][k-1] = mx_R[i-1][k-1] - help;
					}
				}
			}
		}

	// Löse Lz = b durch Vorwärtssubstitution
	for(i = 1; i <= n; i++) {
		help = 0.0;
		for(j = 1; j < i; j++) {
			help += mx_L[i-1][j-1] * z[j-1];
			}
		z[i-1] = (b[i-1] - help) / mx_L[i-1][i-1];
		}

	// Löse Rx = z durch Rückwärtssubstitution
	for(i = n; i >= 1; i--) {
		help = 0.0;
		for(j = i+1; j <= n; j++)
			help += mx_R[i-1][j-1] * x[j-1];
		x[i-1] = (z[i-1] - help) / mx_R[i-1][i-1];
		}

	printf("Lösungsvektor x:\n");
	showVector(x, n);

	kappa = infinityNorm(hilbert, n, n) * infinityNorm(hilbertInverse, n, n);

	printf("Kondition: %e\n\n", kappa);



	}
	/*
	Ende Schleifendurchlauf für n
	*/


	return 0;

}