diff options
author | Nicholas Noll <nbnoll@eml.cc> | 2020-05-08 21:33:24 -0700 |
---|---|---|
committer | Nicholas Noll <nbnoll@eml.cc> | 2020-05-08 21:33:24 -0700 |
commit | 04c688f125069b65517b00660c31c81e210ddf3a (patch) | |
tree | 5a92bd41181ae1c8c586f7da701c4c1115bd5dd5 /sys/libmath/linalg.c | |
parent | 327ca20a2a89d2408b53ff7854982560304cb76c (diff) |
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.
Diffstat (limited to 'sys/libmath/linalg.c')
-rw-r--r-- | sys/libmath/linalg.c | 57 |
1 files changed, 57 insertions, 0 deletions
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; +} |