Matrix multiplication is now ready

6 files changed, 523 insertions, 247 deletions

diff --git a/src/matrixOperations/Makefile.am b/src/matrixOperations/Makefile.am
index 1cc9266c..7a0efe1f 100644
--- a/src/matrixOperations/Makefile.am
+++ b/src/matrixOperations/Makefile.am
@@ -43,7 +43,7 @@ testMatrixSubtraction_SOURCES = testMatrixSubtraction.c
 testMatrixSubtraction_LDADD = $(check_LDADD)
 testMatrixSubtraction_CFLAGS = $(check_INCLUDES)
 
-testMatrixMultiplication_SOURCES = testMatrixMultiplication.c
+testMatrixMultiplication_SOURCES = testFloatMatrixMultiplication.c testDoubleMatrixMultiplication.c
 testMatrixMultiplication_LDADD = $(check_LDADD)
 testMatrixMultiplication_CFLAGS = $(check_INCLUDES)
 
diff --git a/src/matrixOperations/Makefile.in b/src/matrixOperations/Makefile.in
index 73f3b5f6..911cdd4a 100644
--- a/src/matrixOperations/Makefile.in
+++ b/src/matrixOperations/Makefile.in
@@ -76,8 +76,8 @@ testMatrixAddition_LINK = $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) \
 	$(LIBTOOLFLAGS) --mode=link $(CCLD) \
 	$(testMatrixAddition_CFLAGS) $(CFLAGS) $(AM_LDFLAGS) \
 	$(LDFLAGS) -o $@
-am_testMatrixMultiplication_OBJECTS =  \
-	testMatrixMultiplication-testMatrixMultiplication.$(OBJEXT)
+am_testMatrixMultiplication_OBJECTS = testMatrixMultiplication-testFloatMatrixMultiplication.$(OBJEXT) \
+	testMatrixMultiplication-testDoubleMatrixMultiplication.$(OBJEXT)
 testMatrixMultiplication_OBJECTS =  \
 	$(am_testMatrixMultiplication_OBJECTS)
 testMatrixMultiplication_DEPENDENCIES = $(check_LDADD)
@@ -250,7 +250,7 @@ testMatrixAddition_CFLAGS = $(check_INCLUDES)
 testMatrixSubtraction_SOURCES = testMatrixSubtraction.c
 testMatrixSubtraction_LDADD = $(check_LDADD)
 testMatrixSubtraction_CFLAGS = $(check_INCLUDES)
-testMatrixMultiplication_SOURCES = testMatrixMultiplication.c
+testMatrixMultiplication_SOURCES = testFloatMatrixMultiplication.c testDoubleMatrixMultiplication.c
 testMatrixMultiplication_LDADD = $(check_LDADD)
 testMatrixMultiplication_CFLAGS = $(check_INCLUDES)
 ADDITIONHEAD = matrixAddition.h
@@ -348,7 +348,8 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libMatrixOperations_la-matrixMultiplication.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libMatrixOperations_la-matrixSubtraction.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/testMatrixAddition-testMatrixAddition.Po@am__quote@
-@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/testMatrixMultiplication-testMatrixMultiplication.Po@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/testMatrixMultiplication-testDoubleMatrixMultiplication.Po@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/testMatrixMultiplication-testFloatMatrixMultiplication.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/testMatrixSubtraction-testMatrixSubtraction.Po@am__quote@
 
 .c.o:
@@ -407,19 +408,33 @@ testMatrixAddition-testMatrixAddition.obj: testMatrixAddition.c
 @AMDEP_TRUE@@am__fastdepCC_FALSE@	DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCC_FALSE@	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(testMatrixAddition_CFLAGS) $(CFLAGS) -c -o testMatrixAddition-testMatrixAddition.obj `if test -f 'testMatrixAddition.c'; then $(CYGPATH_W) 'testMatrixAddition.c'; else $(CYGPATH_W) '$(srcdir)/testMatrixAddition.c'; fi`
 
-testMatrixMultiplication-testMatrixMultiplication.o: testMatrixMultiplication.c
-@am__fastdepCC_TRUE@	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(testMatrixMultiplication_CFLAGS) $(CFLAGS) -MT testMatrixMultiplication-testMatrixMultiplication.o -MD -MP -MF $(DEPDIR)/testMatrixMultiplication-testMatrixMultiplication.Tpo -c -o testMatrixMultiplication-testMatrixMultiplication.o `test -f 'testMatrixMultiplication.c' || echo '$(srcdir)/'`testMatrixMultiplication.c
-@am__fastdepCC_TRUE@	mv -f $(DEPDIR)/testMatrixMultiplication-testMatrixMultiplication.Tpo $(DEPDIR)/testMatrixMultiplication-testMatrixMultiplication.Po
-@AMDEP_TRUE@@am__fastdepCC_FALSE@	source='testMatrixMultiplication.c' object='testMatrixMultiplication-testMatrixMultiplication.o' libtool=no @AMDEPBACKSLASH@
+testMatrixMultiplication-testFloatMatrixMultiplication.o: testFloatMatrixMultiplication.c
+@am__fastdepCC_TRUE@	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(testMatrixMultiplication_CFLAGS) $(CFLAGS) -MT testMatrixMultiplication-testFloatMatrixMultiplication.o -MD -MP -MF $(DEPDIR)/testMatrixMultiplication-testFloatMatrixMultiplication.Tpo -c -o testMatrixMultiplication-testFloatMatrixMultiplication.o `test -f 'testFloatMatrixMultiplication.c' || echo '$(srcdir)/'`testFloatMatrixMultiplication.c
+@am__fastdepCC_TRUE@	mv -f $(DEPDIR)/testMatrixMultiplication-testFloatMatrixMultiplication.Tpo $(DEPDIR)/testMatrixMultiplication-testFloatMatrixMultiplication.Po
+@AMDEP_TRUE@@am__fastdepCC_FALSE@	source='testFloatMatrixMultiplication.c' object='testMatrixMultiplication-testFloatMatrixMultiplication.o' libtool=no @AMDEPBACKSLASH@
 @AMDEP_TRUE@@am__fastdepCC_FALSE@	DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
