#include <string>
#include <iostream>
#include <iomanip>

#include "JLang/JException.hh"
#include "JMath/JMatrixND.hh"

#include "Jeep/JTimer.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


/**
 * \file
 *
 * Example program to test matrix operations.
 * \author mdejong
 */
int main(int argc, char**argv)
{
  using namespace std;

  int       N;
  int       numberOfEvents;
  int       debug;

  try {

    JParser<> zap("Example program to test matrix operations.");

    zap['N'] = make_field(N);
    zap['n'] = make_field(numberOfEvents);
    zap['d'] = make_field(debug)            = 2;

    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }

  using namespace JPP;
  
 
  JMatrixND A(N);

  const double u = 2.0;
  
  for (int i = 0; i != N; ++i) {

    A(i,i) = u;

    for (int j = 0; j != i; ++j) {
      A(j,i) = A(i,j) = 0.0;
    }
  }

  {
    JTimer t1("standard");
    
    for (int i = 0; i != numberOfEvents; ++i) {

      t1.start();

      JMatrixND B(N);

      B.mul(A,A);

      JMatrixND C(N);

      C.mul(B,A);
      
      JMatrixND D(N);

      D.mul(C,A);
      
      JMatrixND E(N);

      E.mul(D,A);
      
      t1.stop();

      E.div(u*u*u*u*u);

      if (!E.isIdentity()) {
	ERROR("Matrix product error" << endl);
      }
    }

    t1.print(cout);
  }
  {
    JTimer t1("calculator");
    
    JMatrixND B(N);
      
    for (int i = 0; i != numberOfEvents; ++i) {

      t1.start();

      B = A * A * A * A * A;

      t1.stop();

      B.div(u*u*u*u*u);

      if (!B.isIdentity()) {
	ERROR("Matrix product error" << endl);
      }
    }

    t1.print(cout);
  }
}