diff --git a/Projekt1_Knuettel_Daniel.c b/Projekt1_Knuettel_Daniel.c index ac34271..0b23bce 100644 --- a/Projekt1_Knuettel_Daniel.c +++ b/Projekt1_Knuettel_Daniel.c @@ -33,6 +33,11 @@ #include #include +#include + +#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) +{ +