/*!
  \file   matmul.c
  \brief  Matrix multiplication implementation.

  \author Nikos Pitsianis
  \author Dimitris Floros
  \author Christos Choutouridis 8997 <cchoutou@ece.auth.gr>
  \date   2020-05-05
*/
#include "matmul.h"

/*!
 * Square Matrix multiplication - ijk
 * \param	C	pointer to output matrix
 * \param	A	pointer to input matrix A
 * \param	B	pointer to input matrix B
 * \param	n	Size of matrices (both sizes)
 * \return 	   none
 */
void matrixMult_ijk(float * const C, float const * const A, float const * const B, int const n) {
   for (int i = 0; i < n; ++i)
      for (int j = 0; j < n; ++j) {
         int k =0;
         C[ sub2ind(i,j,n) ] = A[ sub2ind(i,k,n) ] * B[ sub2ind(k,j,n) ];
         for (k = 1; k < n; ++k)
            C[ sub2ind(i,j,n) ] += A[ sub2ind(i,k,n) ] * B[ sub2ind(k,j,n) ];
      }
}

/*!
 * Square Matrix multiplication - ikj
 * \param   C  pointer to output matrix
 * \param   A  pointer to input matrix A
 * \param   B  pointer to input matrix B
 * \param   n  Size of matrices (both sizes)
 * \return     none
 */
void matrixMult_ikj(float * const C, float const * const A, float const * const B, int const n) {
   for (int i = 0; i < n; ++i)
      for (int k = 0; k < n; ++k) {
         if (!k) {
            for (int j = 0; j < n; ++j)
               C[ sub2ind(i,j,n) ] = A[ sub2ind(i,k,n) ] * B[ sub2ind(k,j,n) ];
         } else {
            for (int j = 0; j < n; ++j)
               C[ sub2ind(i,j,n) ] += A[ sub2ind(i,k,n) ] * B[ sub2ind(k,j,n) ];
         }
      }
}

/*!
 * Square Matrix multiplication - jik
 * \param   C  pointer to output matrix
 * \param   A  pointer to input matrix A
 * \param   B  pointer to input matrix B
 * \param   n  Size of matrices (both sizes)
 * \return     none
 */
void matrixMult_jik(float * const C, float const * const A, float const * const B, int const n) {
   for (int j = 0; j < n; j++)
      for (int i = 0; i < n; i++) {
         int k =0;
         C[ sub2ind(i,j,n) ] = A[ sub2ind(i,k,n) ] * B[ sub2ind(k,j,n) ];
         for (k = 1; k < n; k++)
            C[ sub2ind(i,j,n) ] += A[ sub2ind(i,k,n) ] * B[ sub2ind(k,j,n) ];
      }
}

/*!
 * Square Matrix multiplication - jki
 * \param   C  pointer to output matrix
 * \param   A  pointer to input matrix A
 * \param   B  pointer to input matrix B
 * \param   n  Size of matrices (both sizes)
 * \return     none
 */
void matrixMult_jki(float * const C, float const * const A, float const * const B, int const n) {
   for (int j = 0; j < n; ++j)
      for (int k = 0; k < n; ++k) {
         if (!k) {
            for (int i = 0; i < n; ++i)
               C[ sub2ind(i,j,n) ] = A[ sub2ind(i,k,n) ] * B[ sub2ind(k,j,n) ];
         } else {
            for (int i = 0; i < n; ++i)
               C[ sub2ind(i,j,n) ] += A[ sub2ind(i,k,n) ] * B[ sub2ind(k,j,n) ];
         }
      }
}

/*!
 * Square Matrix multiplication - kij
 * \param   C  pointer to output matrix
 * \param   A  pointer to input matrix A
 * \param   B  pointer to input matrix B
 * \param   n  Size of matrices (both sizes)
 * \return     none
 */
void matrixMult_kij(float * const C, float const * const A, float const * const B, int const n) {
   for (int k = 0; k < n; ++k) {
      if (!k) {
         for (int i = 0; i < n; ++i)
            for (int j = 0; j < n; ++j)
               C[ sub2ind(i,j,n) ] = A[ sub2ind(i,k,n) ] * B[ sub2ind(k,j,n) ];
      } else {
         for (int i = 0; i < n; ++i)
            for (int j = 0; j < n; ++j)
               C[ sub2ind(i,j,n) ] += A[ sub2ind(i,k,n) ] * B[ sub2ind(k,j,n) ];
      }
   }
}

/*!
 * Square Matrix multiplication - kji
 * \param   C  pointer to output matrix
 * \param   A  pointer to input matrix A
 * \param   B  pointer to input matrix B
 * \param   n  Size of matrices (both sizes)
 * \return     none
 * xxx
 */
void matrixMult_kji(float * const C, float const * const A, float const * const B, int const n) {
   for (int k = 0; k < n; ++k) {
      if (!k) {
         for (int j = 0; j < n; ++j)
            for (int i = 0; i < n; ++i)
               C[ sub2ind(i,j,n) ] = A[ sub2ind(i,k,n) ] * B[ sub2ind(k,j,n) ];
      } else {
         for (int j = 0; j < n; ++j)
            for (int i = 0; i < n; ++i)
               C[ sub2ind(i,j,n) ] += A[ sub2ind(i,k,n) ] * B[ sub2ind(k,j,n) ];
      }
   }
}


/*!
 * Initialize matrix with random indices and return the matrix pointer.
 *
 * \param   n  The size of the matrix (both of them)
 * \return     Pointer to allocated and initialized matrix
 */
float* matrixInit(int const n) {

   float *M = (float *) malloc( n*n*sizeof(float) );

   for (int i = 0; i < n; i++)               /* rows */
      for (int j = 0; j < n; j++)            /* cols */
         M[ sub2ind(i,j,n) ] = (float)rand()/(float)(RAND_MAX);

   return M;
}