lu factorization

5 files changed, 151 insertions, 52 deletions

diff --git a/inc/lizfcm.h b/inc/lizfcm.h
index 2dc9e42..99acbf2 100644
--- a/inc/lizfcm.h
+++ b/inc/lizfcm.h
@@ -12,22 +12,26 @@ extern double forward_derivative_at(double (*f)(double), double a, double h);
 extern double backward_derivative_at(double (*f)(double), double a, double h);
 
 extern double sum_v(Array_double *v);
+extern Array_double *scale_v(Array_double *v1, double m);
 extern Array_double *add_v(Array_double *v1, Array_double *v2);
 extern Array_double *minus_v(Array_double *v1, Array_double *v2);
 extern double dot_v(Array_double *v1, Array_double *v2);
 extern double l2_norm(Array_double *v);
 extern double l1_norm(Array_double *v);
 extern double linf_norm(Array_double *v);
-
 extern double l2_distance(Array_double *v1, Array_double *v2);
 extern double l1_distance(Array_double *v1, Array_double *v2);
 extern double linf_distance(Array_double *v1, Array_double *v2);
-
+extern Array_double *copy_vector(Array_double *v1);
+extern void free_vector(Array_double *v);
 extern void format_vector_into(Array_double *v, char *s);
+extern Line *least_squares_lin_reg(Array_double *x, Array_double *y);
 
-extern void format_matrix_into(Matrix_double *m, char *s);
 extern void put_identity_diagonal(Matrix_double *m);
+extern Matrix_double *copy_matrix(Matrix_double *m);
+extern void free_matrix(Matrix_double *m);
+extern Matrix_double **put_lu_decomp(Matrix_double *m);
 
-extern Line *least_squares_lin_reg(Array_double *x, Array_double *y);
+extern void format_matrix_into(Matrix_double *m, char *s);
 
 #endif // LIZFCM_H
diff --git a/lib/lizfcm.a b/lib/lizfcm.a
index 7a21e87..969c894 100644
--- a/lib/lizfcm.a
+++ b/lib/lizfcm.a
diff --git a/src/matrix.c b/src/matrix.c
index 115e0a3..a23aad4 100644
--- a/src/matrix.c
+++ b/src/matrix.c
@@ -9,14 +9,74 @@ void put_identity_diagonal(Matrix_double *m) {
     m->data[y]->data[y] = 1.0;
 }
 
+Matrix_double *copy_matrix(Matrix_double *m) {
+  Matrix_double *copy = InitMatrixWithSize(double, m->rows, m->cols, 0.0);
+  for (size_t y = 0; y < copy->rows; y++) {
+    free_vector(copy->data[y]);
+    copy->data[y] = copy_vector(m->data[y]);
+  }
+  return copy;
+}
+
+Matrix_double **put_lu_decomp(Matrix_double *m) {
+  assert(m->cols == m->rows);
+
+  Matrix_double *u = copy_matrix(m);
+  Matrix_double *l = InitMatrixWithSize(double, m->rows, m->cols, 0.0);
+  put_identity_diagonal(l);
+
+  Matrix_double **u_l = malloc(sizeof(Matrix_double *) * 2);
+
+  for (size_t y = 0; y < m->rows; y++) {
+    if (u->data[y]->data[y] == 0) {
+      printf("ERROR: a pivot is zero in given matrix\n");
+      exit(-1);
+    }
+  }
+
+  if (u && l) {
+    for (size_t x = 0; x < m->cols; x++) {
+      for (size_t y = x + 1; y < m->rows; y++) {
+        double denom = u->data[x]->data[x];
+
+        if (denom == 0) {
+          printf("ERROR: non-factorable matrix\n");
+          exit(-1);
+        }
+
+        double factor = -(u->data[y]->data[x] / denom);
+
+        Array_double *scaled = scale_v(u->data[x], factor);
+        Array_double *added = add_v(scaled, u->data[y]);
+        free_vector(scaled);
+        free_vector(u->data[y]);
+
+        u->data[y] = added;
+        l->data[y]->data[x] = -factor;
+      }
+    }
+  }
+
+  u_l[0] = u;
+  u_l[1] = l;
+  return u_l;
+}
+
+void free_matrix(Matrix_double *m) {
+  for (size_t y = 0; y < m->rows; ++y)
+    free_vector(m->data[y]);
+  free(m);
+}
+
 void format_matrix_into(Matrix_double *m, char *s) {
   sprintf(s, "");
   if (m->rows == 0)
     sprintf(s, "empty");
 
   for (size_t y = 0; y < m->rows; ++y) {
-    char row_s[256];
+    char *row_s = malloc(sizeof(char) * 256);
     format_vector_into(m->data[y], row_s);
     sprintf(s, "%s %s \n", s, row_s);
+    free(row_s);
   }
 }
