/****************************************************************************
 * Title  : Multiplication of matrices of arbitrary dimension                *
 * Author : Ben Carpenter (www.bencarpenter.co.uk/software)                  *
 * Date   : 1st December 2006                                                *
 ****************************************************************************/
#include <stdio.h>
#include <stdlib.h>

void c_mat( double ***mat, int m, int n);
void f_mat( double ***mat, int m);
void p_mat( double **mat, int m, int n);
void r_mat( double **mat, int m, int n);

int main()
{
  int i,j,k;
  double **mat1,**mat2,**ans;
  
  char carryon = 'y';
  
  int r1,c1; /* A of r1xc1 */
  int r2,c2; /* B of r2xc2 */
  int r3,c3; /* C of r3xc3 => r1xc2 */
  
  while(carryon=='y') {
    printf("\n  Matrix Multiplication");
    printf("\n  ---------------------\n");
  
    int dims_ok = 0;
    do {
      printf("\n  First matrix dimensions:\n");
      printf("    Number of Rows:    ");
      scanf("%d",&r1);
      printf("    Number of Columns: ");
      scanf("%d",&c1);
      
      printf("\n  Second matrix dimensions:\n");
      printf("    Number of Rows:    ");
      scanf("%d",&r2);
      printf("    Number of Columns: ");
      scanf("%d",&c2);
      
      if(r2==c1) dims_ok = 1;
      else {
        printf("\n  Number of columns in matrix 1 must be the same as");
        printf("\n  the number of rows in matrix 2! Please try again . . .\n");
      }
    } while(dims_ok==0);
    
    /* Dimensions of the answer matrix */
    r3 = r1;
    c3 = c2;
    
    /* Allocate memory for all three matrices */
    c_mat( &mat1, r1, c1 );
    c_mat( &mat2, r2, c2 );
    c_mat( &ans, r3, c3 );
    
    /* Read the matrices to be multiplied */
    printf("\n  First matrix values:\n");
    r_mat( mat1, r1, c1 );
    printf("\n  Second matrix values:\n");
    r_mat( mat2, r2, c2 );
      
    /* Multiply the matrices */
    for(i=0; i<r3; i++)
    {
      for(j=0; j<c3; j++)
      {
        ans[i][j] = 0;
        
        for(k=0; k<r2; k++) ans[i][j] += mat1[i][k]*mat2[k][j] ;
      }
    }
    
    /* Show the answer matrix */
    printf("\n  The answer matrix is:\n");
    p_mat( ans, r3, c3 ); 
    
    /* Clear the memory */
    f_mat( &mat1, r1 );
    f_mat( &mat2, r2 );
    f_mat( &ans, r3 );
    
    fflush(stdin);
    printf("\n  Press `y' to continue, or `n' to exit: ");
    scanf("%c", &carryon);
  }
  return 0;
}

/* Print a matrix of order (rows x cols) */
void p_mat( double **mat, int m, int n)
{
  int i,j;
  printf("\n  %c", (char)218);
  for(i=0; i<n; i++) printf("      ");
  printf(" %c", (char)191);
  
  for(i=0; i<m; i++)
  {
    printf("\n  %c", (char)179);
    for(j=0; j<n; j++) printf("%6g", mat[i][j]);
    printf(" %c", (char)179);
  }
  printf("\n  %c", (char)192);
  for(i=0; i<n; i++) printf("      ");
  printf(" %c\n", (char)217);
}

/* Read a matrix of order (rows x cols) */
void r_mat( double **mat, int rows, int cols)
{
  int i,j;
  
  for(i=0; i<rows; i++)
  {
    printf("  Row %d :\n",i+1);
    for(j=0; j<cols; j++)
    {
      printf("    Column %d : ",j+1);
      scanf("%lf",&mat[i][j]);
    }
  }
  /* Show the matrix just read */
  p_mat(mat, rows, cols);
}

/* Allocate memory for a matrix of given order */
void c_mat(double ***mat, int m, int n)
{
  double **array;
  int i;
  
  array = (double**)malloc(m*sizeof(double*));
  for(i=0; i<m; i++) array[i] = (double*)malloc(n*sizeof(double));
  
  *mat = array;
}

/* free memory of matrix allocated in c_mat() */
void f_mat(double ***mat, int m)
{
  double **array;
  int i;
  
  array = *mat;

  for(i=0; i<m; i++) free(array[i]);
  free(array);
  
  *mat = NULL;
}