#include <stdio.h>
#include <malloc.h>
#include <math.h>

#define RANDOM   ( 2.0*((((double) random()) / LONG) - 0.5) )

/* convenient constants */
static int i_one = 1;
static double d_one = 1.0;
static double d_minus_one = -1.0;
static double d_zero = 0.0;

/* Level 1 BLAS */
#define DNRM2 dnrm2_
extern double DNRM2(int *,double *,int *);

#define DAXPY daxpy_
extern void DAXPY(int *,double *,double *,int *,double *,int *);

/* Level 2 BLAS */
#define DGER dger_
extern void DGER(int *,int *,double *,double *,int *,double *,int *,
                 double *,int *);

/* Level 3 BLAS */
#define DGEMM dgemm_
extern void DGEMM(char *, char *, int *, int *, int *,
						double *, double *, int *, double *, int
						*, double *, double *, int *);

/*--------------------------------------------------------------------*/

int main(void) {
  int i, m, n, k, sz;
  double LONG = pow(2.0, 31.0) - 1.0;
  double *a, *b, *c, *d, err;

  /* multiply a 10x5 and a 5x20 matrix with dgemm */
  m = 10;
  k = 5;
  n = 20;
  sz = m*n;

  if ((a = (double *) calloc(m*k, sizeof(double))) ==
		NULL) {
	 fprintf(stderr, "Out of memory\n");
	 exit(1);
  }
  if ((b = (double *) calloc(k*n, sizeof(double))) ==
		NULL) {
	 fprintf(stderr, "Out of memory\n");
	 exit(1);
  }
  if ((c = (double *) calloc(sz, sizeof(double))) ==
		NULL) {
	 fprintf(stderr, "Out of memory\n");
	 exit(1);
  }
  if ((d = (double *) calloc(sz, sizeof(double))) ==
		NULL) {
	 fprintf(stderr, "Out of memory\n");
	 exit(1);
  }

  for (i=0; i<m*k; i++) a[i] = RANDOM;
  for (i=0; i<k*n; i++) b[i] = RANDOM;
  for (i=0; i<sz; i++) c[i] = d[i] = 0.0;

  /* compute c = a*b as sum of outer products */
  for (i=0; i<k; i++) DGER(&m, &n, &d_one, a+(i*m), &i_one,
									b+i, &k, c, &m);
  /* compute d = a*b with dgemm */
  DGEMM("N", "N", &m, &n, &k, &d_one, a, &m, b, &k,
		  &d_zero, d, &m);
	 
  /* print out error */
  DAXPY(&sz, &d_minus_one, c, &i_one, d, &i_one);
  err = DNRM2(&sz, d, &i_one);
  printf("Error = %24.16f\n", err);

  free(a);
  free(b);
  free(c);
  free(d);

  return(0);

}