Added Exercise 1

Projekt1_Knuettel_Daniel.c

@@ -51,7 +51,7 @@
  * \param m Length of A.
  * \return Status.
  * */
-int householder(double ** A, double * alpha, int n, int m);
+int householder(double ** A, double * alpha, int m, int n);
 
 /*
  * \brief Computes the solution of and n times n system using a QR partition.
@@ -79,7 +79,7 @@ int qtb(double ** A, double * alpha, int n, double * b);
 
 /*
  *
- * \brief Uses backwards substitution to solve R = Q^Tb.
+ * \brief Uses backwards substitution to solve R = Q^Tb. b := Q^Tb has already been calculated.
  *
  * \param A Matrix that has been modified by the function householder.
  * \param alpha Diagonal elements of R, produced by the function householder.
@@ -97,11 +97,11 @@ int rw_subs(double ** A, double * alpha, int n, double * b);
  *
  * \param stream The output stream.
  * \param matrix The matrix to print.
- * \param n Length of matrix[i].
+ * \param n Length of A[i].
- * \param m Length of matrix.
+ * \param m Length of A.
  *
  * */
-int fprintm(FILE * stream, double ** matrix, int n, int m);
+int fprintm(FILE * stream, double ** matrix, int m, int n);
 /*
  * \brief Print the vector in a nicely formatted way to stream.
  *
@@ -112,17 +112,113 @@ int fprintm(FILE * stream, double ** matrix, int n, int m);
  * */
 int fprintv(FILE * stream, double * vector, int n);
 
+/*
+ * \brief Extract R from A and alpha, allocates R dynamically.
+ *
+ * \param A Marix computed by householder.
+ * \param alpha Vector computed by householder.
+ * \param m Length of A.
+ * \param n Length of A[i].
+ * \param R Pointer to the result (unallocated).
+ * */
+int unwind_R(double ** A, double * alpha, int m, int n, double *** R);
+/*
+ * \brief Extract v from A and alpha, allocates R dynamically.
+ *
+ * \param A Marix computed by householder.
+ * \param alpha Vector computed by householder.
+ * \param m Length of A.
+ * \param n Length of A[i].
+ * \param v Pointer to the results (unallocated).
+ * */
+int unwind_v(double ** A, int m, int n, double *** v);
+
 #define printm(matrix, n, m) fprintm(stdout, matrix, n, m)
 #define printv(vector, n) fprintv(stdout, vector, n)
 
 int main(void)
 {
 	int status = 0;
+	int i, j;
+	int m, n;
+
+	// PART a: QR-Partition.
+	
+	double A_stack_1[15] = {2, 2.23606797749979, 0
+				, 2, 3, -1
+				, 2, 3, -1
+				, 1, 0, -1
+				, 0, 2.6457513110645907, 5.669467095138408
+				};
+	m = 5;
+	n = 3;
+	double ** A_1 = malloc(sizeof(double *) * m);
+	for(i = 0; i < m; i++)
+	{
+		A_1[i] = malloc(sizeof(double) * n);
+		for(j = 0; j < n; j++)
+		{
+			A_1[i][j] = A_stack_1[n*i + j];
+		}
+	}
+
+	double * alpha_1 = malloc(sizeof(double) * m);
+
+	status = householder(A_1, alpha_1, m, n);
+	if(status)
+	{
+		fprintf(stderr, "failed to compute QR decomposition\n");
+		goto exit;
+	}
+
+	double ** R_1, ** vs_1;
+
+	status = unwind_R(A_1, alpha_1, m, n, &R_1);
+	if(status)
+	{
+		fprintf(stderr, "failed to unwind R\n");
+		goto exit;
+	}
+
+	status = unwind_v(A_1, m, n, &vs_1);
+	if(status)
+	{
+		fprintf(stderr, "failed to unwind v\n");
+		for(i = 0; i < m; i++)
+		{
+			free(R_1[i]);
+		}
+		free(R_1);
+		goto exit;
+	}
+
+	printm(R_1, m, n);
+	for(i = 0; i < n; i++)
+	{
+		printv(vs_1[i], m);
+		free(vs_1[i]);
+	}
+
+	free(vs_1);
+	for(i = 0; i < m; i++)
+	{
+		free(R_1[i]);
+	}
+	free(R_1);
+
+
+	
+
+exit:
+	if(status)
+	{
+		fprintf(stderr, "something went wrong\n");
+	}
 	return status;
 }
 
 
