Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JMatrix4S.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <iomanip>
4 
5 #include "TRandom3.h"
6 #include "TMatrixTSym.h"
7 
8 #include "JLang/JException.hh"
9 #include "JMath/JMatrix4S.hh"
10 
11 #include "Jeep/JParser.hh"
12 #include "Jeep/JMessage.hh"
13 #include "Jeep/JTimer.hh"
14 
15 
16 /**
17  * \file
18  *
19  * Example program to test inversion of symmetric matrix.
20  * \author mdejong
21  */
22 int main(int argc, char**argv)
23 {
24  using namespace std;
25 
26  double precision;
27  int numberOfEvents;
28  int debug;
29 
30  try {
31 
32  JParser<> zap("Example program to test inversion of symmetric matrix.");
33 
34  zap['e'] = make_field(precision) = 1.0e-6;
35  zap['n'] = make_field(numberOfEvents) = 1000;
36  zap['d'] = make_field(debug) = 3;
37 
38  zap(argc, argv);
39  }
40  catch(const exception &error) {
41  FATAL(error.what() << endl);
42  }
43 
44  using namespace JPP;
45 
46  ASSERT(numberOfEvents > 0);
47 
48  gRandom->SetSeed(0);
49 
50  const Double_t xmin = -1.0;
51  const Double_t xmax = +1.0;
52 
53  JTimer t1("Jpp");
54  JTimer t2("ROOT");
55 
56  for (int i = 0; i != numberOfEvents; ++i) {
57 
58  if (i%100 == 0) { STATUS(setw(8) << i << '\r'); DEBUG(endl); }
59 
60  JMatrix4D A(gRandom->Uniform(xmin,xmax), 0.0, 0.0, 0.0,
61  gRandom->Uniform(xmin,xmax), gRandom->Uniform(xmin,xmax), 0.0, 0.0,
62  gRandom->Uniform(xmin,xmax), gRandom->Uniform(xmin,xmax), gRandom->Uniform(xmin,xmax), 0.0,
63  gRandom->Uniform(xmin,xmax), gRandom->Uniform(xmin,xmax), gRandom->Uniform(xmin,xmax), gRandom->Uniform(xmin,xmax));
64 
65  A.a01 = A.a10;
66  A.a02 = A.a20;
67  A.a03 = A.a30;
68  A.a12 = A.a21;
69  A.a13 = A.a31;
70  A.a23 = A.a32;
71 
72 
73  TMatrixTSym<double> R(4);
74 
75  R(0,0) = A.a00; R(0,1) = A.a01; R(0,2) = A.a02; R(0,3) = A.a03;
76  R(1,0) = A.a10; R(1,1) = A.a11; R(1,2) = A.a12; R(1,3) = A.a13;
77  R(2,0) = A.a20; R(2,1) = A.a21; R(2,2) = A.a22; R(2,3) = A.a23;
78  R(3,0) = A.a30; R(3,1) = A.a31; R(3,2) = A.a32; R(3,3) = A.a33;
79 
80 
81  DEBUG("Matrix A" << endl);
82  DEBUG(A << endl);
83 
84  DEBUG("Determinant A " << A.getDeterminant() << endl);
85 
86  try {
87 
88  JMatrix4S B(A);
89 
90  t1.start();
91 
92  B.invert();
93 
94  t1.stop();
95 
96  t2.start();
97 
98  R.Invert();
99 
100  t2.stop();
101 
102  DEBUG("Matrix A^-1" << endl);
103  DEBUG(B << endl);
104 
105  DEBUG("Determinant A^-1 = " << B.getDeterminant() << endl);
106 
107  JMatrix4D C(A * B);
108 
109  DEBUG("Matrix A x A^-1" << endl);
110  DEBUG(C << endl);
111 
112  DEBUG("Determinant (A x A^-1) = " << C.getDeterminant() << endl);
113  DEBUG("Determinant A x Determinant A^-1 = " << A.getDeterminant() * B.getDeterminant() << endl);
114  DEBUG("A x A^-1 = I ? " << C.isIdentity(precision) << endl);
115 
116  if (!C.isIdentity(precision)) {
117  ERROR("Matrix A x A^-1 /= I" << endl);
118  }
119 
120  JMatrix4D D = C - JMatrix4D::getIdentity();
121 
122  DEBUG("Matrix D = C - I" << endl);
123  DEBUG(D << endl);
124 
125  ASSERT(C.isIdentity(precision));
126  }
127  catch (const JException& error) {
128  FATAL(error << endl);
129  }
130  }
131  STATUS(endl);
132 
133  if (numberOfEvents > 0) {
134 
135  const double factor = 1.0; // / numberOfEvents;
136 
137  t1.print(cout, factor, micro_t);
138  t2.print(cout, factor, micro_t);
139  }
140 
141  return 0;
142 }
Utility class to parse command line options.
Definition: JParser.hh:1410
Exceptions.
micro
Definition: JScale.hh:29
#define STATUS(A)
Definition: JMessage.hh:61
#define ASSERT(A)
Assert macro.
Definition: JMessage.hh:72
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1836
#define ERROR(A)
Definition: JMessage.hh:64
int debug
debug level
Definition: JSirene.cc:59
General purpose messaging.
#define FATAL(A)
Definition: JMessage.hh:65
Utility class to parse command line options.
static const double C
Speed of light in vacuum [m/ns].
Definition: JConstants.hh:22
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:60
int main(int argc, char *argv[])
Definition: Main.cpp:15