Added macros to iterate over matrix

This avoids having to write double for-loops again and again.

- 'matrix_loop' just replaces the loops and gives access to the indices i
and j
- 'matrix_foreach' initializes a double pointer to the current element
- 'matrix_foreach_idx' is combination of the two above

All for-loops that simply iterate over a matrix have been replaced.

2 files changed, 49 insertions, 70 deletions

diff --git a/matrix.c b/matrix.c
index 9992bfb..d0c2cd4 100644
--- a/matrix.c
+++ b/matrix.c
@@ -37,12 +37,12 @@ void matrix_freen(Matrix **mats)
 
 void matrix_print(const Matrix *mat)
 {
-    for (size_t i = 0; i < mat->m; ++i)
+    matrix_foreach_idx(mat, it, i, j)
     {
-        for (size_t j = 0; j < mat->n; ++j)
-            printf("%lf ", matrix_at(mat, i, j));
-        printf("\n");
+        printf("%lf ", *it);
+        if (j == mat->n - 1 && i != mat->m - 1) printf("\n");
     }
+    printf("\n");
 }
 
 /* Expected format of a 3x2 matrix [3 2;1 2;0 0]
@@ -101,23 +101,14 @@ Matrix *matrix_id(size_t n)
 Matrix *matrix_const(size_t m, size_t n, double x)
 {
     Matrix *mat = matrix_alloc(m, n);
-    for (size_t i = 0; i < mat->m; ++i)
-    {
-        for (size_t j = 0; j < mat->n; ++j)
-            matrix_at(mat, i, j) = x;
-    }
+    matrix_foreach(mat, it) *it = x;
     return mat;
 }
 
 Matrix *matrix_copy (const Matrix *mat)
 {
     Matrix *tmp = matrix_alloc(mat->m, mat->n);
-
-    for (size_t i = 0; i < mat->m; ++i)
-    {
-        for (size_t j = 0; j < mat->n; ++j)
-            matrix_at(tmp, i, j) = matrix_at(mat, i, j);
-    }
+    matrix_foreach_idx(tmp, it, i, j) *it = matrix_at(mat, i, j);
     return tmp;
 }
 
@@ -183,25 +174,17 @@ Matrix *matrix_transpose(const Matrix *mat)
     assert(mat != NULL && "Matrix points to null");
     assert(mat->xs != NULL && "Matrix elements point to null");
 
-    Matrix *t = matrix_alloc(mat->n, mat->m);
-    for (size_t i = 0; i < mat->n; ++i)
-    {
-        for (size_t j = 0; j < mat->m; ++j)
-            matrix_at(t, i, j) = matrix_at(mat, j, i);
-    }
-    return t;
+    Matrix *T = matrix_alloc(mat->n, mat->m);
+    matrix_foreach_idx(T, t, i, j) *t = matrix_at(mat, j, i);
+    return T;
 }
 
 Matrix *matrix_add(const Matrix *A, const Matrix *B)
 {
     assert(A->m == B->m && A->n == B->n && "Dimension mismatch");
-    Matrix *C = matrix_alloc(A->m, A->n);
 
-    for (size_t i = 0; i < A->m; ++i)
-    {
-        for (size_t j = 0; j < A->n; ++j)
-            matrix_at(C, i, j) = matrix_at(A, i, j) + matrix_at(B, i, j);
-    }
+    Matrix *C = matrix_alloc(A->m, A->n);
+    matrix_foreach_idx(C, c, i, j) *c = matrix_at(A, i, j) + matrix_at(B, i, j);
 
     return C;
 }
@@ -209,19 +192,19 @@ Matrix *matrix_add(const Matrix *A, const Matrix *B)
 Matrix *matrix_scale(const Matrix *A, double x)
 {
     Matrix *tmp = matrix_alloc(A->m, A->n);
-
-    for (size_t i = 0; i < A->m; ++i)
-    {
-        for (size_t j = 0; j < A->n; ++j)
-            matrix_at(tmp, i, j) = x * matrix_at(A, i, j);
-    }
+    matrix_foreach_idx(tmp, it, i, j) *it = matrix_at(A, i, j) * x;
 
     return tmp;
 }
 
 Matrix *matrix_sub(const Matrix *A, const Matrix *B)
 {
-    return matrix_add(A, matrix_scale(B, -1));
+    assert(A->m == B->m && A->n == B->n && "Dimension mismatch");
+
+    Matrix *C = matrix_alloc(A->m, A->n);
+    matrix_foreach_idx(C, c, i, j) *c = matrix_at(A, i, j) - matrix_at(B, i, j);
+
+    return C;
 }
 
 Matrix *matrix_mult(const Matrix *A, const Matrix *B)
@@ -229,13 +212,9 @@ Matrix *matrix_mult(const Matrix *A, const Matrix *B)
     assert(A->n == B->m && "Dimension mismatch");
     Matrix *C = matrix_const(A->m, B->n, 0);
 
-    for (size_t i = 0; i < C->m; ++i)
-    {
-        for (size_t j = 0; j < C->n; ++j)
-        {
-            for (size_t k = 0; k < A->n; ++k)
-                matrix_at(C, i, j) += matrix_at(A, i, k) * matrix_at(B, k, j);
-        }
+    matrix_foreach_idx(C, c, i, j) {
+        for (size_t k = 0; k < A->n; ++k)
+            *c += matrix_at(A, i, k) * matrix_at(B, k, j);
     }
 
     return C;
@@ -252,31 +231,28 @@ Matrix *matrix_rand(size_t m, size_t n, int bound_l, int bound_u, MatrixType typ
 
     Matrix *mat = matrix_const(m, n, 0);
 
-    for (size_t i = 0; i < mat->m; ++i)
+    matrix_foreach_idx(mat, it, i, j)
     {
-        for (size_t j = 0; j < mat->n; ++j)
+        switch (type)
         {
-            switch (type)
-            {
-                case MATRIX_DIAG: 
-                    if (i != j) continue;
-                    break;
-                case MATRIX_TRIAG_UPPER:
-                    if (i > j) continue;
-                    break;
-                case MATRIX_TRIAG_LOWER:
-                    if (i < j) continue;
-                    break;
-                case MATRIX_TRIAG_SUPPER:
-                    if (i >= j) continue;
-                    break;
-                case MATRIX_TRIAG_SLOWER:
-                    if (i <= j) continue;
-                    break;
-                case MATRIX_NONE: break;
-            }
-            matrix_at(mat, i, j) = int_rand(bound_l, bound_u);
+            case MATRIX_DIAG:
+                if (i != j) continue;
+                break;
+            case MATRIX_TRIAG_UPPER:
+                if (i > j) continue;
+                break;
+            case MATRIX_TRIAG_LOWER:
+                if (i < j) continue;
+                break;
+            case MATRIX_TRIAG_SUPPER:
+                if (i >= j) continue;
+                break;
+            case MATRIX_TRIAG_SLOWER:
+                if (i <= j) continue;
+                break;
+            case MATRIX_NONE: break;
         }
+        *it = int_rand(bound_l, bound_u);
     }
 
     return mat;
@@ -290,11 +266,8 @@ int matrix_is_square(const Matrix *mat)
 int matrix_eq(const Matrix *A, const Matrix *B, double tol)
 {
     assert(A->m == B->m && A->n == B->n && "Dimension mismatch");
-    for (size_t i = 0; i < A->m; ++i)
-    {
-        for (size_t j = 0; j < A->n; ++j)
-            if (ABS(matrix_at(A, i, j) - matrix_at(B, i, j)) > tol) return FALSE;
-    }
+    matrix_foreach_idx(A, a, i, j)
+            if (ABS(*a - matrix_at(B, i, j)) > tol) return FALSE;
     return TRUE;
 }
 
@@ -319,7 +292,6 @@ Matrix **matrix_LR(const Matrix *A, const Matrix *b)
         {
             for (size_t j = 0; j < A->n; ++j)
                 matrix_at(A_curr, i, j) = matrix_at(A_prev, i, j);
-
         }
 
         assert(matrix_at(A_prev, k - 1, k - 1) != 0 && "Pivoting strategy needed");
diff --git a/matrix.h b/matrix.h
index 5ee5dc3..dc09ef6 100644
--- a/matrix.h
+++ b/matrix.h
@@ -4,6 +4,13 @@
 #include <stddef.h>
 
 #define matrix_at(mat, i, j) ((mat)->xs[(i) * (mat)->n + (j)])
+
+#define matrix_loop(mat, i, j) for (size_t i = 0; i < (mat)->m; i++) for (size_t j = 0; j < (mat)->m; j++)
+
+#define matrix_foreach(mat, it) for (double *it = (mat)->xs; it < (mat)->xs + ((mat)->m * (mat)->n); ++it)
+
+#define matrix_foreach_idx(mat, it, i, j) double *it = (mat)->xs; for (size_t i = 0; i < (mat)->m; i++) for (size_t j = 0; j < (mat)->n; j++, ++it)
+
 #define TRUE 1
 #define FALSE 0