1 files changed, 108 insertions, 0 deletions
diff --git a/linalg.h b/linalg.h
new file mode 100644
index 0000000..5613100
--- /dev/null
+++ b/linalg.h
@@ -0,0 +1,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 *);