Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
TMatrixTSymCramerInv.cxx
Go to the documentation of this file.
1// @(#)root/base:$Id$
2// Authors: Fons Rademakers, Eddy Offermann Oct 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 TMatrixTSymCramerInv
13 \ingroup Matrix
14
15 TMatrixTSymCramerInv
16
17 Encapsulate templates of Cramer Inversion routines.
18
19 The 4x4, 5x5 and 6x6 are adapted from routines written by
20 Mark Fischler and Steven Haywood as part of the CLHEP package
21
22 Although for sizes <= 6x6 the Cramer Inversion has a gain in speed
23 compared to factorization schemes (like LU) , one pays a price in
24 accuracy .
25
26 For 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
35 however Cramer algorithm is about 10 (!) times faster
36*/
37
39
40#include "TError.h"
41
42#if !defined(R__SOLARIS) && !defined(R__ACC) && !defined(R__FBSD)
44#endif
45
46////////////////////////////////////////////////////////////////////////////////
47
48template<class Element>
50{
51 if (m.GetNrows() != 2) {
52 Error("Inv2x2","matrix should be square 2x2");
53 return kFALSE;
54 }
55
56 Element *pM = m.GetMatrixArray();
57
58 const Double_t det = pM[0] * pM[3] - pM[1] * pM[1];
59
60 if (determ)
61 *determ = det;
62
63 if ( det == 0 ) {
64 Error("Inv2x2","matrix is singular");
65 return kFALSE;
66 }
67
68 const Double_t tmp1 = pM[3] / det;
69 pM[3] = pM[0] / det;
70 pM[2] = pM[1] = -pM[1] / det;
71 pM[0] = tmp1;
72
73 return kTRUE;
74}
75
76////////////////////////////////////////////////////////////////////////////////
77
78template<class Element>
80{
81 if (m.GetNrows() != 3) {
82 Error("Inv3x3","matrix should be square 3x3");
83 return kFALSE;
84 }
85
86 Element *pM = m.GetMatrixArray();
87
88 const Double_t c00 = pM[4] * pM[8] - pM[5] * pM[5];
89 const Double_t c01 = pM[5] * pM[2] - pM[1] * pM[8];
90 const Double_t c02 = pM[1] * pM[5] - pM[4] * pM[2];
91 const Double_t c11 = pM[8] * pM[0] - pM[2] * pM[2];
92 const Double_t c12 = pM[2] * pM[1] - pM[5] * pM[0];
93 const Double_t c22 = pM[0] * pM[4] - pM[1] * pM[1];
94
95 const Double_t t0 = TMath::Abs(pM[0]);
96 const Double_t t1 = TMath::Abs(pM[1]);
97 const Double_t t2 = TMath::Abs(pM[2]);
98
101
102 if (t0 >= t1) {
103 if (t2 >= t0) {
104 tmp = pM[2];
105 det = c12*c01-c11*c02;
106 } else {
107 tmp = pM[0];
108 det = c11*c22-c12*c12;
109 }
110 } else if (t2 >= t1) {
111 tmp = pM[2];
112 det = c12*c01-c11*c02;
113 } else {
114 tmp = pM[1];
115 det = c02*c12-c01*c22;
116 }
117
118 if ( det == 0 || tmp == 0) {
119 Error("Inv3x3","matrix is singular");
120 return kFALSE;
121 }
122
123 Double_t s = tmp/det;
124 if (determ)
125 *determ = 1./s;
126
127 pM[0] = s*c00;
128 pM[1] = s*c01;
129 pM[2] = s*c02;
130 pM[3] = pM[1];
131 pM[4] = s*c11;
132 pM[5] = s*c12;
133 pM[6] = pM[2];
134 pM[7] = pM[5];
135 pM[8] = s*c22;
136
137 return kTRUE;
138}
139
140// SFij are indices for a 4x4 symmetric matrix.
141
142#define SF00 0
143#define SF01 1
144#define SF02 2
145#define SF03 3
146
147#define SF10 1
148#define SF11 5
149#define SF12 6
150#define SF13 7
151
152#define SF20 2
153#define SF21 6
154#define SF22 10
155#define SF23 11
156
157#define SF30 3
158#define SF31 7
159#define SF32 11
160#define SF33 15
161
162////////////////////////////////////////////////////////////////////////////////
163
164template<class Element>
166{
167 if (m.GetNrows() != 4) {
168 Error("Inv4x4","matrix should be square 4x4");
169 return kFALSE;
170 }
171
172 Element *pM = m.GetMatrixArray();
173
174 // Find all NECESSARY 2x2 dets: (14 of them)
175
176 const Double_t mDet2_12_01 = pM[SF10]*pM[SF21] - pM[SF11]*pM[SF20];
177 const Double_t mDet2_12_02 = pM[SF10]*pM[SF22] - pM[SF12]*pM[SF20];
178 const Double_t mDet2_12_12 = pM[SF11]*pM[SF22] - pM[SF12]*pM[SF21];
179 const Double_t mDet2_13_01 = pM[SF10]*pM[SF31] - pM[SF11]*pM[SF30];
180 const Double_t mDet2_13_02 = pM[SF10]*pM[SF32] - pM[SF12]*pM[SF30];
181 const Double_t mDet2_13_03 = pM[SF10]*pM[SF33] - pM[SF13]*pM[SF30];
182 const Double_t mDet2_13_12 = pM[SF11]*pM[SF32] - pM[SF12]*pM[SF31];
183 const Double_t mDet2_13_13 = pM[SF11]*pM[SF33] - pM[SF13]*pM[SF31];
184 const Double_t mDet2_23_01 = pM[SF20]*pM[SF31] - pM[SF21]*pM[SF30];
185 const Double_t mDet2_23_02 = pM[SF20]*pM[SF32] - pM[SF22]*pM[SF30];
186 const Double_t mDet2_23_03 = pM[SF20]*pM[SF33] - pM[SF23]*pM[SF30];
187 const Double_t mDet2_23_12 = pM[SF21]*pM[SF32] - pM[SF22]*pM[SF31];
188 const Double_t mDet2_23_13 = pM[SF21]*pM[SF33] - pM[SF23]*pM[SF31];
189 const Double_t mDet2_23_23 = pM[SF22]*pM[SF33] - pM[SF23]*pM[SF32];
190
191 // SFind all NECESSSFRY 3x3 dets: (10 of them)
192
213
214 // Find the 4x4 det:
215
218
219 if (determ)
220 *determ = det;
221
222 if ( det == 0 ) {
223 Error("Inv4x4","matrix is singular");
224 return kFALSE;
225 }
226
227 const Double_t oneOverDet = 1.0/det;
229
234
238
241
243
244 for (Int_t irow = 0; irow < 4; irow++) {
245 const Int_t rowOff1 = irow*4;
246 for (Int_t icol = 0; icol < irow; icol++) {
247 const Int_t rowOff2 = icol*4;
249 }
250 }
251
252 return kTRUE;
253}
254
255// Mij are indices for a 5x5 matrix.
256
257#define SM00 0
258#define SM01 1
259#define SM02 2
260#define SM03 3
261#define SM04 4
262
263#define SM10 1
264#define SM11 6
265#define SM12 7
266#define SM13 8
267#define SM14 9
268
269#define SM20 2
270#define SM21 7
271#define SM22 12
272#define SM23 13
273#define SM24 14
274
275#define SM30 3
276#define SM31 8
277#define SM32 13
278#define SM33 18
279#define SM34 19
280
281#define SM40 4
282#define SM41 9
283#define SM42 14
284#define SM43 19
285#define SM44 24
286
287////////////////////////////////////////////////////////////////////////////////
288
289template<class Element>
291{
292 if (m.GetNrows() != 5) {
293 Error("Inv5x5","matrix should be square 5x5");
294 return kFALSE;
295 }
296
297 Element *pM = m.GetMatrixArray();
298
299 // Find all NECESSARY 2x2 dets: (25 of them)
300
301 const Double_t mDet2_23_01 = pM[SM20]*pM[SM31] - pM[SM21]*pM[SM30];
302 const Double_t mDet2_23_02 = pM[SM20]*pM[SM32] - pM[SM22]*pM[SM30];
303 const Double_t mDet2_23_03 = pM[SM20]*pM[SM33] - pM[SM23]*pM[SM30];
304 const Double_t mDet2_23_12 = pM[SM21]*pM[SM32] - pM[SM22]*pM[SM31];
305 const Double_t mDet2_23_13 = pM[SM21]*pM[SM33] - pM[SM23]*pM[SM31];
306 const Double_t mDet2_23_23 = pM[SM22]*pM[SM33] - pM[SM23]*pM[SM32];
307 const Double_t mDet2_24_01 = pM[SM20]*pM[SM41] - pM[SM21]*pM[SM40];
308 const Double_t mDet2_24_02 = pM[SM20]*pM[SM42] - pM[SM22]*pM[SM40];
309 const Double_t mDet2_24_03 = pM[SM20]*pM[SM43] - pM[SM23]*pM[SM40];
310 const Double_t mDet2_24_04 = pM[SM20]*pM[SM44] - pM[SM24]*pM[SM40];
311 const Double_t mDet2_24_12 = pM[SM21]*pM[SM42] - pM[SM22]*pM[SM41];
312 const Double_t mDet2_24_13 = pM[SM21]*pM[SM43] - pM[SM23]*pM[SM41];
313 const Double_t mDet2_24_14 = pM[SM21]*pM[SM44] - pM[SM24]*pM[SM41];
314 const Double_t mDet2_24_23 = pM[SM22]*pM[SM43] - pM[SM23]*pM[SM42];
315 const Double_t mDet2_24_24 = pM[SM22]*pM[SM44] - pM[SM24]*pM[SM42];
316 const Double_t mDet2_34_01 = pM[SM30]*pM[SM41] - pM[SM31]*pM[SM40];
317 const Double_t mDet2_34_02 = pM[SM30]*pM[SM42] - pM[SM32]*pM[SM40];
318 const Double_t mDet2_34_03 = pM[SM30]*pM[SM43] - pM[SM33]*pM[SM40];
319 const Double_t mDet2_34_04 = pM[SM30]*pM[SM44] - pM[SM34]*pM[SM40];
320 const Double_t mDet2_34_12 = pM[SM31]*pM[SM42] - pM[SM32]*pM[SM41];
321 const Double_t mDet2_34_13 = pM[SM31]*pM[SM43] - pM[SM33]*pM[SM41];
322 const Double_t mDet2_34_14 = pM[SM31]*pM[SM44] - pM[SM34]*pM[SM41];
323 const Double_t mDet2_34_23 = pM[SM32]*pM[SM43] - pM[SM33]*pM[SM42];
324 const Double_t mDet2_34_24 = pM[SM32]*pM[SM44] - pM[SM34]*pM[SM42];
325 const Double_t mDet2_34_34 = pM[SM33]*pM[SM44] - pM[SM34]*pM[SM43];
326
327 // Find all NECESSARY 3x3 dets: (30 of them)
328
359
360 // Find all NECESSARY 4x4 dets: (15 of them)
361
392
393 // Find the 5x5 det:
394
397 if (determ)
398 *determ = det;
399
400 if ( det == 0 ) {
401 Error("Inv5x5","matrix is singular");
402 return kFALSE;
403 }
404
405 const Double_t oneOverDet = 1.0/det;
407
413
418
422
425
427
428 for (Int_t irow = 0; irow < 5; irow++) {
429 const Int_t rowOff1 = irow*5;
430 for (Int_t icol = 0; icol < irow; icol++) {
431 const Int_t rowOff2 = icol*5;
433 }
434 }
435
436 return kTRUE;
437}
438
439// Aij are indices for a 6x6 symmetric matrix.
440
441#define SA00 0
442#define SA01 1
443#define SA02 2
444#define SA03 3
445#define SA04 4
446#define SA05 5
447
448#define SA10 1
449#define SA11 7
450#define SA12 8
451#define SA13 9
452#define SA14 10
453#define SA15 11
454
455#define SA20 2
456#define SA21 8
457#define SA22 14
458#define SA23 15
459#define SA24 16
460#define SA25 17
461
462#define SA30 3
463#define SA31 9
464#define SA32 15
465#define SA33 21
466#define SA34 22
467#define SA35 23
468
469#define SA40 4
470#define SA41 10
471#define SA42 16
472#define SA43 22
473#define SA44 28
474#define SA45 29
475
476#define SA50 5
477#define SA51 11
478#define SA52 17
479#define SA53 23
480#define SA54 29
481#define SA55 35
482
483////////////////////////////////////////////////////////////////////////////////
484
485template<class Element>
487{
488 if (m.GetNrows() != 6 || m.GetNcols() != 6 || m.GetRowLwb() != m.GetColLwb()) {
489 Error("Inv6x6","matrix should be square 6x6");
490 return kFALSE;
491 }
492
493 Element *pM = m.GetMatrixArray();
494
495 // Find all NECESSSARY 2x2 dets: (39 of them)
496
497 const Double_t mDet2_34_01 = pM[SA30]*pM[SA41] - pM[SA31]*pM[SA40];
498 const Double_t mDet2_34_02 = pM[SA30]*pM[SA42] - pM[SA32]*pM[SA40];
499 const Double_t mDet2_34_03 = pM[SA30]*pM[SA43] - pM[SA33]*pM[SA40];
500 const Double_t mDet2_34_04 = pM[SA30]*pM[SA44] - pM[SA34]*pM[SA40];
501 const Double_t mDet2_34_12 = pM[SA31]*pM[SA42] - pM[SA32]*pM[SA41];
502 const Double_t mDet2_34_13 = pM[SA31]*pM[SA43] - pM[SA33]*pM[SA41];
503 const Double_t mDet2_34_14 = pM[SA31]*pM[SA44] - pM[SA34]*pM[SA41];
504 const Double_t mDet2_34_23 = pM[SA32]*pM[SA43] - pM[SA33]*pM[SA42];
505 const Double_t mDet2_34_24 = pM[SA32]*pM[SA44] - pM[SA34]*pM[SA42];
506 const Double_t mDet2_34_34 = pM[SA33]*pM[SA44] - pM[SA34]*pM[SA43];
507 const Double_t mDet2_35_01 = pM[SA30]*pM[SA51] - pM[SA31]*pM[SA50];
508 const Double_t mDet2_35_02 = pM[SA30]*pM[SA52] - pM[SA32]*pM[SA50];
509 const Double_t mDet2_35_03 = pM[SA30]*pM[SA53] - pM[SA33]*pM[SA50];
510 const Double_t mDet2_35_04 = pM[SA30]*pM[SA54] - pM[SA34]*pM[SA50];
511 const Double_t mDet2_35_05 = pM[SA30]*pM[SA55] - pM[SA35]*pM[SA50];
512 const Double_t mDet2_35_12 = pM[SA31]*pM[SA52] - pM[SA32]*pM[SA51];
513 const Double_t mDet2_35_13 = pM[SA31]*pM[SA53] - pM[SA33]*pM[SA51];
514 const Double_t mDet2_35_14 = pM[SA31]*pM[SA54] - pM[SA34]*pM[SA51];
515 const Double_t mDet2_35_15 = pM[SA31]*pM[SA55] - pM[SA35]*pM[SA51];
516 const Double_t mDet2_35_23 = pM[SA32]*pM[SA53] - pM[SA33]*pM[SA52];
517 const Double_t mDet2_35_24 = pM[SA32]*pM[SA54] - pM[SA34]*pM[SA52];
518 const Double_t mDet2_35_25 = pM[SA32]*pM[SA55] - pM[SA35]*pM[SA52];
519 const Double_t mDet2_35_34 = pM[SA33]*pM[SA54] - pM[SA34]*pM[SA53];
520 const Double_t mDet2_35_35 = pM[SA33]*pM[SA55] - pM[SA35]*pM[SA53];
521 const Double_t mDet2_45_01 = pM[SA40]*pM[SA51] - pM[SA41]*pM[SA50];
522 const Double_t mDet2_45_02 = pM[SA40]*pM[SA52] - pM[SA42]*pM[SA50];
523 const Double_t mDet2_45_03 = pM[SA40]*pM[SA53] - pM[SA43]*pM[SA50];
524 const Double_t mDet2_45_04 = pM[SA40]*pM[SA54] - pM[SA44]*pM[SA50];
525 const Double_t mDet2_45_05 = pM[SA40]*pM[SA55] - pM[SA45]*pM[SA50];
526 const Double_t mDet2_45_12 = pM[SA41]*pM[SA52] - pM[SA42]*pM[SA51];
527 const Double_t mDet2_45_13 = pM[SA41]*pM[SA53] - pM[SA43]*pM[SA51];
528 const Double_t mDet2_45_14 = pM[SA41]*pM[SA54] - pM[SA44]*pM[SA51];
529 const Double_t mDet2_45_15 = pM[SA41]*pM[SA55] - pM[SA45]*pM[SA51];
530 const Double_t mDet2_45_23 = pM[SA42]*pM[SA53] - pM[SA43]*pM[SA52];
531 const Double_t mDet2_45_24 = pM[SA42]*pM[SA54] - pM[SA44]*pM[SA52];
532 const Double_t mDet2_45_25 = pM[SA42]*pM[SA55] - pM[SA45]*pM[SA52];
533 const Double_t mDet2_45_34 = pM[SA43]*pM[SA54] - pM[SA44]*pM[SA53];
534 const Double_t mDet2_45_35 = pM[SA43]*pM[SA55] - pM[SA45]*pM[SA53];
535 const Double_t mDet2_45_45 = pM[SA44]*pM[SA55] - pM[SA45]*pM[SA54];
536
537 // Find all NECESSSARY 3x3 dets: (65 of them)
538
604
605 // Find all NECESSSARY 4x4 dets: (55 of them)
606
717
718 // Find all NECESSSARY 5x5 dets: (19 of them)
719
762
763 // Find the determinant
764
767
768 if (determ)
769 *determ = det;
770
771 if ( det == 0 ) {
772 Error("Inv6x6","matrix is singular");
773 return kFALSE;
774 }
775
776 const Double_t oneOverDet = 1.0/det;
778
785
791
796
800
803
805
806 for (Int_t irow = 0; irow < 6; irow++) {
807 const Int_t rowOff1 = irow*6;
808 for (Int_t icol = 0; icol < irow; icol++) {
809 const Int_t rowOff2 = icol*6;
811 }
812 }
813
814 return kTRUE;
815}
816
817#include "TMatrixFSymfwd.h"
818
819template Bool_t TMatrixTSymCramerInv::Inv2x2<Float_t>(TMatrixFSym&,Double_t*);
820template Bool_t TMatrixTSymCramerInv::Inv3x3<Float_t>(TMatrixFSym&,Double_t*);
821template Bool_t TMatrixTSymCramerInv::Inv4x4<Float_t>(TMatrixFSym&,Double_t*);
822template Bool_t TMatrixTSymCramerInv::Inv5x5<Float_t>(TMatrixFSym&,Double_t*);
823template Bool_t TMatrixTSymCramerInv::Inv6x6<Float_t>(TMatrixFSym&,Double_t*);
824
825#include "TMatrixDSymfwd.h"
826
827template Bool_t TMatrixTSymCramerInv::Inv2x2<Double_t>(TMatrixDSym&,Double_t*);
828template Bool_t TMatrixTSymCramerInv::Inv3x3<Double_t>(TMatrixDSym&,Double_t*);
829template Bool_t TMatrixTSymCramerInv::Inv4x4<Double_t>(TMatrixDSym&,Double_t*);
830template Bool_t TMatrixTSymCramerInv::Inv5x5<Double_t>(TMatrixDSym&,Double_t*);
831template Bool_t TMatrixTSymCramerInv::Inv6x6<Double_t>(TMatrixDSym&,Double_t*);
bool Bool_t
Definition RtypesCore.h:63
int Int_t
Definition RtypesCore.h:45
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:379
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 SM33
#define SA15
#define SA25
#define SA31
#define SA53
#define SM12
#define SM13
#define SM40
#define SA45
#define SA05
#define SF03
#define SA13
#define SA43
#define SM14
#define SA54
#define SA01
#define SF02
#define SA14
#define SF31
#define SA21
#define SF00
#define SF13
#define SA24
#define SA02
#define SF22
#define SA52
#define SF32
#define SM23
#define SF10
#define SM22
#define SA22
#define SF33
#define SM10
#define SM00
#define SF21
#define SA33
#define SA40
#define SM30
#define SF12
#define SA41
#define SM01
#define SF23
#define SM02
#define SF20
#define SM41
#define SA04
#define SA20
#define SM21
#define SA00
#define SA23
#define SM32
#define SM11
#define SA30
#define SM20
#define SM44
#define SA11
#define SA03
#define SA10
#define SM04
#define SF30
#define SF11
#define SA35
#define SM43
#define SA12
#define SA34
#define SM31
#define SM03
#define SM34
#define SM42
#define SA51
#define SA32
#define SA42
#define SA55
#define SA44
#define SA50
#define SF01
#define SM24
TMatrixTSym.
Definition TMatrixTSym.h:36
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:123
Bool_t Inv6x6(TMatrixTSym< Element > &m, Double_t *determ)
Bool_t Inv3x3(TMatrixTSym< Element > &m, Double_t *determ)
Bool_t Inv4x4(TMatrixTSym< Element > &m, Double_t *determ)
Bool_t Inv2x2(TMatrixTSym< Element > &m, Double_t *determ)
Bool_t Inv5x5(TMatrixTSym< Element > &m, Double_t *determ)
TMarker m
Definition textangle.C:8
auto * t1
Definition textangle.C:20