numpy performance reached on reasonable matrix size

.gitignore

@@ -1,6 +1,3 @@
 .swp
-test
-test2
-test_loader
-time
 .vscode
+*.txt

README.md

@@ -1,2 +1,4 @@
 # 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!
+
+The performance branch focuses only on the core functionalities of lina and aims to produce faster and reliable routines.

benchmarks/bench_dot.c

@@ -0,0 +1,181 @@
+#include <time.h>
+#include <stdint.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <assert.h>
+#include "../src/lina.h"
+
+#define A_ROWS 960llu
+#define A_COLS 960llu
+
+#define B_ROWS 960llu
+#define B_COLS 960llu
+
+int saveMatrixToStream(FILE *fp, double *A, int width, int height, char **error);
+static uint64_t nanos();
+
+int main()
+{
+    uint64_t ops = A_ROWS*B_COLS*2*A_COLS;
+    uint64_t start,stop,lina_dot_time, lina_dot1_time, lina_dot2_time, lina_dot3_time, lina_dot4_time;
+    
+    double *A = (double *)aligned_alloc(32,sizeof(double)*A_ROWS*A_COLS);
+    double *B = (double *)aligned_alloc(32,sizeof(double)*B_ROWS*B_COLS);
+
+    double *C1 = (double *)aligned_alloc(32,sizeof(double)*A_ROWS*B_COLS);
+    double *C2 = (double *)aligned_alloc(32,sizeof(double)*A_ROWS*B_COLS);
+    double *C3 = (double *)aligned_alloc(32,sizeof(double)*A_ROWS*B_COLS);
+    double *C4 = (double *)aligned_alloc(32,sizeof(double)*A_ROWS*B_COLS);
+    double *C5 = (double *)aligned_alloc(32,sizeof(double)*A_ROWS*B_COLS);
+    
+    for (int i = 0; i < A_ROWS*A_COLS; i++)
+        A[i] = (double)(rand()%2);
+    for (int i = 0; i < B_ROWS*B_COLS; i++)
+        B[i] = (double)(rand()%2);
+    for (int i = 0; i < A_ROWS*B_COLS; i++)
+    {
+        C1[i] = (double)(rand()%2);
+        C2[i] = (double)(rand()%2);
+        C3[i] = (double)(rand()%2);
+        C4[i] = (double)(rand()%2);
+        C5[i] = (double)(rand()%2);
+    }
+
+    start = nanos();
+    lina_dot(A,B,C1,A_ROWS,A_COLS,B_COLS);
+    stop = nanos();
+
+    lina_dot_time = stop-start;
+    
+    start = nanos();
+    lina_dot1(A,B,C2,A_ROWS,A_COLS,B_COLS);
+    stop = nanos();
+
+    lina_dot1_time = stop-start;
+
+    start = nanos();
+    lina_dot2(A,B,C3,A_ROWS,A_COLS,B_COLS);
+    stop = nanos();
+
+    lina_dot2_time = stop-start;
+
+    start = nanos();
+    lina_dot3(A,B,C4,A_ROWS,A_COLS,B_COLS);
+    stop = nanos();
+
+    lina_dot3_time = stop-start;
+
+    start = nanos();
+    lina_dot4(A,B,C5,A_ROWS,A_COLS,B_COLS);
+    stop = nanos();
+
+    lina_dot4_time = stop-start;
+
+    if(!memcmp(C1,C2,sizeof(double)*A_ROWS*B_COLS)
+        && !memcmp(C2,C3,sizeof(double)*A_ROWS*B_COLS)
+        && !memcmp(C3,C4,sizeof(double)*A_ROWS*B_COLS)
+        && !memcmp(C4,C5,sizeof(double)*A_ROWS*B_COLS))
+    {
+        printf( "lina_dot : %f GFLOPS\n"
+                "lina_dot1: %f GFLOPS\n"
+                "lina_dot2: %f GFLOPS\n"
+                "lina_dot3: %f GFLOPS\n"
+                "lina_dot4: %f GFLOPS\n",   (double)ops/lina_dot_time,
+                                            (double)ops/lina_dot1_time,
+                                            (double)ops/lina_dot2_time,
+                                            (double)ops/lina_dot3_time,
+                                            (double)ops/lina_dot4_time);
+        FILE *fp = fopen("lina_dots_success.txt", "w");
+        if (!fp)
+            return -1;
+        
+        saveMatrixToStream(fp,C1,A_ROWS,A_COLS,NULL);
+        fprintf(fp,"\nFINE\n");
+        saveMatrixToStream(fp,C2,A_ROWS,A_COLS,NULL);
+        fprintf(fp,"\nFINE\n");
+        saveMatrixToStream(fp,C3,A_ROWS,A_COLS,NULL);
+        fprintf(fp,"\nFINE\n");
+        saveMatrixToStream(fp,C4,A_ROWS,A_COLS,NULL);
+        fprintf(fp,"\nFINE\n");
+        saveMatrixToStream(fp,C5,A_ROWS,A_COLS,NULL);
+
+        fclose(fp);
+    }
+    else
+    {
+        printf("ERRORE: i prodotti matriciali sono diversi!\n");
+        FILE *fp = fopen("lina_dots_error.txt", "w");
+        if (!fp)
+            return -1;
+        
+        saveMatrixToStream(fp,C1,A_ROWS,A_COLS,NULL);
+        fprintf(fp,"\nFINE\n");
+        saveMatrixToStream(fp,C2,A_ROWS,A_COLS,NULL);
+        fprintf(fp,"\nFINE\n");
+        saveMatrixToStream(fp,C3,A_ROWS,A_COLS,NULL);
+        fprintf(fp,"\nFINE\n");
+        saveMatrixToStream(fp,C4,A_ROWS,A_COLS,NULL);
+        fprintf(fp,"\nFINE\n");
+        saveMatrixToStream(fp,C5,A_ROWS,A_COLS,NULL);
+
+        fclose(fp);
+    }
+    
+    free(A);
+    free(B);
+    free(C1);
+    free(C2);
+    free(C3);
+    free(C4);
+    free(C5);
+}
+
+static uint64_t nanos()
+{
+    struct timespec time;
+    clock_gettime(CLOCK_MONOTONIC, &time);
+    return (uint64_t)time.tv_sec*1000000000 + (uint64_t)time.tv_nsec;
+}
+
+int 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, "%.1f ", A[i*width + j]);
+
+        fprintf(fp, "%.1f,\n", A[i*width + width-1]);
+    }
+
+    for (int j = 0; j < width-1; j++)
+        fprintf(fp, "%.1f ", A[(height-1)*width + j]);
+
+    fprintf(fp, "%.1f", A[(height-1)*width + width-1]);
+        
+    putc(']',fp);
+
+    return 0;
+}