-@am__fastdepCC_FALSE@	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(testMatrixMultiplication_CFLAGS) $(CFLAGS) -c -o testMatrixMultiplication-testMatrixMultiplication.o `test -f 'testMatrixMultiplication.c' || echo '$(srcdir)/'`testMatrixMultiplication.c
+@am__fastdepCC_FALSE@	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(testMatrixMultiplication_CFLAGS) $(CFLAGS) -c -o testMatrixMultiplication-testFloatMatrixMultiplication.o `test -f 'testFloatMatrixMultiplication.c' || echo '$(srcdir)/'`testFloatMatrixMultiplication.c
 
-testMatrixMultiplication-testMatrixMultiplication.obj: testMatrixMultiplication.c
-@am__fastdepCC_TRUE@	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(testMatrixMultiplication_CFLAGS) $(CFLAGS) -MT testMatrixMultiplication-testMatrixMultiplication.obj -MD -MP -MF $(DEPDIR)/testMatrixMultiplication-testMatrixMultiplication.Tpo -c -o testMatrixMultiplication-testMatrixMultiplication.obj `if test -f 'testMatrixMultiplication.c'; then $(CYGPATH_W) 'testMatrixMultiplication.c'; else $(CYGPATH_W) '$(srcdir)/testMatrixMultiplication.c'; fi`
-@am__fastdepCC_TRUE@	mv -f $(DEPDIR)/testMatrixMultiplication-testMatrixMultiplication.Tpo $(DEPDIR)/testMatrixMultiplication-testMatrixMultiplication.Po
-@AMDEP_TRUE@@am__fastdepCC_FALSE@	source='testMatrixMultiplication.c' object='testMatrixMultiplication-testMatrixMultiplication.obj' libtool=no @AMDEPBACKSLASH@
+testMatrixMultiplication-testFloatMatrixMultiplication.obj: testFloatMatrixMultiplication.c
+@am__fastdepCC_TRUE@	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(testMatrixMultiplication_CFLAGS) $(CFLAGS) -MT testMatrixMultiplication-testFloatMatrixMultiplication.obj -MD -MP -MF $(DEPDIR)/testMatrixMultiplication-testFloatMatrixMultiplication.Tpo -c -o testMatrixMultiplication-testFloatMatrixMultiplication.obj `if test -f 'testFloatMatrixMultiplication.c'; then $(CYGPATH_W) 'testFloatMatrixMultiplication.c'; else $(CYGPATH_W) '$(srcdir)/testFloatMatrixMultiplication.c'; fi`
+@am__fastdepCC_TRUE@	mv -f $(DEPDIR)/testMatrixMultiplication-testFloatMatrixMultiplication.Tpo $(DEPDIR)/testMatrixMultiplication-testFloatMatrixMultiplication.Po
+@AMDEP_TRUE@@am__fastdepCC_FALSE@	source='testFloatMatrixMultiplication.c' object='testMatrixMultiplication-testFloatMatrixMultiplication.obj' libtool=no @AMDEPBACKSLASH@
 @AMDEP_TRUE@@am__fastdepCC_FALSE@	DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
-@am__fastdepCC_FALSE@	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(testMatrixMultiplication_CFLAGS) $(CFLAGS) -c -o testMatrixMultiplication-testMatrixMultiplication.obj `if test -f 'testMatrixMultiplication.c'; then $(CYGPATH_W) 'testMatrixMultiplication.c'; else $(CYGPATH_W) '$(srcdir)/testMatrixMultiplication.c'; fi`
+@am__fastdepCC_FALSE@	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(testMatrixMultiplication_CFLAGS) $(CFLAGS) -c -o testMatrixMultiplication-testFloatMatrixMultiplication.obj `if test -f 'testFloatMatrixMultiplication.c'; then $(CYGPATH_W) 'testFloatMatrixMultiplication.c'; else $(CYGPATH_W) '$(srcdir)/testFloatMatrixMultiplication.c'; fi`
+
+testMatrixMultiplication-testDoubleMatrixMultiplication.o: testDoubleMatrixMultiplication.c
+@am__fastdepCC_TRUE@	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(testMatrixMultiplication_CFLAGS) $(CFLAGS) -MT testMatrixMultiplication-testDoubleMatrixMultiplication.o -MD -MP -MF $(DEPDIR)/testMatrixMultiplication-testDoubleMatrixMultiplication.Tpo -c -o testMatrixMultiplication-testDoubleMatrixMultiplication.o `test -f 'testDoubleMatrixMultiplication.c' || echo '$(srcdir)/'`testDoubleMatrixMultiplication.c
+@am__fastdepCC_TRUE@	mv -f $(DEPDIR)/testMatrixMultiplication-testDoubleMatrixMultiplication.Tpo $(DEPDIR)/testMatrixMultiplication-testDoubleMatrixMultiplication.Po
+@AMDEP_TRUE@@am__fastdepCC_FALSE@	source='testDoubleMatrixMultiplication.c' object='testMatrixMultiplication-testDoubleMatrixMultiplication.o' libtool=no @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCC_FALSE@	DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCC_FALSE@	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(testMatrixMultiplication_CFLAGS) $(CFLAGS) -c -o testMatrixMultiplication-testDoubleMatrixMultiplication.o `test -f 'testDoubleMatrixMultiplication.c' || echo '$(srcdir)/'`testDoubleMatrixMultiplication.c
+
+testMatrixMultiplication-testDoubleMatrixMultiplication.obj: testDoubleMatrixMultiplication.c
+@am__fastdepCC_TRUE@	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(testMatrixMultiplication_CFLAGS) $(CFLAGS) -MT testMatrixMultiplication-testDoubleMatrixMultiplication.obj -MD -MP -MF $(DEPDIR)/testMatrixMultiplication-testDoubleMatrixMultiplication.Tpo -c -o testMatrixMultiplication-testDoubleMatrixMultiplication.obj `if test -f 'testDoubleMatrixMultiplication.c'; then $(CYGPATH_W) 'testDoubleMatrixMultiplication.c'; else $(CYGPATH_W) '$(srcdir)/testDoubleMatrixMultiplication.c'; fi`
+@am__fastdepCC_TRUE@	mv -f $(DEPDIR)/testMatrixMultiplication-testDoubleMatrixMultiplication.Tpo $(DEPDIR)/testMatrixMultiplication-testDoubleMatrixMultiplication.Po
+@AMDEP_TRUE@@am__fastdepCC_FALSE@	source='testDoubleMatrixMultiplication.c' object='testMatrixMultiplication-testDoubleMatrixMultiplication.obj' libtool=no @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCC_FALSE@	DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCC_FALSE@	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(testMatrixMultiplication_CFLAGS) $(CFLAGS) -c -o testMatrixMultiplication-testDoubleMatrixMultiplication.obj `if test -f 'testDoubleMatrixMultiplication.c'; then $(CYGPATH_W) 'testDoubleMatrixMultiplication.c'; else $(CYGPATH_W) '$(srcdir)/testDoubleMatrixMultiplication.c'; fi`
 
 testMatrixSubtraction-testMatrixSubtraction.o: testMatrixSubtraction.c
 @am__fastdepCC_TRUE@	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(testMatrixSubtraction_CFLAGS) $(CFLAGS) -MT testMatrixSubtraction-testMatrixSubtraction.o -MD -MP -MF $(DEPDIR)/testMatrixSubtraction-testMatrixSubtraction.Tpo -c -o testMatrixSubtraction-testMatrixSubtraction.o `test -f 'testMatrixSubtraction.c' || echo '$(srcdir)/'`testMatrixSubtraction.c
