/*

3. 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 'norm'
// liefert die 1-Norm eines Vektors
double norm(double *v, int n) {

	int i;
	double v_norm = 0.0;
	for(i = 0; i < n; i++)
		v_norm += double_abs(v[i]);

	return v_norm;
}


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


// Funktion 'makeA'
// erstellt Tridiagonalmatrix A
double **makeA(int n) {

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

	return mx_A;
}


// Funktion 'makeC'
// erstellt Vektor c = (1/h^2)*(-a,0,...,0,-b)T
double *makeC(int n, double a, double b) {

	int i;
	double *c, h;

	c = (double *) calloc(n, sizeof(double));

	h = 1.0/((double)n + 1.0);

	c[0] = (1.0/(h*h))*(-a);
	c[n-1] = (1.0/(h*h))*(-b);

	for(i = 1; i < n-1; i++)
		c[i] = 0.0;

	return c;
}


// Funktion 'makeOmega'
// erstellt die transponierte Rotationsmatrix mit Einträgen in Zeile k und Spalte l
// (wird benötigt zum Lösen des Gleichungssystem dF(y)*delta_y = -F(y))
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 plus_mx
// liefert komponentenweise Addition zweier Matrizen
double **plus_mx(double **mxA, double **mxB, int n) {

	int i, j;

	double **mxC;

	mxC = (double **) calloc(n, sizeof(double *));

	for(i = 0; i < n; i++) {
		mxC[i] = (double *) calloc(n, sizeof(double));
		for(j = 0; j < n; j++)
			mxC[i][j] = mxA[i][j] + mxB[i][j];
	}

	return mxC;
}

// Funktion minus
// liefert Vektor -u zum Vektor u
double *minus(double *u, int n) {

	int i;
	double *v;

	v = (double *) calloc(n, sizeof(double));

	for(i = 0; i < n; i++)
		v[i] = -u[i];

	return v;
}


// Funktion plus
// liefert komponentenweise Addition zweier Vektoren
double *plus(double *u, double *v, int n) {

	int i;
	double *w;

	w = (double *) calloc(n, sizeof(double));

	for(i = 0; i < n; i++)
		w[i] = u[i] + v[i];

	return w;
}


// Funktion f_A
// liefert Produkt Ay
double *f_A(double **A, double *y, int n) {

	int i, j;
	double *Ay, help;
	Ay = (double *) calloc(n, sizeof(double));
	for(i = 0; i < n; i++) {
		help = 0.0;
		for(j = 0; j < n; j++)
			help += A[i][j] * y[j];
		Ay[i] = help;
		}

	return Ay;
}


// Funktion f_G
// liefert Funktion G(y) (hier: G(y) = exp(y))
double *f_G(double *y, int n) {

	int i;
	double *Gy;

	Gy = (double *) calloc(n, sizeof(double));

	for(i = 0; i < n; i++)
		Gy[i] = exp(y[i]);

	return Gy;
}


// Funktion f_F
// liefert Funktion F(y) = Ay + G(y) + c
double *f_F(double **A, double *y, double *c, int n) {

	return plus(plus(f_A(A, y, n), f_G(y, n), n), c, n);
}


// Funktion f_dF
// liefert Ableitung der Funktion F(y); dF(y) = A + dG(y)
double **f_dF(double **A, double *y, int n) {

	int i, j;
	double **dG;

	dG = (double **) calloc(n, sizeof(double *));

	// Ableitung dG(y)
	for(i = 0; i < n; i++) {
		dG[i] = (double *) calloc(n, sizeof(double));
		for(j = 0; j < n; j++)
			dG[i][j] = (i == j) ? exp(y[i]) : 0.0;
	}

	return plus_mx(A, dG, n);
}



// Funktion linearSolve
// liefert Lösungsvektor x des Gleichungssystems Ax = b (A nxn-Matrix mit Rang n)
// mit Hilfe von Givens-Rotationen aus Programmieraufgabe 2
double *linearSolve(double **A, double *b, int n) {

	int i, j;
	double **mx_R, **mx_Omega_ij, *x, tau, s, c, help;

	mx_R = (double **) calloc(n, sizeof(double *));
	x = (double *) calloc(n, sizeof(double));

	// Matrix A kopieren
	for(i = 0; i < n; i++) {
		mx_R[i] = (double *) calloc(n, sizeof(double));
		for(j = 0; j < n; j++)
			mx_R[i][j] = A[i][j];
	}

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


	return x;
}


// Funktion 'newtonIteration'
// liefert Element delta_y der Gleichung dF(y)*delta_y = -F(y)
double *newtonIteration(double **dFy, double *Fy, int n) {

	return linearSolve(dFy, minus(Fy, n), n);
}


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

	int std_n = 8, i;
	int n = std_n;

	double **mxA, **dF, *delta_y, *y, *c, a = 0.0, b = 1.0;


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

			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);


	mxA = makeA(n);
	c = makeC(n, a, b);
	dF = f_dF(mxA, y, n);

	// Startwert: y = 0
	y = (double *) calloc(n, sizeof(double));
	for(i = 0; i < n; i++)
		y[i] = 0.0;


	// Newton-Iteration wiederholen, bis Norm der Newton-Korrektur kleiner als
	// 10^-6 mal die Norm der neuen Näherungslösung
	do {
		delta_y = newtonIteration(dF, f_F(mxA, y, c, n), n);
		y = plus(y, delta_y, n);
	} while(norm(delta_y, n) >= 10e-6*norm(y, n));

	printf("Naeherungsloesung: y =\n");
	showVector(y, n);
	printf("Probe: F(y) = \n");
	showVector(f_F(mxA, y, c, n), n);

	return 0;

}