benchmarks/benchmarks.sh

@@ -0,0 +1,3 @@
+gcc bench_dot.c ../src/lina.c -O3 -march=native -ffast-math -funroll-loops -o bench_dot
+./bench_dot
+python3 py_dot.py

benchmarks/py_dot.py

@@ -0,0 +1,22 @@
+import os
+os.environ['OMP_NUM_THREADS'] = '1'
+
+import numpy as np
+import time
+
+N = 1024
+
+if __name__ == "__main__":
+    A = np.random.randn(N,N).astype(np.float64)
+    B = np.random.randn(N,N).astype(np.float64)
+
+    start = time.monotonic()
+    C = A @ B
+    stop = time.monotonic()
+    
+    s = stop-start
+    
+    ops = 2*N*N*N
+    
+    print(f"NUMPY: {ops/s * 1e-9} GFLOPS\n")
+    

benchmarks/time.c

@@ -1,87 +0,0 @@
-#include <time.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <assert.h>
-#include <string.h>
-#include "lina.h"
-
-/* This program compares the lina_transpose
-** implementation against the naive implementation.
-** Build it with:
-**   $ gcc time.c lina.c -o time -Wall -Wextra -O3
-*/
-
-#define check assert
-
-static void naive_transpose(double *A, double *B, int m, int n)
-{
-    assert(m > 0 && n > 0);
-    assert(A != NULL && B != NULL);
-
-    double *support;
-
-    if(A == B)
-    	{
-    		support = malloc(sizeof(*support) * m * n);
-    		
-    		check(support != NULL);
-
-		    memcpy(support, A, sizeof(*support) * m * n);
-    	}
-    else
-    	{
-    		support = A;
-    	}
-
-    for(int i = 0; i < n; i++)
-        
-        for(int j = 0; j < m; j++)
-
-            B[j*n + i] = support[i*m + j];
-
-    if(support != A)
-	    free(support);
-}
-
-// Wrap transposing functions and return their
-// execution time.
-static double time_transposition(void (*callback)(double*, double*, int, int), double *A, double *B, int m, int n)
-{
-	clock_t begin = clock();
-
-	callback(A, B, m, n);
-
-	clock_t end = clock();
-	
-	return (double) (end - begin) / CLOCKS_PER_SEC;
-}
-
-int main()
-{
-	int m = 1000;
-	int n = 100000;
-
-	double *big = malloc(sizeof(double) * m * n);
-	check(big != NULL);
-
-	memset(big, 0, sizeof(double) * m * n);
-
-	printf("lina_transpose took %gms (in-place)\n", 
-		1000 * time_transposition(lina_transpose, big, big, m, n));
-
-	printf("naive_transpose took %gms (in-place)\n", 
-		1000 * time_transposition(naive_transpose, big, big, m, n));
-
-	double *big2 = malloc(sizeof(double) * m * n);
-	check(big2 != NULL);
-
-	printf("lina_transpose took %gms\n", 
-		1000 * time_transposition(lina_transpose, big, big2, m, n));
-	
-	printf("naive_transpose took %gms\n", 
-		1000 * time_transposition(naive_transpose, big, big2, m, n));
-
-	free(big);
-	free(big2);
-	return 0;
-}

build.sh

@@ -1,2 +0,0 @@
-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

@@ -1,12 +1,14 @@
-#include <stddef.h>
 #include <assert.h>
 #include <stdlib.h>
 #include <string.h>
-#include <stdio.h>
 #include <errno.h>
 #include <ctype.h>
-#include <math.h>
 #include "lina.h"
