#include <stdlib.h>
#include "toolsmith.h"
#include "Vector.h"
#include "Matrix.h"


/************************************************************
 * Matrix class member functions                            *
 ************************************************************/

void Matrix::set(int m, int n)
{
   int i;

   if(m < 0 || n < 0)
      error("Matrix::set(int, int)", "Matrix dimensions must be non-negative", ERREXIT);

   M = m;
   N = n;
   p = new Vector [m];
   for(i = 0; i < m; i++) {
      p[i].length(n);
   }
}

void Matrix::clear()
{
   if(M > 0) delete [] p;
}

Matrix::Matrix()
{
   M = N = 0;
}

Matrix::Matrix(int m, int n)
{
   set(m, n);
}

Matrix::Matrix(const Matrix &X)
{
   int i;

   M = X.M;
   N = X.N;
   p = new Vector [M];
   for(i = 0; i < M; i++) {
      p[i] = X.p[i];
   }
}

Matrix::~Matrix()
{
   clear();
}

Vector Matrix::dim() const
{
   Vector w(2);
   
   w[0] = M;
   w[1] = N;
   
   return w;
}

Matrix Matrix::dim(int m, int n)
{
   clear();
   set(m, n);

   return *this;
}

Matrix Matrix::operator [](const Vector &v) const
{
   int i, I, j, m = v.length();
   Matrix Y(N, m);

   for(i = 0; i < m; i++) {
      I = (int)v[i];
      for(j = 0; j < N; j++) {
         Y.p[j][i] = p[I][j];
      }
   }

   return Y;
}

Matrix::operator Vector() const
{
   int i, j, k = 0;
   Vector v(M * N);

   for(j = 0; j < N; j++) {
      for(i = 0; i < M; i++) {
         v[k++] = p[i][j];
      }
   }

   return v;
}

Matrix Matrix::operator !() const
{
   int i;
   Matrix Y(M, N);
   Vector *ap = Y.p, *bp = p;

   for(i = 0; i < M; i++) {
      *ap++ = !(*bp++);
   }

   return Y;
}

Matrix Matrix::operator -() const
{
   int i;
   Matrix Y(M, N);
   Vector *ap = Y.p, *bp = p;

   for(i = 0; i < M; i++) {
      *ap++ = -(*bp++);
   }

   return Y;
}

Matrix operator %(const Vector &u, const Vector &v)
{
   int i, n = u.length();
   double sum = 0.0;
   Matrix Y(1, 1);

   if(n != v.length())
      error("Matrix operator % ", "non-conformable vectors", ERRNCM);

   for(i = 0; i < n; i++) {
      sum += u[i] * v[i];
   }
   Y.p[0][0] = sum;

   return Y;
}

Matrix operator %(const Vector &v, const Matrix &X)
{
   int i, j;
   double sum;
   Matrix Y(1, X.N);

   if(v.length() != X.M)
      error("Matrix operator % ", "non-conformable matrices", ERRNCM);

   for(j = 0; j < X.N; j++) {
      sum = 0.0;
      for(i = 0; i < X.M; i++) {
         sum += v[i] * X.p[i][j];
      }
      Y.p[0][j] = sum;
   }

   return Y;
}

Matrix Matrix::operator %(const Vector &v) const
{
   int i, j;
   double sum;
   Matrix Y(M, 1);

   if(N != v.length())
      error("Matrix::operator % ", "non-conformable matrices", ERRNCM);

   for(i = 0; i < M; i++) {
      sum = 0.0;
      for(j = 0; j < N; j++) {
         sum += p[i][j] * v[j];
      }
      Y.p[i][0] = sum;
   }

   return Y;
}

Matrix Matrix::operator %(const Matrix &X) const
{
   int i, j, k;
   double sum;
   Matrix Y(M, X.N);

   if(N != X.M)
      error("Matrix::operator % ", "non-conformable matrices", ERRNCM);

   for(i = 0; i < M; i++) {
      for(j = 0; j < X.N; j++) {
         sum = 0.0;
         for(k = 0; k < N; k++) {
            sum += p[i][k] * X.p[k][j];
         }
         Y.p[i][j] = sum;
      }
   }

   return Y;
}


#define arithmetic(op)                                                         \
Matrix operator op (double x, const Matrix &X)                                 \
{                                                                              \
   int i;                                                                      \
   Matrix Y(X.M, X.N);                                                         \
   Vector *ap = Y.p, *bp = X.p;                                                \
                                                                               \
   for(i = 0; i < X.M; i++) {                                                  \
      *ap++ = x op *bp++;                                                      \
   }                                                                           \
                                                                               \
   return Y;                                                                   \
}                                                                              \
                                                                               \
