Logo ROOT   6.10/09
Reference Guide
testInversion.cxx
Go to the documentation of this file.
1 // stress inversion of matrices with varios methods
2 #include "Math/SMatrix.h"
3 #include "TMatrixTSym.h"
4 #include "TDecompChol.h"
5 #include "TDecompBK.h"
6 
7 #include "TRandom.h"
8 #include <vector>
9 #include <string>
10 #include <iostream>
11 #include <cmath>
12 #include <limits>
13 
14 #include "TStopwatch.h"
15 
16 // matrix size
17 #ifndef N
18 #define N 5
19 #endif
20 
21 bool doSelfTest = true;
22 
23 //timer
24 namespace test {
25 
26 
27 #ifdef REPORT_TIME
28  void reportTime( std::string s, double time);
29 #endif
30 
31  void printTime(TStopwatch & time, std::string s) {
32  int pr = std::cout.precision(8);
33  std::cout << s << "\t" << " time = " << time.RealTime() << "\t(sec)\t"
34  // << time.CpuTime()
35  << std::endl;
36  std::cout.precision(pr);
37  }
38 
39 
40 
41  class Timer {
42 
43  public:
44 
45  Timer(const std::string & s = "") : fName(s), fTime(0)
46  {
47  fWatch.Start();
48  }
49  Timer(double & t, const std::string & s = "") : fName(s), fTime(&t)
50  {
51  fWatch.Start();
52  }
53 
54  ~Timer() {
55  fWatch.Stop();
56  printTime(fWatch,fName);
57 #ifdef REPORT_TIME
58  // report time
59  reportTime(fName, fWatch.RealTime() );
60 #endif
61  if (fTime) *fTime += fWatch.RealTime();
62  }
63 
64 
65  private:
66 
67  std::string fName;
68  double * fTime;
69  TStopwatch fWatch;
70 
71  };
72 }
73 
74 using namespace ROOT::Math;
75 
76 
77 
79 
80 // create matrix
81 template<class M>
82 M * createMatrix() {
83  return new M();
84 }
85 // specialized for TMatrix
86 template<>
87 TMatrixTSym<double> * createMatrix<TMatrixTSym<double> >() {
88  return new TMatrixTSym<double>(N);
89 }
90 
91 //print matrix
92 template<class M>
93 void printMatrix(const M & m) {
94  std::cout << m << std::endl;
95 }
96 template<>
97 void printMatrix<TMatrixTSym<double> >(const TMatrixTSym<double> & m ) {
98  m.Print();
99 }
100 
101 // generate matrices
102 template<class M>
103 void genMatrix(M & m ) {
104  TRandom & r = *gRandom;
105  // generate first diagonal elemets
106  for (int i = 0; i < N; ++i) {
107  double maxVal = i*10000/(N-1) + 1; // max condition is 10^4
108  m(i,i) = r.Uniform(0, maxVal);
109  }
110  for (int i = 0; i < N; ++i) {
111  for (int j = 0; j < i; ++j) {
112  double v = 0.3*std::sqrt( m(i,i) * m(j,j) ); // this makes the matrix pos defined
113  m(i,j) = r.Uniform(0, v);
114  m(j,i) = m(i,j); // needed for TMatrix
115  }
116  }
117 }
118 
119 // generate all matrices
120 template<class M>
121 void generate(std::vector<M*> & v) {
122  int n = v.size();
123  gRandom->SetSeed(111);
124  for (int i = 0; i < n; ++i) {
125  v[i] = createMatrix<M>();
126  genMatrix(*v[i] );
127  }
128 }
129 
130 
131 struct Choleski {};
132 struct BK {};
133 struct QR {};
134 struct Cramer {};
135 struct Default {};
136 
137 template <class M, class Type>
138 struct TestInverter {
139  static bool Inv ( const M & , M & ) { return false;}
140  static bool Inv2 ( M & ) { return false;}
141 };
142 
143 template <>
144 struct TestInverter<SymMatrix, Choleski> {
145  static bool Inv ( const SymMatrix & m, SymMatrix & result ) {
146  int ifail = 0;
147  result = m.InverseChol(ifail);
148  return ifail == 0;
149  }
150  static bool Inv2 ( SymMatrix & m ) {
151  return m.InvertChol();
152  }
153 };
154 
155 template <>
156 struct TestInverter<SymMatrix, BK> {
157  static bool Inv ( const SymMatrix & m, SymMatrix & result ) {
158  int ifail = 0;
159  result = m.Inverse(ifail);
160  return ifail==0;
161  }
162  static bool Inv2 ( SymMatrix & m ) {
163  return m.Invert();
164  }
165 };
166 
167 template <>
168 struct TestInverter<SymMatrix, Cramer> {
169  static bool Inv ( const SymMatrix & m, SymMatrix & result ) {
170  int ifail = 0;
171  result = m.InverseFast(ifail);
172  return ifail==0;
173  }
174  static bool Inv2 ( SymMatrix & m ) {
175  return m.InvertFast();
176  }
177 };
178 
179 #ifdef LATER
180 template <>
181 struct TestInverter<SymMatrix, QR> {
182  static bool Inv ( const SymMatrix & m, SymMatrix & result ) {
183  ROOT::Math::QRDecomposition<double> d;
184  int ifail = 0;
185  result = m.InverseFast(ifail);
186  return ifail==0;
187  }
188 };
189 #endif
190 
191 //TMatrix functions
192 
193 template <>
194 struct TestInverter<TMatrixDSym, Default> {
195  static bool Inv ( const TMatrixDSym & m, TMatrixDSym & result ) {
196  result = m;
197  result.Invert();
198  return true;
199  }
200  static bool Inv2 ( TMatrixDSym & m ) {
201  m.Invert();
202  return true;
203  }
204 };
205 
206 template <>
207 struct TestInverter<TMatrixDSym, Cramer> {
208  static bool Inv ( const TMatrixDSym & m, TMatrixDSym & result ) {
209  result = m;
210  result.InvertFast();
211  return true;
212  }
213  static bool Inv2 ( TMatrixDSym & m ) {
214  m.InvertFast();
215  return true;
216  }
217 };
218 
219 template <>
220 struct TestInverter<TMatrixDSym, Choleski> {
221  static bool Inv ( const TMatrixDSym & m, TMatrixDSym & result ) {
222  TDecompChol chol(m);
223  if (!chol.Decompose() ) return false;
224  chol.Invert(result);
225  return true;
226  }
227 };
228 
229 template <>
230 struct TestInverter<TMatrixDSym, BK> {
231  static bool Inv ( const TMatrixDSym & m, TMatrixDSym & result ) {
232  TDecompBK d(m);
233  if (!d.Decompose() ) return false;
234  d.Invert(result);
235  return true;
236  }
237 };
238 
239 
240 template<class M, class T>
241 double invert( const std::vector<M* > & matlist, double & time,std::string s) {
242  M result = *(matlist.front());
243  test::Timer t(time,s);
244  int nloop = matlist.size();
245  double sum = 0;
246  for (int l = 0; l < nloop; l++)
247  {
248  const M & m = *(matlist[l]);
249  bool ok = TestInverter<M,T>::Inv(m,result);
250  if (!ok) {
251  std::cout << "inv failed for matrix " << l << std::endl;
252  printMatrix<M>( m);
253  return -1;
254  }
255  sum += result(0,1);
256  }
257  return sum;
258 }
259 
260 // invert without copying the matrices (une INv2)
261 
262 template<class M, class T>
263 double invert2( const std::vector<M* > & matlist, double & time,std::string s) {
264 
265  // copy vector of matrices
266  int nloop = matlist.size();
267  std::vector<M *> vmat(nloop);
268  for (int l = 0; l < nloop; l++)
269  {
270  vmat[l] = new M( *matlist[l] );
271  }
272 
273  test::Timer t(time,s);
274  double sum = 0;
275  for (int l = 0; l < nloop; l++)
276  {
277  M & m = *(vmat[l]);
278  bool ok = TestInverter<M,T>::Inv2(m);
279  if (!ok) {
280  std::cout << "inv failed for matrix " << l << std::endl;
281  printMatrix<M>( m);
282  return -1;
283  }
284  sum += m(0,1);
285  }
286  return sum;
287 }
288 
289 bool equal(double d1, double d2, double stol = 10000) {
290  std::cout.precision(12); // tolerance is 1E-12
291  double eps = stol * std::numeric_limits<double>::epsilon();
292  if ( std::abs(d1) > 0 && std::abs(d2) > 0 )
293  return ( std::abs( d1-d2) < eps * std::max(std::abs(d1), std::abs(d2) ) );
294  else if ( d1 == 0 )
295  return std::abs(d2) < eps;
296  else // d2 = 0
297  return std::abs(d1) < eps;
298 }
299 
300 // test matrices symmetric and positive defines
301 bool stressSymPosInversion(int n, bool selftest ) {
302 
303  // test smatrix
304 
305  std::vector<SymMatrix *> v1(n);
306  generate(v1);
307  std::vector<TMatrixDSym *> v2(n);
308  generate(v2);
309 
310 
311  bool iret = true;
312  double time = 0;
313  double s1 = invert<SymMatrix, Choleski> (v1, time,"SMatrix Chol");
314  double s2 = invert<SymMatrix, BK> (v1, time,"SMatrix BK");
315  double s3 = invert<SymMatrix, Cramer> (v1, time,"SMatrix Cram");
316  bool ok = ( equal(s1,s2) && equal(s1,s3) );
317  if (!ok) {
318  std::cout << "result SMatrix choleski " << s1 << " BK " << s2 << " cramer " << s3 << std::endl;
319  std::cerr <<"Error: inversion test for SMatrix FAILED ! " << std::endl;
320  }
321  iret &= ok;
322  std::cout << std::endl;
323 
324  double m1 = invert<TMatrixDSym, Choleski> (v2, time,"TMatrix Chol");
325  double m2 = invert<TMatrixDSym, BK> (v2, time,"TMatrix BK");
326  double m3 = invert<TMatrixDSym, Cramer> (v2, time,"TMatrix Cram");
327  double m4 = invert<TMatrixDSym, Default> (v2, time,"TMatrix Def");
328 
329  ok = ( equal(m1,m2) && equal(m1,m3) && equal(m1,m4) );
330  if (!ok) {
331  std::cout << "result TMatrix choleski " << m1 << " BK " << m2
332  << " cramer " << m3 << " default " << m4 << std::endl;
333  std::cerr <<"Error: inversion test for TMatrix FAILED ! " << std::endl;
334  }
335  iret &= ok;
336  std::cout << std::endl;
337 
338 
339  // test using self inversion
340  if (selftest) {
341  std::cout << "\n - self inversion test \n";
342  double s11 = invert2<SymMatrix, Choleski> (v1, time,"SMatrix Chol");
343  double s12 = invert2<SymMatrix, BK> (v1, time,"SMatrix BK");
344  double s13 = invert2<SymMatrix, Cramer> (v1, time,"SMatrix Cram");
345  ok = ( equal(s11,s12) && equal(s11,s13) );
346  if (!ok) {
347  std::cout << "result SMatrix choleski " << s11 << " BK " << s12 << " cramer " << s13 << std::endl;
348  std::cerr <<"Error: self inversion test for SMatrix FAILED ! " << std::endl;
349  }
350  iret &= ok;
351  std::cout << std::endl;
352 
353  double m13 = invert2<TMatrixDSym, Cramer> (v2, time,"TMatrix Cram");
354  double m14 = invert2<TMatrixDSym, Default> (v2, time,"TMatrix Def");
355  ok = ( equal(m13,m14) );
356  if (!ok) {
357  std::cout << "result TMatrix cramer " << m13 << " default " << m14 << std::endl;
358  std::cerr <<"Error: self inversion test for TMatrix FAILED ! " << std::endl;
359  }
360  iret &= ok;
361  std::cout << std::endl;
362  }
363 
364  return iret;
365 }
366 
367 int testInversion(int n = 100000) {
368  std::cout << "Test Inversion for matrix with N = " << N << std::endl;
369  bool ok = stressSymPosInversion(n, doSelfTest);
370  std::cerr << "Test inversion of positive defined matrix ....... ";
371  if (ok) std::cerr << "OK \n";
372  else std::cerr << "FAILED \n";
373  return (ok) ? 0 : -1;
374 }
375 
376 int main() {
377  return testInversion();
378 }
SMatrix< T, D1, D2, R > InverseChol(int &ifail) const
Invert of a symmetric positive defined Matrix using Choleski decomposition.
Definition: SMatrix.icc:446
static long int sum(long int i)
Definition: Factory.cxx:2162
Double_t RealTime()
Stop the stopwatch (if it is running) and return the realtime (in seconds) passed between the start a...
Definition: TStopwatch.cxx:110
virtual void Print(Option_t *option="") const
Dump this marker with its attributes.
Definition: TMarker.cxx:281
bool equal(double d1, double d2, double stol=10000)
SMatrix< double, N, N, MatRepSym< double, N > > SymMatrix
Definition: test.py:1
bool doSelfTest
bool Invert()
Invert a square Matrix ( this method changes the current matrix).
Definition: SMatrix.icc:406
int main()
bool InvertFast()
Fast Invertion of a square Matrix ( this method changes the current matrix).
Definition: SMatrix.icc:423
bool InvertChol()
Invertion of a symmetric positive defined Matrix using Choleski decomposition.
Definition: SMatrix.icc:440
SMatrix< T, D1, D2, R > InverseFast(int &ifail) const
Invert a square Matrix and returns a new matrix.
Definition: SMatrix.icc:430
double invert(const std::vector< M * > &matlist, double &time, std::string s)
The Bunch-Kaufman diagonal pivoting method decomposes a real symmetric matrix A using.
Definition: TDecompBK.h:26
bool stressSymPosInversion(int n, bool selftest)
double sqrt(double)
SMatrix: a generic fixed size D1 x D2 Matrix class.
This is the base class for the ROOT Random number generators.
Definition: TRandom.h:27
virtual void SetSeed(ULong_t seed=0)
Set the random generator seed.
Definition: TRandom.cxx:568
Bool_t Invert(TMatrixDSym &inv)
For a symmetric matrix A(m,m), its inverse A_inv(m,m) is returned .
TMatrixTSym.
M * createMatrix()
double invert2(const std::vector< M * > &matlist, double &time, std::string s)
void genMatrix(M &m)
Cholesky Decomposition class.
Definition: TDecompChol.h:24
Bool_t Invert(TMatrixDSym &inv)
For a symmetric matrix A(m,m), its inverse A_inv(m,m) is returned .
Definition: TDecompBK.cxx:610
RooCmdArg Timer(Bool_t flag=kTRUE)
TRandom2 r(17)
SVector< double, 2 > v
Definition: Dict.h:5
void printMatrix(const M &m)
TMarker * m
Definition: textangle.C:8
TLine * l
Definition: textangle.C:4
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
SMatrix< T, D1, D2, R > Inverse(int &ifail) const
Invert a square Matrix and returns a new matrix.
Definition: SMatrix.icc:413
REAL epsilon
Definition: triangle.c:617
void printTime(TStopwatch &time, std::string s)
TMatrixTSym< Element > & InvertFast(Double_t *det=0)
Invert the matrix and calculate its determinant.
virtual Double_t Uniform(Double_t x1=1)
Returns a uniform deviate on the interval (0, x1).
Definition: TRandom.cxx:606
virtual Bool_t Decompose()
Matrix A is decomposed in component U so that A = U^T * U If the decomposition succeeds, bit kDecomposed is set , otherwise kSingular.
#define N
TMatrixTSym< Element > & Invert(Double_t *det=0)
Invert the matrix and calculate its determinant Notice that the LU decomposition is used instead of B...
void generate(std::vector< M *> &v)
double result[121]
int testInversion(int n=100000)
virtual Bool_t Decompose()
Matrix A is decomposed in components U and D so that A = U*D*U^T If the decomposition succeeds...
Definition: TDecompBK.cxx:131
const Int_t n
Definition: legend1.C:16
Stopwatch class.
Definition: TStopwatch.h:28