+#include <immintrin.h>
+#include <stdint.h>
+
+static void
+dot_kernel_6x8(double *A_sub, double *B_sub, double *C_sub, int x, int y, int c_min, int c_max, int n, int l);
 
 /* Function: lina_dot
 **
@@ -43,7 +45,7 @@ void lina_dot(double *A, double *B, double *C, int m, int n, int l)
             double sum = 0;
 
             // Iteration over the single B column 
-            // for executing the product of sum
+            // for executing the sum of product
                     
             for(int j=0; j < n; j++)
                 sum += A[i * n + j] * B[j * l + k];
@@ -53,6 +55,370 @@ void lina_dot(double *A, double *B, double *C, int m, int n, int l)
     }
 }
 
+/* Function: lina_dot1
+**
+**   Evaluates the dot product C = A * B. The A,B
+**   matrices are, respectively, mxn and nxl, which
+**   means C is mxl. The resulting C matrix is stored
+**   in a memory  region specified by the caller. 
+**
+** Variant 1 of lina_dot:
+**   The idea of this variant is that inverting the order
+**   of the first and the third loop cicle we can avoid the
+**   rolling sum and so breaking the depencency chain 
+**   among subsequent add thus increasing the IPC.
+**
+** Notes:
+**
+**   - A,B must be provided as contiguous memory regions
+**     represented in row-major order. Also, C is stored
+**     that way too.
+**
+**   - The C pointer CAN'T refer to the same memory region 
+**     of either A or B.
+**
+**   - m,n,l must be greater than 0.
+**
+**   - This function can never fail.
+*/
+void lina_dot1(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);
+
+    // Since the C matrix can contain any value,
+    // this first pass is done to overwrite the values
+
+    // Iteration over A's rows
+    for(int i = 0; i < m; i++) {
+        // Iteration over B's columns
+        for(int k = 0; k < l; k++)
+            C[i * l + k] = A[i * n] * B[k];
+    }
+
+    // Iteration over the single B column 
+    // for executing the sum of product
+    for(int j=1; j < n; j++)
+    {
+        // Iteration over A's rows
+        for(int i = 0; i < m; i++) {
+            // Iteration over B's columns
+            for(int k = 0; k < l; k++)
+                C[i * l + k] += A[i * n + j] * B[j * l + k];
+        }
+    }
+}
+
+/* Function: lina_dot2
+**
+**   Evaluates the dot product C = A * B. The A,B
+**   matrices are, respectively, mxn and nxl, which
+**   means C is mxl. The resulting C matrix is stored
+**   in a memory  region specified by the caller. 
+**
+** Variant 2 of lina_dot:
+**   Other than inverting the order of the first and the
+**   third loop cicle this version does the dot product in block
+**   of 32x32 values. Doing so the number of cache misses decreases.
+**
+** Notes:
+**
+**   - A,B must be provided as contiguous memory regions
+**     represented in row-major order. Also, C is stored
+**     that way too.
+**
+**   - The C pointer CAN'T refer to the same memory region 
+**     of either A or B.
+**
+**   - m,n,l must be greater than 0.
+**
+**   - This function can never fail.
+*/
+void lina_dot2(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);
+
+    // This size is based on experimental results
+    #define BLOCKSIZE 32
+
+    const int br_max = (m & ~(BLOCKSIZE - 1));
+    const int bc_max = (l & ~(BLOCKSIZE - 1));
+
+    // Dealing with the squared submatrix of C
+    for (int br = 0; br < br_max; br += BLOCKSIZE)
+    {
+        for (int bc = 0; bc < bc_max; bc += BLOCKSIZE)
+        {
+            double block[BLOCKSIZE*BLOCKSIZE];
+            
+            // 1. Compute block
+
+            // Iteration over A's rows
+            for(int i = br; i < br+BLOCKSIZE; i++) {
+
+                // Iteration over B's columns
+                for(int k = bc; k < bc+BLOCKSIZE; k++)
+                    block[(i-br)*BLOCKSIZE + (k-bc)] = A[i * n] * B[k];
+            }
+
+            // Iteration over the single B column 
+            // for executing the sum of product
+            for(int j=1; j < n; j++)
+            {
+                // Iteration over A's rows
+                for(int i = br; i < br+BLOCKSIZE; i++) {
+
+                    // Iteration over B's columns
+                    for(int k = bc; k < bc+BLOCKSIZE; k++)
+                        block[(i-br)*BLOCKSIZE + (k-bc)] += A[i * n + j] * B[j * l + k];
+                }
+            }
+
+            // 2. Copy block to C
+            for (int i = 0; i < BLOCKSIZE; i++)
+                memcpy(&C[(i+br)*l + bc],&block[i*BLOCKSIZE], sizeof(double)*BLOCKSIZE);
+        }
+    }
+
+    // Dealing with the last rows and cols
+    
+    // Last rows
+    // Iteration over A's rows
+    for(int i = br_max; i < m; i++) {
+        // Iteration over B's columns
+        for(int k = 0; k < l; k++)
+            C[i*l + k] = A[i * n ] * B[k];
+    }
+
+    // Last cols
+    // Iteration over A's rows
+    for (int i = 0; i < br_max; i++)
+    {
+        // Iteration over B's columns
+        for(int k = bc_max; k < l; k++)
+            C[i*l + k] = A[i * n] * B[k];
+    }
+
+    // Iteration over the single B column 
+    // for executing the product of sum
+    for(int j=1; j < n; j++)
+    {
+        // Iteration over A's rows
+        for(int i = br_max; i < m; i++) {
+            // Iteration over B's columns
+            for(int k = 0; k < l; k++)
+                C[i*l + k] += A[i * n + j] * B[j * l + k];
+        }
+
+        // Iteration over A's rows
+        for (int i = 0; i < br_max; i++)
+        {
+            // Iteration over B's columns
+            for(int k = bc_max; k < l; k++)
+                C[i*l + k] += A[i * n + j] * B[j * l + k];
+        }
+    }
+}
+
+/* Function: lina_dot3
+**
+**   Evaluates the dot product C = A * B. The A,B
+**   matrices are, respectively, mxn and nxl, which
+**   means C is mxl. The resulting C matrix is stored
+**   in a memory  region specified by the caller. 
+**
+** Variant 3 of lina_dot:
+**   This include the changes of lina_dot2 but uses
+**   simd instructions to compute products and sums.
+**
+** Notes:
+**
+**   - A,B must be provided as contiguous memory regions
+**     represented in row-major order. Also, C is stored
+**     that way too.
+**
+**   - The C pointer CAN'T refer to the same memory region 
+**     of either A or B.
+**
+**   - m,n,l must be greater than 0.
+**
+**   - This function can never fail.
+*/
+void lina_dot3(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);
+
+    // This size is based on experimental results
+    #define BLOCK_ROWS 6
+    #define BLOCK_COLS 8
+
+    const int br_max = (m & ~(BLOCK_ROWS - 1));
+    const int bc_max = (l & ~(BLOCK_COLS - 1));
+
+    __m256d *Bm = (__m256d *)B;
+    __m256d *Cm = (__m256d *)C;
+
+    // problema: B non è allineato a 32 byte, cosa che pare essere il problema
+
+    // Dealing with the squared submatrix of C
+    for (int br = 0; br < br_max; br += BLOCK_ROWS)
+    {
+        for (int bc = 0; bc < bc_max; bc += BLOCK_COLS)
+        {
+            __m256d mblock[BLOCK_ROWS][BLOCK_COLS/4] = {0};
+            
+            // 1. Compute block
+
+            for(int j=0; j < n; j++)
+            {
+                for(int i = 0; i < BLOCK_ROWS; i++)
+                {
+                    __m256d A_brdcst = _mm256_broadcast_sd(&A[(i+br) * n + j]);
+                    for(int k = 0; k < BLOCK_COLS/4; k++)
+                    {
+                            mblock[i][k] = _mm256_fmadd_pd(A_brdcst, Bm[(j * l + bc)/4 + k], mblock[i][k]);
+                    }
+                }
+                // Iteration over A's rows
+            }
+
+            // 2. Copy block to C
+            for (int i = 0; i < BLOCK_ROWS; i++)
+                for (int j = 0; j < BLOCK_COLS/4; j++)
+                   Cm[((i+br)*l + bc)/4 + j] = mblock[i][j];
+        }
+    }
+
+    // Dealing with the last rows and cols
+    //printf("br_max: %d\nbc_max: %d\n",br_max,bc_max);
+    // Last rows
+    // Iteration over A's rows
+    for(int i = br_max; i < m; i++) {
+        // Iteration over B's columns
+        for(int k = 0; k < l; k++)
+            C[i*l + k] = A[i * n ] * B[k];
+    }
+
+    // Last cols
+    // Iteration over A's rows
+    for (int i = 0; i < br_max; i++)
+    {
+        // Iteration over B's columns
+        for(int k = bc_max; k < l; k++)
+            C[i*l + k] = A[i * n] * B[k];
+    }
+
+    // Iteration over the single B column 
+    // for executing the product of sum
+    for(int j=1; j < n; j++)
+    {
+        // Iteration over A's rows
+        for(int i = br_max; i < m; i++) {
+            // Iteration over B's columns
+            for(int k = 0; k < l; k++)
+                C[i*l + k] += A[i * n + j] * B[j * l + k];
+        }
+
+        // Iteration over A's rows
+        for (int i = 0; i < br_max; i++)
+        {
+            // Iteration over B's columns
+            for(int k = bc_max; k < l; k++)
+                C[i*l + k] += A[i * n + j] * B[j * l + k];
+        }
+    }
+}
+
+/* Function: lina_dot4
+**
+**   Evaluates the dot product C = A * B. The A,B
+**   matrices are, respectively, mxn and nxl, which
+**   means C is mxl. The resulting C matrix is stored
+**   in a memory  region specified by the caller. 
+**
+** Variant 4 of lina_dot:
+**   This include the changes of lina_dot3 but uses the
+**   micro kernel subroutine
+**
+** Notes:
+**
+**   - A,B must be provided as contiguous memory regions
+**     represented in row-major order. Also, C is stored
+**     that way too.
+**
+**   - The C pointer CAN'T refer to the same memory region 
+**     of either A or B.
+**
+**   - m,n,l must be greater than 0.
+**
+**   - This function can never fail.
+*/
+void lina_dot4(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);
+    // A_sub, B_sub and C_sub must be 32 byte aligned
+    assert(!((uintptr_t)A & 31llu) && !((uintptr_t)B & 31llu) && !((uintptr_t)C & 31llu));
+
+    #define KERNEL_ROW 6
+    #define KERNEL_COLS 8
+
+    const int br_max = (m & ~(KERNEL_ROW - 1));
+    const int bc_max = (l & ~(KERNEL_COLS - 1));
+
+    for (int br = 0; br < br_max; br += KERNEL_ROW)
+    {
+        for (int bc = 0; bc < bc_max; bc += KERNEL_COLS)
+        {
+            dot_kernel_6x8(A, B, C, br, bc, 0, n, n, l);
+        }
+    }
+}
+
+/*
+*   
+*   Computes C_sub += A_sub * B_sub where:
+*       - C_sub = C[x:x+6][y:y+8]
+*       - A_sub = A[x:x+6][c_min:c_max]
+*       - B_sub = B[c_min:c_max][y:y+8]
+*       - n is the number of columns of A
+*       - l the number of columns of B
+*/
+static void
+dot_kernel_6x8(double *A_sub, double *B_sub, double *C_sub, int x, int y, int c_min, int c_max, int n, int l)
+{
+    // A_sub, B_sub and C_sub must be 32 byte aligned
+    // assert is done in the main lina_dot function
+    //assert(!((uintptr_t)A_sub & 31llu) && !((uintptr_t)B_sub & 31llu) && !((uintptr_t)C_sub & 31llu));
+
+    // This structure should use 12 YMM registers
+
+    __m256d *Bm_sub = (__m256d *)B_sub;
+    __m256d *Cm_sub = (__m256d *)C_sub;
+    __m256d acc[6][2] = {0};
+
+
+    for (int k = c_min; k < c_max; k++)
+    {
+        for (int i = 0; i < 6; i++)
+        {
+            __m256d A_brdcst = _mm256_broadcast_sd(&A_sub[(x + i)*n + k]);
+            for (int j = 0; j < 2; j++)
+                acc[i][j] = _mm256_fmadd_pd(A_brdcst,Bm_sub[(k*l + y)/4 + j],acc[i][j]);
+        }
+    }
+
+    for (int i = 0; i < 6; i++)
+        for (int j = 0; j < 2; j++)
+            Cm_sub[((x+i)*l + y)/4 + j] = acc[i][j];
+}
+
 /* Function: lina_add
 **
 **   Evaluates the matrix addition C = A + B. The result
@@ -639,406 +1005,7 @@ void lina_conv(double *A, double *B, double *C,
         }
 }
 
-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;
+    // To be done
 }

src/lina.h

@@ -1,18 +1,15 @@
-#include <complex.h>
 #include <stdbool.h>
+#include <stdio.h>
 
 void lina_dot(double *A, double *B, double *C, int m, int n, int l);
+void lina_dot1(double *A, double *B, double *C, int m, int n, int l);
+void lina_dot2(double *A, double *B, double *C, int m, int n, int l);
+void lina_dot3(double *A, double *B, double *C, int m, int n, int l);
+void lina_dot4(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_conv(double *A, double *B, double *C, int Aw, int Ah, int Bw, int Bh);
 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

@@ -1,142 +0,0 @@
-#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

@@ -1,135 +0,0 @@
-#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;
-}

tests/README.md

@@ -1,16 +0,0 @@
-## Description
-Here is developed the testing unit for all the lina functions that need numerical testing.
-
-## Usage
-For each function in the lina library named in the form of _lina_something()_, here is defined a folder named _something_. In each folder there are many tests, each one identified by a ti.txt file, for i=1,...,n.
-Each test file is defined as follows: The first matrix/matrices are the inputs of the function (depending on the function, for example: lina_add() has two inputs A,B and one output C=A+B. A and B have to be the first two matrices in the test file), the last matrix/matrices are the output of the function and after there are input scalar values (ordered in the same order of the function under test) of the function represented as a 1x1 matrix.
-
-For example, a scale test file, that is a test for the lina_scale() function is defined as follows:
-
-[1 1 1,1 1 1,1 1 1]    
-[2 2 2,2 2 2,2 2 2]   
-[2]
-
-Where the first matrix is the input, the second the output and the last is the scalar value.
-
-By default, executing the test file will generate all the testing and provide the results on the stdout.

tests/dot/t1.txt

@@ -1,11 +0,0 @@
-[1 0 0, 
- 0 1 0, 
- 0 0 1]
- 
-[1 0 0,
- 0 1 0,
- 0 0 1]
- 
-[1 0 0,
- 0 1 0,
- 0 0 1]

tests/scale/t1.txt

@@ -1,9 +0,0 @@
-[1 1 1, 
- 1 1 1, 
- 1 1 1]
- 
-[2 2 2,
- 2 2 2,
- 2 2 2]
- 
-[2]

tests/test.c

@@ -1,553 +0,0 @@
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <assert.h>
-#include <dirent.h>
-#include <sys/types.h>
-#include <errno.h>
-#include "lina.h"
-
-#define check assert
-
-//Print the matrix A with size m by n
-static void pmatrix(FILE *fp, double *A, int m, int n);
-
-
-typedef struct dot_test{
-
-	double *A;
-	double *B;
-	double *C;
-	int m;
-	int n;
-	int l;
-
-}dot_test;
-
-typedef struct add_test{
-
-	double *A;
-	double *B;
-	double *C;
-	int m;
-	int n;
-
-}add_test;
-
-typedef struct scale_test{
-
-	double *A;
-	double *B;
-	double s;
-	int m;
-	int n;
-
-}scale_test;
-
-typedef struct transpose_test{
-
-	double *A;
-	double *B;
-	int m;
-	int n;
-
-}transpose_test;
-
-#define PATH "./tests/"
-
-int main()
-{	
-	//Defining pointers to the test structures
-	add_test *add_tests;
-	dot_test *dot_tests;
-	scale_test *scale_tests;
-	transpose_test *transpose_tests;
-
-	//Number of tests for each lina functions
-	int n_dot_tests, n_add_tests, n_scale_tests, n_transpose_tests;
-	
-	//Opening dir stream
-	DIR *dir = opendir(PATH);
-	struct dirent *ep;
-	check(dir != NULL);
-
-	//Loading all the tests from files
-	while (ep = readdir(dir))
-	{
-		if(ep->d_type != DT_DIR || !strcmp(ep->d_name, ".") || !strcmp(ep->d_name, ".."))
-			continue;
-		
-		if(!strcmp(ep->d_name,"add"))
-		{
-
-			//Implement loading lina_add test matrices
-			char sub_path[256] = PATH;
-			strcat(sub_path,ep->d_name);
-			strcat(sub_path,"/");
-			
-			DIR *sub_dir = opendir(sub_path);
-			check(sub_dir != NULL);
-
-			struct dirent *sub_ep;
-
-			unsigned int count = 0;
-
-			while (sub_ep = readdir(sub_dir))
-			{
-				if(sub_ep->d_type == DT_DIR)
-					continue;
-				count += 1;
-			}
-
-			closedir(sub_dir);
-
-			sub_dir = opendir(sub_path);
-			check(sub_dir != NULL);
-			
-			add_tests = malloc(sizeof(add_test)*count);
-			n_add_tests = count;
-			int i = 0;
-
-			while (sub_ep = readdir(sub_dir))
-			{
-				if(sub_ep->d_type == DT_DIR)
-					continue;
-
-				char file_pos[256];
-				strcat(file_pos,sub_path);
-				strcat(file_pos,sub_ep->d_name);
-
-				FILE *fp;
-				fp = fopen(file_pos,"r");
-				check(fp != NULL);
-
-				int m,n;
-				char *error;
-
-				add_tests[i].A = lina_loadMatrixFromStream(fp,&n,&m,&error);
-				check(add_tests[i].A  != NULL);
-
-				add_tests[i].B = lina_loadMatrixFromStream(fp,&n,&m,&error);
-				check(add_tests[i].B  != NULL);
-
-				add_tests[i].C = lina_loadMatrixFromStream(fp,&n,&m,&error);
-				check(add_tests[i].C  != NULL);
-
-				add_tests[i].m = m;
-				add_tests[i].n = n;
-
-				i += 1;
-				fclose(fp);
-
-			}
-
-			closedir(sub_dir);
-
-		}
-		else if (!strcmp(ep->d_name,"dot"))
-		{
-			//Implement loading lina_dot test matrices
-			char sub_path[256] = PATH;
-			strcat(sub_path,ep->d_name);
-			strcat(sub_path,"/");
-			
-			DIR *sub_dir = opendir(sub_path);
-			check(sub_dir != NULL);
-
-			struct dirent *sub_ep;
-
-			unsigned int count = 0;
-
-			while (sub_ep = readdir(sub_dir))
-			{
-				if(sub_ep->d_type == DT_DIR)
-					continue;
-				count += 1;
-			}
-
-			closedir(sub_dir);
-
-			sub_dir = opendir(sub_path);
-			check(sub_dir != NULL);
-			dot_tests = malloc(sizeof(dot_test)*count);
-			n_dot_tests = count;
-
-			int i = 0;
-
-			while (sub_ep = readdir(sub_dir))
-			{
-				if(sub_ep->d_type == DT_DIR)
-					continue;
-
-				char file_pos[256];
-				strcpy(file_pos,sub_path);
-				strcat(file_pos,sub_ep->d_name);
-
-				FILE *fp;
-				fp = fopen(file_pos,"r");
-				check(fp != NULL);
-
-				int m,n,l;
-				char *error;
-
-				dot_tests[i].A = lina_loadMatrixFromStream(fp,&n,&m,&error);
-				check(dot_tests[i].A  != NULL);
-
-				dot_tests[i].B = lina_loadMatrixFromStream(fp,&l,&n,&error);
-				check(dot_tests[i].B  != NULL);
-
-				dot_tests[i].C = lina_loadMatrixFromStream(fp,&l,&m,&error);
-				check(dot_tests[i].C  != NULL);
-
-				dot_tests[i].m = m;
-				dot_tests[i].n = n;
-				dot_tests[i].l = l;
-
-				i += 1;
-				fclose(fp);
-
-			}
-
-			closedir(sub_dir);
-
-
-		}
-		else if (!strcmp(ep->d_name,"scale"))
-		{
-			//Implement loading lina_scale test matrices
-			char sub_path[256] = PATH;
-			strcat(sub_path,ep->d_name);
-			strcat(sub_path,"/");
-			
-			DIR *sub_dir = opendir(sub_path);
-			check(sub_dir != NULL);
-
-			struct dirent *sub_ep;
-
-			unsigned int count = 0;
-
-			while (sub_ep = readdir(sub_dir))
-			{
-				if(sub_ep->d_type == DT_DIR)
-					continue;
-				count += 1;
-			}
-
-			closedir(sub_dir);
-
-			sub_dir = opendir(sub_path);
-			check(sub_dir != NULL);
-
-			scale_tests = malloc(sizeof(scale_test)*count);
-			
-			n_scale_tests = count;
-
-			int i = 0;
-
-			while (sub_ep = readdir(sub_dir))
-			{
-				if(sub_ep->d_type == DT_DIR)
-					continue;
-
-				char file_pos[256];
-				strcpy(file_pos,sub_path);
-
-				strcat(file_pos,sub_ep->d_name);
-
-				FILE *fp;
-				fp = fopen(file_pos,"r");
-				check(fp != NULL);
-
-				int m,n;
-				int useless1,useless2;
-				double *scale;
-				char *error;
-
-				scale_tests[i].A = lina_loadMatrixFromStream(fp,&n,&m,&error);
-				check(scale_tests[i].A  != NULL);
-
-				scale_tests[i].B = lina_loadMatrixFromStream(fp,&n,&m,&error);
-				check(scale_tests[i].B  != NULL);
-
-				scale = lina_loadMatrixFromStream(fp,&useless1,&useless2,&error);
-				check(scale != NULL);
-
-				scale_tests[i].m = m;
-				scale_tests[i].n = n;
-				scale_tests[i].s = scale[0];
-				free(scale);
-
-
-				i += 1;
-				fclose(fp);
-
-			}
-
-			closedir(sub_dir);
-
-
-		}
-		else if (!strcmp(ep->d_name,"transpose"))
-		{
-			//Implement loading lina_transpose test matrices
-			char sub_path[256] = PATH;
-			strcat(sub_path,ep->d_name);
-			strcat(sub_path,"/");
-			
-			DIR *sub_dir = opendir(sub_path);
-			check(sub_dir != NULL);
-
-			struct dirent *sub_ep;
-
-			unsigned int count = 0;
-
-			while (sub_ep = readdir(sub_dir))
-			{
-				if(sub_ep->d_type == DT_DIR)
-					continue;
-				count += 1;
-			}
-
-			closedir(sub_dir);
-
-			sub_dir = opendir(sub_path);
-			check(sub_dir != NULL);
-			
-			transpose_tests = malloc(sizeof(transpose_test)*count);
-			n_transpose_tests = count;
-
-			int i = 0;
-
-			while (sub_ep = readdir(sub_dir))
-			{
-				if(sub_ep->d_type == DT_DIR)
-					continue;
-
-				char file_pos[256];
-				strcpy(file_pos,sub_path);
-				strcat(file_pos,sub_ep->d_name);
-
-				FILE *fp;
-				fp = fopen(file_pos,"r");
-				check(fp != NULL);
-
-				int m,n;
-				char *error;
-				
-				transpose_tests[i].A = lina_loadMatrixFromStream(fp,&n,&m,&error);
-				check(transpose_tests[i].A  != NULL);
-
-				transpose_tests[i].B = lina_loadMatrixFromStream(fp,&m,&n,&error);
-				check(transpose_tests[i].B  != NULL);
-
-				transpose_tests[i].m = m;
-				transpose_tests[i].n = n;
-
-				i += 1;
-				fclose(fp);
-
-			}
-
-			closedir(sub_dir);
-
-
-		}
-
-	}
-	
-	closedir(dir);
-
-	//Starting the lina_add tests
-	{
-		int passed_tests = 0;
-		fprintf(stdout,"Starting tests on lina_add():\n");
-		for(int i=0;i<n_add_tests;i++){
-			
-			double *C = (double*) malloc(sizeof(*C)*add_tests[i].m * add_tests[i].n);
-			check(C != NULL);
-
-			lina_add(add_tests[i].A, add_tests[i].B, C,add_tests[i].m,add_tests[i].n);
-			
-			if( !memcmp(add_tests[i].C, C, sizeof(*C)*add_tests[i].m * add_tests[i].n) )
-				passed_tests += 1;
-			else{
-				
-				fprintf(stderr,"----------------------------------------------------\n");
-				fprintf(stderr,"Test on lina_add() failed on the following matrices:\n");
-				pmatrix(stderr,add_tests[i].A, add_tests[i].m, add_tests[i].n);
-				fprintf(stderr,"+\n");
-				pmatrix(stderr,add_tests[i].B, add_tests[i].m, add_tests[i].n);
-				fprintf(stderr,"lina_add() gives following output:\n");
-				pmatrix(stderr,C, add_tests[i].m, add_tests[i].n);
-				fprintf(stderr,"instead of:\n");
-				pmatrix(stderr,add_tests[i].C, add_tests[i].m, add_tests[i].n);
-				fprintf(stderr,"----------------------------------------------------\n");
-
-
-			}
-			free(C);
-
-		}
-		if(n_add_tests != 0)
-			fprintf(stdout, "Test on lina_add() finished: %d out of %d tests were succesfull\n",passed_tests,n_add_tests);
-		else
-			fprintf(stdout, "There are no tests for lina_add() function.\n");
-	}
-
-	//Starting the lina_dot tests
-	{
-		int passed_tests = 0;
-		fprintf(stdout,"\nStarting tests on lina_dot():\n");
-		for(int i=0;i<n_dot_tests;i++){
-			
-			double *C = (double*) malloc(sizeof(*C)*dot_tests[i].m * dot_tests[i].l);
-			check(C != NULL);
-
-			lina_dot(dot_tests[i].A, dot_tests[i].B, C, dot_tests[i].m, dot_tests[i].n, dot_tests[i].l);
-			
-			if( !memcmp(dot_tests[i].C, C, sizeof(*C)*dot_tests[i].m * dot_tests[i].l) )
-				passed_tests += 1;
-			else{
-				
-				fprintf(stderr,"----------------------------------------------------\n");
-				fprintf(stderr,"Test on lina_dot() failed on the following matrices:\n");
-				pmatrix(stderr,dot_tests[i].A, dot_tests[i].m, dot_tests[i].n);
-				fprintf(stderr,"*\n");
-				pmatrix(stderr,dot_tests[i].B, dot_tests[i].n, dot_tests[i].l);
-				fprintf(stderr,"lina_dot() gives following output:\n");
-				pmatrix(stderr,C, dot_tests[i].m, dot_tests[i].l);
-				fprintf(stderr,"instead of:\n");
-				pmatrix(stderr,dot_tests[i].C, dot_tests[i].m, dot_tests[i].l);
-				fprintf(stderr,"----------------------------------------------------\n");
-
-
-			}
-			free(C);
-
-		}
-		if(n_dot_tests != 0)
-			fprintf(stdout, "Test on lina_dot() finished: %d out of %d tests were succesfull\n",passed_tests,n_dot_tests);
-		else
-			fprintf(stdout, "There are no tests for lina_dot() function.\n");
-	}
-
-	//Starting the lina_transpose tests
-	{
-		int passed_tests = 0;
-		fprintf(stdout,"\nStarting tests on lina_transpose():\n");
-		for(int i=0;i<n_transpose_tests;i++){
-			
-			double *C = (double*) malloc(sizeof(*C)*transpose_tests[i].m * transpose_tests[i].n);
-			check(C != NULL);
-
-			lina_transpose(transpose_tests[i].A, C, transpose_tests[i].m, transpose_tests[i].n);
-			
-			if( !memcmp(transpose_tests[i].B, C, sizeof(*C)*transpose_tests[i].m * transpose_tests[i].n) )
-				passed_tests += 1;
-			else{
-				
-				fprintf(stderr,"----------------------------------------------------\n");
-				fprintf(stderr,"Test on lina_transpose() failed on the following matrices:\n");
-				pmatrix(stderr,transpose_tests[i].A, transpose_tests[i].m, transpose_tests[i].n);
-				fprintf(stderr,"lina_transpose() gives following output:\n");
-				pmatrix(stderr,C, transpose_tests[i].n, transpose_tests[i].m);
-				fprintf(stderr,"instead of:\n");
-				pmatrix(stderr,transpose_tests[i].B, transpose_tests[i].n, transpose_tests[i].m);
-				fprintf(stderr,"----------------------------------------------------\n");
-
-
-			}
-			free(C);
-
-		}
-		if(n_transpose_tests != 0)
-			fprintf(stdout, "Test on lina_transpose() finished: %d out of %d tests were succesfull\n",passed_tests,n_transpose_tests);
-		else
-			fprintf(stdout, "There are no tests for lina_transpose() function.\n");
-	}
-
-	//Starting the lina_scale tests
-	{
-		int passed_tests = 0;
-		fprintf(stdout,"\nStarting tests on lina_scale():\n");
-		for(int i=0;i<n_scale_tests;i++){
-			
-			double *C = (double*) malloc(sizeof(*C)*scale_tests[i].m * scale_tests[i].n);
-			check(C != NULL);
-
-			lina_scale(scale_tests[i].A, C, scale_tests[i].s, scale_tests[i].m, scale_tests[i].n);
-			
-			if( !memcmp(scale_tests[i].B, C, sizeof(*C)*scale_tests[i].m * scale_tests[i].n) )
-				passed_tests += 1;
-			else{
-				
-				fprintf(stderr,"----------------------------------------------------\n");
-				fprintf(stderr,"Test on lina_scale() failed on the following matrices:\n");
-				pmatrix(stderr,scale_tests[i].A, scale_tests[i].m, scale_tests[i].n);
-				fprintf(stderr,"lina_scale() gives following output:\n");
-				pmatrix(stderr, C, scale_tests[i].m, scale_tests[i].n);
-				fprintf(stderr,"instead of:\n");
-				pmatrix(stderr,scale_tests[i].B, scale_tests[i].m, scale_tests[i].n);
-				fprintf(stderr,"----------------------------------------------------\n");
-
-
-			}
-			free(C);
-
-		}
-		if(n_scale_tests != 0)
-			fprintf(stdout, "Test on lina_scale() finished: %d out of %d tests were succesfull\n",passed_tests,n_scale_tests);
-		else
-			fprintf(stdout, "There are no tests for lina_scale() function.\n");
-	}
-
-	//Freeing the memory of all the heap variables
-	{
-		//Freeing add_tests memory
-		for(int i=0;i< n_add_tests;i++)
-		{
-			free(add_tests[i].A);
-			free(add_tests[i].B);
-			free(add_tests[i].C);
-		}
-
-		//Freeing dot_tests memory
-		for(int i=0;i< n_dot_tests;i++)
-		{
-			free(dot_tests[i].A);
-			free(dot_tests[i].B);
-			free(dot_tests[i].C);
-		}
-
-		//Freeing scale_tests memory
-		for(int i=0;i< n_scale_tests;i++)
-		{
-			free(scale_tests[i].A);
-			free(scale_tests[i].B);
-		}
-
-		//Freeing transpose_tests memory
-		for(int i=0;i< n_transpose_tests;i++)
-		{
-			free(transpose_tests[i].A);
-			free(transpose_tests[i].B);
-		}
-	}
-
-	return 0;
-}
-
-static void pmatrix(FILE *fp, double *A,int m,int n){
-
-	for(int i = 0; i<m; i++){
-		fprintf(fp, "    | ");
-		for(int j = 0; j< n; j++)
-			fprintf(fp, "%g ",A[i*n + j]);
-		
-		fprintf(fp, "|\n");
-
-	}
-
-	fprintf(fp, "\n");
-
-}

tests/transpose/t1.txt

@@ -1,9 +0,0 @@
-[1 0 0 0, 
- 0 2 0 0,
- 0 0 3 0, 
- 0 0 0 4]
- 
-[1 0 0 0, 
- 0 2 0 0,
- 0 0 3 0, 
- 0 0 0 4]

tests/transpose/t10.txt

@@ -1,6 +0,0 @@
-[1 2, 
- 3 4,
- 5 6]
- 
-[1 3 5, 
- 2 4 6]

tests/transpose/t2.txt

@@ -1,9 +0,0 @@
-[1 2 3 4, 
- 0 0 0 0, 
- 0 0 0 0,
- 0 0 0 0]
- 
-[1 0 0 0, 
- 2 0 0 0, 
- 3 0 0 0,
- 4 0 0 0]

tests/transpose/t3.txt

@@ -1,9 +0,0 @@
-[0 1 0 0, 
- 0 2 0 0, 
- 0 3 0 0,
- 0 4 0 0]
- 
-[0 0 0 0, 
- 1 2 3 4, 
- 0 0 0 0,
- 0 0 0 0]

tests/transpose/t4.txt

@@ -1,9 +0,0 @@
-[0 0 1 0, 
- 0 0 2 0, 
- 0 0 3 0,
- 0 0 4 0]
- 
-[0 0 0 0, 
- 0 0 0 0, 
- 1 2 3 4,
- 0 0 0 0]

tests/transpose/t5.txt

@@ -1,9 +0,0 @@
-[0 0 0 0, 
- 0 0 0 0, 
- 0 0 0 0,
- 1 2 3 4]
- 
-[0 0 0 1, 
- 0 0 0 2, 
- 0 0 0 3,
- 0 0 0 4]

tests/transpose/t6.txt

@@ -1,7 +0,0 @@
-[1 0 0, 
- 0 2 0, 
- 0 0 3]
- 
-[1 0 0, 
- 0 2 0,
- 0 0 3]

tests/transpose/t7.txt

@@ -1,7 +0,0 @@
-[1 2 3, 
- 0 0 0, 
- 0 0 0]
- 
-[1 0 0, 
- 2 0 0,
- 3 0 0]

tests/transpose/t8.txt

@@ -1,7 +0,0 @@
-[0 0 1, 
- 0 0 2, 
- 0 0 3]
- 
-[0 0 0, 
- 0 0 0,
- 1 2 3]

tests/transpose/t9.txt

@@ -1,6 +0,0 @@
-[1 2 3, 
- 4 5 6]
- 
-[1 4, 
- 2 5,
- 3 6]