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
-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!
+
+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 -o test  -Wall -Wextra -g -Isrc/
-gcc tests/test_loader.c src/lina.c -o test_loader -Wall -Wextra -g -Isrc/

src/lina.c

@@ -1,11 +1,14 @@
-#include <stddef.h>
 #include <assert.h>
 #include <stdlib.h>
 #include <string.h>
-#include <stdio.h>
 #include <errno.h>
 #include <ctype.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
 **
@@ -27,30 +30,393 @@
 **
 **   - 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 sum = 0;
+
+            // Iteration over the single B column 
+            // for executing the sum of product
+                    
+            for(int j=0; j < n; j++)
+                sum += A[i * n + j] * B[j * l + k];
+
+            C[i * l + k] = sum;
+        }
+    }
+}
+
+/* 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++)
     {
-                    double pos = 0;
+        // 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++)
-
-                        pos += A[i*n + j] * B[j*l + k];
-
-                    C[i*l + k] = pos;
+            {
+                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
@@ -72,8 +438,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 +488,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 +497,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 +510,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 +537,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 +595,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 +606,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 +616,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 +630,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 +650,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 +722,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 +744,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 +754,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 +764,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 +790,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 +815,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 +831,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 +868,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 +878,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 +890,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 +905,107 @@ 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];
+        }
+}
+
+bool lina_inverse(double *M, double *D, int n)
+{
+    // To be done
+}

src/lina.h

@@ -1,6 +1,15 @@
+#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);
 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);
 double *lina_loadMatrixFromStream(FILE *fp, int *width, int *height, char **error);
+int     lina_saveMatrixToStream(FILE *fp, double *A, int width, int height, char **error);

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]