diff options
Diffstat (limited to 'homeworks/hw-8.org')
| -rw-r--r-- | homeworks/hw-8.org | 288 |
1 files changed, 281 insertions, 7 deletions
diff --git a/homeworks/hw-8.org b/homeworks/hw-8.org index f4f4ebd..92a380f 100644 --- a/homeworks/hw-8.org +++ b/homeworks/hw-8.org @@ -4,11 +4,9 @@ #+LATEX: \setlength\parindent{0pt} #+OPTIONS: toc:nil -TODO: Update LIZFCM org file with jacobi solve - * Question One See ~UTEST(jacobi, solve_jacobi)~ in ~test/jacobi.t.c~ and the entry -~Jacobi -> solve_jacobi~ in the LIZFCM API documentation. +~Jacobi / Gauss-Siedel -> solve_jacobi~ in the LIZFCM API documentation. * Question Two We cannot just perform the Jacobi algorithm on a Leslie matrix since it is obviously not diagonally dominant - which is a requirement. It is @@ -24,14 +22,290 @@ direct solution method. It's simply the nature of the Jacobi algorithm being a convergent solution that determines its accuracy. LU factorization also performs in order $O(n^3)$ runtime for an $n \times n$ -matrix, whereas the Jacobi algorithm runs in order $O(k n^2) = O(n^2)$ but with the -con that $k$ is given by the convergence criteria, which might end up worse in -some cases, than LU. +matrix, whereas the Jacobi algorithm runs in order $O(k n^2) = O(n^2)$ on average +but with the con that $k$ is given by some function on both the convergence criteria and the number of +nonzero entries in the matrix - which might end up worse in some cases than the LU decomp approach. * Question Three See ~UTEST(jacobi, gauss_siedel_solve)~ in ~test/jacobi.t.c~ which runs the same unit test as ~UTEST(jacobi, solve_jacobi)~ but using the -~Jacobi -> gauss_siedel_solve~ method as documented in the LIZFCM API reference. +~Jacobi / Gauss-Siedel -> gauss_siedel_solve~ method as documented in the LIZFCM API reference. * Question Four, Five +We produce the following operation counts (by hackily adding the operation count as the last element +to the solution vector) and errors - the sum of each vector elements' absolute value away from 1.0 +using the proceeding patch and unit test. + +| N | JAC opr | JAC err | GS opr | GS err | LU opr | LU err | +| 5 | 1622 | 0.001244 | 577 | 0.000098 | 430 | 0.000000 | +| 6 | 2812 | 0.001205 | 775 | 0.000080 | 681 | 0.000000 | +| 7 | 5396 | 0.001187 | 860 | 0.000178 | 1015 | 0.000000 | +| 8 | 5618 | 0.001468 | 1255 | 0.000121 | 1444 | 0.000000 | +| 9 | 7534 | 0.001638 | 1754 | 0.000091 | 1980 | 0.000000 | +| 10 | 10342 | 0.001425 | 1847 | 0.000435 | 2635 | 0.000000 | +| 11 | 12870 | 0.001595 | 2185 | 0.000368 | 3421 | 0.000000 | +| 12 | 17511 | 0.001860 | 2912 | 0.000322 | 4350 | 0.000000 | +| 13 | 16226 | 0.001631 | 3362 | 0.000270 | 5434 | 0.000000 | +| 14 | 34333 | 0.001976 | 3844 | 0.000121 | 6685 | 0.000000 | +| 15 | 38474 | 0.001922 | 4358 | 0.000311 | 8115 | 0.000000 | +| 16 | 40405 | 0.002061 | 4904 | 0.000204 | 9736 | 0.000000 | +| 17 | 58518 | 0.002125 | 5482 | 0.000311 | 11560 | 0.000000 | +| 18 | 68079 | 0.002114 | 6092 | 0.000279 | 13599 | 0.000000 | +| 19 | 95802 | 0.002159 | 6734 | 0.000335 | 15865 | 0.000000 | +| 20 | 85696 | 0.002141 | 7408 | 0.000289 | 18370 | 0.000000 | +| 21 | 89026 | 0.002316 | 8114 | 0.000393 | 21126 | 0.000000 | +| 22 | 101537 | 0.002344 | 8852 | 0.000414 | 24145 | 0.000000 | +| 23 | 148040 | 0.002323 | 9622 | 0.000230 | 27439 | 0.000000 | +| 24 | 137605 | 0.002348 | 10424 | 0.000213 | 31020 | 0.000000 | +| 25 | 169374 | 0.002409 | 11258 | 0.000894 | 34900 | 0.000000 | +| 26 | 215166 | 0.002502 | 12124 | 0.000564 | 39091 | 0.000000 | +| 27 | 175476 | 0.002616 | 13022 | 0.000535 | 43605 | 0.000000 | +| 28 | 268454 | 0.002651 | 13952 | 0.000690 | 48454 | 0.000000 | +| 29 | 267034 | 0.002697 | 14914 | 0.000675 | 53650 | 0.000000 | +| 30 | 277193 | 0.002686 | 15908 | 0.000542 | 59205 | 0.000000 | +| 31 | 336792 | 0.002736 | 16934 | 0.000390 | 65131 | 0.000000 | +| 32 | 293958 | 0.002741 | 17992 | 0.000660 | 71440 | 0.000000 | +| 33 | 323638 | 0.002893 | 19082 | 0.001072 | 78144 | 0.000000 | +| 34 | 375104 | 0.003001 | 20204 | 0.001018 | 85255 | 0.000000 | +| 35 | 436092 | 0.003004 | 21358 | 0.000912 | 92785 | 0.000000 | +| 36 | 538143 | 0.003005 | 22544 | 0.000954 | 100746 | 0.000000 | +| 37 | 511886 | 0.003029 | 23762 | 0.000462 | 109150 | 0.000000 | +| 38 | 551332 | 0.003070 | 25012 | 0.000996 | 118009 | 0.000000 | +| 39 | 592750 | 0.003110 | 26294 | 0.000989 | 127335 | 0.000000 | +| 40 | 704208 | 0.003165 | 27608 | 0.000583 | 137140 | 0.000000 | + +#+BEGIN_SRC +diff --git a/src/matrix.c b/src/matrix.c +index 901a426..af5529f 100644 +--- a/src/matrix.c ++++ b/src/matrix.c +@@ -144,20 +144,54 @@ Array_double *solve_matrix_lu_bsubst(Matrix_double *m, Array_double *b) { + assert(b->size == m->rows); + assert(m->rows == m->cols); + ++ double opr = 0; ++ ++ opr += b->size; + Array_double *x = copy_vector(b); ++ ++ size_t n = m->rows; ++ opr += n * n; // (u copy) ++ opr += n * n; // l_empty ++ opr += n * n + n; // copy + put_identity_diagonal ++ opr += n; // pivot check ++ opr += m->cols; ++ for (size_t x = 0; x < m->cols; x++) { ++ opr += (m->rows - (x + 1)); ++ for (size_t y = x + 1; y < m->rows; y++) { ++ opr += 1; ++ opr += 2; // -factor ++ opr += 4 * n; // scale, add_v, free_vector ++ opr += 1; // -factor ++ } ++ } ++ opr += n; + Matrix_double **u_l = lu_decomp(m); ++ + Matrix_double *u = u_l[0]; + Matrix_double *l = u_l[1]; + ++ opr += n; ++ for (int64_t row = n - 1; row >= 0; row--) { ++ opr += 2 * (n - row); ++ opr += 1; ++ } + Array_double *b_fsub = fsubst(l, b); ++ ++ opr += n; ++ for (size_t x = 0; x < n; x++) { ++ opr += 2 * (x + 1); ++ opr += 1; // /= l->data[row]->data[row] ++ } + x = bsubst(u, b_fsub); +- free_vector(b_fsub); + ++ free_vector(b_fsub); + free_matrix(u); + free_matrix(l); + free(u_l); + +- return x; ++ Array_double *copy = add_element(x, opr); ++ free_vector(x); ++ return copy; + } + + Matrix_double *gaussian_elimination(Matrix_double *m) { +@@ -231,18 +265,36 @@ Array_double *jacobi_solve(Matrix_double *m, Array_double *b, + assert(b->size == m->cols); + size_t iter = max_iterations; + ++ double opr = 0; ++ ++ opr += 2 * b->size; // to initialize two vectors with the same dim of b twice + Array_double *x_k = InitArrayWithSize(double, b->size, 0.0); + Array_double *x_k_1 = + InitArrayWithSize(double, b->size, rand_from(0.1, 10.0)); + ++ // add since these wouldn't be accounter for after the loop ++ opr += 1; // iter decrement ++ opr += ++ 3 * x_k_1->size; // 1 to perform x_k_1, x_k and 2 to perform ||x_k_1||_2 + while ((--iter) > 0 && l2_distance(x_k_1, x_k) > l2_convergence_tolerance) { ++ opr += 1; // iter decrement ++ opr += ++ 3 * x_k_1->size; // 1 to perform x_k_1, x_k and 2 to perform ||x_k_1||_2 ++ ++ opr += m->rows; // row for add oprs + for (size_t i = 0; i < m->rows; i++) { + double delta = 0.0; ++ ++ opr += m->cols; + for (size_t j = 0; j < m->cols; j++) { + if (i == j) + continue; ++ ++ opr += 1; + delta += m->data[i]->data[j] * x_k->data[j]; + } ++ ++ opr += 2; + x_k_1->data[i] = (b->data[i] - delta) / m->data[i]->data[i]; + } + +@@ -251,8 +303,9 @@ Array_double *jacobi_solve(Matrix_double *m, Array_double *b, + x_k_1 = tmp; + } + +- free_vector(x_k); +- return x_k_1; ++ Array_double *copy = add_element(x_k_1, opr); ++ free_vector(x_k_1); ++ return copy; + } + + Array_double *gauss_siedel_solve(Matrix_double *m, Array_double *b, +@@ -262,30 +315,48 @@ Array_double *gauss_siedel_solve(Matrix_double *m, Array_double *b, + assert(b->size == m->cols); + size_t iter = max_iterations; + ++ double opr = 0; ++ ++ opr += 2 * b->size; // to initialize two vectors with the same dim of b twice + Array_double *x_k = InitArrayWithSize(double, b->size, 0.0); + Array_double *x_k_1 = + InitArrayWithSize(double, b->size, rand_from(0.1, 10.0)); + + while ((--iter) > 0) { ++ opr += 1; // iter decrement ++ ++ opr += x_k->size; // copy oprs + for (size_t i = 0; i < x_k->size; i++) + x_k->data[i] = x_k_1->data[i]; + ++ opr += m->rows; // row for add oprs + for (size_t i = 0; i < m->rows; i++) { + double delta = 0.0; ++ ++ opr += m->cols; + for (size_t j = 0; j < m->cols; j++) { + if (i == j) + continue; ++ ++ opr += 1; + delta += m->data[i]->data[j] * x_k_1->data[j]; + } ++ ++ opr += 2; + x_k_1->data[i] = (b->data[i] - delta) / m->data[i]->data[i]; + } + ++ opr += ++ 3 * x_k_1->size; // 1 to perform x_k_1, x_k and 2 to perform ||x_k_1||_2 + if (l2_distance(x_k_1, x_k) <= l2_convergence_tolerance) + break; + } + + free_vector(x_k); +- return x_k_1; ++ ++ Array_double *copy = add_element(x_k_1, opr); ++ free_vector(x_k_1); ++ return copy; + } +#+END_SRC + + +And this unit test: +#+BEGIN_SRC c +UTEST(hw_8, p4_5) { + printf("| N | JAC opr | JAC err | GS opr | GS err | LU opr | LU err | \n"); + + for (size_t i = 5; i < 100; i++) { + Matrix_double *m = generate_ddm(i); + double oprs[3] = {0.0, 0.0, 0.0}; + double errs[3] = {0.0, 0.0, 0.0}; + + Array_double *b_1 = InitArrayWithSize(double, m->rows, 1.0); + Array_double *b = m_dot_v(m, b_1); + double tolerance = 0.001; + size_t max_iter = 400; + + // JACOBI + { + Array_double *solution_with_opr_count = + jacobi_solve(m, b, tolerance, max_iter); + Array_double *solution = slice_element(solution_with_opr_count, + solution_with_opr_count->size - 1); + + for (size_t i = 0; i < solution->size; i++) + errs[0] += fabs(solution->data[i] - 1.0); + + oprs[0] = + solution_with_opr_count->data[solution_with_opr_count->size - 1]; + + free_vector(solution); + free_vector(solution_with_opr_count); + } + + // GAUSS-SIEDEL + { + Array_double *solution_with_opr_count = + gauss_siedel_solve(m, b, tolerance, max_iter); + Array_double *solution = slice_element(solution_with_opr_count, + solution_with_opr_count->size - 1); + + for (size_t i = 0; i < solution->size; i++) + errs[1] += fabs(solution->data[i] - 1.0); + + oprs[1] = + solution_with_opr_count->data[solution_with_opr_count->size - 1]; + + free_vector(solution); + free_vector(solution_with_opr_count); + } + + // LU-BSUBST + { + Array_double *solution_with_opr_count = solve_matrix_lu_bsubst(m, b); + Array_double *solution = slice_element(solution_with_opr_count, + solution_with_opr_count->size - 1); + + for (size_t i = 0; i < solution->size; i++) + errs[2] += fabs(solution->data[i] - 1.0); + + oprs[2] = + solution_with_opr_count->data[solution_with_opr_count->size - 1]; + + free_vector(solution); + free_vector(solution_with_opr_count); + } + free_matrix(m); + free_vector(b_1); + free_vector(b); + printf("| %zu | %f | %f | %f | %f | %f | %f | \n", i, oprs[0], errs[0], + oprs[1], errs[1], oprs[2], errs[2]); + } +} +#+END_SRC |