bug fix in matrix inversion routine

README.md

@@ -1,2 +1,2 @@
-# Lina
-Lina (***Lin**ear **A**lgebra*) is a C library that implements common linear algebra operations.
+# Lina, the nice-to-read linear algebra toolkit!
+Lina (***Lin**ear **A**lgebra*) is a C library that implements common linear algebra operations with the aim to be nice to read!

build.sh

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

src/lina.h

@@ -1,6 +1,18 @@
+#include <complex.h>
+#include <stdbool.h>
 
-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);
-void    lina_scale(double *A, double *B, double k, int m, int n);
-void    lina_transpose(double *A, double *B, int m, int n);
-double *lina_loadMatrixFromStream(FILE *fp, int *width, int *height, char **error);
+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);
+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);
+bool lina_inverse(double *M, double *D, int n);
+void lina_conv(double *A, double *B, double *C, int Aw, int Ah, int Bw, int Bh);
+bool lina_eig(double *M, double complex *E, int n);
+void lina_reallyP(int *P, double *P2, int n);
+int  lina_decompLUP(double *A, double *L, double *U, int *P, int n);
+void lina_decompQR(double *A, double *Q, double *R, int n);
+void lina_orthoNormGramSchmidt(double *A, double *Q, int n);
+
+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 
+calculate_norm(vector_t V)
+{
+    double sum_of_squares = scalar_product(V, V);
+    return sqrt(sum_of_squares);
+}
+
+static double 
+normalize_inplace(vector_t V)
+{
+    double norm = calculate_norm(V);
+    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,135 @@
+#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 complex *V, int n, FILE *stream)
+{
+    fprintf(stream, "[ ");
+    for (int i = 0; i < n; i++)
+        fprintf(stderr, "(%2.2f + i%2.2f) ", creal(V[i]), cimag(V[i]));
+    fprintf(stream, "]\n");
+}
+
+int main(void)
+{
+    
+    double M[25] = {
+        1, 2, 3, 4, 5,
+        5, 1, 2, 3, 4,
+        4, 5, 1, 2, 3,
+        3, 4, 5, 1, 2,
+        2, 3, 4, 5, 1,
+    };
+
+    fprintf(stderr, "# --- M --- #\n");
+    print_square_matrix(M,  5, stderr);
+    
+
+    /*
+    double  L[25];
+    double  U[25];
+    int     P[5];
+    double P2[25];
+    lina_decompLUP(M, L, U, P, 5);
+    lina_reallyP(P, P2, 5);
+
+    fprintf(stderr, "# --- L --- #\n");
+    print_square_matrix(L, 5, stderr);
+
+    fprintf(stderr, "# --- U --- #\n");
+    print_square_matrix(U, 5, stderr);
+
+    fprintf(stderr, "# --- P2 --- #\n");
+    print_square_matrix(P2, 5, stderr);
+
+    double PA[25];
+    lina_dot(P2, M, PA, 5, 5, 5);
+    fprintf(stderr, "# --- PA --- #\n");
+    print_square_matrix(PA,  5, stderr);
+
+    double LU[25];
+    lina_dot(L, U, LU, 5, 5, 5);
+    fprintf(stderr, "# --- LU --- #\n");
+    print_square_matrix(LU, 5, stderr);
+
+    double det;
+    lina_det(M, 5, &det);
+    fprintf(stderr, "det=%2.2f\n", det);
+
+    fprintf(stderr, "# --- eig(M) --- #\n");
+    double complex E[5];
+    lina_eig(M, E, 5);
+    print_vector(E, 5, stderr);
+    */
+
+    
+    double invM[25];
+    lina_inverse(M, invM, 5);
+
+    double expI[25];
+    lina_dot(M, invM, expI, 5, 5, 5);
+
+    fprintf(stderr, "# --- inv(M) --- #\n");
+    print_square_matrix(invM, 5, stderr);
+
+    fprintf(stderr, "# --- I? --- #\n");
+    print_square_matrix(expI, 5, stderr);
+    
+
+    /*
+    double M[16] = {
+        1, 5, 4, 2,
+        2, 1, 5, 3,
+        4, 3, 2, 5,
+        5, 4, 3, 1,
+    };
+
+    fprintf(stderr, "# --- M --- #\n");
+    print_square_matrix(M, 4, stderr);
+
+    double L[16];
+    double U[16];
+    int P[4];
+    lina_decompLUP(M, L, U, P, 4);
+    fprintf(stderr, "# --- L,U,P --- #\n");
+    print_square_matrix(L, 4, stderr);
+    print_square_matrix(U, 4, stderr);
+
+    fprintf(stderr, "[ ");
+    for(int i = 0; i < 4; i++)
+        fprintf(stderr, "%d ", P[i]);
+    fprintf(stderr, "]\n");
+
+    double RP[16];
+    double PM[16];
+    lina_reallyP(P, RP, 4);
+    lina_dot(RP, M, PM, 4, 4, 4);
+
+    double LU[16];
+    lina_dot(L, U, LU, 4, 4, 4);
+
+    fprintf(stderr, "# --- P,PM,LU --- #\n");
+    print_square_matrix(RP, 4, stderr);
+    print_square_matrix(PM, 4, stderr);
+    print_square_matrix(LU, 4, stderr);
+
+    double det;
+    lina_det(M, 4, &det);
+
+    fprintf(stderr, "det(M) = %2.2f\n", det);
+    */
+    return 0;
+}