Jpp  18.2.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JMatrixNS.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/JVectorND.hh"
10 #include "JMath/JMatrixNS.hh"
11 #include "JMath/JMathTestkit.hh"
12 
13 #include "Jeep/JTimer.hh"
14 #include "Jeep/JParser.hh"
15 #include "Jeep/JMessage.hh"
16 
17 
18 /**
19  * \file
20  *
21  * Example program to test inversion of symmetric matrix.
22  * \author mdejong
23  */
24 int main(int argc, char**argv)
25 {
26  using namespace std;
27  using namespace JPP;
28 
29  int N;
30  double precision;
31  size_t numberOfInversions;
32  UInt_t seed;
33  int debug;
34 
35  try {
36 
37  JParser<> zap("Example program to test inversion of symmetric matrix.");
38 
39  zap['N'] = make_field(N);
40  zap['e'] = make_field(precision) = 1.0e-6;
41  zap['n'] = make_field(numberOfInversions) = 1;
42  zap['S'] = make_field(seed) = 0;
43  zap['d'] = make_field(debug) = 3;
44 
45  zap(argc, argv);
46  }
47  catch(const exception &error) {
48  FATAL(error.what() << endl);
49  }
50 
51 
52  gRandom->SetSeed(seed);
53 
54  const Double_t xmin = -1.0;
55  const Double_t xmax = +1.0;
56 
57  JMatrixNS A(N);
58 
59  for (int i = 0; i != N; ++i) {
60 
61  A(i,i) = gRandom->Uniform(xmin,xmax);
62 
63  for (int j = 0; j != i; ++j) {
64  A(j,i) = A(i,j) = gRandom->Uniform(xmin,xmax);
65  }
66  }
67 
68 
69  try {
70 
71  DEBUG("Matrix A" << endl);
72  DEBUG(A << endl);
73 
74  JMatrixNS B(A);
75 
76  JTimer timer;
77 
78  timer.start();
79 
80  for (size_t i = 2*(numberOfInversions/2) + 1; i != 0; --i) {
81  B.invert();
82  }
83 
84  timer.stop();
85 
86  NOTICE("Elapsed time (Jpp) [us] " << setw(8) << timer.usec_wall << endl);
87  DEBUG("Matrix A^-1" << endl);
88  DEBUG(B << endl);
89 
90  {
91  TMatrixTSym<double> R(B.size());
92 
93  for (size_t row = 0; row != B.size(); ++row) {
94  for (size_t col = 0; col != B.size(); ++col) {
95  R(row, col) = B(row, col);
96  }
97  }
98 
99  timer.start();
100 
101  for (size_t i = 2*(numberOfInversions/2) + 1; i != 0; --i) {
102  R.InvertFast();
103  }
104 
105  timer.stop();
106 
107  NOTICE("Elapsed time (ROOT) [us] " << setw(8) << timer.usec_wall << endl);
108  }
109 
110  JMatrixNS C(A);
111 
112  C.mul(B);
113 
114  DEBUG("Matrix A x A^-1" << endl);
115  DEBUG(C << endl);
116 
117  NOTICE("A x A^-1 = I ? " << C.isIdentity(precision) << endl);
118 
119  if (!C.isIdentity(precision)) {
120  ERROR("Matrix A x A^-1 /= I" << endl);
121  }
122 
123  ASSERT(C.isIdentity(precision));
124 
125 
126  JVectorND X(N);
127 
128  for (int i = 0; i != N; ++i) {
129  X[i] = gRandom->Uniform(xmin,xmax);
130  }
131 
132  JVectorND Y(N);
133 
134  for (int i = 0; i != N; ++i) {
135  for (int j = 0; j != N; ++j) {
136  Y[i] += B(i,j) * X[j];
137  }
138  }
139 
140  C = A;
141 
142  C.solve(X);
143 
144  for (int i = 0; i != N; ++i) {
145 
146  DEBUG(setw(4) << i << ' ' << FIXED(9,5) << X[i] << ' ' << FIXED(9,5) << Y[i] << endl);
147 
148  ASSERT(fabs(X[i] - Y[i]) <= precision);
149  }
150 
151 
152  const double big = 1.0e5;
153 
154  for (int i = 0; i != N; ++i) {
155 
156  A(i,i) += big;
157 
158  timer.reset();
159  timer.start();
160 
161  B.update(i, big);
162 
163  timer.stop();
164 
165  NOTICE("Elapsed time [us] " << setw(8) << timer.usec_wall << endl);
166  NOTICE("Pivot " << setw(2) << i << '+' << big << ' ' << (JMatrixNS(A).mul(B).isIdentity(precision) ? "okay" : "error") << endl);
167 
168  ASSERT(JMatrixNS(A).mul(B).isIdentity(precision));
169 
170  A(i,i) -= big;
171 
172  B.update(i, -big);
173  }
174 
175  {
176  A.reset();
177 
178  for (int i = 0; i != N; ++i) {
179  for (int j = 0; j != N; ++j) {
180  ASSERT(A(i,j) == 0, "Test reset (" << i << "," << j << ")");
181  }
182  }
183  }
184  }
185  catch (const JException& error) {
186  FATAL(error << endl);
187  }
188 
189  return 0;
190 }
const double xmax
Definition: JQuadrature.cc:24
Utility class to parse command line options.
Definition: JParser.hh:1514
Exceptions.
int main(int argc, char *argv[])
Definition: Main.cc:15
then fatal Wrong number of arguments fi set_variable STRING $argv[1] set_variable DETECTORXY_TXT $WORKDIR $DETECTORXY_TXT tail read X Y CHI2 RMS printf optimum n $X $Y $CHI2 $RMS awk v Y
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
static const double C
Physics constants.
#define ASSERT(A,...)
Assert macro.
Definition: JMessage.hh:90
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1989
#define NOTICE(A)
Definition: JMessage.hh:64
#define ERROR(A)
Definition: JMessage.hh:66
General purpose messaging.
#define FATAL(A)
Definition: JMessage.hh:67
then JCookie sh JDataQuality D $DETECTOR_ID R
Definition: JDataQuality.sh:41
const double xmin
Definition: JQuadrature.cc:23
then usage $script< input file >[option[primary[working directory]]] nWhere option can be N
Definition: JMuonPostfit.sh:40
Utility class to parse command line options.
no fit printf nominal n $STRING awk v X
int j
Definition: JPolint.hh:792
source $JPP_DIR setenv csh $JPP_DIR &dev null eval JShellParser o a A
int debug
debug level
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62