#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;
}