Matrices Direct Product C GSL routine

This is my amateur routine of Matrices Direct Product or Kronecker Product using GSL. There is no advantage in doing this with GSL library, this is just a exercise.

#include <stdio.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_math.h>

double *produtodireto (int a,int b, double A[a*a],double B[b*b])
     {
 double *AxB;
 AxB = malloc(a*a*b*b*(sizeof *AxB));
 int i, j, k,l;
 gsl_matrix *M = gsl_matrix_alloc (a,a);
 gsl_matrix *L = gsl_matrix_alloc (b,b);
 for (i = 0; i < a; i++)
  {
  for (j = 0; j < a; j++)
   {
             gsl_matrix_set (M, i, j, A[i+a*j]);
   }
  }
 for (i = 0; i < b; i++)
  {
  for (j = 0; j < b; j++)
   {
      gsl_matrix_set (L, i, j, B[i+b*j]);
   }
  }
 gsl_matrix *R = gsl_matrix_alloc (a*b,a*b);
 for (i = 0; i < a; i++)
  {
  for (j = 0; j < a; j++)
   { 
   for (k = 0; k < b; k++)
    {
    for (l = 0; l < b; l++)
     {
     gsl_matrix_set (R, b*i+k-b+2,b*j+l-b+2, gsl_matrix_get (M, i, j)*gsl_matrix_get (L, k, l));
     }
    }
   }
  }
for (i = 0; i < a*b; i++)  
 {/* OUT OF RANGE ERRORRRRR */
 for (j = 0; j < a*b; j++)
  {
  AxB[i+a*b*j]=gsl_matrix_get (R, i, j);
  }
 }

gsl_matrix_free (M);
gsl_matrix_free (L);
gsl_matrix_free (R);
       return AxB;
     }