#include "common.h"
#include "vector.h"

#if PROTOTYPE_ALLOWED
void make_identity_matrix(float *imatrix, int n)
#else
void make_identity_matrix(imatrix, n)
float *imatrix;                        /* must be dimension [n][n] */
int n;
#endif
{
     int i,j;

     for (i=0; i<=n-1; i++)
     {
          for (j=0; j<=n-1; j++)  (*(imatrix + n * i + j)) = 0.0;
          (*(imatrix + n * i + i)) = 1.0;
     }
}

#if PROTOTYPE_ALLOWED
void matrix_mult(float *c, float *a, float *b, int m, int n, int p) /* c=a*b */
#else
void matrix_mult(c, a, b, m, n, p)    /*  c=a*b  */
float *c,       /* m by p */
      *a,       /* m by n */
      *b;       /* n by p */
int   m, n, p;
#endif
{
   int i, j, k;   /* cik = aij bjk */
   double x;

   for (i=0; i<m; i++)
   {
      for (k=0; k<p; k++)
      {
         for (j=0, x=0.0; j<n; j++)  x += *(a+i*n+j) * *(b+j*p+k);
         *(c+i*p+k) = x;
      }
   }
   return;
}