added QR and LU decomposition

build.sh

@@ -1,2 +1,2 @@
-gcc tests/test.c  src/lina.c -o test  -Wall -Wextra -g -Isrc/
+gcc tests/test.c  src/lina.c -o test  -Wall -Wextra -g -Isrc/ -lm
-gcc tests/test_loader.c src/lina.c -o test_loader -Wall -Wextra -g -Isrc/
+gcc tests/test_loader.c src/lina.c -o test_loader -Wall -Wextra -g -Isrc/ -lm

src/lina.c

@@ -5,6 +5,7 @@
 #include <stdio.h>
 #include <errno.h>
 #include <ctype.h>
+#include <math.h>
 #include "lina.h"
 
 /* Function: lina_dot
@@ -677,3 +678,152 @@ void lina_conv(double *A, double *B, double *C,
                         C[j * Cw + i] += A[(i - Bw/2 + u) * Aw + (i - Bh/2 + v)] * B[v * Bw + u];
             }
 }
+
+void lina_decompLU(double *A, double *L, double *U, int n)
+{
+    assert(n > 0);
+    assert(A != L && A != U && L != U);
+
+    // TODO: Handle the case when A can not be
+    //       decomposed.
+
+    memset(L, 0, sizeof(double) * n * n);
+    memset(U, 0, sizeof(double) * n * n);
+
+    /*
+    // Zero-out the lower half of L and the upper
+    // half of U.
+    for (int i = 0; i < n; i++)
+        for (int j = i+1; j < n; j++) 
+            {
+                L[j * n + i] = 0;
+                U[i * n + j] = 0;
+            }
+    */
+
+    for (int i = 0; i < n; i++)
+        {
+            for (int k = i; k < n; k++)
+                {
+                    int sum = 0; // L[i,j] * U[j,k]
+                    for (int j = 0; j < i; j++)
+                        sum += L[i * n + j] * U[j * n + k];
+
+                    U[i * n + k] = A[i * n + k] - sum;
+                }
+
+            for (int k = i; k < n; k++)
+                {
+                    if (i == k)
+                        L[i * n + i] = 1;
+                    else 
+                        {
+                            int sum = 0;
+                            for (int j = 0; j < i; j++)
+                                sum += L[k * n + j] * U[j * n + i];
+
+                            L[k * n + i] = (A[k * n + i] - sum) / U[i * n + i];
+                        }
+                }
+        }
+}
+
+bool lina_det(double *A, int n, double *det)
+{
+    // Allocate the space for the L,U matrices.
+    // I can't think of a version of this algorithm
+    // where a temporary buffer isn't necessary.
+    double *T = malloc(sizeof(double) * n * n * 2);
+    if (T == NULL)
+        return false;
+
+    // Do the decomposition
+    double *L = T;
+    double *U = T + (n * n);
+    lina_decompLU(A, L, U, n);
+
+    // Knowing that
+    //
+    //   A = LU
+    //
+    // then
+    //
+    //   det(A) = det(LU) = det(L)det(U)
+    //
+    // Since L and U are triangular, their 
+    // determinant is the product of their 
+    // diagonals, so the product of the 
+    // determinants is the product of both 
+    // the diagonals.
+
+    double prod = 1;
+    for (int i = 0; i < n; i++)
+        prod *= L[i * n + i] * U[i * n + i];
+    
+    if (det)
+        *det = prod;
+
+    free(T);
+    return true;
+}
+
+/* Checks that [A] is kind of upper triangular.
+**
+*/
+static bool isUpperTriangularEnough(double *A, int n, double eps)
+{
+    assert(A != NULL && n > 0 && epd > 0);
+
+    // Check that the lower triangular portion (without
+    // considering the diagonal) is zero.
+    for (int i = 0; i < n; i++)
+        for (int j = 0; j < i-1; i++)
+            if (A[i * n + j] > eps)
+                return false;
+
+    // Now check that the diagonal is also zero. Though
+    // since we are using the real version of the QR
+    // algorithm, only real eigenvalues can be found.
+    // Any comples eigenvalues will manifest as 2x2 blocks
+    // on the diagonal, so we need to allow such blocks.
+    // To do this, a non-zero block is allowed if it's
+    // not following another non-zero block.
+    //
+    // An important thing to note is that 2x2 matrices
+    // will always be considered upper triangular by this
+    // function, so the caller must manage this case.
+    bool flag = false;
+    for (int i = 0; i < n-1; i++) {
+        if (fabs(A[(i + 1) * n + i]) > eps) {
+            if (flag)
+                return false;
+            flag = true;
+        } else
+            flag = false;
+    }
+    return true;
+}
+
+bool lina_eig(double *M, double *E, int n)
+{
+    double *A = malloc(sizeof(double) * n * n * 3);
+    if (A == NULL)
+        return false;
+    memcpy(A, M, sizeof(double) * n * n);
+
+    double *Q = A + n * n;
+    double *R = Q + n * n;
+
+    do {
+        for (int i = 0; i < 100; i++) {
+            lina_decompQR(A, Q, R, n); // A(n) = QR
+            lina_dot(R, Q, A, n, n, n); // A(n+1) = RQ
+        }
+    } while (!isUpperTriangularEnough(A, n, 0.1));
+
+    for (int i = 0; i < n; i++)
+        E[i] = A[i * n + i];
+
+    free(A);
+    return true;
+}

