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;
}
