#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include <stdio.h>
#include <gsl/gsl_blas.h>


int main(int argc, char *argv[]){
  double *a,*b,*c;
  gsl_matrix_view A ,B,C;
  int i, j, ndim1,ndim2;
  FILE *fp1,*fp2; 
  
  if(argc < 3) {
    fprintf(stderr,"./compute_multi <A.mat > <B.mat>\n"); 
    return(-1); 
  }
  
  fp1=fopen(argv[1],"r");
  fp2=fopen(argv[2],"r");
  fscanf(fp1,"%d\n",&ndim1);
  fscanf(fp2,"%d\n",&ndim2);
  
  a=(double  *)malloc(ndim1*ndim1*sizeof(double));
  b=(double  *)malloc(ndim1*ndim1*sizeof(double));
  c=(double  *)malloc(ndim1*ndim1*sizeof(double));

  if(ndim1 !=ndim2){
    fprintf(stderr,"Error ! matrix not compatible !\n");
    return(-1); 
  }
  
  for(i=0; i < ndim1;i++){ 
    for(j=0; j < ndim1; j++){
      fscanf(fp1," %lf ",&a[j+ndim1*i]);
      fscanf(fp2," %lf ",&b[j+ndim1*i]);
    }
    fscanf(fp1,"\n"); 
    fscanf(fp2,"\n"); 
  }
  
  A = gsl_matrix_view_array(a, ndim1, ndim1);
  B = gsl_matrix_view_array(b, ndim1, ndim1);
  C = gsl_matrix_view_array(c, ndim1, ndim1);
  
  /* Compute C = A B */
  
  gsl_blas_dgemm (CblasNoTrans,CblasNoTrans,1.0,&A.matrix, &B.matrix,0.0,&C.matrix);
  

 fprintf(stderr,"--------------A-------------------------\n");
 for(i=0; i < ndim1;i++){
    for(j=0; j < ndim1; j++){
      fprintf(stderr," %2.6f ",a[j+ndim1*i]);
    }
    fprintf(stdout,"\n");
  }
 fprintf(stderr,"--------------B-------------------------\n");

 for(i=0; i < ndim1;i++){
    for(j=0; j < ndim1; j++){
      fprintf(stderr," %2.6f ",b[j+ndim1*i]);
    }
    fprintf(stdout,"\n");
  }
 fprintf(stderr,"--------------C=AXB----------------------\n");
 
  fprintf(stdout,"%d\n",ndim1);
  for(i=0; i < ndim1;i++){ 
    for(j=0; j < ndim1; j++){
      fprintf(stdout," %2.6f ",c[j+ndim1*i]);
    }
    fprintf(stdout,"\n"); 
  }
  

 
 return 0;  
}