Matrix transpose doesn't edit memory representation

- transpose swaps dimensions
- transpose sets variable to indicate it has been transposed
- matrix_at acts accordingly
- matrix_foreach_idx now uses matrix_at to set the iterator to allow
  iterating over a transposed matrix

2 files changed, 24 insertions, 9 deletions

diff --git a/matrix.c b/matrix.c
index 646be42..222ad6d 100644
--- a/matrix.c
+++ b/matrix.c
@@ -20,6 +20,7 @@ Matrix *matrix_alloc(size_t m, size_t n)
     mat->xs = xs;
     mat->m = m;
     mat->n = n;
+    mat->T = false;
     return mat;
 }
 
@@ -196,14 +197,19 @@ double matrix_trace(const Matrix *mat)
 
 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);
-    matrix_foreach_idx(T, t, i, j) *t = matrix_at(mat, j, i);
+    Matrix *T = matrix_copy(mat);
+    matrix_transpose1(T);
     return T;
 }
 
+void matrix_transpose1(Matrix *mat)
+{
+    size_t m = mat->m;
+    mat->m = mat->n;
+    mat->n = m;
+    mat->T = !mat->T;
+}
+
 Matrix *matrix_add(const Matrix *A, const Matrix *B)
 {
     assert(A->m == B->m && A->n == B->n && "Dimension mismatch");
diff --git a/matrix.h b/matrix.h
index 7501092..f0603fd 100644
--- a/matrix.h
+++ b/matrix.h
@@ -4,13 +4,20 @@
 #include <stddef.h>
 #include <stdbool.h>
 
-#define matrix_at(mat, i, j) ((mat)->xs[(i) * (mat)->n + (j)])
+#define matrix_at(mat, i, j) \
+    (mat)->xs[(mat)->T ? (j) * (mat)->m + (i) : (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)->n; j++)
+#define matrix_loop(mat, i, j) \
+    for (size_t i = 0; i < (mat)->m; ++i) \
+        for (size_t j = 0; j < (mat)->n; ++j)
 
-#define matrix_foreach(mat, it) for (double *it = (mat)->xs; it < (mat)->xs + ((mat)->m * (mat)->n); ++it)
+#define matrix_foreach(mat, it) \
+    matrix_foreach_idx(mat, it, _dummy1, _dummy2)
 
-#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 matrix_foreach_idx(mat, it, i, j) \
+    for (size_t i = 0, _d = 1; _d; --_d) \
+        for (double _dd = 0, *it = &_dd; i < (mat)->m; ++i) \
+            for (size_t j = 0; j < (mat)->n && (it = &matrix_at((mat), i, j), 1); ++j)
 
 #define matrix_is_colvec(mat) ((mat)->n == 1)
 
@@ -40,6 +47,7 @@ typedef struct Matrix {
     double *xs;
     size_t m;
     size_t n;
+    bool T;
 } Matrix;
 
 Matrix *matrix_alloc(size_t m, size_t n);
@@ -68,6 +76,7 @@ void matrix_swap1(Matrix *mat, MatrixSwapType t, size_t i, size_t j);
 double matrix_trace(const Matrix *mat);
 
 Matrix *matrix_transpose(const Matrix *mat);
+void matrix_transpose1(Matrix *mat);
 
 Matrix *matrix_add(const Matrix *A, const Matrix *B);
 void matrix_add1(Matrix *A, const Matrix *B);