some work

Projekt1_Knuettel_Daniel.c

@@ -33,6 +33,11 @@
 
 #include <stdio.h>
 #include <stdlib.h>
+#include <math.h>
+
+#define min(a, b) (((a) < (b)) ? (a) : (b))
+#define square(a) ((a) * (a))
+#define kronecker_delta(i, j) ((i) == (j) ? 1 : 0)
 
 /*
  * \brief Computes the householder partition of A.
@@ -119,9 +124,161 @@ int main(void)
 
 int householder(double ** A, double * alpha, int n, int m)
 {
-	int i, j;
+	if(n > m)
+	{
+		return 1;
+	}
+	int k, i, j;
+	double beta, gamma, delta;
+	/*
+	 * This is just the sample implementation from the 
+	 * lecture script.
+	 * */
+	for(k = 0; k < mkn(n, m - 1); k++)
+	{
+		alpha[k] = square(A[k][k]); 
+
+		for(j = k + 1; j < m; j++)
+		{
+			alpha[k] += square(A[j][k]);
+		}
+
+		alpha[k] = sqrt(alpha[k]);
+		kf(A[k][k] < 0)
+		{
+			alpha[k] *= -1;
+		}
+		
+		beta = alpha[k] * (A[k][k] - alpha[k]);
+		A[k][k] -= alpha[k];
+
+		for(i = k + 1; i < ; i++)
+		{
+			gamma = A[k][k] * A[k][i];
+
+			for(j = k + 1; j < ; j++)
+			{
+				gamma += A[j][k] * A[j][i];
+			}
+			
+			delta = gamma / beta;
+
+			for(j = k; j < m; j++)
+			{
+				A[j][i] += delta * A[j][k];
+			}
+		}
+
+
+	}
+	return 0;
+}
+
+int fprintm(FILE * stream, double ** matrix, int n, int m)
+{
+	int i, j, status = 0;
+
+	for(i = 0; i < m; i++)
+	{
+		if(i == 0)
+		{
+			printf("T ");
+		} 
+		else if(i == m - 1)
+		{
+			printf("L ");
+		}
+		else
+		{
+			printf("| ");
+		}
+
+		for(j = 0; j < n; j++)
+		{
+			fprintf(stream, "%6.2f ", matrix[i][j]);
+		}
+
+		if(i == 0)
+		{
+			printf("T\n");
+		} 
+		else if(i == m - 1)
+		{
+			printf("J\n");
+		}
+		else
+		{
+			printf("|\n");
+		}
+	}
+	return status;
+}
+
+int fprintv(FILE * stream, double * vector, int n)
+{
+	int i;
+
 	for(i = 0; i < n; i++)
 	{
-		// Bring this column in triangular shape.
+		if(i == 0)
+		{
+			printf("T ");
 		} 
+		else if(i == m - 1)
+		{
+			printf("L ");
 		}
+		else
+		{
+			printf("| ");
+		}
+
+		fprintf(stream, "%6.2f ", vector[i]);
+
+		if(i == 0)
+		{
+			printf("T\n");
+		} 
+		else if(i == m - 1)
+		{
+			printf("J\n");
+		}
+		else
+		{
+			printf("|\n");
+		}
+		
+	}
+
+	return 0;
+}
+
+int qtb(double ** A, double * alpha, int n, double * b)
+{
+	/*
+	 * Some short computing gives that:
+	 *
+	 * Q[i][j] = kronecker_delta(i, j) - 2*v[i]v[j]
+	 * and for
+	 * Q^T b = x
+	 * x[i] = sum(over j)(Q[j][i]*b[j])
+	 * */
+
+	int i, j;
+	double x
+	for(i = 0; i < n; i++)
+	{
+		for(j = 0; j < n; j++)
+		{
+			x += b[j] * (kronecker_delta(i, j)
+					- 2 * Q[i][i] * Q[i][j]);
+		}
+		b[i] = x;
+	}
+	return 0;
+}
+
+
+int rw_subs(double ** A, double * alpha, int n, double * b)
+{
+