lizfcm/homeworks/hw-8.org

Homework 8

Question One

See UTEST(jacobi, solve_jacobi) in test/jacobi.t.c and the entry 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 certainly not always the case, but, if a Leslie matrix $L$ is invertible, we can first perform gaussian elimination on $L$ augmented with $n_{k+1}$ to obtain $n_k$ with the Jacobi method. See UTEST(jacobi, leslie_solve) in test/jacobi.t.c for an example wherein this method is tested on a Leslie matrix to recompute a given initial population distribution.

In terms of accuracy, an LU factorization and back substitution approach will always be as correct as possible within the limits of computation; it's a 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)$ 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 -> 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

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;
 }

And this unit test:

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]);
  }
}