diff --git a/src/matrixOperations/matrixMultiplication.c b/src/matrixOperations/matrixMultiplication.c
index e8e66e60..dbed4330 100644
--- a/src/matrixOperations/matrixMultiplication.c
+++ b/src/matrixOperations/matrixMultiplication.c
@@ -10,8 +10,9 @@
  *
  */
 
-#include <assert.h>
+#ifndef WITHOUT_BLAS
 #include "blas.h"
+#endif
 #include "matrixMultiplication.h"
 
 /*
@@ -32,15 +33,19 @@ void	smulma(float *in1, int lines1, int columns1,
   int k = 0;
   float accu = 0;
 
-#define in1(a, b)	in1[a+b*lines1]
-#define in2(c, d)	in2[c+d*lines2]
+  /*
+  **  How to convert 2 index matrixes to one.
+  ** #define in1(a, b)	in1[a+b*lines1]
+  ** #define in2(c, d)	in2[c+d*lines2]
+  */
 
   for (i = 0 ; i < lines1 * columns2 ; ++i)
     {
       accu = 0;
       for (k = 0; k < columns1 ; ++k)
 	{
-	  accu += in1(i % lines1, k) * in2(k, i / lines1);
+	  accu += in1[i % lines1 + k * lines1]
+	    * in2[k + (i / lines1) * lines2];
 	}
       out[i] = accu;
     }
@@ -61,12 +66,41 @@ void	dmulma(double *in1, int lines1, int columns1,
 	       double *in2, int lines2, int columns2,
 	       double *out)
 {
+#ifndef WITHOUT_BLAS
+  /*
+  ** USES BLAS DGEMM FUNCTION.
+  */
   double One		= 1;
   double Zero		= 0;
 
   /* Cr <-  1*Ar*Br + 0*Cr */
   dgemm_("N","N", &columns2, &columns2, &columns1, &One,
 	 in1 , &lines1, in2, &lines2, &Zero, out, &columns2);
+#else
+  /*
+  ** DO NOT USE ANY BLAS FUNCTION.
+  */
+  int i = 0;
+  int k = 0;
+  double accu = 0;
+
+  /*
+  **  How to convert 2 index matrixes to one.
+  ** #define in1(a, b)	in1[a+b*lines1]
+  ** #define in2(c, d)	in2[c+d*lines2]
+  */
+
+  for (i = 0 ; i < lines1 * columns2 ; ++i)
+    {
+      accu = 0;
+      for (k = 0; k < columns1 ; ++k)
+	{
+	  accu += in1[i % lines1 + k * lines1]
+	    * in2[k + (i / lines1) * lines2];
+	}
+      out[i] = accu;
+    }
+#endif
 }
 
 /*
@@ -92,7 +126,9 @@ void	cmulma(floatComplex *in1, int lines1, int columns1,
       accu = FloatComplex(0,0);
       for (k = 0; k < columns1 ; ++k)
 	{
-	  cadds(accu, ctimess(in1(i % lines1, k) , in2(k, i / lines1)));
+	  accu = cadds(accu,
+		ctimess(in1[i % lines1 + k *lines1] ,
+			in2[k + (i / lines1) *lines2] ));
 	}
       out[i] = accu;
     }
@@ -112,6 +148,11 @@ void	zmulma(doubleComplex *in1, int lines1, int columns1,
 	       doubleComplex *in2, int lines2, int columns2,
 	       doubleComplex *out)
 {
+#ifndef WITHOUT_BLAS
+  /*
+  ** USES BLAS DGEMM FUNCTION.
+  */
+  int i = 0;
   double One		= 1;
   double MinusOne	= -1;
   double Zero		= 0;
@@ -141,13 +182,40 @@ void	zmulma(doubleComplex *in1, int lines1, int columns1,
   dgemm_("N","N", &lines1, &columns2, &columns1, &One,
 	 in1Real, &lines1, in2Imag, &lines2, &Zero, ImagOut, &lines1);
 
-  /*Ci <-  1*Ai*Br + 1*Ci */
+  /* Ci <-  1*Ai*Br + 1*Ci */
   dgemm_("N","N", &lines1, &columns2, &columns1, &One,
 	 in1Imag, &lines1, in2Real, &lines2, &One, ImagOut, &lines1);
 
-  out = DoubleComplexMatrix(RealOut, ImagOut, lines1 * columns2);
+  /* Now fill output matrix */
+  for(i = 0 ; i < lines1 * columns2 ; ++i)
+    {
+      out[i] = DoubleComplex(RealOut[i], ImagOut[i]);
+    }
+  /* FREE allocated variables */
   free(in1Real);
   free(in2Real);
   free(in1Imag);
-  free(in1Imag);
+  free(in2Imag);
+  free(RealOut);
+  free(ImagOut);
+#else
+  /*
+  ** DO NOT USE ANY BLAS FUNCTION.
+  */
+  int i = 0;
+  int k = 0;
+  doubleComplex accu = DoubleComplex(0, 0);
+
+  for (i = 0 ; i < lines1 * columns2 ; ++i)
+    {
+      accu = DoubleComplex(0,0);
+      for (k = 0; k < columns1 ; ++k)
+	{
+	  accu = zadds(accu,
+		ztimess(in1[i % lines1 + k *lines1] ,
+			in2[k + (i / lines1) *lines2] ));
+	}
+      out[i] = accu;
+    }
+#endif
 }
