#ifndef MATRIX_H
#define MATRIX_H

#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 matrix_is_colvec(mat) ((mat)->n == 1)

#define matrix_is_rowvec(mat) ((mat)->m == 1)

#define matrix_is_vec(mat) (matrix_is_colvec(mat) || matrix_is_rowvec(mat))

#define TRUE 1
#define FALSE 0

typedef enum {
    MATRIX_DIAG,
    MATRIX_TRIAG_UPPER,
    MATRIX_TRIAG_LOWER,
    MATRIX_TRIAG_SUPPER, /* S for strict */
    MATRIX_TRIAG_SLOWER,
    MATRIX_NONE
} MatrixType;

typedef enum {
    MATRIX_SWAP_ROWS,
    MATRIX_SWAP_COLS
} MatrixSwapType;

typedef struct {
    double *xs;
    size_t m;
    size_t n;
} Matrix;

Matrix *matrix_alloc(size_t m, size_t n);
void matrix_free(Matrix *mat);
void matrix_freen(Matrix **mats);
void matrix_print(const Matrix *mat);
Matrix *matrix_from_str(char *str);
Matrix *matrix_from_arr(double arr[], size_t m, size_t n);
char *matrix_to_str(const Matrix *mat);
Matrix *matrix_id(size_t n);
Matrix *matrix_const(size_t m, size_t n, double x);
Matrix *matrix_copy (const Matrix *mat);
Matrix *matrix_swap(const Matrix *mat, MatrixSwapType t, size_t i, size_t j);
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);
Matrix *matrix_add(const Matrix *A, const Matrix *B);
Matrix *matrix_scale(const Matrix *A, double x);
Matrix *matrix_sub(const Matrix *A, const Matrix *B);
Matrix *matrix_mult(const Matrix *A, const Matrix *B);
Matrix *matrix_rand(size_t m, size_t n, int bound_l, int bound_u, MatrixType type);
int matrix_is_square(const Matrix *mat);
int matrix_eq(const Matrix *A, const Matrix *B, double tol);
double matrix_norm_frob(const Matrix *mat);
Matrix **matrix_LR(const Matrix *A, const Matrix *b);
Matrix *matrix_forwardel(const Matrix *L, const Matrix *b);
Matrix *matrix_backsubst(const Matrix *R, const Matrix *y);

#endif