-int householder(double ** A, double * alpha, int n, int m)
+int householder(double ** A, double * alpha, int m, int n)
 {
 	if(n > m)
 	{
@@ -134,7 +230,7 @@ int householder(double ** A, double * alpha, int n, int m)
 	 * This is just the sample implementation from the 
 	 * lecture script.
 	 * */
-	for(k = 0; k < mkn(n, m - 1); k++)
+	for(k = 0; k < min(n, m - 1); k++)
 	{
 		alpha[k] = square(A[k][k]); 
 
@@ -144,7 +240,7 @@ int householder(double ** A, double * alpha, int n, int m)
 		}
 
 		alpha[k] = sqrt(alpha[k]);
-		kf(A[k][k] < 0)
+		if(A[k][k] < 0)
 		{
 			alpha[k] *= -1;
 		}
@@ -152,11 +248,11 @@ int householder(double ** A, double * alpha, int n, int m)
 		beta = alpha[k] * (A[k][k] - alpha[k]);
 		A[k][k] -= alpha[k];
 
-		for(i = k + 1; i < ; i++)
+		for(i = k + 1; i < n ; i++)
 		{
 			gamma = A[k][k] * A[k][i];
 
-			for(j = k + 1; j < ; j++)
+			for(j = k + 1; j < m; j++)
 			{
 				gamma += A[j][k] * A[j][i];
 			}
@@ -174,7 +270,7 @@ int householder(double ** A, double * alpha, int n, int m)
 	return 0;
 }
 
-int fprintm(FILE * stream, double ** matrix, int n, int m)
+int fprintm(FILE * stream, double ** matrix, int m, int n)
 {
 	int i, j, status = 0;
 
@@ -182,15 +278,15 @@ int fprintm(FILE * stream, double ** matrix, int n, int m)
 	{
 		if(i == 0)
 		{
-			printf("T ");
+			fprintf(stream, "T ");
 		} 
 		else if(i == m - 1)
 		{
-			printf("L ");
+			fprintf(stream, "L ");
 		}
 		else
 		{
-			printf("| ");
+			fprintf(stream, "| ");
 		}
 
 		for(j = 0; j < n; j++)
@@ -200,17 +296,18 @@ int fprintm(FILE * stream, double ** matrix, int n, int m)
 
 		if(i == 0)
 		{
-			printf("T\n");
+			fprintf(stream, "T\n");
 		} 
 		else if(i == m - 1)
 		{
-			printf("J\n");
+			fprintf(stream, "J\n");
 		}
 		else
 		{
-			printf("|\n");
+			fprintf(stream, "|\n");
 		}
 	}
+	fflush(stream);
 	return status;
 }
 
@@ -222,34 +319,35 @@ int fprintv(FILE * stream, double * vector, int n)
 	{
 		if(i == 0)
 		{
-			printf("T ");
+			fprintf(stream, "T ");
 		} 
-		else if(i == m - 1)
+		else if(i == n - 1)
 		{
-			printf("L ");
+			fprintf(stream, "L ");
 		}
 		else
 		{
-			printf("| ");
+			fprintf(stream, "| ");
 		}
 
 		fprintf(stream, "%6.2f ", vector[i]);
 
 		if(i == 0)
 		{
-			printf("T\n");
+			fprintf(stream, "T\n");
 		} 
-		else if(i == m - 1)
+		else if(i == n - 1)
 		{
-			printf("J\n");
+			fprintf(stream, "J\n");
 		}
 		else
 		{
-			printf("|\n");
+			fprintf(stream, "|\n");
 		}
 		
 	}
 
+	fflush(stream);
 	return 0;
 }
 
@@ -265,13 +363,13 @@ int qtb(double ** A, double * alpha, int n, double * b)
 	 * */
 
 	int i, j;
-	double x
+	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]);
+					- 2 * A[i][i] * A[i][j]);
 		}
 		b[i] = x;
 	}