diff --git a/src/vector.c b/src/vector.c
index 61692e1..e6e2544 100644
--- a/src/vector.c
+++ b/src/vector.c
@@ -41,12 +41,19 @@ double sum_v(Array_double *v) {
   return sum;
 }
 
+Array_double *scale_v(Array_double *v, double m) {
+  Array_double *copy = copy_vector(v);
+  for (size_t i = 0; i < v->size; i++)
+    copy->data[i] *= m;
+  return copy;
+}
+
 Array_double *add_v(Array_double *v1, Array_double *v2) {
   assert(v1->size == v2->size);
 
-  Array_double *sum = InitArrayWithSize(double, v1->size, 0);
+  Array_double *sum = copy_vector(v1);
   for (size_t i = 0; i < v1->size; i++)
-    sum->data[i] = v1->data[i] + v2->data[i];
+    sum->data[i] += v2->data[i];
   return sum;
 }
 
@@ -80,6 +87,18 @@ double linf_distance(Array_double *v1, Array_double *v2) {
   return dist;
 }
 
+Array_double *copy_vector(Array_double *v) {
+  Array_double *copy = InitArrayWithSize(double, v->size, 0.0);
+  for (size_t i = 0; i < copy->size; ++i)
+    copy->data[i] = v->data[i];
+  return copy;
+}
+
+void free_vector(Array_double *v) {
+  free(v->data);
+  free(v);
+}
+
 void format_vector_into(Array_double *v, char *s) {
   sprintf(s, "");
   if (v->size == 0)
diff --git a/test/main.c b/test/main.c
index 798bfbf..3551af2 100644
--- a/test/main.c
+++ b/test/main.c
@@ -5,56 +5,72 @@
 double f(double x) { return (x - 1) / (x + 1); }
 
 int main() {
-  printf("smaceps(): %.10e\n", smaceps());
-  printf("dmaceps(): %.10e\n", dmaceps());
-
-  Array_double *v = InitArray(double, {3, 1, -4, 1, 5, -9, 3});
-  Array_double *w = InitArray(double, {-2, 7, 1, -8, -2, 8, 5});
-
-  char v_s[256];
-  char w_s[256];
-  format_vector_into(v, v_s);
-  format_vector_into(w, w_s);
-
-  printf("v: %s\n", v_s);
-  printf("w: %s\n", w_s);
-  printf("l1_norm(v): %f\n", l1_norm(v));
-  printf("l2_norm(v): %f\n", l2_norm(v));
-  printf("linf_norm(v): %f\n", linf_norm(v));
-
-  printf("l1_dist(v, w): %f\n", l1_distance(v, w));
-  printf("l2_dist(v, w): %f\n", l2_distance(v, w));
-  printf("linf_dist(v, w): %f\n", linf_distance(v, w));
-
-  double h = 0.001;
-  printf("approx f'(1) w/ c.d.: %f\n", central_derivative_at(&f, 1, h));
-  printf("approx f'(1) w/ fw.d.: %f\n", forward_derivative_at(&f, 1, h));
-  printf("approx f'(1) w/ bw.d.: %f\n", backward_derivative_at(&f, 1, h));
-
-  v = InitArray(double, {1, 2, 3, 4, 5});
-  w = InitArray(double, {2, 3, 4, 5, 6});
-  format_vector_into(v, v_s);
-  format_vector_into(w, w_s);
-  printf("v: %s\n", v_s);
-  printf("w: %s\n", w_s);
-
-  Line *line = least_squares_lin_reg(v, w);
-  printf("least_squares_lin_reg(v, w): (%f)x + %f\n", line->m, line->a);
-
-  v = InitArray(double, {1, 2, 3, 4, 5, 6, 7});
-  w = InitArray(double, {0.5, 3, 2, 3.5, 5, 6, 7.5});
-  format_vector_into(v, v_s);
-  format_vector_into(w, w_s);
-  printf("v: %s\n", v_s);
-  printf("w: %s\n", w_s);
-  line = least_squares_lin_reg(v, w);
-  printf("least_squares_lin_reg(v, w): (%f)x + %f\n", line->m, line->a);
+  //  printf("smaceps(): %.10e\n", smaceps());
+  //  printf("dmaceps(): %.10e\n", dmaceps());
+  //
+  //  Array_double *v = InitArray(double, {3, 1, -4, 1, 5, -9, 3});
+  //  Array_double *w = InitArray(double, {-2, 7, 1, -8, -2, 8, 5});
+  //
+  //  char v_s[256];
+  //  char w_s[256];
+  //  format_vector_into(v, v_s);
+  //  format_vector_into(w, w_s);
+  //
+  //  printf("v: %s\n", v_s);
+  //  printf("w: %s\n", w_s);
+  //  printf("l1_norm(v): %f\n", l1_norm(v));
+  //  printf("l2_norm(v): %f\n", l2_norm(v));
+  //  printf("linf_norm(v): %f\n", linf_norm(v));
+  //
+  //  printf("l1_dist(v, w): %f\n", l1_distance(v, w));
+  //  printf("l2_dist(v, w): %f\n", l2_distance(v, w));
+  //  printf("linf_dist(v, w): %f\n", linf_distance(v, w));
+  //
+  //  double h = 0.001;
+  //  printf("approx f'(1) w/ c.d.: %f\n", central_derivative_at(&f, 1, h));
+  //  printf("approx f'(1) w/ fw.d.: %f\n", forward_derivative_at(&f, 1, h));
+  //  printf("approx f'(1) w/ bw.d.: %f\n", backward_derivative_at(&f, 1, h));
+  //
+  //  v = InitArray(double, {1, 2, 3, 4, 5});
+  //  w = InitArray(double, {2, 3, 4, 5, 6});
+  //  format_vector_into(v, v_s);
+  //  format_vector_into(w, w_s);
+  //  printf("v: %s\n", v_s);
+  //  printf("w: %s\n", w_s);
+  //
+  //  Line *line = least_squares_lin_reg(v, w);
+  //  printf("least_squares_lin_reg(v, w): (%f)x + %f\n", line->m, line->a);
+  //
+  //  v = InitArray(double, {1, 2, 3, 4, 5, 6, 7});
+  //  w = InitArray(double, {0.5, 3, 2, 3.5, 5, 6, 7.5});
+  //
+  //  format_vector_into(v, v_s);
+  //  format_vector_into(w, w_s);
+  //  printf("v: %s\n", v_s);
+  //  printf("w: %s\n", w_s);
+  //  line = least_squares_lin_reg(v, w);
+  //  printf("least_squares_lin_reg(v, w): (%f)x + %f\n", line->m, line->a);
 
   char m_s[256];
   Matrix_double *m = InitMatrixWithSize(double, 8, 8, 0.0);
   put_identity_diagonal(m);
+  for (int i = 0; i < 8; i++) {
+    for (int j = 0; j < 8; j++) {
+      m->data[i]->data[j] = (i + 1.0) + j * 3 + (rand() % 12);
+    }
+  }
+
   format_matrix_into(m, m_s);
-  printf("%s", m_s);
+  printf("%s\n", m_s);
+
+  Matrix_double **u_l = put_lu_decomp(m);
+  Matrix_double *u = u_l[0];
+  Matrix_double *l = u_l[1];
+
+  format_matrix_into(u, m_s);
+  printf("%s\n", m_s);
+  format_matrix_into(l, m_s);
+  printf("%s\n", m_s);
 
   return 0;
 }