diff options
Diffstat (limited to 'sys/libmath/linalg.c')
-rw-r--r-- | sys/libmath/linalg.c | 63 |
1 files changed, 0 insertions, 63 deletions
diff --git a/sys/libmath/linalg.c b/sys/libmath/linalg.c deleted file mode 100644 index 8551ff1..0000000 --- a/sys/libmath/linalg.c +++ /dev/null @@ -1,63 +0,0 @@ -#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; -} |