/*
   file: dgemm.cc
   purpose: packages dgemm.f into a simple use c++ class
   does: matrix multiply
*/

#ifndef maxdim
#define maxdim 130
#endif

extern "C" 
{
/* pound-include headers for pure C functions that are incompatible with c++	  
 go here. eg.:
   pound-include "pure.h"						  */

/* you can also directly include headers or pure c code within the extern {..}
   brackets--if your compiler is lenient enough, as is done here:	  */
void dgemm(char*,char*,int*,int*,int*,double*,double[maxdim][maxdim],int*,
           double[maxdim][maxdim],int*,double*,double[maxdim][maxdim],int*); 	
}

class dgemmSimple	// uses dgemm to multiply square matrixes
{

 public:
 
  dgemmSimple(int dim_);
  int dim;
  void mul(double A[maxdim][maxdim],double B[maxdim][maxdim],
           double C[maxdim][maxdim], double alph=1.);
  
// dgemm arguments:  
  char transa[1];
  char transb[1];
  int m;
  int n;
  int k;
  double alpha;
//double A[][]
  int lda;
//double B[][]
  int ldb;
  double gamma;
//double C[][]
  int ldc;

};
  
 
 
dgemmSimple::dgemmSimple(int d) 
{
  dim=m=n=k=d;
  alpha=1.; gamma=0.;
  ldc=ldb=lda=maxdim;
  transa[0]=transb[0]='n';
}


void dgemmSimple::mul(double A[][maxdim], double B[][maxdim],
                           double C[][maxdim], double alph)
{
 alpha=alph;
 dgemm(transa,transb,&m,&n,&k,&alpha,A,&lda,B,&ldb,&gamma,C,&ldc); 	
}
 
 
 

/* TEST PROG: *************************************************************
*
* main() 
{
int i,j;
dgemmSimple dgemmy(2);
double A[130][130],B[130][130],C[130][130];
for (i=0;i<maxdim;i++) 
  for (j=0;j<maxdim;j++) A[i][j]=B[i][j]=0.;
A[0][0]=1.;
A[1][1]=1.;
B[0][1]=1.;
B[1][0]=1.;
dgemmy.mul(A,B,C);
printf("\n%g   %g   %g   %g\n",C[0][0],C[0][1],C[1][0],C[1][1]);
}
*
***************************************************************************/