/* Copyright (c) 2001-2004 by SoftIntegration, Inc. All Rights Reserved */
/* arrmul.c 
   this case is for the C function which takes
   arguments of double ** (two dimensions), 
   rather than double * (one dimension)
*/
#include <ch.h>
#include <stdio.h>

int arrmul(double **pa, double **pb, double **pc, int n, int m, int r);
EXPORTCH int arrmul_chdl(void *varg) {
  ChInterp_t interp;
  ChVaList_t ap;
  int i;
  double *a, *b, *c;
  double **pa, **pb, **pc; /* pointer to temporary memory */ 
  int ext_a_0, ext_a_1;
  int ext_b_0, ext_b_1;
  int ext_c_0, ext_c_1;
  int retval;

  Ch_VaStart(interp, ap, varg);

  if(!Ch_VaArrayType(interp, ap)) {
      printf("The 1st argument should be an array.\n");
      exit(-1); 
  }
  if(Ch_VaArrayDim(interp, ap) != 2) {        
      printf("Wrong dimension of the 1st argument.\n");
      exit(-1); 
  }
  ext_a_0 = Ch_VaArrayExtent(interp, ap, 0);
  ext_a_1 = Ch_VaArrayExtent(interp, ap, 1);
  a = Ch_VaArg(interp, ap, double*);  

  if(!Ch_VaArrayType(interp, ap)) {
      printf("The 2nd argument should be an array.\n");
      exit(-1); 
  }
  if(Ch_VaArrayDim(interp, ap) != 2) {        
      printf("Wrong dimension of the 2nd argument.\n");
      exit(-1); 
  }
  ext_b_0 = Ch_VaArrayExtent(interp, ap, 0);
  ext_b_1 = Ch_VaArrayExtent(interp, ap, 1);
  b = Ch_VaArg(interp, ap, double*);

  if(!Ch_VaArrayType(interp, ap)) {
      printf("The 3rd argument should be an array.\n");
      exit(-1); 
  }
  if(Ch_VaArrayDim(interp, ap) != 2) { 
      printf("Wrong dimension of the 3rd argument.\n");
      exit(-1); 
  }
  ext_c_0 = Ch_VaArrayExtent(interp, ap, 0);
  ext_c_1 = Ch_VaArrayExtent(interp, ap, 1);
  c = Ch_VaArg(interp, ap, double*);

  /* check if extents of these three arrays match for array multiplication */ 
  if(ext_a_1 == ext_b_0 && ext_a_0 == ext_c_0 && ext_b_1 == ext_c_1) {
      if( !(pa = (double **)malloc(ext_a_0 * sizeof(double*))) ||
          !(pb = (double **)malloc(ext_b_0 * sizeof(double*))) ||
          !(pc = (double **)malloc(ext_c_0 * sizeof(double*))) ) {
         printf("Can't allocate memory.\n");
         exit(-1);
      } 
      for(i = 0; i <= ext_a_0-1; i++) {
          pa[i] = &a[i*ext_a_1]; 
          pc[i] = &c[i*ext_c_1];
      }
      for(i = 0; i <= ext_b_0-1; i++) {
          pb[i] = &b[i*ext_b_1];
      }
      retval = arrmul(pa, pb, pc, ext_a_0, ext_a_1, ext_b_1);
  }
  else {
      printf("The extents of arrays don't match.\n");
      exit(-1);    
  }
  free(pa); free(pb); free(pc);
  Ch_VaEnd(interp, ap);
  return retval;
}

/* c[n][r] = a[n][m] * b[m][r] */
int arrmul(double **pa, double **pb, double **pc, int n, int m, int r) {
  int i, j, k;

  for(i = 0; i <= n-1; i++) {
      for(j = 0; j <= r-1; j++) {
          pc[i][j] = 0; 
          for(k = 0; k <= m-1; k++)
              pc[i][j] += pa[i][k] * pb[k][j];
      }
  }
  return 0;
}