src/lina.h

@@ -1,7 +1,10 @@
+#include <stdbool.h>
 
 /* ---- Operations ---- */
 void   lina_dot(double *A, double *B, double *C, int m, int n, int l);
 void   lina_add(double *A, double *B, double *C, int m, int n);
+double lina_mul(double *v, int n);
+bool   lina_det(double *A, int n, double *det);
 void   lina_scale(double *A, double *B, double k, int m, int n);
 void   lina_transpose(double *A, double *B, int m, int n);
 void   lina_conv(double *A, double *B, double *C, 
@@ -11,6 +14,12 @@ void lina_dot2(double *A, double *B, double *C,
                int As, int Bs, int Cs, 
                int m, int n, int l);
 
+bool lina_eig(double *M, double *E, int n);
+void lina_decompLU(double *A, double *L, double *U, int n);
+void lina_decompQR(double *A, double *Q, double *R, int n);
+void lina_orthoNormGramSchmidt(double *A, double *Q, int n);
+
+
 /* ---- Utilities ---- */
 double *lina_loadMatrixFromStream(FILE *fp, int *width, int *height, char **error);
 int     lina_saveMatrixToStream(FILE *fp, double *A, int width, int height, char **error);

src/qr.c

@@ -0,0 +1,142 @@
+#include <math.h>
+#include <assert.h>
+
+typedef struct {
+    double *items;
+    int size;
+} square_matrix_t;
+
+typedef struct {
+    double *items;
+    int stride;
+    int size;
+} vector_t;
+
+typedef struct {
+    vector_t base;
+    double scale;
+} scaled_vector_t;
+
+static square_matrix_t 
+square_matrix_from_raw(double *M, int n)
+{
+    return (square_matrix_t) {.items=M, .size=n};
+}
+
+static vector_t 
+get_column_of_square_matrix(square_matrix_t M, int i)
+{
+    assert(i > -1 && i < M.size);
+
+    return (vector_t) { 
+        .items  = M.items + i, 
+        .stride = M.size, 
+        .size   = M.size
+    };
+}
+
+static void 
+copy_vector(vector_t V, vector_t S)
+{
+    assert(V.size == S.size);
+    for (int i = 0; i < V.size; i++)
+        V.items[V.stride * i] = S.items[S.stride * i];
+}
+
+static void 
+subtract_vector_inplace(vector_t V, scaled_vector_t S)
+{
+    assert(V.size == S.base.size);
+
+    for (int i = 0; i < V.size; i++)
+        V.items[V.stride * i] -= S.scale * S.base.items[S.base.stride * i];
+}
+
+static void
+scale_vector_inplace(vector_t V, double a)
+{
+    for (int i = 0; i < V.size; i++)
+        V.items[V.stride * i] *= a;
+}
+
+static scaled_vector_t
+scale_vector_lazily(vector_t V, double a)
+{
+    return (scaled_vector_t) {.base=V, .scale=a};
+}
+
+static double 
+scalar_product(vector_t V, vector_t U)
+{
+    assert(V.size == U.size);
+
+    double scale = 0;
+    for (int i = 0; i < V.size; i++)
+        scale += V.items[i * V.stride] * U.items[i * U.stride];
+    return scale;
+}
+
+static double 
+normalize_inplace(vector_t V)
+{
+    // Calculate the sum of the component's squares
+    double sum = scalar_product(V, V);
+
+    // Calculate the norm and scale the column
+    // only if the norm isn't zero.
+    double norm = sqrt(sum);
+    
+    if (norm != 0)
+        scale_vector_inplace(V, 1/norm);
+    
+    return norm;
+}
+
+static scaled_vector_t
+project(vector_t V, vector_t U)
+{
+    double scale_vu = scalar_product(V, U);
+    double scale_uu = scalar_product(U, U);
+    double ratio = scale_vu / scale_uu;
+    return scale_vector_lazily(U, ratio);
+}
+
+/** Gram-Schmidt orthonormalization
+ **/
+void lina_orthoNormGramSchmidt(double *A, double *Q, int n)
+{
+    square_matrix_t A2 = square_matrix_from_raw(A, n);
+    square_matrix_t Q2 = square_matrix_from_raw(Q, n);
+
+    for (int i = 0; i < n; i++) {
+
+        vector_t Qi = get_column_of_square_matrix(Q2, i);
+        vector_t Ai = get_column_of_square_matrix(A2, i);
+        copy_vector(Qi, Ai);
+
+        for (int j = 0; j < i; j++) {
+            vector_t Qj = get_column_of_square_matrix(Q2, j);
+            subtract_vector_inplace(Qi, project(Ai, Qj));
+        }
+
+        normalize_inplace(Qi);
+        // TODO: Handle case of zero norm 
+    }
+}
+
+void lina_decompQR(double *A, double *Q, double *R, int n)
+{
+    lina_orthoNormGramSchmidt(A, Q, n);
+
+    square_matrix_t A2 = square_matrix_from_raw(A, n);
+    square_matrix_t Q2 = square_matrix_from_raw(Q, n);
+
+    // Now calculate R by multiplying Q^t and A
+    for(int i = 0; i < n; i++) { // Iterate over each column i of Q..
+        for(int j = 0; j < n; j++) { // ..and over each column j of A
+            vector_t Qi = get_column_of_square_matrix(Q2, i);
+            vector_t Aj = get_column_of_square_matrix(A2, j);
+            R[i * n + j] = scalar_product(Qi, Aj);
+        }
+    }
+}

test.c

@@ -0,0 +1,58 @@
+#include <stdio.h>
+#include "src/lina.h"
+
+void print_square_matrix(double *M, int n, FILE *stream)
+{
+    for (int i = 0; i < n; i++)
+        {
+            fprintf(stream, "| ");
+            for (int j = 0; j < n; j++)
+                {
+                    fprintf(stderr, "%2.2f ", M[i * n + j]);
+                }
+            fprintf(stream, "|\n");
+        }
+    fprintf(stream, "\n");
+}
+
+void print_vector(double *V, int n, FILE *stream)
+{
+    fprintf(stream, "[ ");
+    for (int i = 0; i < n; i++)
+        fprintf(stderr, "%2.2f ", V[i]);
+    fprintf(stream, "]\n");
+}
+
+int main(void)
+{
+    double A[4] = {1, 2, 3, 4};
+    
+    double L[4];
+    double U[4];
+    double LU[4];
+    lina_decompLU(A, L, U, 2);
+    lina_dot(L, U, LU, 2, 2, 2);
+    print_square_matrix(A,  2, stderr);
+    print_square_matrix(L,  2, stderr);
+    print_square_matrix(U,  2, stderr);
+    print_square_matrix(LU, 2, stderr);
+
+    double det;
+    lina_det(A, 2, &det);
+    fprintf(stderr, "det=%2.2f\n", det);
+
+    double Q[4];
+    double R[4];
+    double QR[4];
+    lina_decompQR(A, Q, R, 2);
+    lina_dot(Q, R, QR, 2, 2, 2);
+    print_square_matrix(Q,  2, stderr);
+    print_square_matrix(R,  2, stderr);
+    print_square_matrix(QR, 2, stderr);
+
+    double E[4];
+    lina_eig(A, E, 2);
+    print_vector(E, 2, stderr);
+
+    return 0;
+}