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.c

@@ -5,6 +5,7 @@
 #include <stdio.h>
 #include <errno.h>
 #include <ctype.h>
+#include <math.h>
 #include "lina.h"
 
 /* Function: lina_dot
@@ -27,28 +28,27 @@
 **
 **   - This function can never fail.
 */
-void lina_dot(double *A, double *B, double *C, int m, int n, int l){
-
+void lina_dot(double *A, double *B, double *C, int m, int n, int l)
+{
     assert(m > 0 && n > 0 && l > 0);
     assert(A != NULL && B != NULL && C != NULL);
     assert(A != C && B != C);
 
     // Iteration over A's rows
-    for(int i = 0; i < m; i++)
-        {
+    for(int i = 0; i < m; i++) {
+
         // Iteration over B's columns
-            for(int k = 0; k < l; k++)
-                {
-                    double pos = 0;
+        for(int k = 0; k < l; k++) {
+
+            double sum = 0;
 
             // Iteration over the single B column 
             // for executing the product of sum
                     
             for(int j=0; j < n; j++)
+                sum += A[i * n + j] * B[j * l + k];
 
-                        pos += A[i*n + j] * B[j*l + k];
-
-                    C[i*l + k] = pos;
+            C[i * l + k] = sum;
         }
     }
 }
@@ -72,8 +72,8 @@ void lina_dot(double *A, double *B, double *C, int m, int n, int l){
 **
 **   - This function can never fail.
 */
-void lina_add(double *A, double *B, double *C, int m, int n){
-
+void lina_add(double *A, double *B, double *C, int m, int n)
+{
     assert(m > 0 && n > 0);
     assert(A != NULL && B != NULL && C != NULL);
     
@@ -122,8 +122,7 @@ void lina_transpose(double *A, double *B, int m, int n)
     assert(m > 0 && n > 0);
     assert(A != NULL && B != NULL);
 
-    if(m == 1 || n == 1)
-        {
+    if(m == 1 || n == 1) {
         // For a matrix with height or width of 1
         // row-major and column-major order coincide,
         // so the stransposition doesn't change the
@@ -132,9 +131,9 @@ void lina_transpose(double *A, double *B, int m, int n)
 
             if(A != B) // Does the copy or the branch cost more?
                 memcpy(B, A, sizeof(A[0]) * m * n);
-        }
-    else if(m == n)
-        {
+    
+    } else if(m == n) {
+
         // Iterate over the upper triangular portion of
         // the matrix and switch each element with the
         // corresponding one in the lower triangular portion.
@@ -145,15 +144,13 @@ void lina_transpose(double *A, double *B, int m, int n)
         //       is avoided.
 
         for(int i = 0; i < n; i += 1)
-                for(int j = 0; j < i+1; j += 1)
-                    {
+            for(int j = 0; j < i+1; j += 1) {
                 double temp = A[i*n + j];
                 B[i*n + j] = A[j*n + i];
                 B[j*n + i] = temp;
             }
-        }
-    else
-        {
+
+    } else {
         // Not only the matrix needs to be transposed
         // assuming the destination matrix is the same
         // as the source matrix, but the memory representation
@@ -174,8 +171,7 @@ void lina_transpose(double *A, double *B, int m, int n)
         double item = A[1];
         int    next = m;
 
-            while(next != 1)
-                {
+        while(next != 1) {
             double temp = A[next];
             B[next] = item;
             item = temp;
@@ -233,12 +229,9 @@ static int scanValue(FILE *fp, char *buffer, int max_length, char first, char *f
     // Scan the integer portion of
     // the numeric value and copy it
     // into the buffer.
-    do
-        {
-            if(n == max_length)
-                {
-                    // ERROR: Internal buffer is too small to hold
-                    //        the representation of this item.
+    do {
+    
+        if(n == max_length) {
             *error = "Internal buffer is too small to hold "
                      "the representation of a numeric value";
             return 0;
@@ -247,8 +240,8 @@ static int scanValue(FILE *fp, char *buffer, int max_length, char first, char *f
         buffer[n++] = c;
             
         c = getc(fp);
-        }
-    while(c != EOF && isdigit(c));
+    
+    } while(c != EOF && isdigit(c));
 
     // Did the integer part end with
     // a dot?
@@ -257,10 +250,8 @@ static int scanValue(FILE *fp, char *buffer, int max_length, char first, char *f
     // Now scan and copy the decimal 
     // part of the numeric value if 
     // a dot was found.
-    if(dot)
-        {
-            if(n == max_length)
-                {
+    if(dot) {
+        if(n == max_length) {
             // ERROR: Internal buffer is too small to hold
             //        the representation of this item.
             //        (The dot doesn't fit.)
@@ -273,18 +264,15 @@ static int scanValue(FILE *fp, char *buffer, int max_length, char first, char *f
 
         c = getc(fp);
 
-            if(!isdigit(c))
-                {
+        if(!isdigit(c)) {
             // ERROR: Got something other than a 
             //        digit after the dot.
             *error = "Got something other than a digit after the dot.";
             return 0;
         }
             
-            do
-                {
-                    if(n == max_length)
-                        {
+        do {
+            if(n == max_length) {
                 // ERROR: Internal buffer is too small 
                 //        to hold the representation of
                 //        this item.
@@ -296,8 +284,7 @@ static int scanValue(FILE *fp, char *buffer, int max_length, char first, char *f
             buffer[n++] = c;
 
             c = getc(fp);
-                }
-            while(c != EOF && isdigit(c));
+        } while(c != EOF && isdigit(c));
     }
 
     buffer[n] = '\0';
@@ -369,16 +356,14 @@ double *lina_loadMatrixFromStream(FILE *fp, int *width, int *height, char **erro
     while(c != EOF && isspace(c))
         c = getc(fp);
 
-    if(c == EOF)
-        {
+    if(c == EOF) {
         // ERROR: Stream ended before a matrix was
         //        found.
         *error = "Stream ended before a matrix was found";
         return NULL;
     }
 
-    if(c != '[')
-        {
+    if(c != '[') {
         // ERROR: Was expected a '[' as the first 
         //        character of a matrix, but got 
         //        something else instead.
@@ -393,8 +378,7 @@ double *lina_loadMatrixFromStream(FILE *fp, int *width, int *height, char **erro
     while(c != EOF && isspace(c))
         c = getc(fp);
 
-    if(c == EOF)
-        {
+    if(c == EOF) {
         // ERROR: Stream ended where a numeric value
         //        was expected. 
         *error = "Stream ended where a numeric value "
@@ -404,8 +388,7 @@ double *lina_loadMatrixFromStream(FILE *fp, int *width, int *height, char **erro
 
     double *matrix = malloc(sizeof(matrix[0]) * 64);
 
-    if(matrix == NULL)
-        {
+    if(matrix == NULL) {
         // ERROR: Insufficient memory.
         *error = "Insufficient memory";
         return NULL;
@@ -415,10 +398,8 @@ double *lina_loadMatrixFromStream(FILE *fp, int *width, int *height, char **erro
         w = -1, i = 0, j = 0;
 
     if(c != ']')
-        while(1)
-            {
-                if(!isdigit(c))
-                    {
+        while(1) {
+            if(!isdigit(c)) {
                 // ERROR: Got something other than a digit 
                 //        where a numeric value was expected.
                 *error = "Got something other than a numeric "
@@ -443,14 +424,12 @@ double *lina_loadMatrixFromStream(FILE *fp, int *width, int *height, char **erro
             assert(res == 1 || res == -1);
 
             // Make sure the matrix has enough space.
-                if(size == capacity)
-                    {
+            if(size == capacity) {
                 int new_capacity = capacity * 2;
 
                 double *temp = realloc(matrix, sizeof(double) * new_capacity);
 
-                        if(temp == NULL)
-                            {
+                if(temp == NULL) {
                     // ERROR: Insufficient memory.
                     *error = "Insufficient memory";
                     free(matrix);
@@ -470,8 +449,7 @@ double *lina_loadMatrixFromStream(FILE *fp, int *width, int *height, char **erro
             else
                 casted = strtod(buffer, NULL);
 
-                if(errno)
-                    {
+            if(errno) {
                 // ERROR: Failed to convert a numeric value
                 //        from it's string form to a numeric 
                 //        variable.
@@ -487,23 +465,20 @@ double *lina_loadMatrixFromStream(FILE *fp, int *width, int *height, char **erro
             while(c != EOF && isspace(c))
                 c = getc(fp);
 
-                if(c == ']' || c == ',')
-                    {
+            if(c == ']' || c == ',') {
                 // The matrix's row just ended.
 
                 if(w == -1)
                     // This was the first row.
                     w = i;
-                        else
-                            {
+                else {
                     // This wasn't the first row,
                     // so it's possible that it's
                     // length is different from the
                     // previous ones.
                     assert(w > -1);
 
-                                if(i != w)
-                                    {
+                    if(i != w) {
                         // ERROR: The j-th row has the wrong
                         //        number of elements.
                         if(i < w)
@@ -527,8 +502,7 @@ double *lina_loadMatrixFromStream(FILE *fp, int *width, int *height, char **erro
                     c = getc(fp);
             }
 
-                if(c == EOF)
-                    {
+            if(c == EOF) {
                 // ERROR: Stream ended inside a matrix, where 
                 //        either ',', ']' or a numeric value was
                 //        expected.
@@ -538,8 +512,7 @@ double *lina_loadMatrixFromStream(FILE *fp, int *width, int *height, char **erro
             }
         }
 
-    if(size == 0)
-        {
+    if(size == 0) {
         free(matrix);
         *error = "Empty matrix";
         return NULL;
@@ -551,8 +524,8 @@ double *lina_loadMatrixFromStream(FILE *fp, int *width, int *height, char **erro
     // build the matrix.
     int fragm_threshold = 30; // (It's a percentage)
     
-    if(100.0 * size/capacity < fragm_threshold)
-        {
+    if(100.0 * size/capacity < fragm_threshold) {
+
         int new_capacity = (size == 0) ? 1 : size;
 
         double *temp = realloc(matrix, new_capacity * sizeof(double));
@@ -566,3 +539,506 @@ double *lina_loadMatrixFromStream(FILE *fp, int *width, int *height, char **erro
 
     return matrix;
 }
+
+/* Function: lina_saveMatrixToStream
+**
+**   Save to the stream [fp] a matrix [A] encoding it as an
+**   ASCII sequence in the form:
+**
+**     [a b c .. , d e f .. , ..]
+**
+**   For instance, the 4x4 identity matrix will
+**   be encoded as:
+**
+**      [1 0 0 0, 0 1 0 0, 0 0 1 0, 0 0 0 1]
+**
+**   Since the matrix is in row-major order, the caller must
+**   specify the collumns and the rows of the matrix
+**   through [width] and [height] input arguments.
+**
+**   If an error occurres, a negative integer is returned 
+**   and a human-readable description of what happened 
+**   is returned through the [error] pointer.
+**
+** Notes:
+**   - It can be called multiple times on a stream to write
+**     more than one matrix on it.
+**
+**   - The [error] pointer is optional (it can be NULL).
+**
+**   - If the stream [fp] is NULL, then [stdout] is used.
+*/
+int lina_saveMatrixToStream(FILE *fp, double *A, int width, int height, char **error)
+{
+    assert(A != NULL);
+
+    char *dummy;
+    if (error == NULL)
+        error = &dummy;
+    else
+        *error = NULL;
+    
+    if (width < 1) {
+        *error = "The provided width is less than one";
+        return -1;
+    }
+    
+    if (height < 1) {
+        *error = "The provided height is less than one";
+        return -1;
+    }
+
+    if (fp == NULL)
+        fp = stdout;
+
+    putc('[',fp);
+
+    for (int i = 0; i < height-1; i++) {
+        for (int j = 0; j < width-1; j++)
+            fprintf(fp, "%f ", A[i*width + j]);
+
+        fprintf(fp, "%f, ", A[i*width + width-1]);
+    }
+
+    for (int j = 0; j < width-1; j++)
+        fprintf(fp, "%f ", A[(height-1)*width + j]);
+
+    fprintf(fp, "%f", A[(height-1)*width + width-1]);
+        
+    putc(']',fp);
+
+    return 0;
+}
+
+void lina_conv(double *A, double *B, double *C, 
+               int Aw, int Ah, int Bw, int Bh)
+{
+    assert(A != NULL && B != NULL && C != NULL);
+    assert(A != B && B != C && C != A);
+    assert(Aw > 0 && Ah > 0 && Bw > 0 && Bh > 0);
+    assert((Bw & 1) && (Bh & 1)); // B must have odd height and width.
+
+    // NOTE: The output C matrix is smaller than 
+    //       A proportionally to B's size.
+    
+    int Cw = Aw - Bw + 1;
+    int Ch = Ah - Bh + 1;
+    assert(Cw > 0 && Ch > 0);
+
+    // Iterate over each pixel of the result matrix..
+    for(int j = 0; j < Ch; j += 1)
+        for(int i = 0; i < Cw; i += 1) {
+            // ..and calculate it's value as
+            // the scalar product between the
+            // mask B and a portion of A.
+
+            C[j * Cw + i] = 0;
+            for(int v = 0; v < Bh; v += 1)
+                for(int u = 0; u < Bw; u += 1)
+                    C[j * Cw + i] += A[(i - Bw/2 + u) * Aw + (i - Bh/2 + v)] * B[v * Bw + u];
+        }
+}
+
+void lina_reallyP(int *P, double *P2, int n)
+{
+    memset(P2, 0, sizeof(double) * n * n);
+
+    for (int i = 0; i < n; i++)
+        P2[i * n + P[i]] = 1;
+}
+
+int lina_decompLUP(double *A, double *L, 
+                   double *U, int    *P, 
+                   int n)
+{
+    assert(n > 0);
+    assert(A != L && A != U && L != U);
+
+    for (int i = 0; i < n; i++)
+        P[i] = i;
+
+    int swaps = 0;
+    for (int i = 0; i < n; i++) {
+
+        int v = P[i];
+        double max_v = A[v * n + i];
+        int    max_i = i;
+        
+        for (int j = i+1; j < n; j++) {
+            int u = P[j];
+            double abs = fabs(A[u * n + j]);
+            if (abs > max_v) {
+                max_v = abs;
+                max_i = j;
+            }
+        }
+
+        if (max_i != i) {
+
+            // Swap rows
+            int temp = P[i];
+            P[i] = P[max_i];
+            P[max_i] = temp;
+
+            swaps++;
+        }
+    }
+
+    for (int i = 0; i < n; i++)
+        for (int j = 0; j < n; j++)
+            U[i * n + j] = A[P[i] * n + j];
+
+    memset(L, 0, sizeof(double) * n * n);
+    for (int i = 0; i < n; i++)
+        L[i * n + i] = 1;
+
+    for (int i = 0; i < n; i++)
+        for (int j = i+1; j < n; j++) {
+            double u = U[i * n + i];
+            L[j * n + i] = U[j * n + i] / u;
+            for (int k = 0; k < n; k++)
+                U[j * n + k] -= L[j * n + i] * U[i * n + k];
+        }
+
+    return swaps;
+}
+
+static void 
+printSquareMatrix(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");
+}
+
+/* Function: lina_det
+**
+**   Calculates the determinant of the n by n matrix A
+**   and returns it throught the output parameter [det].
+**
+**   If not enough memory is available, false is returned,
+**   else true is returned.
+**
+** Notes:
+**   - The output parameter [det] is optional. (you can
+**     ignore the result by passing NULL).
+*/
+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 + sizeof(int) * n);
+    if (T == NULL)
+        return false;
+
+    // Do the decomposition
+    double *L = T;
+    double *U = L + (n * n);
+    int    *P = (int*) (U + (n * n));
+    
+    int swaps = lina_decompLUP(A, L, U, P, n);
+    if (swaps < 0) {
+        free(T);
+        return false;
+    }
+
+    // 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++) {
+        double l = L[i * n + i];
+        double u = U[i * n + i];
+        prod *= l * u;
+    }
+
+    if (swaps & 1)
+        prod = -prod;
+
+    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 && eps > 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; j++)
+            if (A[i * n + j] > eps)
+                return false;
+
+    // Now check that the subdiagonal 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;
+    }
+
+    // NOTE: Ideas were taken from [https://math.stackexchange.com/questions/4352389/exact-stop-condition-for-qr-algorithm]
+    return true;
+}
+
+/* Function: lina_eig
+**
+**   Calculates the eigenvalues of the n by n matrix M
+**   using the (unshifted) QR algorithm and stores them 
+**   in the E vector.
+**
+**   If not enough memory is available, this function
+**   aborts returning false. If all went well, true is
+**   returned.
+**
+** Algorithm:
+**
+**   The algorithm works by decomposing the M matrix into
+**   the product of two matrices Q and R, such that Q is
+**   orthonormal and R is upper triangular:
+**
+**       M = QR
+**
+**   Q and R are then multiplied in inverse order to obtain
+**   a new matrix M1, which is then decomposed in two new
+**   matrices Q1,R1. The algorithm is iterated n times until 
+**   the matrix Mn is upper triangular:
+**
+**       M = QR  ->  RQ = M(1)
+**
+**       M(1) = Q(1)R(1)  ->  R(1)Q(1) = M(2)
+**
+**       M(2) = Q(2)R(2)  ->  R(2)Q(2) = M(3)
+**
+**              ...
+**
+**       M(n-1) = Q(n-1)R(n-1)  ->  R(n-1)Q(n-1) = M(n)
+**
+**       M(n) <--- Triangular!
+**
+**   The eigenvalues of M(n) are the same as M. Being upper 
+**   triangular, M(n) has its eigenvalues on its diagonal,
+**   so we just need to scan the diagonal and store it into
+**   the E vector. If the original matrix has complex roots,
+**   the M(n) sequence will converge to a matrix with a 
+**   non-zero 2x2 block on the diagonal for each pair of
+**   complex roots. If that's the case, these blocks must
+**   be unpacked into the complex values using the quadratic
+**   formula.
+**   
+*/
+bool lina_eig(double *M, double complex *E, int n)
+{
+    // Allocate space for three matrices n by n
+    double *T = malloc(sizeof(double) * n * n * 3);
+    if (T == NULL)
+        return false;
+
+    double *A = T;
+    double *Q = A + n * n;
+    double *R = Q + n * n;
+    memcpy(A, M, sizeof(double) * n * n);
+
+    // At least 100 iterations are done. This is because
+    // the QR algorithm doesn't allow complex eigenvalues,
+    // so the A matrix may converge to a matrix with 2x2
+    // blocks on the diagonal. In general, the algorithm
+    // must iterate until the end result is triangular,
+    // but that may never be the case, so we end when the
+    // result matrix is "kind of triangular" (triangular
+    // with 2x2 blocks on the diagonal). But by using this
+    // rule, a 2x2 matrix will be considered as tringular
+    // from the start, which is not right! That's why we
+    // do at least 100 warm-up iterations.
+    double eps = 0.1;
+    int batch = 100;
+    do {
+        for (int i = 0; i < batch; 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, eps));
+
+    // Now we export the diagonal of the iteration result
+    // also looking out for 2x2 diagonal blocks, in which
+    // case we need to unpack their complex eigenvalues
+    for (int i = 0; i < n; i++) {
+
+        // The current diagonal entry is A[i*n + i],
+        // so if this is the first entry of a 2x2 block,
+        // its lower entry A[(i+1)*n + i] will be non-zero
+        if (i+1 < n && fabs(A[(i+1) * n + i]) > eps) {
+
+            // It's a 2x2 block. Unpack the complex eigenvalues
+            // using the quadratic formula. (Is there a better
+            // way?)
+            
+            double a = A[(i+0) * n + (i+0)];
+            double b = A[(i+0) * n + (i+1)];
+            double c = A[(i+1) * n + (i+0)];
+            double d = A[(i+1) * n + (i+1)];
+            
+            // Given the block is:
+            //
+            //   | a  b |
+            //   | c  d |
+            //
+            // Then the eigenvalues are the roots of:
+            //
+            //   det(| a-y  b   |) = (a-y)(d-y) - bc = y^2 - (a + d)y + (ad - bc)
+            //       | c    d-y |
+            //
+            // For simplicity:
+            //
+            //   D = (a + d)^2 - 4(ad - bc)
+            //
+            // so that
+            //
+            //   y1, y2 = (a + d)/2 +/- 1/2 sqrt{D}
+            //
+            // y1 and y2 are one the conjugate of the other. Their
+            // real part is
+            // 
+            //   Re{y1, y2} = (a+d)/2
+            //
+            // While their immaginary part (in absolute value) is
+            //
+            //   Imm{y1, y2} = 1/2 sqrt{-D}
+
+            double D = (a+d)*(a+d) - 4*(a*d - b*c);
+            assert(D < 0);
+
+            double re = 0.5 * (a+d);
+            double im = 0.5 * sqrt(-D);
+
+            double complex y1 = re + im * I;
+            double complex y2 = re - im * I;
+
+            // Now place the results into the output vector
+            // and tell the loop to skip one iteration
+            E[i]   = y1; 
+            E[i+1] = y2;
+            i++;
+
+        } else
+            E[i] = A[i * n + i];
+    }
+
+    free(T);
+    return true;
+}
+
+/* Create the n-1 by n-1 matrix D obtained by
+** removing the [del_col] column and [del_row]
+** frow the n by n matrix M.
+*/
+static void 
+copyMatrixWithoutRowAndCol(double *M, double *D, int n, 
+                           int del_col, int del_row)
+{
+    // Copy the upper-left portion of matrix M
+    // that comes before the deleted column and
+    // row.
+    for (int i = 0; i < del_row; i++)
+        for (int j = 0; j < del_col; j++)
+            D[i * (n-1) + j] = M[i * n + j];
+
+    // Copy the lower left portion that comes
+    // after both the deleted column and row.
+    for (int i = del_row+1; i < n; i++)
+        for (int j = del_col+1; j < n; j++)
+            D[(i-1) * (n-1) + (j-1)] = M[i * n + j];
+
+    // Copy the bottom portion that comes after
+    // the deleted row but before the deleted column.
+    for (int i = del_row+1; i < n; i++)
+        for (int j = 0; j < del_col; j++)
+            D[(i-1) * (n-1) + j] = M[i * n + j];
+
+    // Copy the right portion that comes after
+    // the deleted column but before the deleted row.
+    for (int i = 0; i < del_row; i++)
+        for (int j = del_col+1; j < n; j++)
+            D[i * (n-1) + (j-1)] = M[i * n + j];
+}
+
+bool lina_inverse(double *M, double *D, int n)
+{
+    double det;
+    if (!lina_det(M, n, &det))
+        return false;
+
+    if (det == 0)
+        return false; // The matrix can't be inverted
+
+    double *T = malloc(sizeof(double) * ((n-1) * (n-1) + n * n));
+    if (T == NULL)
+        return false;
+
+    double *M_t = T + (n-1) * (n-1);    
+    lina_transpose(M, M_t, n, n);
+
+    for (int i = 0; i < n; i++)
+        for (int j = 0; j < n; j++) {
+            
+            copyMatrixWithoutRowAndCol(M_t, T, n, j, i);
+
+            double det2;
+            if (!lina_det(T, n-1, &det2)) {
+                free(T);
+                return false;
+            }
+
+            // If the determinant of M isn't zero,
+            // neither is this!
+            assert(det2 != 0);
+
+            bool i_is_odd = i & 1;
+            bool j_is_odd = j & 1;
+            int sign = (i_is_odd == j_is_odd) ? 1 : -1;
+
+            D[i * n + j] = sign * det2 / det;
+        }
+
+    free(T);
+    return true;
+}

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);
+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;
+}