diff --git a/homeworks/hw-4.org b/homeworks/hw-4.org
new file mode 100644
index 0000000..ba82878
--- /dev/null
+++ b/homeworks/hw-4.org
@@ -0,0 +1,10 @@
+* My Basic Routines
+* Linear Solution Routines
+** Gaussian Elimination with Backwards Substitution
+** LU factorization
+* Matrix Routines
+** dotproduct
+** crossproduct
+** matrixvector
+** matrixmatrix
+
diff --git a/inc/lizfcm.h b/inc/lizfcm.h
index b390394..2dc9e42 100644
--- a/inc/lizfcm.h
+++ b/inc/lizfcm.h
@@ -25,6 +25,9 @@ extern double linf_distance(Array_double *v1, Array_double *v2);
 
 extern void format_vector_into(Array_double *v, char *s);
 
+extern void format_matrix_into(Matrix_double *m, char *s);
+extern void put_identity_diagonal(Matrix_double *m);
+
 extern Line *least_squares_lin_reg(Array_double *x, Array_double *y);
 
 #endif // LIZFCM_H
diff --git a/inc/macros.h b/inc/macros.h
index 841084d..f6e531e 100644
--- a/inc/macros.h
+++ b/inc/macros.h
@@ -11,15 +11,20 @@
     size_t size;                                                               \
   } Array_##TYPE
 
+#define DEFINE_MATRIX(TYPE)                                                    \
+  typedef struct {                                                             \
+    Array_##TYPE **data;                                                       \
+    size_t cols;                                                               \
+    size_t rows;                                                               \
+  } Matrix_##TYPE
+
 #define InitArray(TYPE, ...)                                                   \
   ({                                                                           \
     TYPE temp[] = __VA_ARGS__;                                                 \
     Array_##TYPE *arr = malloc(sizeof(Array_##TYPE));                          \
     arr->size = sizeof(temp) / sizeof(temp[0]);                                \
     arr->data = malloc(arr->size * sizeof(TYPE));                              \
-    if (arr->data) {                                                           \
-      memcpy(arr->data, temp, arr->size * sizeof(TYPE));                       \
-    }                                                                          \
+    memcpy(arr->data, temp, arr->size * sizeof(TYPE));                         \
     arr;                                                                       \
   })
 
@@ -28,13 +33,22 @@
     Array_##TYPE *arr = malloc(sizeof(Array_##TYPE));                          \
     arr->size = SIZE;                                                          \
     arr->data = malloc(arr->size * sizeof(TYPE));                              \
-    if (arr->data) {                                                           \
-      for (size_t i = 0; i < arr->size; i++)                                   \
-        arr->data[i] = INIT_VALUE;                                             \
-    }                                                                          \
+    for (size_t i = 0; i < arr->size; i++)                                     \
+      arr->data[i] = INIT_VALUE;                                               \
     arr;                                                                       \
   })
 
+#define InitMatrixWithSize(TYPE, ROWS, COLS, INIT_VALUE)                       \
+  ({                                                                           \
+    Matrix_##TYPE *matrix = malloc(sizeof(Matrix_##TYPE));                     \
+    matrix->rows = ROWS;                                                       \
+    matrix->cols = COLS;                                                       \
+    matrix->data = malloc(matrix->rows * sizeof(Array_##TYPE *));              \
+    for (size_t y = 0; y < matrix->rows; y++)                                  \
+      matrix->data[y] = InitArrayWithSize(TYPE, COLS, INIT_VALUE);             \
+    matrix;                                                                    \
+  })
+
 #define c_max(x, y) (((x) >= (y)) ? (x) : (y))
 #define c_min(x, y) (((x) <= (y)) ? (x) : (y))
 
diff --git a/inc/types.h b/inc/types.h
index 4af17db..6df04c2 100644
--- a/inc/types.h
+++ b/inc/types.h
@@ -10,6 +10,12 @@ DEFINE_ARRAY(int32_t);
 DEFINE_ARRAY(float);
 DEFINE_ARRAY(double);
 
+DEFINE_MATRIX(int);
+DEFINE_MATRIX(uint32_t);
+DEFINE_MATRIX(int32_t);
+DEFINE_MATRIX(float);
+DEFINE_MATRIX(double);
+
 typedef struct Line {
   double m;
   double a;
diff --git a/lib/lizfcm.a b/lib/lizfcm.a
index fe9a6ca..7a21e87 100644
Binary files a/lib/lizfcm.a and b/lib/lizfcm.a differ
diff --git a/src/matrix.c b/src/matrix.c
new file mode 100644
index 0000000..115e0a3
--- /dev/null
+++ b/src/matrix.c
@@ -0,0 +1,22 @@
+#include "lizfcm.h"
+#include <assert.h>
+#include <stdio.h>
+
+void put_identity_diagonal(Matrix_double *m) {
+  assert(m->rows == m->cols);
+
+  for (size_t y = 0; y < m->rows; ++y)
+    m->data[y]->data[y] = 1.0;
+}
+
+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];
+    format_vector_into(m->data[y], row_s);
+    sprintf(s, "%s %s \n", s, row_s);
+  }
+}
diff --git a/test/main.c b/test/main.c
index b0faf97..798bfbf 100644
--- a/test/main.c
+++ b/test/main.c
@@ -50,5 +50,11 @@ int main() {
   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);
+  format_matrix_into(m, m_s);
+  printf("%s", m_s);
+
   return 0;
 }