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
106
107
108
|
#pragma src "/usr/doc/src/linalg"
#pragma lib "liblinalg.a"
typedef float real;
extern const float meps;
/* Cubes, matrics, and Vector data types */
/* Morton ordering */
#define evenbits 0xAAAAAAAA
#define oddbits 0x55555555
//#define threebits 0x49249249
#define oinc2(i) (((i) - oddbits) & oddbits)
#define einc2(i) (((i) - evenbits) & evenbits)
#define odec2(i) (((i) - 1) & oddbits)
#define edec2(i) (((i) - 1) & evenbits)
#define tr(i, j) (((i) << 1) | ((j) >> 1))
ushort undilate2(uint);
uint dilate2(ushort);
/* Iterators */
#define miter(m, i) for((i) = 0; (i) < (m)->len; (i)++) if(mbndok((m), (i)))
#define viter(v, i) for((i) = 0; (i) < (v)->len; (i)++)
typedef struct block block;
typedef struct cube cube;
typedef struct matrix matrix;
typedef struct vector vector;
#define C(x, i, j, k) ((x)->data[(k) * (x)->stride + (dilate2(i) | (dilate2(j) << 1))])
#define CV(x, i, k) ((x)->data[(k) * (x)->stride + (i)])
#define M(x, i, j) ((x)->data[dilate2(i) | (dilate2(j) << 1)])
#define V(x, i) ((x)->data[(i)])
#define mbndok(a, i) (((i) & oddbits) < (a)->mr && ((i) & evenbits) < (a)->mc)
struct block{
real *data;
ulong sz;
uchar allocated:1;
};
struct vector {
block;
uint len;
};
struct matrix{
vector;
ushort nr, nc;
uint mr, mc;
uint levels, offset;
};
struct cube{
matrix;
ushort ns;
uint stride;
};
cube *ncube(cube *, real *, ushort, ushort, ushort);
matrix *nmatrix(matrix *, real *, ushort, ushort);
vector *nvector(vector *, real *, uint);
void fblock(block *);
block *zblock(block *);
block *cblock(block *, block *);
matrix *qdot(matrix *, matrix *, matrix *, int);
matrix *mdot(matrix *, matrix *, matrix *, int);
matrix *madd(matrix *, matrix *, matrix *);
matrix *smprod(real, matrix *, matrix *);
real dot(vector *, vector *);
vector *sprod(real, vector *, vector *);
vector *vadd(vector *, vector *, vector *);
vector *vprod(vector *, vector *, vector *);
vector *vsub(vector *, vector *, vector *);
vector *vdot(matrix *, vector *, vector *, int);
void vnorm(uint, ...);
void vsort(vector *);
uint *ivsort(vector *);
extern uint MAXITERS;
matrix *eigen(matrix *, matrix *, vector *);
matrix *jacobi(matrix *, matrix *, vector *);
matrix *tridiag(matrix *, matrix *, matrix *);
matrix *nystrom(matrix *, matrix *, vector *);
matrix *qr(matrix *, matrix *, matrix *);
matrix *gs(matrix *);
real entropy(matrix *);
matrix *mahalanobis(matrix *, matrix *);
matrix *msqrt(matrix *, matrix *);
matrix *minv(matrix *, matrix *, real);
vector *solve(matrix *, vector *, vector *);
void svd(matrix *, matrix *, vector *, matrix *);
matrix *col(matrix *, matrix *, uint);
matrix *row(matrix *, matrix *, uint);
matrix *slice(cube *, matrix *, uint);
vector *drop(matrix *, vector *);
matrix *raise(vector *, matrix *);
matrix *trans(matrix *, matrix *);
matrix *diag(vector *, matrix *);
vector *msum(matrix *, vector *, uint);
vector *getdiag(matrix *, vector *);
|
