JMatrixND.cc

Go to the documentation of this file.

#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 = i; 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);
  }
}