Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
TMatrixTCramerInv.cxx
Go to the documentation of this file.
1// @(#)root/base:$Id$
2// Authors: Fons Rademakers, Eddy Offermann Jan 2004
3
4/*************************************************************************
5 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6 * All rights reserved. *
7 * *
8 * For the licensing terms see $ROOTSYS/LICENSE. *
9 * For the list of contributors see $ROOTSYS/README/CREDITS. *
10 *************************************************************************/
11
12/** \class TMatrixTCramerInv
13 \ingroup Matrix
14
15TMatrixTCramerInv
16
17Encapsulate templates of Cramer Inversion routines.
18
19The 4x4, 5x5 and 6x6 are adapted from routines written by
20Mark Fischler and Steven Haywood as part of the CLHEP package
21
22Although for sizes <= 6x6 the Cramer Inversion has a gain in speed
23compared to factorization schemes (like LU) , one pays a price in
24accuracy .
25
26For Example:
27~~~
28 H * H^-1 = U, where H is a 5x5 Hilbert matrix
29 U is a 5x5 Unity matrix
30
31 LU : |U_jk| < 10e-13 for j!=k
32 Cramer: |U_jk| < 10e-7 for j!=k
33~~~
34
35however Cramer algorithm is about 10 (!) times faster
36*/
37
38#include "TMatrixTCramerInv.h"
39#include "TError.h"
40
41#if !defined(R__SOLARIS) && !defined(R__ACC) && !defined(R__FBSD)
43#endif
44
45////////////////////////////////////////////////////////////////////////////////
46
47template<class Element>
49{
50 if (m.GetNrows() != 2 || m.GetNcols() != 2 || m.GetRowLwb() != m.GetColLwb()) {
51 Error("Inv2x2","matrix should be square 2x2");
52 return kFALSE;
53 }
54
55 Element *pM = m.GetMatrixArray();
56
57 const Double_t det = pM[0] * pM[3] - pM[2] * pM[1];
58
59 if (determ)
60 *determ = det;
61
62 const Double_t s = 1./det;
63 if ( det == 0 ) {
64 Error("Inv2x2","matrix is singular");
65 return kFALSE;
66 }
67
68 const Double_t tmp = s*pM[3];
69 pM[1] *= -s;
70 pM[2] *= -s;
71 pM[3] = s*pM[0];
72 pM[0] = tmp;
73
74 return kTRUE;
75}
76
77////////////////////////////////////////////////////////////////////////////////
78
79template<class Element>
81{
82 if (m.GetNrows() != 3 || m.GetNcols() != 3 || m.GetRowLwb() != m.GetColLwb()) {
83 Error("Inv3x3","matrix should be square 3x3");
84 return kFALSE;
85 }
86
87 Element *pM = m.GetMatrixArray();
88
89 const Double_t c00 = pM[4] * pM[8] - pM[5] * pM[7];
90 const Double_t c01 = pM[5] * pM[6] - pM[3] * pM[8];
91 const Double_t c02 = pM[3] * pM[7] - pM[4] * pM[6];
92 const Double_t c10 = pM[7] * pM[2] - pM[8] * pM[1];
93 const Double_t c11 = pM[8] * pM[0] - pM[6] * pM[2];
94 const Double_t c12 = pM[6] * pM[1] - pM[7] * pM[0];
95 const Double_t c20 = pM[1] * pM[5] - pM[2] * pM[4];
96 const Double_t c21 = pM[2] * pM[3] - pM[0] * pM[5];
97 const Double_t c22 = pM[0] * pM[4] - pM[1] * pM[3];
98
99 const Double_t t0 = TMath::Abs(pM[0]);
100 const Double_t t1 = TMath::Abs(pM[3]);
101 const Double_t t2 = TMath::Abs(pM[6]);
104 if (t0 >= t1) {
105 if (t2 >= t0) {
106 tmp = pM[6];
107 det = c12*c01-c11*c02;
108 } else {
109 tmp = pM[0];
110 det = c11*c22-c12*c21;
111 }
112 } else if (t2 >= t1) {
113 tmp = pM[6];
114 det = c12*c01-c11*c02;
115 } else {
116 tmp = pM[3];
117 det = c02*c21-c01*c22;
118 }
119
120 if ( det == 0 || tmp == 0) {
121 Error("Inv3x3","matrix is singular");
122 return kFALSE;
123 }
124
125 const Double_t s = tmp/det;
126 if (determ)
127 *determ = 1./s;
128
129 pM[0] = s*c00;
130 pM[1] = s*c10;
131 pM[2] = s*c20;
132 pM[3] = s*c01;
133 pM[4] = s*c11;
134 pM[5] = s*c21;
135 pM[6] = s*c02;
136 pM[7] = s*c12;
137 pM[8] = s*c22;
138
139 return kTRUE;
140}
141
142// GFij are indices for a 4x4 matrix.
143
144#define GF00 0
145#define GF01 1
146#define GF02 2
147#define GF03 3
148
149#define GF10 4
150#define GF11 5
151#define GF12 6
152#define GF13 7
153
154#define GF20 8
155#define GF21 9
156#define GF22 10
157#define GF23 11
158
159#define GF30 12
160#define GF31 13
161#define GF32 14
162#define GF33 15
163
164////////////////////////////////////////////////////////////////////////////////
165
166template<class Element>
168{
169 if (m.GetNrows() != 4 || m.GetNcols() != 4 || m.GetRowLwb() != m.GetColLwb()) {
170 Error("Inv4x4","matrix should be square 4x4");
171 return kFALSE;
172 }
173
174 Element *pM = m.GetMatrixArray();
175
176 // Find all NECESSARY 2x2 dets: (18 of them)
177
178 const Double_t det2_12_01 = pM[GF10]*pM[GF21] - pM[GF11]*pM[GF20];
179 const Double_t det2_12_02 = pM[GF10]*pM[GF22] - pM[GF12]*pM[GF20];
180 const Double_t det2_12_03 = pM[GF10]*pM[GF23] - pM[GF13]*pM[GF20];
181 const Double_t det2_12_13 = pM[GF11]*pM[GF23] - pM[GF13]*pM[GF21];
182 const Double_t det2_12_23 = pM[GF12]*pM[GF23] - pM[GF13]*pM[GF22];
183 const Double_t det2_12_12 = pM[GF11]*pM[GF22] - pM[GF12]*pM[GF21];
184 const Double_t det2_13_01 = pM[GF10]*pM[GF31] - pM[GF11]*pM[GF30];
185 const Double_t det2_13_02 = pM[GF10]*pM[GF32] - pM[GF12]*pM[GF30];
186 const Double_t det2_13_03 = pM[GF10]*pM[GF33] - pM[GF13]*pM[GF30];
187 const Double_t det2_13_12 = pM[GF11]*pM[GF32] - pM[GF12]*pM[GF31];
188 const Double_t det2_13_13 = pM[GF11]*pM[GF33] - pM[GF13]*pM[GF31];
189 const Double_t det2_13_23 = pM[GF12]*pM[GF33] - pM[GF13]*pM[GF32];
190 const Double_t det2_23_01 = pM[GF20]*pM[GF31] - pM[GF21]*pM[GF30];
191 const Double_t det2_23_02 = pM[GF20]*pM[GF32] - pM[GF22]*pM[GF30];
192 const Double_t det2_23_03 = pM[GF20]*pM[GF33] - pM[GF23]*pM[GF30];
193 const Double_t det2_23_12 = pM[GF21]*pM[GF32] - pM[GF22]*pM[GF31];
194 const Double_t det2_23_13 = pM[GF21]*pM[GF33] - pM[GF23]*pM[GF31];
195 const Double_t det2_23_23 = pM[GF22]*pM[GF33] - pM[GF23]*pM[GF32];
196
197 // Find all NECESSARY 3x3 dets: (16 of them)
198
200 + pM[GF02]*det2_12_01;
202 + pM[GF03]*det2_12_01;
204 + pM[GF03]*det2_12_02;
206 + pM[GF03]*det2_12_12;
208 + pM[GF02]*det2_13_01;
210 + pM[GF03]*det2_13_01;
212 + pM[GF03]*det2_13_02;
214 + pM[GF03]*det2_13_12;
216 + pM[GF02]*det2_23_01;
218 + pM[GF03]*det2_23_01;
220 + pM[GF03]*det2_23_02;
222 + pM[GF03]*det2_23_12;
224 + pM[GF12]*det2_23_01;
226 + pM[GF13]*det2_23_01;
228 + pM[GF13]*det2_23_02;
230 + pM[GF13]*det2_23_12;
231
232 // Find the 4x4 det:
233
236 if (determ)
237 *determ = det;
238
239 if ( det == 0 ) {
240 Error("Inv4x4","matrix is singular");
241 return kFALSE;
242 }
243
244 const Double_t oneOverDet = 1.0/det;
246
251
256
261
266
267 return kTRUE;
268}
269
270// GMij are indices for a 5x5 matrix.
271
272#define GM00 0
273#define GM01 1
274#define GM02 2
275#define GM03 3
276#define GM04 4
277
278#define GM10 5
279#define GM11 6
280#define GM12 7
281#define GM13 8
282#define GM14 9
283
284#define GM20 10
285#define GM21 11
286#define GM22 12
287#define GM23 13
288#define GM24 14
289
290#define GM30 15
291#define GM31 16
292#define GM32 17
293#define GM33 18
294#define GM34 19
295
296#define GM40 20
297#define GM41 21
298#define GM42 22
299#define GM43 23
300#define GM44 24
301
302////////////////////////////////////////////////////////////////////////////////
303
304template<class Element>
306{
307 if (m.GetNrows() != 5 || m.GetNcols() != 5 || m.GetRowLwb() != m.GetColLwb()) {
308 Error("Inv5x5","matrix should be square 5x5");
309 return kFALSE;
310 }
311
312 Element *pM = m.GetMatrixArray();
313
314 // Find all NECESSARY 2x2 dets: (30 of them)
315
316 const Double_t det2_23_01 = pM[GM20]*pM[GM31] - pM[GM21]*pM[GM30];
317 const Double_t det2_23_02 = pM[GM20]*pM[GM32] - pM[GM22]*pM[GM30];
318 const Double_t det2_23_03 = pM[GM20]*pM[GM33] - pM[GM23]*pM[GM30];
319 const Double_t det2_23_04 = pM[GM20]*pM[GM34] - pM[GM24]*pM[GM30];
320 const Double_t det2_23_12 = pM[GM21]*pM[GM32] - pM[GM22]*pM[GM31];
321 const Double_t det2_23_13 = pM[GM21]*pM[GM33] - pM[GM23]*pM[GM31];
322 const Double_t det2_23_14 = pM[GM21]*pM[GM34] - pM[GM24]*pM[GM31];
323 const Double_t det2_23_23 = pM[GM22]*pM[GM33] - pM[GM23]*pM[GM32];
324 const Double_t det2_23_24 = pM[GM22]*pM[GM34] - pM[GM24]*pM[GM32];
325 const Double_t det2_23_34 = pM[GM23]*pM[GM34] - pM[GM24]*pM[GM33];
326 const Double_t det2_24_01 = pM[GM20]*pM[GM41] - pM[GM21]*pM[GM40];
327 const Double_t det2_24_02 = pM[GM20]*pM[GM42] - pM[GM22]*pM[GM40];
328 const Double_t det2_24_03 = pM[GM20]*pM[GM43] - pM[GM23]*pM[GM40];
329 const Double_t det2_24_04 = pM[GM20]*pM[GM44] - pM[GM24]*pM[GM40];
330 const Double_t det2_24_12 = pM[GM21]*pM[GM42] - pM[GM22]*pM[GM41];
331 const Double_t det2_24_13 = pM[GM21]*pM[GM43] - pM[GM23]*pM[GM41];
332 const Double_t det2_24_14 = pM[GM21]*pM[GM44] - pM[GM24]*pM[GM41];
333 const Double_t det2_24_23 = pM[GM22]*pM[GM43] - pM[GM23]*pM[GM42];
334 const Double_t det2_24_24 = pM[GM22]*pM[GM44] - pM[GM24]*pM[GM42];
335 const Double_t det2_24_34 = pM[GM23]*pM[GM44] - pM[GM24]*pM[GM43];
336 const Double_t det2_34_01 = pM[GM30]*pM[GM41] - pM[GM31]*pM[GM40];
337 const Double_t det2_34_02 = pM[GM30]*pM[GM42] - pM[GM32]*pM[GM40];
338 const Double_t det2_34_03 = pM[GM30]*pM[GM43] - pM[GM33]*pM[GM40];
339 const Double_t det2_34_04 = pM[GM30]*pM[GM44] - pM[GM34]*pM[GM40];
340 const Double_t det2_34_12 = pM[GM31]*pM[GM42] - pM[GM32]*pM[GM41];
341 const Double_t det2_34_13 = pM[GM31]*pM[GM43] - pM[GM33]*pM[GM41];
342 const Double_t det2_34_14 = pM[GM31]*pM[GM44] - pM[GM34]*pM[GM41];
343 const Double_t det2_34_23 = pM[GM32]*pM[GM43] - pM[GM33]*pM[GM42];
344 const Double_t det2_34_24 = pM[GM32]*pM[GM44] - pM[GM34]*pM[GM42];
345 const Double_t det2_34_34 = pM[GM33]*pM[GM44] - pM[GM34]*pM[GM43];
346
347 // Find all NECESSARY 3x3 dets: (40 of them)
348
389
390 // Find all NECESSARY 4x4 dets: (25 of them)
391
442
443 // Find the 5x5 det:
444
447 if (determ)
448 *determ = det;
449
450 if ( det == 0 ) {
451 Error("Inv5x5","matrix is singular");
452 return kFALSE;
453 }
454
455 const Double_t oneOverDet = 1.0/det;
457
463
469
475
481
487
488 return kTRUE;
489}
490
491// Aij are indices for a 6x6 matrix.
492
493#define GA00 0
494#define GA01 1
495#define GA02 2
496#define GA03 3
497#define GA04 4
498#define GA05 5
499
500#define GA10 6
501#define GA11 7
502#define GA12 8
503#define GA13 9
504#define GA14 10
505#define GA15 11
506
507#define GA20 12
508#define GA21 13
509#define GA22 14
510#define GA23 15
511#define GA24 16
512#define GA25 17
513
514#define GA30 18
515#define GA31 19
516#define GA32 20
517#define GA33 21
518#define GA34 22
519#define GA35 23
520
521#define GA40 24
522#define GA41 25
523#define GA42 26
524#define GA43 27
525#define GA44 28
526#define GA45 29
527
528#define GA50 30
529#define GA51 31
530#define GA52 32
531#define GA53 33
532#define GA54 34
533#define GA55 35
534
535////////////////////////////////////////////////////////////////////////////////
536
537template<class Element>
539{
540 if (m.GetNrows() != 6 || m.GetNcols() != 6 || m.GetRowLwb() != m.GetColLwb()) {
541 Error("Inv6x6","matrix should be square 6x6");
542 return kFALSE;
543 }
544
545 Element *pM = m.GetMatrixArray();
546
547 // Find all NECESSGARY 2x2 dets: (45 of them)
548
549 const Double_t det2_34_01 = pM[GA30]*pM[GA41] - pM[GA31]*pM[GA40];
550 const Double_t det2_34_02 = pM[GA30]*pM[GA42] - pM[GA32]*pM[GA40];
551 const Double_t det2_34_03 = pM[GA30]*pM[GA43] - pM[GA33]*pM[GA40];
552 const Double_t det2_34_04 = pM[GA30]*pM[GA44] - pM[GA34]*pM[GA40];
553 const Double_t det2_34_05 = pM[GA30]*pM[GA45] - pM[GA35]*pM[GA40];
554 const Double_t det2_34_12 = pM[GA31]*pM[GA42] - pM[GA32]*pM[GA41];
555 const Double_t det2_34_13 = pM[GA31]*pM[GA43] - pM[GA33]*pM[GA41];
556 const Double_t det2_34_14 = pM[GA31]*pM[GA44] - pM[GA34]*pM[GA41];
557 const Double_t det2_34_15 = pM[GA31]*pM[GA45] - pM[GA35]*pM[GA41];
558 const Double_t det2_34_23 = pM[GA32]*pM[GA43] - pM[GA33]*pM[GA42];
559 const Double_t det2_34_24 = pM[GA32]*pM[GA44] - pM[GA34]*pM[GA42];
560 const Double_t det2_34_25 = pM[GA32]*pM[GA45] - pM[GA35]*pM[GA42];
561 const Double_t det2_34_34 = pM[GA33]*pM[GA44] - pM[GA34]*pM[GA43];
562 const Double_t det2_34_35 = pM[GA33]*pM[GA45] - pM[GA35]*pM[GA43];
563 const Double_t det2_34_45 = pM[GA34]*pM[GA45] - pM[GA35]*pM[GA44];
564 const Double_t det2_35_01 = pM[GA30]*pM[GA51] - pM[GA31]*pM[GA50];
565 const Double_t det2_35_02 = pM[GA30]*pM[GA52] - pM[GA32]*pM[GA50];
566 const Double_t det2_35_03 = pM[GA30]*pM[GA53] - pM[GA33]*pM[GA50];
567 const Double_t det2_35_04 = pM[GA30]*pM[GA54] - pM[GA34]*pM[GA50];
568 const Double_t det2_35_05 = pM[GA30]*pM[GA55] - pM[GA35]*pM[GA50];
569 const Double_t det2_35_12 = pM[GA31]*pM[GA52] - pM[GA32]*pM[GA51];
570 const Double_t det2_35_13 = pM[GA31]*pM[GA53] - pM[GA33]*pM[GA51];
571 const Double_t det2_35_14 = pM[GA31]*pM[GA54] - pM[GA34]*pM[GA51];
572 const Double_t det2_35_15 = pM[GA31]*pM[GA55] - pM[GA35]*pM[GA51];
573 const Double_t det2_35_23 = pM[GA32]*pM[GA53] - pM[GA33]*pM[GA52];
574 const Double_t det2_35_24 = pM[GA32]*pM[GA54] - pM[GA34]*pM[GA52];
575 const Double_t det2_35_25 = pM[GA32]*pM[GA55] - pM[GA35]*pM[GA52];
576 const Double_t det2_35_34 = pM[GA33]*pM[GA54] - pM[GA34]*pM[GA53];
577 const Double_t det2_35_35 = pM[GA33]*pM[GA55] - pM[GA35]*pM[GA53];
578 const Double_t det2_35_45 = pM[GA34]*pM[GA55] - pM[GA35]*pM[GA54];
579 const Double_t det2_45_01 = pM[GA40]*pM[GA51] - pM[GA41]*pM[GA50];
580 const Double_t det2_45_02 = pM[GA40]*pM[GA52] - pM[GA42]*pM[GA50];
581 const Double_t det2_45_03 = pM[GA40]*pM[GA53] - pM[GA43]*pM[GA50];
582 const Double_t det2_45_04 = pM[GA40]*pM[GA54] - pM[GA44]*pM[GA50];
583 const Double_t det2_45_05 = pM[GA40]*pM[GA55] - pM[GA45]*pM[GA50];
584 const Double_t det2_45_12 = pM[GA41]*pM[GA52] - pM[GA42]*pM[GA51];
585 const Double_t det2_45_13 = pM[GA41]*pM[GA53] - pM[GA43]*pM[GA51];
586 const Double_t det2_45_14 = pM[GA41]*pM[GA54] - pM[GA44]*pM[GA51];
587 const Double_t det2_45_15 = pM[GA41]*pM[GA55] - pM[GA45]*pM[GA51];
588 const Double_t det2_45_23 = pM[GA42]*pM[GA53] - pM[GA43]*pM[GA52];
589 const Double_t det2_45_24 = pM[GA42]*pM[GA54] - pM[GA44]*pM[GA52];
590 const Double_t det2_45_25 = pM[GA42]*pM[GA55] - pM[GA45]*pM[GA52];
591 const Double_t det2_45_34 = pM[GA43]*pM[GA54] - pM[GA44]*pM[GA53];
592 const Double_t det2_45_35 = pM[GA43]*pM[GA55] - pM[GA45]*pM[GA53];
593 const Double_t det2_45_45 = pM[GA44]*pM[GA55] - pM[GA45]*pM[GA54];
594
595 // Find all NECESSGARY 3x3 dets: (80 of them)
596
677
678 // Find all NECESSGARY 4x4 dets: (75 of them)
679
830
831 // Find all NECESSGARY 5x5 dets: (36 of them)
832
905
906 // Find the determinant
907
910 if (determ)
911 *determ = det;
912
913 if ( det == 0 ) {
914 Error("Inv6x6","matrix is singular");
915 return kFALSE;
916 }
917
918 const Double_t oneOverDet = 1.0/det;
920
927
934
941
948
955
962
963 return kTRUE;
964}
965
966#include "TMatrixFfwd.h"
967
968template Bool_t TMatrixTCramerInv::Inv2x2<Float_t>(TMatrixF&,Double_t*);
969template Bool_t TMatrixTCramerInv::Inv3x3<Float_t>(TMatrixF&,Double_t*);
970template Bool_t TMatrixTCramerInv::Inv4x4<Float_t>(TMatrixF&,Double_t*);
971template Bool_t TMatrixTCramerInv::Inv5x5<Float_t>(TMatrixF&,Double_t*);
972template Bool_t TMatrixTCramerInv::Inv6x6<Float_t>(TMatrixF&,Double_t*);
973
974#include "TMatrixDfwd.h"
975
976template Bool_t TMatrixTCramerInv::Inv2x2<Double_t>(TMatrixD&,Double_t*);
977template Bool_t TMatrixTCramerInv::Inv3x3<Double_t>(TMatrixD&,Double_t*);
978template Bool_t TMatrixTCramerInv::Inv4x4<Double_t>(TMatrixD&,Double_t*);
979template Bool_t TMatrixTCramerInv::Inv5x5<Double_t>(TMatrixD&,Double_t*);
980template Bool_t TMatrixTCramerInv::Inv6x6<Double_t>(TMatrixD&,Double_t*);
bool Bool_t
Definition RtypesCore.h:63
constexpr Bool_t kFALSE
Definition RtypesCore.h:94
double Double_t
Definition RtypesCore.h:59
constexpr Bool_t kTRUE
Definition RtypesCore.h:93
#define NamespaceImp(name)
Definition Rtypes.h:398
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition TError.cxx:185
#define GA22
#define GF30
#define GM21
#define GA15
#define GM30
#define GA45
#define GF13
#define GA00
#define GA43
#define GA25
#define GF32
#define GM22
#define GM14
#define GM40
#define GA03
#define GF10
#define GA41
#define GM02
#define GF22
#define GF02
#define GA23
#define GM42
#define GA40
#define GF03
#define GA05
#define GM11
#define GA32
#define GA34
#define GA33
#define GA30
#define GM31
#define GA52
#define GM03
#define GF21
#define GM01
#define GA42
#define GM33
#define GA35
#define GA13
#define GA21
#define GF11
#define GM41
#define GF23
#define GA20
#define GM32
#define GM10
#define GF31
#define GM00
#define GA12
#define GM04
#define GF20
#define GA55
#define GM34
#define GM20
#define GA02
#define GA50
#define GA44
#define GA51
#define GA11
#define GM44
#define GF01
#define GM13
#define GF00
#define GF33
#define GM12
#define GA04
#define GA54
#define GM43
#define GA14
#define GA10
#define GA01
#define GM23
#define GM24
#define GA53
#define GA31
#define GA24
#define GF12
TMatrixT.
Definition TMatrixT.h:40
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:123
Bool_t Inv6x6(TMatrixT< Element > &m, Double_t *determ)
Bool_t Inv3x3(TMatrixT< Element > &m, Double_t *determ)
Bool_t Inv4x4(TMatrixT< Element > &m, Double_t *determ)
Bool_t Inv2x2(TMatrixT< Element > &m, Double_t *determ)
Bool_t Inv5x5(TMatrixT< Element > &m, Double_t *determ)
TMarker m
Definition textangle.C:8
auto * t1
Definition textangle.C:20