diff --git a/include/matio.h b/include/matio.h
new file mode 100644
index 0000000..fb908df
--- /dev/null
+++ b/include/matio.h
@@ -0,0 +1,14 @@
+#include <gsl/gsl_matrix.h>
+
+#ifndef MATIO_H
+#define MATIO_H
+
+void
+print_matrix
+( gsl_matrix *matrix );
+
+void
+print_flatten
+( gsl_matrix *matrix );
+
+#endif
diff --git a/src/matio.c b/src/matio.c
new file mode 100644
index 0000000..1d76fdd
--- /dev/null
+++ b/src/matio.c
@@ -0,0 +1,38 @@
+#include <gsl/gsl_matrix.h>
+
+#include "matio.h"
+
+void print_matrix ( gsl_matrix *matrix ) {
+  printf("\n");
+  for ( int i = 0 ; i < matrix->size1 ; i++ ) {
+    for ( int j = 0 ; j < matrix->size2 ; j++ ) {
+      printf(" %.6e ", gsl_matrix_get(matrix, i, j));
+    }
+    printf("\n");
+  }
+  printf("\n");
+}
+
+void print_flatten ( gsl_matrix *matrix ) {
+  printf("\n");
+  if ( matrix->size1 == 1 ) {
+    for ( int i = 0 ; i < matrix->size2 ; i++ ) {
+      if ( i == matrix->size2-1 ) {
+        printf("%g\n", gsl_matrix_get(matrix, 0, i));
+        return;
+      }
+      printf("%g, ", gsl_matrix_get(matrix, 0, i));
+    }
+    printf("\n");
+    return;
+  }
+
+  for ( int i = 0 ; i < matrix->size1 ; i++ ) {
+    if ( i == matrix->size1-1 ) {
+      printf("%g\n", gsl_matrix_get(matrix, i, 0));
+      return;
+    }
+    printf("%g, ", gsl_matrix_get(matrix, i, 0));
+  }
+  printf("\n");
+}