diff options
Diffstat (limited to 'matrix.c')
| -rw-r--r-- | matrix.c | 112 |
1 files changed, 42 insertions, 70 deletions
@@ -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"); |