diff --git a/Projekt1_Knuettel_Daniel.c b/Projekt1_Knuettel_Daniel.c index 0b23bce..750c1cf 100644 --- a/Projekt1_Knuettel_Daniel.c +++ b/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 m Length of matrix. + * \param n Length of A[i]. + * \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; + +}