Matrix operator op (const Vector &v, const Matrix &X)                          \
{                                                                              \
   int i, j, k = 0, l = v.length();                                            \
   Matrix Y(X.M, X.N);                                                         \
                                                                               \
   if(l == 0 || l > (X.M * X.N))                                               \
      error("Matrix operator " #op, "non-conformable matrices", ERRNCM);       \
                                                                               \
   for(j = 0; j < X.N; j++) {                                                  \
      for(i = 0; i < X.M; i++) {                                               \
         Y.p[i][j] = v[k] op X.p[i][j];                                        \
         k = (k + 1) % l;                                                      \
      }                                                                        \
   }                                                                           \
                                                                               \
   return Y;                                                                   \
}                                                                              \
                                                                               \
Matrix Matrix::operator op (double x) const                                    \
{                                                                              \
   int i;                                                                      \
   Matrix Y(M, N);                                                             \
   Vector *ap = Y.p, *bp = p;                                                  \
                                                                               \
   for(i = 0; i < M; i++) {                                                    \
      *ap++ = *bp++ op x;                                                      \
   }                                                                           \
                                                                               \
   return Y;                                                                   \
}                                                                              \
                                                                               \
Matrix Matrix::operator op (const Vector &v) const                             \
{                                                                              \
   int i, j, k = 0, l = v.length();                                            \
   Matrix Y(M, N);                                                             \
                                                                               \
   if(l == 0 || l > (M * N))                                                   \
      error("Matrix::operator " #op, "non-conformable matrices", ERRNCM);      \
                                                                               \
   for(j = 0; j < N; j++) {                                                    \
      for(i = 0; i < M; i++) {                                                 \
         Y.p[i][j] = p[i][j] op v[k];                                          \
         k = (k + 1) % l;                                                      \
      }                                                                        \
   }                                                                           \
                                                                               \
   return Y;                                                                   \
}                                                                              \
                                                                               \
Matrix Matrix::operator op (const Matrix &X) const                             \
{                                                                              \
   int i;                                                                      \
   Matrix Y(M, N);                                                             \
   Vector *ap = Y.p, *bp = p, *cp = X.p;                                       \
                                                                               \
   if(M != X.M || N != X.N)                                                    \
      error("Matrix::operator " #op, "non-conformable matrices", ERRNCM);      \
                                                                               \
   for(i = 0; i < M; i++) {                                                    \
      *ap++ = *bp++ op *cp++;                                                  \
   }                                                                           \
                                                                               \
   return Y;                                                                   \
}                                                                              \

arithmetic(*);
arithmetic(/);

arithmetic(+);
arithmetic(-);

arithmetic(<);
arithmetic(<=);
arithmetic(>);
arithmetic(>=);

arithmetic(==);
arithmetic(!=);

arithmetic(&&);

arithmetic(||);

#undef arithmetic


Matrix Matrix::operator =(double x)
{
   int i;
   Vector *ap = p;

   for(i = 0; i < M; i++) {
      *ap++ = x;
   }

   return *this;
}

Matrix Matrix::operator =(const Vector &v)
{
   int i, j, k = 0, l = v.length();

   if(l == 0 || l > (M * N))
      error("Matrix::operator = ", "non-conformable matrices", ERRNCM);

   for(j = 0; j < N; j++) {
      for(i = 0; i < M; i++) {
         p[i][j] = v[k];
         k = (k + 1) % l;
      }
   }

   return *this;
}

Matrix Matrix::operator =(const Matrix &X)
{
   int i;
   Vector *ap, *bp = X.p;

   if(M != X.M || N != X.N) {
      clear();
      M = X.M;
      N = X.N;
      p = new Vector [M];
   }
   ap = p;
   for(i = 0; i < M; i++) {
      *ap++ = *bp++;
   }

   return *this;
}


#define assignment(op)                                                         \
Matrix Matrix::operator op (double x)                                          \
{                                                                              \
   int i;                                                                      \
   Vector *ap = p;                                                             \
                                                                               \
   for(i = 0; i < M; i++) {                                                    \
      *ap++ op x;                                                              \
   }                                                                           \
                                                                               \
   return *this;                                                               \
}                                                                              \
                                                                               \
Matrix Matrix::operator op (const Vector &v)                                   \
{                                                                              \
   int i, j, k = 0, l = v.length();                                            \
                                                                               \
   if(l == 0 || l > (M * N))                                                   \
      error("Matrix::operator " #op, "non-conformable matrices", ERRNCM);      \
                                                                               \
   for(j = 0; j < N; j++) {                                                    \
      for(i = 0; i < M; i++) {                                                 \
         p[i][j] op v[k];                                                      \
         k = (k + 1) % l;                                                      \
      }                                                                        \
   }                                                                           \
                                                                               \
   return *this;                                                               \
}                                                                              \
                                                                               \
Matrix Matrix::operator op (const Matrix &X)                                   \
{                                                                              \
   int i;                                                                      \
   Vector *ap = p, *bp = X.p;                                                  \
                                                                               \
   if(M != X.M || N != X.N)                                                    \
      error("Matrix::operator " #op, "non-conformable matrices", ERRNCM);      \
                                                                               \
   for(i = 0; i < M; i++) {                                                    \
      *ap++ op *bp++;                                                          \
   }                                                                           \
                                                                               \
   return *this;                                                               \
}                                                                              \

assignment(+=);
assignment(-=);
assignment(*=);
assignment(/=);

#undef assignment


Matrix Matrix::operator %=(const Vector &v)
{
   if(v.length() != 1 || (M * N) != 1)
      error("Matrix::operator %=", "non-conformable matrices", ERRNCM);

   p[0][0] *= v[0];

   return *this;
}

Matrix Matrix::operator %=(const Matrix &X)
{
   int i, j, k;
   double sum;
   Vector v(N);

   if(N != X.M || X.N != N)
      error("Matrix::operator %=", "non-conformable matrices", ERRNCM);

   for(i = 0; i < M; i++) {
      v = p[i];
      for(j = 0; j < N; j++) {
         sum = 0.0;
         for(k = 0; k < N; k++) {
            sum += v[k] * X.p[k][j];
         }
         p[i][j] = sum;
      }
   }

   return *this;
}


ostream &operator <<(ostream &os, const Matrix &X)
{
    int i;
    Vector *ap = X.p;
    
    if(X.M > 0) {
       os << *ap++;
       for(i = 1; i < X.M; i++) {
          os << endl << *ap++;
       }
    } else {
       os << "NULL";
    }

    return os;
}