1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
|
#ifndef MATRIX_H
#define MATRIX_H
#include <stddef.h>
#include <stdbool.h>
#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_foreach(mat, it) \
matrix_foreach_idx(mat, it, _dummy1, _dummy2)
#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)
#define matrix_is_rowvec(mat) ((mat)->m == 1)
#define matrix_is_vec(mat) (matrix_is_colvec(mat) || matrix_is_rowvec(mat))
#define matrix_is_square(mat) ((mat)->m == (mat)->n)
typedef struct Vector Vector;
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 Matrix {
double *xs;
size_t m;
size_t n;
bool T;
} Matrix;
Matrix *matrix_alloc(size_t m, size_t n);
void matrix_free(Matrix *mat);
void matrix_free_many(size_t argcount, ...);
void matrix_free_many1(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);
void matrix_const1(Matrix *mat, double x);
Matrix *matrix_copy (const Matrix *mat);
void matrix_copy1(Matrix *A, const Matrix *B);
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);
void matrix_transpose1(Matrix *mat);
Matrix *matrix_add(const Matrix *A, const Matrix *B);
void matrix_add1(Matrix *A, const Matrix *B);
Matrix *matrix_scale(const Matrix *A, double x);
void matrix_scale1(Matrix *A, double x);
Matrix *matrix_sub(const Matrix *A, const Matrix *B);
void matrix_sub1(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);
void matrix_rand1(Matrix *mat, int bound_l, int bound_u, MatrixType type);
bool matrix_eq(const Matrix *A, const Matrix *B, double tol);
double matrix_norm_frob(const Matrix *mat);
Matrix **matrix_LR(const Matrix *A, const Vector *b);
Vector *matrix_forwardel(const Matrix *L, const Vector *b);
Vector *matrix_backsubst(const Matrix *R, const Vector *y);
#endif
|
