/*
 *  multiply two m x m square matrices: c = a * b
 *  using blas subroutine dgemm
 */

/* 
 *  DGEMM  performs one of the matrix-matrix operations
 *
 *     C := alpha*op( A )*op( B ) + beta*C,
 *
 *  where  op( X ) is one of
 *
 *     op( X ) = X   or   op( X ) = X',
 *
 *  alpha and beta are scalars, and A, B and C are matrices, with op( A )
 *  an m by k matrix,  op( B )  a  k by n matrix and  C an m by n matrix.
 *
 *  dgemm (TRANA, TRANB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)

 *  Arguments
 *  ==========
 *
 *  TRANA - CHARACTER*1.
 *           On entry, TRANA specifies the form of op( A ) to be used in
 *           the matrix multiplication as follows:
 *
 *              TRANA = 'N' or 'n',  op( A ) = A.
 *
 *              TRANA = 'T' or 't',  op( A ) = A'.
 *
 *              TRANA = 'C' or 'c',  op( A ) = A'.
 *
 *           Unchanged on exit.
 *
 *  TRANB - CHARACTER*1.
 *           On entry, TRANB specifies the form of op( B ) to be used in
 *           the matrix multiplication as follows:
 *
 *              TRANB = 'N' or 'n',  op( B ) = B.
 *
 *              TRANB = 'T' or 't',  op( B ) = B'.
 *
 *              TRANB = 'C' or 'c',  op( B ) = B'.
 *
 *           Unchanged on exit.
 *
 *  M      - INTEGER.
 *           On entry,  M  specifies  the number  of rows  of the  matrix
 *           op( A )  and of the  matrix  C.  M  must  be at least  zero.
 *           Unchanged on exit.
 *
 *  N      - INTEGER.
 *           On entry,  N  specifies the number  of columns of the matrix
 *           op( B ) and the number of columns of the matrix C. N must be
 *           at least zero.
 *           Unchanged on exit.
 *
 *  K      - INTEGER.
 *           On entry,  K  specifies  the number of columns of the matrix
 *           op( A ) and the number of rows of the matrix op( B ). K must
 *           be at least  zero.
 *           Unchanged on exit.
 *
 *  ALPHA  - DOUBLE PRECISION.
 *           On entry, ALPHA specifies the scalar alpha.
 *           Unchanged on exit.
 *
 *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
 *           k  when  TRANA = 'N' or 'n',  and is  m  otherwise.
 *           Before entry with  TRANA = 'N' or 'n',  the leading  m by k
 *           part of the array  A  must contain the matrix  A,  otherwise
 *           the leading  k by m  part of the array  A  must contain  the
 *           matrix A.
 *           Unchanged on exit.
 *
 *  LDA    - INTEGER.
 *           On entry, LDA specifies the first dimension of A as declared
 *           in the calling (sub) program. When  TRANA = 'N' or 'n' then
 *           LDA must be at least  max( 1, m ), otherwise  LDA must be at
 *           least  max( 1, k ).
 *           Unchanged on exit.
 *
 *  B      - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is
 *           n  when  TRANB = 'N' or 'n',  and is  k  otherwise.
 *           Before entry with  TRANB = 'N' or 'n',  the leading  k by n
 *           part of the array  B  must contain the matrix  B,  otherwise
 *           the leading  n by k  part of the array  B  must contain  the
 *           matrix B.
 *           Unchanged on exit.
 *
 *  LDB    - INTEGER.
 *           On entry, LDB specifies the first dimension of B as declared
 *           in the calling (sub) program. When  TRANB = 'N' or 'n' then
 *           LDB must be at least  max( 1, k ), otherwise  LDB must be at
 *           least  max( 1, n ).
 *           Unchanged on exit.
 *
 *  BETA   - DOUBLE PRECISION.
 *           On entry,  BETA  specifies the scalar  beta.  When  BETA  is
 *           supplied as zero then C need not be set on input.
 *           Unchanged on exit.
 *
 *  C      - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
 *           Before entry, the leading  m by n  part of the array  C must
 *           contain the matrix  C,  except when  beta  is zero, in which
 *           case C need not be set on entry.
 *           On exit, the array  C  is overwritten by the  m by n  matrix
 *           ( alpha*op( A )*op( B ) + beta*C ).
 *
 *  LDC    - INTEGER.
 *           On entry, LDC specifies the first dimension of C as declared
 *           in  the  calling  (sub)  program.   LDC  must  be  at  least
 *           max( 1, m ).
 *           Unchanged on exit.
 *
 *  Level 3 Blas routine.
 */

void matmul (int loops, int ldim, int m, double *a, double *b, double *c);

extern void dgemm_ (char *, char *, int *, int *, int *, double *, double *,
    int *, double *, int *, double *, double *, int *);

void matmul (int loops, int ldim, int m, double *a, double *b, double *c)
{
    double alpha = 1.0;
    double beta = 0.0;
    char trans = 'N';
    int l;

    for (l = 0; l < loops; l++)
    {
        dgemm_ (&trans, &trans, &m, &m, &m, &alpha, a, &ldim, b, &ldim, &beta,
            c, &ldim);
    }
}