From 04c688f125069b65517b00660c31c81e210ddf3a Mon Sep 17 00:00:00 2001
From: Nicholas Noll <nbnoll@eml.cc>
Date: Fri, 8 May 2020 21:33:24 -0700
Subject: Adding strided computation to blas kernels.

I started implementing LQ factorization and immediately realized I needed strided views.
For simplicity, I will just implement them in the most portable, C native way (no vectorization).
Speed can come later.
---
 sys/libmath/linalg.c | 57 ++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 57 insertions(+)

(limited to 'sys/libmath/linalg.c')

diff --git a/sys/libmath/linalg.c b/sys/libmath/linalg.c
index 57f799b..5a73527 100644
--- a/sys/libmath/linalg.c
+++ b/sys/libmath/linalg.c
@@ -2,4 +2,61 @@
 #include <libn.h>
 #include <libmath.h>
 
+// -----------------------------------------------------------------------
+// Vector
 
+void
+linalg·normalize(math·Vector vec)
+{
+    double norm;
+
+    norm = blas·norm(vec.len, vec.data);
+    blas·scale(vec.len, 1/norm, vec.data);
+}
+// 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·norm(len, row);
+        blas·axpy(len, 1.0, row, w.data);
+        mag  = blas·norm(len, w.data);
+        blas·scale(len, 1/mag, w.data);
+
+        blas·copy(len - m.dim[0], w.data, m.data + i);
+    }
+
+    return err·nil;
+}
-- 
cgit v1.2.1