@@ -282,3 +380,156 @@ int qtb(double ** A, double * alpha, int n, double * b)
 int rw_subs(double ** A, double * alpha, int n, double * b)
 {
 
+	/*
+	 * In the lecture section 1 about the gaussian algorithm 
+	 * we have given the following formula: 
+	 *
+	 * x[j] = 1/A[j][j] * (b[j] - sum(from j + 1 over k to n)(a[j][k]*x[k]))
+	 *
+	 * So loop from k = n to 1 to avoid recursion and 
+	 * substitute the content of b with the solution.
+	 * */
+	int j, k;
+	double tmp;
+	// Note: indices start at 0.
+	for(j = n - 1; j >= 0; j--)
+	{
+		// Note: we start from k + 1, so there are no
+		// diagonal elements of A contained.
+		for(k = j + 1; k < n; k++)
+		{
+			tmp += b[k] * A[j][k];
+		}
+
+		// A[j][j] = alpha[j].
+		b[j] = (b[j] - tmp) / alpha[j];
+
+	}
+	return 0;
+}
+
+int solve_QR(double ** A, double * alpha, int n, double * b)
+{
+	if(qtb(A, alpha, n, b))
+	{
+		return 1;
+	}
+	if(rw_subs(A, alpha, n, b))
+	{
+		return 1;
+	}
+	return 0;
+}
+
+int unwind_R(double ** A, double * alpha, int m, int n, double *** R)
+{
+	double ** result = malloc(sizeof(double *) * m);
+	int i, j;
+	int have_to_free_all = 0;
+	if(!result)
+	{
+		return 1;
+	}
+
+	// Allocate the memory.
+	for(i = 0; i < m; i++)
+	{
+		result[i] = malloc(sizeof(double) * n);
+		if(!result[i])
+		{
+			have_to_free_all = 1;
+			break;
+		}
+	}
+
+	// Out of Memory.
+	if(have_to_free_all)
+	{
+		for(j = 0; j < i; j++)
+		{
+			free(result[i]);
+		}
+		free(result);
+		return 1;
+	}
+
+	// Alloc is OK.
+	// Calculate the result.
+	for(i = 0; i < m; i++)
+	{
+		for(j = 0; j < n; j++)
+		{
+			if(i > j)
+			{
+				result[i][j] = 0;
+			}
+			if(i == j)
+			{
+				result[i][j] = alpha[j];
+			}
+			if(i < j)
+			{
+				result[i][j] = A[i][j];
+			}
+		}
+	}
+	*R = result;
+	return 0;
+
+}
+
+
+int unwind_v(double ** A, int m, int n, double *** v)
+{
+	// We have n vectors.
+	double ** result = malloc(sizeof(double *) * n);
+	int i, j;
+	int have_to_free_all = 0;
+	if(!result)
+	{
+		return 1;
+	}
+
+	// Allocate the memory.
+	for(i = 0; i < n; i++)
+	{
+		// They are m long.
+		result[i] = malloc(sizeof(double) * m);
+		if(!result[i])
+		{
+			have_to_free_all = 1;
+			break;
+		}
+	}
+
+	// Out of Memory.
+	if(have_to_free_all)
+	{
+		for(j = 0; j < i; j++)
+		{
+			free(result[i]);
+		}
+		free(result);
+		return 1;
+	}
+
+	// Alloc is OK.
+	// Calculate the result.
+	for(i = 0; i < m; i++)
+	{
+		for(j = 0; j < n; j++)
+		{
+			if(i > j)
+			{
+				result[j][i] = 0;
+			}
+			else
+			{
+				result[j][i] = A[j][i];
+			}
+		}
+	}
+	*v = result;
+	return 0;
+
+}