diff --git a/src/matrixOperations/testDoubleMatrixMultiplication.c b/src/matrixOperations/testDoubleMatrixMultiplication.c
new file mode 100644
index 00000000..e9b71a87
--- /dev/null
+++ b/src/matrixOperations/testDoubleMatrixMultiplication.c
@@ -0,0 +1,211 @@
+/*
+ *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+ *  Copyright (C) 2008-2008 - INRIA - Bruno JOFRET
+ *
+ *  This file must be used under the terms of the CeCILL.
+ *  This source file is licensed as described in the file COPYING, which
+ *  you should have received as part of this distribution.  The terms
+ *  are also available at
+ *  http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+ *
+ */
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <assert.h>
+#include <math.h>
+#include "matrixMultiplication.h"
+
+int testFloatMultiplication(void);
+
+static void zmulmaTest(void) {
+  double realM1[4] = {1.0, 2.0, 3.0, 4.0};
+  double imagM1[4] = {1.0, 2.0, 3.0, 4.0};
+  double realM3[6] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
+  double imagM3[6] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
+
+  doubleComplex *M1;
+  doubleComplex *M2;
+  doubleComplex M1_mul_M2[4];
+  doubleComplex *M3;
+  doubleComplex *M4;
+  doubleComplex M3_mul_M4[4];
+  doubleComplex miscM3_mul_M4[9];
+
+  int i = 0;
+
+  printf("\n>>>> Matrix Complex Double Multiplication Tests\n");
+
+  M1 = DoubleComplexMatrix(realM1, imagM1, 4);
+  M2 = DoubleComplexMatrix(realM1, imagM1, 4);
+
+  /*
+    [ 1+1.%i 3+3.%i] * [ 1+1.%i 3+3.%i] = [ 14.%i 30.%i ]
+    [ 2+2.%i 4+4.%i]   [ 2+2.%i 4+4.%i]   [ 20.%i 44.%i ]
+  */
+  zmulma(M1, 2, 2, M2, 2, 2, M1_mul_M2);
+  for (i = 0; i < 4; ++i)
+    {
+      printf("M1_mul_M2[%d] = %e + %e i\n", i, zreals(M1_mul_M2[i]), zimags(M1_mul_M2[i]));
+    }
+  for (i = 0; i < 4; ++i)
+    {
+      assert(zreals(M1_mul_M2[i]) == 0.0);
+    }
+  assert(zimags(M1_mul_M2[0]) == 14.0);
+  assert(zimags(M1_mul_M2[1]) == 20.0);
+  assert(zimags(M1_mul_M2[2]) == 30.0);
+  assert(zimags(M1_mul_M2[3]) == 44.0);
+
+  M3 = DoubleComplexMatrix(realM3, imagM3, 6);
+  M4 = DoubleComplexMatrix(realM3, imagM3, 6);
+
+  /*
+    [ 1+1.%i 3+3.%i 5+5.%i ] * [ 1+1.%i 4+4.%i ] = [ 44.%i  98.%i ]
+    [ 2+2.%i 4+4.%i 6+6.%i ]   [ 2+2.%i 5+5.%i ]   [ 56.%i 128.%i ]
+			       [ 3+3.%i 6+6.%i ]
+   */
+  zmulma(M3, 2, 3, M4, 3, 2, M3_mul_M4);
+  for (i = 0; i < 4; ++i)
+    {
+      printf("M3_mul_M4[%d] = %e + %e i\n", i, zreals(M3_mul_M4[i]), zimags(M3_mul_M4[i]));
+    }
+  for (i = 0; i < 4; ++i)
+    {
+      assert(zreals(M3_mul_M4[i]) == 0.0);
+    }
+  assert(zimags(M3_mul_M4[0]) == 44.0);
+  assert(zimags(M3_mul_M4[1]) == 56.0);
+  assert(zimags(M3_mul_M4[2]) == 98.0);
+  assert(zimags(M3_mul_M4[3]) == 128.0);
+
+  /*
+    [ 1+1.%i 4+4.%i ] * [ 1+1.%i 3+3.%i 5+5.%i ] = [ 18.%i 38.%i  58.%i ]
+    [ 2+2.%i 5+5.%i ]   [ 2+2.%i 4+4.%i 6+6.%i ]   [ 24.%i 52.%i  80.%i ]
+    [ 3+3.%i 6+6.%i ]				   [ 30.%i 66.%i 102.%i ]
+   */
+  zmulma(M3, 3, 2, M4, 2, 3, miscM3_mul_M4);
+  for (i = 0; i < 9; ++i)
+    {
+      printf("miscM3_mul_M4[%d] = %e + %e i\n", i, zreals(miscM3_mul_M4[i]), zimags(miscM3_mul_M4[i]));
+    }
+  for (i = 0; i < 9; ++i)
+    {
+      assert(zreals(miscM3_mul_M4[i]) == 0.0);
+    }
+  assert(zimags(miscM3_mul_M4[0]) == 18.0);
+  assert(zimags(miscM3_mul_M4[1]) == 24.0);
+  assert(zimags(miscM3_mul_M4[2]) == 30.0);
+  assert(zimags(miscM3_mul_M4[3]) == 38.0);
+  assert(zimags(miscM3_mul_M4[4]) == 52.0);
+  assert(zimags(miscM3_mul_M4[5]) == 66.0);
+  assert(zimags(miscM3_mul_M4[6]) == 58.0);
+  assert(zimags(miscM3_mul_M4[7]) == 80.0);
+  assert(zimags(miscM3_mul_M4[8]) == 102.0);
+}
+
+static void dmulmaTest(void) {
+  double M1[4] = {1.0, 2.0, 3.0, 4.0};
+  double M2[4] = {1.0, 2.0, 3.0, 4.0};
+  double M1_by_M2[4];
+  double M3[4] = {1.0, 0.0, 1.0, 0.0};
+  double M4[4] = {0.0, 1.0, 0.0, 1.0};
+  double M3_by_M4[4];
+  double M5[4] = {1.0, 0.0, 0.0, 1.0};
+  double M6[4] = {42.0, 51.0, 69.0, 1664.0};
+  double M5_by_M6[4];
+  double M7[6] = {1.0, 4.0, 2.0, 5.0, 3.0, 6.0};
+  double M8[6] = {1.0, 3.0, 5.0, 2.0, 4.0, 6.0};
+  double M7_by_M8[4];
+  double miscM7_by_M8[9];
+
+  int i = 0;
+
+  printf("\n>>>> Matrix Real Double Multiplication Tests\n");
+  /*
+    [ 1 3 ]  * [ 1 3 ] = [ 7 15 ]
+    [ 2 4 ]    [ 2 4 ]   [10 22 ]
+  */
+  dmulma(M1, 2, 2, M2, 2, 2, M1_by_M2);
+  for (i = 0; i < 4; ++i) {
+    printf("M1_by_M2[%d] = %e\n", i, M1_by_M2[i]);
+  }
+  assert(M1_by_M2[0] == 7.0);
+  assert(M1_by_M2[1] == 10.0);
+  assert(M1_by_M2[2] == 15.0);
+  assert(M1_by_M2[3] == 22.0);
+  /*
+    [ 1 1 ]  * [ 0 0 ] = [ 1 1 ]
+    [ 0 0 ]    [ 1 1 ]   [ 0 0 ]
+  */
+  dmulma(M3, 2, 2, M4, 2, 2, M3_by_M4);
+  for (i = 0; i < 4; ++i) {
+    printf("M3_by_M4[%d] = %e\n", i, M3_by_M4[i]);
+  }
+  assert(M3_by_M4[0] == 1.0);
+  assert(M3_by_M4[1] == 0.0);
+  assert(M3_by_M4[2] == 1.0);
+  assert(M3_by_M4[3] == 0.0);
+
+  /*
+    [ 1 0 ]  * [ 42   69 ] = [ 42   69 ]
+    [ 0 1 ]    [ 51 1664 ]   [ 51 1664 ]
+  */
+  dmulma(M5, 2, 2, M6, 2, 2, M5_by_M6);
+  for (i = 0; i < 4; ++i) {
+    printf("M5_by_M6[%d] = %e\n", i, M5_by_M6[i]);
+  }
+  assert(M5_by_M6[0] == 42.0);
+  assert(M5_by_M6[1] == 51.0);
+  assert(M5_by_M6[2] == 69.0);
+  assert(M5_by_M6[3] == 1664.0);
+
+  /*
+    [ 1 2 3 ]  * [ 1 2 ] = [ 22 28 ]
+    [ 4 5 6 ]    [ 3 4 ]   [ 49 64 ]
+		 [ 5 6 ]
+  */
+  dmulma(M7, 2, 3, M8, 3, 2, M7_by_M8);
+  for (i = 0; i < 4; ++i) {
+    printf("M7_by_M8[%d] = %e\n", i, M7_by_M8[i]);
+  }
+  assert(M7_by_M8[0] == 22.0);
+  assert(M7_by_M8[1] == 49.0);
+  assert(M7_by_M8[2] == 28.0);
+  assert(M7_by_M8[3] == 64.0);
+
+  /*
+    [ 1 5 ]  * [ 1 5 4 ] = [ 16 15 34 ]
+    [ 4 3 ]    [ 3 2 6 ]   [ 13 26 34 ]
+    [ 2 6 ]		   [ 20 22 44 ]
+  */
+  dmulma(M7, 3, 2, M8, 2, 3, miscM7_by_M8);
+  for (i = 0; i < 9; ++i) {
+    printf("miscM7_by_M8[%d] = %e\n", i, miscM7_by_M8[i]);
+  }
+  assert(miscM7_by_M8[0] == 16.0);
+  assert(miscM7_by_M8[1] == 13.0);
+  assert(miscM7_by_M8[2] == 20.0);
+  assert(miscM7_by_M8[3] == 15.0);
+  assert(miscM7_by_M8[4] == 26.0);
+  assert(miscM7_by_M8[5] == 22.0);
+  assert(miscM7_by_M8[6] == 34.0);
+  assert(miscM7_by_M8[7] == 34.0);
+  assert(miscM7_by_M8[8] == 44.0);
+}
+
+static int testDoubleMultiplication(void) {
+
+  printf("\n>>>> Matrix Double Multiplication Tests\n");
+  dmulmaTest();
+  zmulmaTest();
+  return 0;
+}
+
+
+
+int main(void) {
+  assert(testFloatMultiplication() == 0);
+  assert(testDoubleMultiplication() == 0);
+  return 0;
+}
diff --git a/src/matrixOperations/testFloatMatrixMultiplication.c b/src/matrixOperations/testFloatMatrixMultiplication.c
new file mode 100644
index 00000000..6091fc66
--- /dev/null
+++ b/src/matrixOperations/testFloatMatrixMultiplication.c
@@ -0,0 +1,206 @@
+/*
+ *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+ *  Copyright (C) 2008-2008 - INRIA - Bruno JOFRET
+ *
+ *  This file must be used under the terms of the CeCILL.
+ *  This source file is licensed as described in the file COPYING, which
+ *  you should have received as part of this distribution.  The terms
+ *  are also available at
+ *  http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+ *
+ */
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <assert.h>
+#include <math.h>
+#include "matrixMultiplication.h"
+
+int testFloatMultiplication(void);
+
+static void cmulmaTest(void) {
+  float realM1[4] = {1.0f, 2.0f, 3.0f, 4.0f};
+  float imagM1[4] = {1.0f, 2.0f, 3.0f, 4.0f};
+  float realM3[6] = {1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f};
+  float imagM3[6] = {1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f};
+
+  floatComplex *M1;
+  floatComplex *M2;
+  floatComplex M1_mul_M2[4];
+  floatComplex *M3;
+  floatComplex *M4;
+  floatComplex M3_mul_M4[4];
+  floatComplex miscM3_mul_M4[9];
+
+  int i = 0;
+
+  printf("\n>>>> Matrix Complex Float Multiplication Tests\n");
+
+  M1 = FloatComplexMatrix(realM1, imagM1, 4);
+  M2 = FloatComplexMatrix(realM1, imagM1, 4);
+
+  /*
+    [ 1+1.%i 3+3.%i] * [ 1+1.%i 3+3.%i] = [ 14.%i 30.%i ]
+    [ 2+2.%i 3+3.%i]   [ 2+2.%i 3+3.%i]   [ 20.%i 44.%i ]
+  */
+  cmulma(M1, 2, 2, M2, 2, 2, M1_mul_M2);
+  for (i = 0; i < 4; ++i)
+    {
+      printf("M1_mul_M2[%d] = %e + %e i\n", i, creals(M1_mul_M2[i]), cimags(M1_mul_M2[i]));
+    }
+  for (i = 0; i < 4; ++i)
+    {
+      assert(creals(M1_mul_M2[i]) == 0.0f);
+    }
+  assert(cimags(M1_mul_M2[0]) == 14.0f);
+  assert(cimags(M1_mul_M2[1]) == 20.0f);
+  assert(cimags(M1_mul_M2[2]) == 30.0f);
+  assert(cimags(M1_mul_M2[3]) == 44.0f);
+
+  M3 = FloatComplexMatrix(realM3, imagM3, 6);
+  M4 = FloatComplexMatrix(realM3, imagM3, 6);
+
+  /*
+    [ 1+1.%i 3+3.%i 5+5.%i ] * [ 1+1.%i 4+4.%i ] = [ 44.%i  98.%i ]
+    [ 2+2.%i 4+4.%i 6+6.%i ]   [ 2+2.%i 5+5.%i ]   [ 56.%i 128.%i ]
+			       [ 3+3.%i 6+6.%i ]
+   */
+  cmulma(M3, 2, 3, M4, 3, 2, M3_mul_M4);
+  for (i = 0; i < 4; ++i)
+    {
+      printf("M3_mul_M4[%d] = %e + %e i\n", i, creals(M3_mul_M4[i]), cimags(M3_mul_M4[i]));
+    }
+  for (i = 0; i < 4; ++i)
+    {
+      assert(creals(M3_mul_M4[i]) == 0.0f);
+    }
+  assert(cimags(M3_mul_M4[0]) == 44.0f);
+  assert(cimags(M3_mul_M4[1]) == 56.0f);
+  assert(cimags(M3_mul_M4[2]) == 98.0f);
+  assert(cimags(M3_mul_M4[3]) == 128.0f);
+
+  /*
+    [ 1+1.%i 4+4.%i ] * [ 1+1.%i 3+3.%i 5+5.%i ] = [ 18.%i 38.%i  58.%i ]
+    [ 2+2.%i 5+5.%i ]   [ 2+2.%i 4+4.%i 6+6.%i ]   [ 24.%i 52.%i  80.%i ]
+    [ 3+3.%i 6+6.%i ]				   [ 30.%i 66.%i 102.%i ]
+   */
+  cmulma(M3, 3, 2, M4, 2, 3, miscM3_mul_M4);
+  for (i = 0; i < 9; ++i)
+    {
+      printf("miscM3_mul_M4[%d] = %e + %e i\n", i, creals(miscM3_mul_M4[i]), cimags(miscM3_mul_M4[i]));
+    }
+  for (i = 0; i < 9; ++i)
+    {
+      assert(creals(miscM3_mul_M4[i]) == 0.0);
+    }
+  assert(cimags(miscM3_mul_M4[0]) == 18.0);
+  assert(cimags(miscM3_mul_M4[1]) == 24.0);
+  assert(cimags(miscM3_mul_M4[2]) == 30.0);
+  assert(cimags(miscM3_mul_M4[3]) == 38.0);
+  assert(cimags(miscM3_mul_M4[4]) == 52.0);
+  assert(cimags(miscM3_mul_M4[5]) == 66.0);
+  assert(cimags(miscM3_mul_M4[6]) == 58.0);
+  assert(cimags(miscM3_mul_M4[7]) == 80.0);
+  assert(cimags(miscM3_mul_M4[8]) == 102.0);
+
+}
+
+static void smulmaTest(void) {
+  float M1[4] = {1.0f, 2.0f, 3.0f, 4.0f};
+  float M2[4] = {1.0f, 2.0f, 3.0f, 4.0f};
+  float M1_by_M2[4];
+  float M3[4] = {1.0f, 0.0f, 1.0f, 0.0f};
+  float M4[4] = {0.0f, 1.0f, 0.0f, 1.0f};
+  float M3_by_M4[4];
+  float M5[4] = {1.0f, 0.0f, 0.0f, 1.0f};
+  float M6[4] = {42.0f, 51.0f, 69.0f, 1664.0f};
+  float M5_by_M6[4];
+  float M7[6] = {1.0f, 4.0f, 2.0f, 5.0f, 3.0f, 6.0f};
+  float M8[6] = {1.0f, 3.0f, 5.0f, 2.0f, 4.0f, 6.0f};
+  float M7_by_M8[4];
+  float miscM7_by_M8[9];
+
+  int i = 0;
+
+  printf("\n>>>> Matrix Float Multiplication Tests\n");
+  /*
+    [ 1 3 ]  * [ 1 3 ] = [ 7 15 ]
+    [ 2 4 ]    [ 2 4 ]   [10 22 ]
+  */
+  smulma(M1, 2, 2, M2, 2, 2, M1_by_M2);
+  for (i = 0; i < 4; ++i) {
+    printf("M1_by_M2[%d] = %e\n", i, M1_by_M2[i]);
+  }
+  assert(M1_by_M2[0] == 7.0f);
+  assert(M1_by_M2[1] == 10.0f);
+  assert(M1_by_M2[2] == 15.0f);
+  assert(M1_by_M2[3] == 22.0f);
+  /*
+    [ 1 1 ]  * [ 0 0 ] = [ 1 1 ]
+    [ 0 0 ]    [ 1 1 ]   [ 0 0 ]
+  */
+  smulma(M3, 2, 2, M4, 2, 2, M3_by_M4);
+  for (i = 0; i < 4; ++i) {
+    printf("M3_by_M4[%d] = %e\n", i, M3_by_M4[i]);
+  }
+  assert(M3_by_M4[0] == 1.0f);
+  assert(M3_by_M4[1] == 0.0f);
+  assert(M3_by_M4[2] == 1.0f);
+  assert(M3_by_M4[3] == 0.0f);
+
+  /*
+    [ 1 0 ]  * [ 42   69 ] = [ 42   69 ]
+    [ 0 1 ]    [ 51 1664 ]   [ 51 1664 ]
+  */
+  smulma(M5, 2, 2, M6, 2, 2, M5_by_M6);
+  for (i = 0; i < 4; ++i) {
+    printf("M5_by_M6[%d] = %e\n", i, M5_by_M6[i]);
+  }
+  assert(M5_by_M6[0] == 42.0f);
+  assert(M5_by_M6[1] == 51.0f);
+  assert(M5_by_M6[2] == 69.0f);
+  assert(M5_by_M6[3] == 1664.0f);
+
+  /*
+    [ 1 2 3 ]  * [ 1 2 ] = [ 22 28 ]
+    [ 4 5 6 ]    [ 3 4 ]   [ 49 64 ]
+    [ 5 6 ]
+  */
+  smulma(M7, 2, 3, M8, 3, 2, M7_by_M8);
+  for (i = 0; i < 4; ++i) {
+    printf("M7_by_M8[%d] = %e\n", i, M7_by_M8[i]);
+  }
+  assert(M7_by_M8[0] == 22.0f);
+  assert(M7_by_M8[1] == 49.0f);
+  assert(M7_by_M8[2] == 28.0f);
+  assert(M7_by_M8[3] == 64.0f);
+
+  /*
+    [ 1 5 ]  * [ 1 5 4 ] = [ 16 15 34 ]
+    [ 4 3 ]    [ 3 2 6 ]   [ 13 26 34 ]
+    [ 2 6 ]		   [ 20 22 44 ]
+  */
+  smulma(M7, 3, 2, M8, 2, 3, miscM7_by_M8);
+  for (i = 0; i < 9; ++i) {
+    printf("miscM7_by_M8[%d] = %e\n", i, miscM7_by_M8[i]);
+  }
+  assert(miscM7_by_M8[0] == 16.0f);
+  assert(miscM7_by_M8[1] == 13.0f);
+  assert(miscM7_by_M8[2] == 20.0f);
+  assert(miscM7_by_M8[3] == 15.0f);
+  assert(miscM7_by_M8[4] == 26.0f);
+  assert(miscM7_by_M8[5] == 22.0f);
+  assert(miscM7_by_M8[6] == 34.0f);
+  assert(miscM7_by_M8[7] == 34.0f);
+  assert(miscM7_by_M8[8] == 44.0f);
+
+}
+
+int testFloatMultiplication(void) {
+
+  printf("\n>>>> Matrix Multiplication Tests\n");
+  smulmaTest();
+  cmulmaTest();
+
+  return 0;
+}
diff --git a/src/matrixOperations/testMatrixMultiplication.c b/src/matrixOperations/testMatrixMultiplication.c
deleted file mode 100644
index 38f59405..00000000
--- a/src/matrixOperations/testMatrixMultiplication.c
+++ /dev/null
@@ -1,224 +0,0 @@
-/*
- *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- *  Copyright (C) 2008-2008 - INRIA - Bruno JOFRET
- *
- *  This file must be used under the terms of the CeCILL.
- *  This source file is licensed as described in the file COPYING, which
- *  you should have received as part of this distribution.  The terms
- *  are also available at
- *  http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
- *
- */
-
-#include <stdlib.h>
-#include <stdio.h>
-#include <assert.h>
-#include <math.h>
-#include "matrixMultiplication.h"
-
-
-#define LOCAL_DEBUG
-
-/*
-  static void zaddmaTest(void) {
-  doubleComplex M1[SIZE];
-  doubleComplex M2[SIZE];
-  doubleComplex M1_mul_M2[SIZE];
-
-  int i = 0;
-
-  srand(1);
-  for (i = 0; i < SIZE; ++i) {
-  M1[i] = DoubleComplex((double) rand(), (double) rand());
-  M2[i] = DoubleComplex((double) rand(), (double) rand());
-  }
-  zaddma(M1, SIZE, M2, SIZE, M1_and_M2);
-  for (i = 0; i < SIZE; ++i) {
-  #ifdef LOCAL_DEBUG
-  printf("M1_and_M2 = %e + %e i\n", zreals(M1_and_M2[i]), zimags(M1_and_M2[i]));
-  #endif
-  assert(zreals(M1_and_M2[i]) -  (zreals(M1[i]) + zreals(M2[i])) == 0);
-  assert(zimags(M1_and_M2[i]) -  (zimags(M1[i]) + zimags(M2[i])) == 0);
-  }
-  }
-
-  static void caddmaTest(void) {
-  floatComplex M1[SIZE];
-  floatComplex M2[SIZE];
-  floatComplex M1_and_M2[SIZE];
-
-  int i = 0;
-
-  srand(1);
-  for (i = 0; i < SIZE; ++i) {
-  M1[i] = FloatComplex((float) rand(), (float) rand());
-  M2[i] = FloatComplex((float) rand(), (float) rand());
-  }
-  caddma(M1, SIZE, M2, SIZE, M1_and_M2);
-  for (i = 0; i < SIZE; ++i) {
-  #ifdef LOCAL_DEBUG
-  printf("M1_and_M2 = %e + %e i\n", creals(M1_and_M2[i]), cimags(M1_and_M2[i]));
-  #endif
-  assert(fabsf(creals(M1_and_M2[i]) -  (creals(M1[i]) + creals(M2[i]))) / creals(M1_and_M2[i]) < 1e-07);
-  assert(fabsf(cimags(M1_and_M2[i]) -  (cimags(M1[i]) + cimags(M2[i]))) / cimags(M1_and_M2[i]) < 1e-07);
-  }
-  }
-
-*/
-static void dmulmaTest(void) {
-  double M1[4] = {1.0, 2.0, 3.0, 4.0};
-  double M2[4] = {1.0, 2.0, 3.0, 4.0};
-  double M1_by_M2[4];
-  double M3[4] = {1.0, 0.0, 1.0, 0.0};
-  double M4[4] = {0.0, 1.0, 0.0, 1.0};
-  double M3_by_M4[4];
-  double M5[4] = {1.0, 0.0, 0.0, 1.0};
-  double M6[4] = {42.0, 51.0, 69.0, 1664.0};
-  double M5_by_M6[4];
-  double M7[6] = {1.0, 4.0, 2.0, 5.0, 3.0, 6.0};
-  double M8[6] = {1.0, 3.0, 5.0, 2.0, 4.0, 6.0};
-  double M7_by_M8[4];
-
- int i = 0;
-
-  printf("\n>>>> Matrix Double Multiplication Tests\n");
-  /*
-    [ 1 3 ]  * [ 1 3 ] = [ 7 15 ]
-    [ 2 4 ]    [ 2 4 ]   [10 22 ]
-  */
-  dmulma(M1, 2, 2, M2, 2, 2, M1_by_M2);
-  for (i = 0; i < 4; ++i) {
-    printf("M1_by_M2[%d] = %e\n", i, M1_by_M2[i]);
-  }
-  assert(M1_by_M2[0] == 7.0);
-  assert(M1_by_M2[1] == 10.0);
-  assert(M1_by_M2[2] == 15.0);
-  assert(M1_by_M2[3] == 22.0);
-  /*
-    [ 1 1 ]  * [ 0 0 ] = [ 1 1 ]
-    [ 0 0 ]    [ 1 1 ]   [ 0 0 ]
-  */
-  dmulma(M3, 2, 2, M4, 2, 2, M3_by_M4);
-  for (i = 0; i < 4; ++i) {
-    printf("M3_by_M4[%d] = %e\n", i, M3_by_M4[i]);
-  }
-  assert(M3_by_M4[0] == 1.0);
-  assert(M3_by_M4[1] == 0.0);
-  assert(M3_by_M4[2] == 1.0);
-  assert(M3_by_M4[3] == 0.0);
-
-  /*
-    [ 1 0 ]  * [ 42   69 ] = [ 42   69 ]
-    [ 0 1 ]    [ 51 1664 ]   [ 51 1664 ]
-  */
-  dmulma(M5, 2, 2, M6, 2, 2, M5_by_M6);
-  for (i = 0; i < 4; ++i) {
-    printf("M5_by_M6[%d] = %e\n", i, M5_by_M6[i]);
-  }
-  assert(M5_by_M6[0] == 42.0);
-  assert(M5_by_M6[1] == 51.0);
-  assert(M5_by_M6[2] == 69.0);
-  assert(M5_by_M6[3] == 1664.0);
-
-  /*
-    [ 1 2 3 ]  * [ 1 2 ] = [ 22 28 ]
-    [ 4 5 6 ]    [ 3 4 ]   [ 49 64 ]
-		 [ 5 6 ]
-  */
-  dmulma(M7, 2, 3, M8, 3, 2, M7_by_M8);
-  for (i = 0; i < 4; ++i) {
-    printf("M7_by_M8[%d] = %e\n", i, M7_by_M8[i]);
-  }
-  assert(M7_by_M8[0] == 22.0);
-  assert(M7_by_M8[1] == 49.0);
-  assert(M7_by_M8[2] == 28.0);
-  assert(M7_by_M8[3] == 64.0);
-}
-
-static void smulmaTest(void) {
-  float M1[4] = {1.0f, 2.0f, 3.0f, 4.0f};
-  float M2[4] = {1.0f, 2.0f, 3.0f, 4.0f};
-  float M1_by_M2[4];
-  float M3[4] = {1.0f, 0.0f, 1.0f, 0.0f};
-  float M4[4] = {0.0f, 1.0f, 0.0f, 1.0f};
-  float M3_by_M4[4];
-  float M5[4] = {1.0f, 0.0f, 0.0f, 1.0f};
-  float M6[4] = {42.0f, 51.0f, 69.0f, 1664.0f};
-  float M5_by_M6[4];
-  float M7[6] = {1.0f, 4.0f, 2.0f, 5.0f, 3.0f, 6.0f};
-  float M8[6] = {1.0f, 3.0f, 5.0f, 2.0f, 4.0f, 6.0f};
-  float M7_by_M8[4];
-  int i = 0;
-
-  printf("\n>>>> Matrix Float Multiplication Tests\n");
-  /*
-    [ 1 3 ]  * [ 1 3 ] = [ 7 15 ]
-    [ 2 4 ]    [ 2 4 ]   [10 22 ]
-  */
-  smulma(M1, 2, 2, M2, 2, 2, M1_by_M2);
-  for (i = 0; i < 4; ++i) {
-    printf("M1_by_M2[%d] = %e\n", i, M1_by_M2[i]);
-  }
-  assert(M1_by_M2[0] == 7.0f);
-  assert(M1_by_M2[1] == 10.0f);
-  assert(M1_by_M2[2] == 15.0f);
-  assert(M1_by_M2[3] == 22.0f);
-  /*
-    [ 1 1 ]  * [ 0 0 ] = [ 1 1 ]
-    [ 0 0 ]    [ 1 1 ]   [ 0 0 ]
-  */
-  smulma(M3, 2, 2, M4, 2, 2, M3_by_M4);
-  for (i = 0; i < 4; ++i) {
-    printf("M3_by_M4[%d] = %e\n", i, M3_by_M4[i]);
-  }
-  assert(M3_by_M4[0] == 1.0f);
-  assert(M3_by_M4[1] == 0.0f);
-  assert(M3_by_M4[2] == 1.0f);
-  assert(M3_by_M4[3] == 0.0f);
-
-  /*
-    [ 1 0 ]  * [ 42   69 ] = [ 42   69 ]
-    [ 0 1 ]    [ 51 1664 ]   [ 51 1664 ]
-  */
-  smulma(M5, 2, 2, M6, 2, 2, M5_by_M6);
-  for (i = 0; i < 4; ++i) {
-    printf("M5_by_M6[%d] = %e\n", i, M5_by_M6[i]);
-  }
-  assert(M5_by_M6[0] == 42.0f);
-  assert(M5_by_M6[1] == 51.0f);
-  assert(M5_by_M6[2] == 69.0f);
-  assert(M5_by_M6[3] == 1664.0f);
-
-  /*
-    [ 1 2 3 ]  * [ 1 2 ] = [ 22 28 ]
-    [ 4 5 6 ]    [ 3 4 ]   [ 49 64 ]
-		 [ 5 6 ]
-  */
-  smulma(M7, 2, 3, M8, 3, 2, M7_by_M8);
-  for (i = 0; i < 4; ++i) {
-    printf("M7_by_M8[%d] = %e\n", i, M7_by_M8[i]);
-  }
-  assert(M7_by_M8[0] == 22.0f);
-  assert(M7_by_M8[1] == 49.0f);
-  assert(M7_by_M8[2] == 28.0f);
-  assert(M7_by_M8[3] == 64.0f);
-}
-
-static int testMultiplication(void) {
-
-  printf("\n>>>> Matrix Multiplication Tests\n");
-  dmulmaTest();
-  smulmaTest();
-  /*
-    cmulmaTest();
-    zmulmaTest();
-  */
-  return 0;
-}
-
-
-
-int main(void) {
-  assert(testMultiplication() == 0);
-  return 0;
-}