1 files changed, 63 insertions, 0 deletions
diff --git a/src/libmath/linalg.c b/src/libmath/linalg.c
new file mode 100644
index 0000000..8551ff1
--- /dev/null
+++ b/src/libmath/linalg.c
@@ -0,0 +1,63 @@
+#include <u.h>
+#include <libn.h>
+#include <libmath.h>
+#include <libmath/blas.h>
+
+// -----------------------------------------------------------------------
+// Vector
+
+void
+linalg·normalize(math·Vector vec)
+{
+    double norm;
+
+    norm = blas·normd(vec.len, vec.data, 1);
+    blas·scaled(vec.len, 1/norm, vec.data, 1);
+}
+// TODO: Write blas wrappers that eat vectors for convenience
+
+// -----------------------------------------------------------------------
+// Matrix
+//
+// NOTE: all matrices are row major oriented
+
+/*
+ * linalg·lq
+ * computes the LQ decomposition of matrix M: M = LQ
+ * L is lower triangular
+ * Q is orthogonal -> transp(Q) * Q = I
+ *
+ * m: matrix to factorize. changes in place
+ *     + lower triangle -> L
+ *     + upper triangle -> all reflection vectors stored in rows
+ * w: working buffer: len = ncols!
+ */
+error
+linalg·lq(math·Matrix m, math·Vector w)
+{
+    int i, j, len;
+    double *row, mag;
+    enum {
+        err·nil,
+        err·baddims,
+    };
+
+    if (m.dim[0] > m.dim[1]) {
+        return err·baddims;
+    }
+
+    for (i = 0; i < m.dim[0]; i++, m.data += m.dim[1]) {
+        row = m.data   + i;
+        len = m.dim[0] - i;
+
+        // TODO: Don't want to compute norm twice!!
+        w.data[0] = math·sgn(row[0]) * blas·normd(len, row, 1);
+        blas·axpyd(len, 1.0, row, 1, w.data, 1);
+        mag  = blas·normd(len, w.data, 1);
+        blas·scaled(len, 1/mag, w.data, 1);
+
+        blas·copyd(len - m.dim[0], w.data, 1, m.data + i, 1);
+    }
+
+    return err·nil;
+}