/*

2. Programmieraufgabe

Peter Hartmann, Matrikelnummer 840151


*/



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


// Funktion 'double_abs'
//	liefert Betrag von Zahl a
double double_abs(double a) {
	return (a < 0) ? -a : a;
	}


// 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 'makeA'
// erstellt Tridiagonalmatrix wie vorgegeben
double **makeA(int rows, int cols) {

	int i, j;
	double **mx_A = (double **) calloc(rows, sizeof(double *));
	for(i = 0; i < rows; i++) {
		mx_A[i] = (double *) calloc(cols, sizeof(double));
		for(j = 0; j < cols; j++) {
			if(i == j)
				mx_A[i][j] = 2.0;
			else
				mx_A[i][j] = (double_abs((double) i-j) == 1.0) ? 1.0 : 0.0;
			}
		}

	return mx_A;
}


// Funktion 'makeB'
// erstellt Vektor b = (3,4,...,4,3)T
double *makeB(int rows) {

	int i;
	double *b = (double *) calloc(rows, sizeof(double));

	for(i = 0; i < rows; i++)
		b[i] = (i == 0 || i == rows-1) ? 3.0 : 4.0;

	return b;
}


// Funktion 'makeOmega'
// erstellt die transponierte Rotationsmatrix mit Einträgen in Zeile k und Spalte l
double **makeOmega(int n, int k, int l, double c, double s) {

	int i, j;
	double **mx = (double **) calloc(n, sizeof(double *));

	for(i = 0; i < n; i++) {
		mx[i] = (double *) calloc(n, sizeof(double));
		for(j = 0; j < n; j++) {
			if(i == k && j == k)
				mx[i][j] = c;
			else if(i == k && j == l)
				mx[i][j] = -s;
			else if(i == l && j == k)
				mx[i][j] = s;
			else if(i == l && j == l)
				mx[i][j] = c;
			else if(i == j)
				mx[i][j] = 1.0;
			else
				mx[i][j] = 0.0;
		}
	}

	return mx;
}


// 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 mxn-Matrix und n-Vektor (weil Vektoren vom Typ double * sind)
double *matrixProductVector(double **mx, double *ve, int m, int n) {

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

	return ve2;
}




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

	int std_n = 9, i, j;
	int n = std_n;
	/* mx_A ist die Tridiagonalmatrix, mx_Omega_ij jeweils Rotationsmatrix
	   für Zeile i und Spalte j, b und x Vektoren des Gleichungssystems Ax = b,
	   tau, c und s sind Hilfswerte zur
	   Berechnung der Rotationsmatrix
	*/
	double **mx_A, **mx_R, **mx_Omega_ij, *b, *x, tau, c, s, help;

	if(argc > 1)
		n = atoi(argv[1]);

	printf("Programm wegen Wurzelfunktion 'sqrt' in Datei 'math.h' bitte mit ");
	printf("Option -lm kompilieren\n");
	printf("Benutzung: <Programm> [n]\n");
	printf("Standardwert: n = %d.\n", std_n);

	if(n == 0 || argc > 2)
		return 1;

	printf("\nDurchlauf mit Wert n = %d:\n", n);

	mx_A = makeA(n, n);
	mx_R = makeA(n, n);
	b = makeB(n);
	x = (double *) calloc(n, sizeof(double));


	// Givens-Rotationen

	// Spalten durchlaufen
	for(i = 0; i < n; i++) {
		// Zeilen in Spalte i rückwärts durchlaufen
		for(j = n-1; j > i; j--) {
			if(mx_R[j][i] != 0.0) {
				if(double_abs(mx_R[j-1][i]) > double_abs(mx_R[j][i])) {
					tau = mx_R[j-1][i] / mx_R[j][i];
					s = 1 / sqrt(1.0+tau*tau);
					c = s*tau;
					}
				else {
					tau = mx_R[j][i] / mx_R[j-1][i];
					c = 1 / sqrt(1.0+tau*tau);
					s = c*tau;
					}
				mx_Omega_ij = makeOmega(n, j, i, c, s);
				mx_R = matrixProduct(mx_Omega_ij, mx_R, n, n);
				b = matrixProductVector(mx_Omega_ij, b, n, n);
				}
		}
	}


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

	printf("\nLoesungsvektor x:\n");
	showVector(x, n);

	printf("Probe: Matrixprodukt Ax:\n");
	showVector(matrixProductVector(mx_A, x, n, n), n);

	return 0;

	}