Jpp
 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 
7 #include "JLang/JException.hh"
8 #include "JMath/JMatrixNS.hh"
9 
10 #include "Jeep/JTimer.hh"
11 #include "Jeep/JParser.hh"
12 #include "Jeep/JMessage.hh"
13 
14 
15 /**
16  * \file
17  *
18  * Example program to test inversion of symmetric matrix.
19  * \author mdejong
20  */
21 int main(int argc, char**argv)
22 {
23  using namespace std;
24 
25  int N;
26  double precision;
27  UInt_t seed;
28  int debug;
29 
30  try {
31 
32  JParser<> zap("Example program to test inversion of symmetric matrix.");
33 
34  zap['N'] = make_field(N);
35  zap['e'] = make_field(precision) = 1.0e-6;
36  zap['S'] = make_field(seed) = 0;
37  zap['d'] = make_field(debug) = 3;
38 
39  zap(argc, argv);
40  }
41  catch(const exception &error) {
42  FATAL(error.what() << endl);
43  }
44 
45  using namespace JPP;
46 
47  gRandom->SetSeed(seed);
48 
49  const Double_t xmin = -1.0;
50  const Double_t xmax = +1.0;
51 
52  JMatrixNS A(N);
53 
54  for (int i = 0; i != N; ++i) {
55 
56  A(i,i) = gRandom->Uniform(xmin,xmax);
57 
58  for (int j = i; j != i; ++j) {
59  A(j,i) = A(i,j) = gRandom->Uniform(xmin,xmax);
60  }
61  }
62 
63 
64  try {
65 
66  DEBUG("Matrix A" << endl);
67  DEBUG(A << endl);
68 
69  NOTICE("Determinant A = " << A.getDeterminant() << endl);
70 
71  JMatrixNS B(A);
72 
73  JTimer timer;
74 
75  timer.start();
76 
77  B.invert();
78 
79  timer.stop();
80 
81  NOTICE("Elapsed time [us] " << setw(8) << timer.usec_wall << endl);
82  DEBUG("Matrix A^-1" << endl);
83  DEBUG(B << endl);
84 
85  NOTICE("Determinant A^-1 = " << B.getDeterminant() << endl);
86 
87  JMatrixNS C(A);
88 
89  C.mul(B);
90 
91  DEBUG("Matrix A x A^-1" << endl);
92  DEBUG(C << endl);
93 
94  NOTICE("Determinant (A x A^-1) = " << C.getDeterminant() << endl);
95  NOTICE("Determinant A x Determinant A^-1 = " << A.getDeterminant() * B.getDeterminant() << endl);
96  NOTICE("A x A^-1 = I ? " << C.isIdentity(precision) << endl);
97 
98  if (!C.isIdentity(precision)) {
99  ERROR("Matrix A x A^-1 /= I" << endl);
100  }
101 
102  ASSERT (C.isIdentity(precision));
103 
104  const double big = 1.0e5;
105 
106  for (int i = 0; i != N; ++i) {
107 
108  A(i,i) += big;
109 
110  timer.reset();
111  timer.start();
112 
113  B.update(i, big);
114 
115  timer.stop();
116 
117  NOTICE("Elapsed time [us] " << setw(8) << timer.usec_wall << endl);
118  NOTICE("Pivot " << setw(2) << i << '+' << big << ' ' << (JMatrixNS(A).mul(B).isIdentity(precision) ? "okay" : "error") << endl);
119 
120  ASSERT(JMatrixNS(A).mul(B).isIdentity(precision));
121 
122  A(i,i) -= big;
123 
124  B.update(i, -big);
125  }
126 
127  }
128  catch (const JException& error) {
129  FATAL(error << endl);
130  }
131 
132  return 0;
133 }
Utility class to parse command line options.
Definition: JParser.hh:1410
Exceptions.
#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 NOTICE(A)
Definition: JMessage.hh:62
#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