root/Math/CramerInversion.icc

0001 // @(#)root/smatrix:$Id$
0002 // Authors: L. Moneta    2005
0003
0004
0005 /**********************************************************************
0006  *                                                                    *
0007  * Copyright (c) 2005 , LCG ROOT MathLib Team                         *
0008  *                                                                    *
0009  *                                                                    *
0010  **********************************************************************/
0011 //
0012 // Cramer optimized inversion for matrices up to size 5x5.
0013 // Code from ROOT TMatrixDCramerInv which originates from CLHEP
0014 // (original author Mark Fischler)
0015 //
0016 // Modified by L. Moneta 22/03/07: specialize only until 5x5 (before was up to 6x6)
0017 //  tests show that on 64 machines (like on slc4) it is faster the general method
0018 //
0019
0020 #ifndef ROOT_Math_CramerInversion_icc
0021 #define ROOT_Math_CramerInversion_icc
0022
0023 #ifndef ROOT_Math_Dinv
0024 #error "Do not use CramerInversion.icc directly. #include \"Math/Dinv.h\" instead."
0025 #endif // ROOT_Math_Dinv
0026
0027 #include <cmath>
0028
0029
0030 namespace ROOT {
0031
0032   namespace Math {
0033
0034
0035
0036 //==============================================================================
0037 // Inversion for 3x3 matrices
0038 //==============================================================================
0039
0040 /**
0041    Inversion for a 3x3 matrix
0042  */
0043 template <class MatrixRep>
0044 bool FastInverter<3>::Dinv(MatrixRep & rhs) {
0045
0046   typedef typename MatrixRep::value_type Scalar;
0047
0048   // check matrix sizes ??
0049
0050   //  Scalar * pM = rhs.Array();
0051
0052   const Scalar c00 = rhs[4] * rhs[8] - rhs[5] * rhs[7];
0053   const Scalar c01 = rhs[5] * rhs[6] - rhs[3] * rhs[8];
0054   const Scalar c02 = rhs[3] * rhs[7] - rhs[4] * rhs[6];
0055   const Scalar c10 = rhs[7] * rhs[2] - rhs[8] * rhs[1];
0056   const Scalar c11 = rhs[8] * rhs[0] - rhs[6] * rhs[2];
0057   const Scalar c12 = rhs[6] * rhs[1] - rhs[7] * rhs[0];
0058   const Scalar c20 = rhs[1] * rhs[5] - rhs[2] * rhs[4];
0059   const Scalar c21 = rhs[2] * rhs[3] - rhs[0] * rhs[5];
0060   const Scalar c22 = rhs[0] * rhs[4] - rhs[1] * rhs[3];
0061
0062   const Scalar t0 = std::abs(rhs[0]);
0063   const Scalar t1 = std::abs(rhs[3]);
0064   const Scalar t2 = std::abs(rhs[6]);
0065   Scalar det;
0066   Scalar tmp;
0067   if (t0 >= t1) {
0068     if (t2 >= t0) {
0069     tmp = rhs[6];
0070     det = c12*c01-c11*c02;
0071     } else {
0072       tmp = rhs[0];
0073       det = c11*c22-c12*c21;
0074     }
0075   } else if (t2 >= t1) {
0076     tmp = rhs[6];
0077     det = c12*c01-c11*c02;
0078   } else {
0079     tmp = rhs[3];
0080     det = c02*c21-c01*c22;
0081   }
0082
0083   if ( det == 0 || tmp == 0) {
0084     return false;
0085   }
0086
0087   const Scalar s = tmp/det;
0088
0089 //   if (determ)
0090 //     *determ = 1./s;
0091
0092   rhs[0] = s*c00;
0093   rhs[1] = s*c10;
0094   rhs[2] = s*c20;
0095   rhs[3] = s*c01;
0096   rhs[4] = s*c11;
0097   rhs[5] = s*c21;
0098   rhs[6] = s*c02;
0099   rhs[7] = s*c12;
0100   rhs[8] = s*c22;
0101
0102   return true;
0103 }
0104
0105
0106 //==============================================================================
0107 // Inversion for 4x4 matrices
0108 //==============================================================================
0109 // Fij are indices for a 4x4 matrix.
0110
0111 #define F00 0
0112 #define F01 1
0113 #define F02 2
0114 #define F03 3
0115
0116 #define F10 4
0117 #define F11 5
0118 #define F12 6
0119 #define F13 7
0120
0121 #define F20 8
0122 #define F21 9
0123 #define F22 10
0124 #define F23 11
0125
0126 #define F30 12
0127 #define F31 13
0128 #define F32 14
0129 #define F33 15
0130
0131 /**
0132    Inversion for a 4x4 matrix
0133  */
0134 template <class MatrixRep>
0135 bool FastInverter<4>::Dinv(MatrixRep & rhs) {
0136
0137   typedef typename MatrixRep::value_type Scalar;
0138
0139   // check matrix sizes ??
0140
0141   //  Scalar * pM = rhs.Array();
0142
0143   // Find all NECESSARY 2x2 dets:  (18 of them)
0144
0145   const Scalar det2_12_01 = rhs[F10]*rhs[F21] - rhs[F11]*rhs[F20];
0146   const Scalar det2_12_02 = rhs[F10]*rhs[F22] - rhs[F12]*rhs[F20];
0147   const Scalar det2_12_03 = rhs[F10]*rhs[F23] - rhs[F13]*rhs[F20];
0148   const Scalar det2_12_13 = rhs[F11]*rhs[F23] - rhs[F13]*rhs[F21];
0149   const Scalar det2_12_23 = rhs[F12]*rhs[F23] - rhs[F13]*rhs[F22];
0150   const Scalar det2_12_12 = rhs[F11]*rhs[F22] - rhs[F12]*rhs[F21];
0151   const Scalar det2_13_01 = rhs[F10]*rhs[F31] - rhs[F11]*rhs[F30];
0152   const Scalar det2_13_02 = rhs[F10]*rhs[F32] - rhs[F12]*rhs[F30];
0153   const Scalar det2_13_03 = rhs[F10]*rhs[F33] - rhs[F13]*rhs[F30];
0154   const Scalar det2_13_12 = rhs[F11]*rhs[F32] - rhs[F12]*rhs[F31];
0155   const Scalar det2_13_13 = rhs[F11]*rhs[F33] - rhs[F13]*rhs[F31];
0156   const Scalar det2_13_23 = rhs[F12]*rhs[F33] - rhs[F13]*rhs[F32];
0157   const Scalar det2_23_01 = rhs[F20]*rhs[F31] - rhs[F21]*rhs[F30];
0158   const Scalar det2_23_02 = rhs[F20]*rhs[F32] - rhs[F22]*rhs[F30];
0159   const Scalar det2_23_03 = rhs[F20]*rhs[F33] - rhs[F23]*rhs[F30];
0160   const Scalar det2_23_12 = rhs[F21]*rhs[F32] - rhs[F22]*rhs[F31];
0161   const Scalar det2_23_13 = rhs[F21]*rhs[F33] - rhs[F23]*rhs[F31];
0162   const Scalar det2_23_23 = rhs[F22]*rhs[F33] - rhs[F23]*rhs[F32];
0163
0164   // Find all NECESSARY 3x3 dets:   (16 of them)
0165
0166   const Scalar det3_012_012 = rhs[F00]*det2_12_12 - rhs[F01]*det2_12_02
0167                                 + rhs[F02]*det2_12_01;
0168   const Scalar det3_012_013 = rhs[F00]*det2_12_13 - rhs[F01]*det2_12_03
0169                                 + rhs[F03]*det2_12_01;
0170   const Scalar det3_012_023 = rhs[F00]*det2_12_23 - rhs[F02]*det2_12_03
0171                                 + rhs[F03]*det2_12_02;
0172   const Scalar det3_012_123 = rhs[F01]*det2_12_23 - rhs[F02]*det2_12_13
0173                                 + rhs[F03]*det2_12_12;
0174   const Scalar det3_013_012 = rhs[F00]*det2_13_12 - rhs[F01]*det2_13_02
0175                                 + rhs[F02]*det2_13_01;
0176   const Scalar det3_013_013 = rhs[F00]*det2_13_13 - rhs[F01]*det2_13_03
0177                                 + rhs[F03]*det2_13_01;
0178   const Scalar det3_013_023 = rhs[F00]*det2_13_23 - rhs[F02]*det2_13_03
0179                                 + rhs[F03]*det2_13_02;
0180   const Scalar det3_013_123 = rhs[F01]*det2_13_23 - rhs[F02]*det2_13_13
0181                                 + rhs[F03]*det2_13_12;
0182   const Scalar det3_023_012 = rhs[F00]*det2_23_12 - rhs[F01]*det2_23_02
0183                                 + rhs[F02]*det2_23_01;
0184   const Scalar det3_023_013 = rhs[F00]*det2_23_13 - rhs[F01]*det2_23_03
0185                                 + rhs[F03]*det2_23_01;
0186   const Scalar det3_023_023 = rhs[F00]*det2_23_23 - rhs[F02]*det2_23_03
0187                                 + rhs[F03]*det2_23_02;
0188   const Scalar det3_023_123 = rhs[F01]*det2_23_23 - rhs[F02]*det2_23_13
0189                                 + rhs[F03]*det2_23_12;
0190   const Scalar det3_123_012 = rhs[F10]*det2_23_12 - rhs[F11]*det2_23_02
0191                                 + rhs[F12]*det2_23_01;
0192   const Scalar det3_123_013 = rhs[F10]*det2_23_13 - rhs[F11]*det2_23_03
0193                                 + rhs[F13]*det2_23_01;
0194   const Scalar det3_123_023 = rhs[F10]*det2_23_23 - rhs[F12]*det2_23_03
0195                                 + rhs[F13]*det2_23_02;
0196   const Scalar det3_123_123 = rhs[F11]*det2_23_23 - rhs[F12]*det2_23_13
0197                                 + rhs[F13]*det2_23_12;
0198
0199   // Find the 4x4 det:
0200
0201   const Scalar det = rhs[F00]*det3_123_123 - rhs[F01]*det3_123_023
0202                        + rhs[F02]*det3_123_013 - rhs[F03]*det3_123_012;
0203
0204 //   if (determ)
0205 //     *determ = det;
0206
0207   if ( det == 0 ) {
0208     return false;
0209   }
0210
0211   // use 1.0f to remove warning C4244 on Windows when using float
0212   const Scalar oneOverDet = 1.0f / det;
0213   const Scalar mn1OverDet = - oneOverDet;
0214
0215   rhs[F00] =  det3_123_123 * oneOverDet;
0216   rhs[F01] =  det3_023_123 * mn1OverDet;
0217   rhs[F02] =  det3_013_123 * oneOverDet;
0218   rhs[F03] =  det3_012_123 * mn1OverDet;
0219
0220   rhs[F10] =  det3_123_023 * mn1OverDet;
0221   rhs[F11] =  det3_023_023 * oneOverDet;
0222   rhs[F12] =  det3_013_023 * mn1OverDet;
0223   rhs[F13] =  det3_012_023 * oneOverDet;
0224
0225   rhs[F20] =  det3_123_013 * oneOverDet;
0226   rhs[F21] =  det3_023_013 * mn1OverDet;
0227   rhs[F22] =  det3_013_013 * oneOverDet;
0228   rhs[F23] =  det3_012_013 * mn1OverDet;
0229
0230   rhs[F30] =  det3_123_012 * mn1OverDet;
0231   rhs[F31] =  det3_023_012 * oneOverDet;
0232   rhs[F32] =  det3_013_012 * mn1OverDet;
0233   rhs[F33] =  det3_012_012 * oneOverDet;
0234
0235   return true;
0236 }
0237
0238 //==============================================================================
0239 // Inversion for 5x5 matrices
0240 //==============================================================================
0241 // Mij are indices for a 5x5 matrix.
0242 #define M00 0
0243 #define M01 1
0244 #define M02 2
0245 #define M03 3
0246 #define M04 4
0247
0248 #define M10 5
0249 #define M11 6
0250 #define M12 7
0251 #define M13 8
0252 #define M14 9
0253
0254 #define M20 10
0255 #define M21 11
0256 #define M22 12
0257 #define M23 13
0258 #define M24 14
0259
0260 #define M30 15
0261 #define M31 16
0262 #define M32 17
0263 #define M33 18
0264 #define M34 19
0265
0266 #define M40 20
0267 #define M41 21
0268 #define M42 22
0269 #define M43 23
0270 #define M44 24
0271
0272
0273 /**
0274    Inversion for a 5x5 matrix
0275  */
0276 template <class MatrixRep>
0277 bool FastInverter<5>::Dinv(MatrixRep & rhs) {
0278
0279   typedef typename MatrixRep::value_type Scalar;
0280
0281   // check matrix sizes ??
0282
0283   //  Scalar * pM = rhs.Array();
0284
0285
0286   // Find all NECESSARY 2x2 dets:  (30 of them)
0287
0288   const Scalar det2_23_01 = rhs[M20]*rhs[M31] - rhs[M21]*rhs[M30];
0289   const Scalar det2_23_02 = rhs[M20]*rhs[M32] - rhs[M22]*rhs[M30];
0290   const Scalar det2_23_03 = rhs[M20]*rhs[M33] - rhs[M23]*rhs[M30];
0291   const Scalar det2_23_04 = rhs[M20]*rhs[M34] - rhs[M24]*rhs[M30];
0292   const Scalar det2_23_12 = rhs[M21]*rhs[M32] - rhs[M22]*rhs[M31];
0293   const Scalar det2_23_13 = rhs[M21]*rhs[M33] - rhs[M23]*rhs[M31];
0294   const Scalar det2_23_14 = rhs[M21]*rhs[M34] - rhs[M24]*rhs[M31];
0295   const Scalar det2_23_23 = rhs[M22]*rhs[M33] - rhs[M23]*rhs[M32];
0296   const Scalar det2_23_24 = rhs[M22]*rhs[M34] - rhs[M24]*rhs[M32];
0297   const Scalar det2_23_34 = rhs[M23]*rhs[M34] - rhs[M24]*rhs[M33];
0298   const Scalar det2_24_01 = rhs[M20]*rhs[M41] - rhs[M21]*rhs[M40];
0299   const Scalar det2_24_02 = rhs[M20]*rhs[M42] - rhs[M22]*rhs[M40];
0300   const Scalar det2_24_03 = rhs[M20]*rhs[M43] - rhs[M23]*rhs[M40];
0301   const Scalar det2_24_04 = rhs[M20]*rhs[M44] - rhs[M24]*rhs[M40];
0302   const Scalar det2_24_12 = rhs[M21]*rhs[M42] - rhs[M22]*rhs[M41];
0303   const Scalar det2_24_13 = rhs[M21]*rhs[M43] - rhs[M23]*rhs[M41];
0304   const Scalar det2_24_14 = rhs[M21]*rhs[M44] - rhs[M24]*rhs[M41];
0305   const Scalar det2_24_23 = rhs[M22]*rhs[M43] - rhs[M23]*rhs[M42];
0306   const Scalar det2_24_24 = rhs[M22]*rhs[M44] - rhs[M24]*rhs[M42];
0307   const Scalar det2_24_34 = rhs[M23]*rhs[M44] - rhs[M24]*rhs[M43];
0308   const Scalar det2_34_01 = rhs[M30]*rhs[M41] - rhs[M31]*rhs[M40];
0309   const Scalar det2_34_02 = rhs[M30]*rhs[M42] - rhs[M32]*rhs[M40];
0310   const Scalar det2_34_03 = rhs[M30]*rhs[M43] - rhs[M33]*rhs[M40];
0311   const Scalar det2_34_04 = rhs[M30]*rhs[M44] - rhs[M34]*rhs[M40];
0312   const Scalar det2_34_12 = rhs[M31]*rhs[M42] - rhs[M32]*rhs[M41];
0313   const Scalar det2_34_13 = rhs[M31]*rhs[M43] - rhs[M33]*rhs[M41];
0314   const Scalar det2_34_14 = rhs[M31]*rhs[M44] - rhs[M34]*rhs[M41];
0315   const Scalar det2_34_23 = rhs[M32]*rhs[M43] - rhs[M33]*rhs[M42];
0316   const Scalar det2_34_24 = rhs[M32]*rhs[M44] - rhs[M34]*rhs[M42];
0317   const Scalar det2_34_34 = rhs[M33]*rhs[M44] - rhs[M34]*rhs[M43];
0318
0319   // Find all NECESSARY 3x3 dets:   (40 of them)
0320
0321   const Scalar det3_123_012 = rhs[M10]*det2_23_12 - rhs[M11]*det2_23_02 + rhs[M12]*det2_23_01;
0322   const Scalar det3_123_013 = rhs[M10]*det2_23_13 - rhs[M11]*det2_23_03 + rhs[M13]*det2_23_01;
0323   const Scalar det3_123_014 = rhs[M10]*det2_23_14 - rhs[M11]*det2_23_04 + rhs[M14]*det2_23_01;
0324   const Scalar det3_123_023 = rhs[M10]*det2_23_23 - rhs[M12]*det2_23_03 + rhs[M13]*det2_23_02;
0325   const Scalar det3_123_024 = rhs[M10]*det2_23_24 - rhs[M12]*det2_23_04 + rhs[M14]*det2_23_02;
0326   const Scalar det3_123_034 = rhs[M10]*det2_23_34 - rhs[M13]*det2_23_04 + rhs[M14]*det2_23_03;
0327   const Scalar det3_123_123 = rhs[M11]*det2_23_23 - rhs[M12]*det2_23_13 + rhs[M13]*det2_23_12;
0328   const Scalar det3_123_124 = rhs[M11]*det2_23_24 - rhs[M12]*det2_23_14 + rhs[M14]*det2_23_12;
0329   const Scalar det3_123_134 = rhs[M11]*det2_23_34 - rhs[M13]*det2_23_14 + rhs[M14]*det2_23_13;
0330   const Scalar det3_123_234 = rhs[M12]*det2_23_34 - rhs[M13]*det2_23_24 + rhs[M14]*det2_23_23;
0331   const Scalar det3_124_012 = rhs[M10]*det2_24_12 - rhs[M11]*det2_24_02 + rhs[M12]*det2_24_01;
0332   const Scalar det3_124_013 = rhs[M10]*det2_24_13 - rhs[M11]*det2_24_03 + rhs[M13]*det2_24_01;
0333   const Scalar det3_124_014 = rhs[M10]*det2_24_14 - rhs[M11]*det2_24_04 + rhs[M14]*det2_24_01;
0334   const Scalar det3_124_023 = rhs[M10]*det2_24_23 - rhs[M12]*det2_24_03 + rhs[M13]*det2_24_02;
0335   const Scalar det3_124_024 = rhs[M10]*det2_24_24 - rhs[M12]*det2_24_04 + rhs[M14]*det2_24_02;
0336   const Scalar det3_124_034 = rhs[M10]*det2_24_34 - rhs[M13]*det2_24_04 + rhs[M14]*det2_24_03;
0337   const Scalar det3_124_123 = rhs[M11]*det2_24_23 - rhs[M12]*det2_24_13 + rhs[M13]*det2_24_12;
0338   const Scalar det3_124_124 = rhs[M11]*det2_24_24 - rhs[M12]*det2_24_14 + rhs[M14]*det2_24_12;
0339   const Scalar det3_124_134 = rhs[M11]*det2_24_34 - rhs[M13]*det2_24_14 + rhs[M14]*det2_24_13;
0340   const Scalar det3_124_234 = rhs[M12]*det2_24_34 - rhs[M13]*det2_24_24 + rhs[M14]*det2_24_23;
0341   const Scalar det3_134_012 = rhs[M10]*det2_34_12 - rhs[M11]*det2_34_02 + rhs[M12]*det2_34_01;
0342   const Scalar det3_134_013 = rhs[M10]*det2_34_13 - rhs[M11]*det2_34_03 + rhs[M13]*det2_34_01;
0343   const Scalar det3_134_014 = rhs[M10]*det2_34_14 - rhs[M11]*det2_34_04 + rhs[M14]*det2_34_01;
0344   const Scalar det3_134_023 = rhs[M10]*det2_34_23 - rhs[M12]*det2_34_03 + rhs[M13]*det2_34_02;
0345   const Scalar det3_134_024 = rhs[M10]*det2_34_24 - rhs[M12]*det2_34_04 + rhs[M14]*det2_34_02;
0346   const Scalar det3_134_034 = rhs[M10]*det2_34_34 - rhs[M13]*det2_34_04 + rhs[M14]*det2_34_03;
0347   const Scalar det3_134_123 = rhs[M11]*det2_34_23 - rhs[M12]*det2_34_13 + rhs[M13]*det2_34_12;
0348   const Scalar det3_134_124 = rhs[M11]*det2_34_24 - rhs[M12]*det2_34_14 + rhs[M14]*det2_34_12;
0349   const Scalar det3_134_134 = rhs[M11]*det2_34_34 - rhs[M13]*det2_34_14 + rhs[M14]*det2_34_13;
0350   const Scalar det3_134_234 = rhs[M12]*det2_34_34 - rhs[M13]*det2_34_24 + rhs[M14]*det2_34_23;
0351   const Scalar det3_234_012 = rhs[M20]*det2_34_12 - rhs[M21]*det2_34_02 + rhs[M22]*det2_34_01;
0352   const Scalar det3_234_013 = rhs[M20]*det2_34_13 - rhs[M21]*det2_34_03 + rhs[M23]*det2_34_01;
0353   const Scalar det3_234_014 = rhs[M20]*det2_34_14 - rhs[M21]*det2_34_04 + rhs[M24]*det2_34_01;
0354   const Scalar det3_234_023 = rhs[M20]*det2_34_23 - rhs[M22]*det2_34_03 + rhs[M23]*det2_34_02;
0355   const Scalar det3_234_024 = rhs[M20]*det2_34_24 - rhs[M22]*det2_34_04 + rhs[M24]*det2_34_02;
0356   const Scalar det3_234_034 = rhs[M20]*det2_34_34 - rhs[M23]*det2_34_04 + rhs[M24]*det2_34_03;
0357   const Scalar det3_234_123 = rhs[M21]*det2_34_23 - rhs[M22]*det2_34_13 + rhs[M23]*det2_34_12;
0358   const Scalar det3_234_124 = rhs[M21]*det2_34_24 - rhs[M22]*det2_34_14 + rhs[M24]*det2_34_12;
0359   const Scalar det3_234_134 = rhs[M21]*det2_34_34 - rhs[M23]*det2_34_14 + rhs[M24]*det2_34_13;
0360   const Scalar det3_234_234 = rhs[M22]*det2_34_34 - rhs[M23]*det2_34_24 + rhs[M24]*det2_34_23;
0361
0362   // Find all NECESSARY 4x4 dets:   (25 of them)
0363
0364   const Scalar det4_0123_0123 = rhs[M00]*det3_123_123 - rhs[M01]*det3_123_023
0365                                   + rhs[M02]*det3_123_013 - rhs[M03]*det3_123_012;
0366   const Scalar det4_0123_0124 = rhs[M00]*det3_123_124 - rhs[M01]*det3_123_024
0367                                   + rhs[M02]*det3_123_014 - rhs[M04]*det3_123_012;
0368   const Scalar det4_0123_0134 = rhs[M00]*det3_123_134 - rhs[M01]*det3_123_034
0369                                   + rhs[M03]*det3_123_014 - rhs[M04]*det3_123_013;
0370   const Scalar det4_0123_0234 = rhs[M00]*det3_123_234 - rhs[M02]*det3_123_034
0371                                   + rhs[M03]*det3_123_024 - rhs[M04]*det3_123_023;
0372   const Scalar det4_0123_1234 = rhs[M01]*det3_123_234 - rhs[M02]*det3_123_134
0373                                   + rhs[M03]*det3_123_124 - rhs[M04]*det3_123_123;
0374   const Scalar det4_0124_0123 = rhs[M00]*det3_124_123 - rhs[M01]*det3_124_023
0375                                   + rhs[M02]*det3_124_013 - rhs[M03]*det3_124_012;
0376   const Scalar det4_0124_0124 = rhs[M00]*det3_124_124 - rhs[M01]*det3_124_024
0377                                   + rhs[M02]*det3_124_014 - rhs[M04]*det3_124_012;
0378   const Scalar det4_0124_0134 = rhs[M00]*det3_124_134 - rhs[M01]*det3_124_034
0379                                   + rhs[M03]*det3_124_014 - rhs[M04]*det3_124_013;
0380   const Scalar det4_0124_0234 = rhs[M00]*det3_124_234 - rhs[M02]*det3_124_034
0381                                   + rhs[M03]*det3_124_024 - rhs[M04]*det3_124_023;
0382   const Scalar det4_0124_1234 = rhs[M01]*det3_124_234 - rhs[M02]*det3_124_134
0383                                   + rhs[M03]*det3_124_124 - rhs[M04]*det3_124_123;
0384   const Scalar det4_0134_0123 = rhs[M00]*det3_134_123 - rhs[M01]*det3_134_023
0385                                   + rhs[M02]*det3_134_013 - rhs[M03]*det3_134_012;
0386   const Scalar det4_0134_0124 = rhs[M00]*det3_134_124 - rhs[M01]*det3_134_024
0387                                   + rhs[M02]*det3_134_014 - rhs[M04]*det3_134_012;
0388   const Scalar det4_0134_0134 = rhs[M00]*det3_134_134 - rhs[M01]*det3_134_034
0389                                   + rhs[M03]*det3_134_014 - rhs[M04]*det3_134_013;
0390   const Scalar det4_0134_0234 = rhs[M00]*det3_134_234 - rhs[M02]*det3_134_034
0391                                   + rhs[M03]*det3_134_024 - rhs[M04]*det3_134_023;
0392   const Scalar det4_0134_1234 = rhs[M01]*det3_134_234 - rhs[M02]*det3_134_134
0393                                   + rhs[M03]*det3_134_124 - rhs[M04]*det3_134_123;
0394   const Scalar det4_0234_0123 = rhs[M00]*det3_234_123 - rhs[M01]*det3_234_023
0395                                   + rhs[M02]*det3_234_013 - rhs[M03]*det3_234_012;
0396   const Scalar det4_0234_0124 = rhs[M00]*det3_234_124 - rhs[M01]*det3_234_024
0397                                   + rhs[M02]*det3_234_014 - rhs[M04]*det3_234_012;
0398   const Scalar det4_0234_0134 = rhs[M00]*det3_234_134 - rhs[M01]*det3_234_034
0399                                   + rhs[M03]*det3_234_014 - rhs[M04]*det3_234_013;
0400   const Scalar det4_0234_0234 = rhs[M00]*det3_234_234 - rhs[M02]*det3_234_034
0401                                   + rhs[M03]*det3_234_024 - rhs[M04]*det3_234_023;
0402   const Scalar det4_0234_1234 = rhs[M01]*det3_234_234 - rhs[M02]*det3_234_134
0403                                   + rhs[M03]*det3_234_124 - rhs[M04]*det3_234_123;
0404   const Scalar det4_1234_0123 = rhs[M10]*det3_234_123 - rhs[M11]*det3_234_023
0405                                   + rhs[M12]*det3_234_013 - rhs[M13]*det3_234_012;
0406   const Scalar det4_1234_0124 = rhs[M10]*det3_234_124 - rhs[M11]*det3_234_024
0407                                   + rhs[M12]*det3_234_014 - rhs[M14]*det3_234_012;
0408   const Scalar det4_1234_0134 = rhs[M10]*det3_234_134 - rhs[M11]*det3_234_034
0409                                   + rhs[M13]*det3_234_014 - rhs[M14]*det3_234_013;
0410   const Scalar det4_1234_0234 = rhs[M10]*det3_234_234 - rhs[M12]*det3_234_034
0411                                   + rhs[M13]*det3_234_024 - rhs[M14]*det3_234_023;
0412   const Scalar det4_1234_1234 = rhs[M11]*det3_234_234 - rhs[M12]*det3_234_134
0413                                   + rhs[M13]*det3_234_124 - rhs[M14]*det3_234_123;
0414
0415   // Find the 5x5 det:
0416
0417   const Scalar det = rhs[M00]*det4_1234_1234 - rhs[M01]*det4_1234_0234 + rhs[M02]*det4_1234_0134
0418                        - rhs[M03]*det4_1234_0124 + rhs[M04]*det4_1234_0123;
0419
0420 //   if (determ)
0421 //     *determ = det;
0422
0423   if ( det == 0 ) {
0424     //Error("Inv5x5","matrix is singular");
0425     //m.Invalidate();
0426     return false;
0427   }
0428
0429   const Scalar oneOverDet = 1.0f / det;
0430   const Scalar mn1OverDet = - oneOverDet;
0431
0432   rhs[M00] =  det4_1234_1234 * oneOverDet;
0433   rhs[M01] =  det4_0234_1234 * mn1OverDet;
0434   rhs[M02] =  det4_0134_1234 * oneOverDet;
0435   rhs[M03] =  det4_0124_1234 * mn1OverDet;
0436   rhs[M04] =  det4_0123_1234 * oneOverDet;
0437
0438   rhs[M10] =  det4_1234_0234 * mn1OverDet;
0439   rhs[M11] =  det4_0234_0234 * oneOverDet;
0440   rhs[M12] =  det4_0134_0234 * mn1OverDet;
0441   rhs[M13] =  det4_0124_0234 * oneOverDet;
0442   rhs[M14] =  det4_0123_0234 * mn1OverDet;
0443
0444   rhs[M20] =  det4_1234_0134 * oneOverDet;
0445   rhs[M21] =  det4_0234_0134 * mn1OverDet;
0446   rhs[M22] =  det4_0134_0134 * oneOverDet;
0447   rhs[M23] =  det4_0124_0134 * mn1OverDet;
0448   rhs[M24] =  det4_0123_0134 * oneOverDet;
0449
0450   rhs[M30] =  det4_1234_0124 * mn1OverDet;
0451   rhs[M31] =  det4_0234_0124 * oneOverDet;
0452   rhs[M32] =  det4_0134_0124 * mn1OverDet;
0453   rhs[M33] =  det4_0124_0124 * oneOverDet;
0454   rhs[M34] =  det4_0123_0124 * mn1OverDet;
0455
0456   rhs[M40] =  det4_1234_0123 * oneOverDet;
0457   rhs[M41] =  det4_0234_0123 * mn1OverDet;
0458   rhs[M42] =  det4_0134_0123 * oneOverDet;
0459   rhs[M43] =  det4_0124_0123 * mn1OverDet;
0460   rhs[M44] =  det4_0123_0123 * oneOverDet;
0461
0462   return true;
0463 }
0464
0465
0466   }  // namespace Math
0467
0468 }  // namespace ROOT
0469
0470
0471
0472
0473 // undef macros to avoid conflicts
0474
0475 // 4x4 indices
0476
0477 #undef F00
0478 #undef F01
0479 #undef F02
0480 #undef F03
0481
0482 #undef F10
0483 #undef F11
0484 #undef F12
0485 #undef F13
0486
0487 #undef F20
0488 #undef F21
0489 #undef F22
0490 #undef F23
0491
0492 #undef F30
0493 #undef F31
0494 #undef F32
0495 #undef F33
0496
0497 // undef 5x5 indices
0498
0499 #undef M00
0500 #undef M01
0501 #undef M02
0502 #undef M03
0503 #undef M04
0504
0505 #undef M10
0506 #undef M11
0507 #undef M12
0508 #undef M13
0509 #undef M14
0510
0511 #undef M20
0512 #undef M21
0513 #undef M22
0514 #undef M23
0515 #undef M24
0516
0517 #undef M30
0518 #undef M31
0519 #undef M32
0520 #undef M33
0521 #undef M34
0522
0523 #undef M40
0524 #undef M41
0525 #undef M42
0526 #undef M43
0527 #undef M44
0528
0529
0530 #endif