better transposition algorithm

.gitignore

@@ -1,3 +1,4 @@
 .swp
 test
+time
 .vscode

build.sh

@@ -1 +1,2 @@
 gcc test.c lina.c -o test -Wall -Wextra -g
+gcc time.c lina.c -o time -Wall -Wextra -O3

lina.c

@@ -12,24 +12,25 @@ void lina_dot(double *A, double *B, double *C, int m, int n, int l){
     assert(A != NULL && B != NULL && C != NULL);
     assert(A != C && B != C);
 
-    //Actual dot
+    // Iteration over A's rows
-    //Iteration over A's rows
+    for(int i = 0; i < m; i++)
-    for(int i=0; i < m;i++){
+        {
-
+            // Iteration over B's columns
-        //Iteration over B's columns
+            for(int k = 0; k < l; k++)
-        for(int k=0; k < l; k++){
+                {
-
                     double pos = 0;
 
-            //Iteration over the single B column for executing the product of sum
+                    // Iteration over the single B column 
+                    // for executing the product of sum
+                    
                     for(int j=0; j < n; j++)
-                pos += A[i * n + j] * B[j * l + k];
 
-            //The usage of a support variable 'pos' is for safety purpose (non-zero C matrix)
+                        pos += A[i*n + j] * B[j*l + k];
-            C[i*l+k] = pos;
 
+                    // The usage of a support variable 'pos'
+                    // is for safety purpose (non-zero C matrix)
+                    C[i*l + k] = pos;
                 }
-        
         }
 
 }
@@ -40,10 +41,9 @@ void lina_add(double *A, double *B, double *C, int m, int n){
     assert(A != NULL && B != NULL && C != NULL);
     
     for(int i = 0; i < m; i++)
-
         for(int j = 0; j < n; j++)
 
-            C[i * n + j] = A[i * n + j] + B[i * n + j];
+            C[i*n + j] = A[i*n + j] + B[i*n + j];
 }
 
 void lina_scale(double *A, double *B, double k, int m, int n){
@@ -55,27 +55,46 @@ void lina_scale(double *A, double *B, double k, int m, int n){
 
         for(int j = 0; j < n; j++)
 
-            B[i * n + j] = k * A[i * n + j];
+            B[i*n + j] = k * A[i*n + j];
 }
 
-void lina_transpose(double *A, double *B, int m, int n){
+void lina_transpose(double *A, double *B, int m, int n)
-
+{
     assert(m > 0 && n > 0);
     assert(A != NULL && B != NULL);
 
-#warning "Try to think better algorithm"
+    if(m == 1 || n == 1)
+        {
+            memcpy(B, A, sizeof(A[0]) * m * n);
+            return;
+        }
 
-    double *support = malloc(sizeof(*support) * m * n);
+    if(m == n)
+        {
+            for(int i = 0; i < n; i += 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
+        {
+            B[0] = A[0];
+            B[m*n-1] = A[m*n-1];
 
-    check(support != NULL);
+            double item = A[1];
+            int    next = m;
 
-    memcpy(support, A, sizeof(*support) * m * n);
+            while(next != 1)
+                {
+                    double temp = A[next];
+                    B[next] = item;
+                    item = temp;
+                    next = (next % n) * m + (next / n);
+                }
 
-    for(int i = 0; i < n; i++)
+            B[1] = item;
-        
+        }
-        for(int j = 0; j < m; j++)
-
-            B[j*n + i] = support[i*m + j];
-
-    free(support);
 }

test.c

@@ -1,8 +1,11 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
+#include <assert.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);
 
@@ -51,57 +54,91 @@ struct {
 };
 
 struct {
-	double A[9], B[9];
+	double *A, *B;
+	int m, n;
 } transp_tests[] = {
 	{
-		.A = {
+		.A = (double[]) {
 			1, 0, 0,
 			0, 2, 0,
 			0, 0, 3,
 		},
-		.B = {
+		.B = (double[]) {
 			1, 0, 0,
 			0, 2, 0,
 			0, 0, 3,
 		},
+		.m = 3,
+		.n = 3,
 	},
 	{
-		.A = {
+		.A = (double[]) {
 			1, 2, 3,
 			0, 0, 0,
 			0, 0, 0,
 		},
-		.B = {
+		.B = (double[]) {
 			1, 0, 0,
 			2, 0, 0,
 			3, 0, 0,
 		},
+		.m = 3,
+		.n = 3,
 	},
 	{
-		.A = {
+		.A = (double[]) {
 			0, 1, 0,
 			0, 2, 0,
 			0, 3, 0,
 		},
-		.B = {
+		.B = (double[]) {
 			0, 0, 0,
 			1, 2, 3,
 			0, 0, 0,
 		},
+		.m = 3,
+		.n = 3,
 	},
 	{
-		.A = {
+		.A = (double[]) {
 			0, 0, 1,
 			0, 0, 2,
 			0, 0, 3,
 		},
-		.B = {
+		.B = (double[]) {
 			0, 0, 0,
 			0, 0, 0,
 			1, 2, 3,
 		},
+		.m = 3,
+		.n = 3,
+	},
+	{
+		.A = (double[]) {
+			1, 2, 3,
+			4, 5, 6,
+		},
+		.B = (double[]) {
+			1, 4,
+			2, 5,
+			3, 6,
+		},
+		.m = 2,
+		.n = 3,
+	},
+	{
+		.A = (double[]) {
+			1, 4,
+			2, 5,
+			3, 6,
+		},
+		.B = (double[]) {
+			1, 2, 3,
+			4, 5, 6,
+		},
+		.m = 3,
+		.n = 2,
 	},
-
 };
 
 int main()
@@ -168,11 +205,16 @@ int main()
 
 		for(int i = 0; i < transp_total; i += 1)
 			{
-				double R[9];
+				int m = transp_tests[i].m;
+				int n = transp_tests[i].n;
 
-				lina_transpose(transp_tests[i].A, R, 3, 3);
+				double R[32];
 
-				if(!memcmp(R, transp_tests[i].B, sizeof(R)))
+				assert(sizeof(R) >= m * n * sizeof(R[0]));
+
+				lina_transpose(transp_tests[i].A, R, m, n);
+
+				if(!memcmp(R, transp_tests[i].B, sizeof(R[0]) * m * n))
 					{
 						fprintf(stderr, "Transposition test %d passed.\n", i);
 						transp_passed += 1;
@@ -180,9 +222,9 @@ int main()
 				else
 					{
 						fprintf(stderr, "Transposition test %d failed:\n  got matrix:\n\n", i);
-						pmatrix(stderr, R, 3, 3);
+						pmatrix(stderr, R, m, n);
 						fprintf(stderr, "  instead of:\n\n");
-						pmatrix(stderr, transp_tests[i].B, 3, 3);
+						pmatrix(stderr, transp_tests[i].B, m, n);
 					}
 			}
 		fprintf(stderr, "\n\t%d transpositions out of %d were succesful.\n\n", transp_passed, transp_total);

time.c

@@ -0,0 +1,62 @@
+#include <time.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <assert.h>
+#include <string.h>
+#include "lina.h"
+
+/* This function compares the lina_transpose
+** implementation against the naive implementation.
+*/
+
+#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 = malloc(sizeof(*support) * m * n);
+
+    check(support != NULL);
+
+    memcpy(support, A, sizeof(*support) * m * n);
+
+    for(int i = 0; i < n; i++)
+        
+        for(int j = 0; j < m; j++)
+
+            B[j*n + i] = support[i*m + j];
+
+    free(support);
+}
+
+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);
+
+	double t1, t2;
+	clock_t begin, end;
+
+	begin = clock();
+	lina_transpose(big, big, m, n);
+	end = clock();
+	t1 = (double) (end - begin) / CLOCKS_PER_SEC;
+
+	begin = clock();
+	naive_transpose(big, big, m, n);
+	end = clock();
+	t2 = (double) (end - begin) / CLOCKS_PER_SEC;
+
+	printf("lina_transpose took %gms\n", t1*1000);
+	printf("naive_transpose took %gms\n", t2*1000);
+
+	free(big);
+	return 0;
+}