Logo ROOT   6.13/01
Reference Guide
Dfactir.h
Go to the documentation of this file.
1 // @(#)root/smatrix:$Id$
2 // Authors: T. Glebe, L. Moneta 2005
3 
4 #ifndef ROOT_Math_Dfactir
5 #define ROOT_Math_Dfactir
6 // ********************************************************************
7 //
8 // source:
9 //
10 // type: source code
11 //
12 // created: 02. Apr 2001
13 //
14 // author: Thorsten Glebe
15 // HERA-B Collaboration
16 // Max-Planck-Institut fuer Kernphysik
17 // Saupfercheckweg 1
18 // 69117 Heidelberg
19 // Germany
20 // E-mail: T.Glebe@mpi-hd.mpg.de
21 //
22 // Description: Determinant of a square matrix, needed for Dfinv()
23 // Code was taken from CERNLIB::kernlib dfact function, translated
24 // from FORTRAN to C++ and optimized.
25 //
26 // changes:
27 // 02 Apr 2001 (TG) creation
28 //
29 // ********************************************************************
30 
31 #include <cmath>
32 
33 namespace ROOT {
34 
35  namespace Math {
36 
37 
38 /** Dfactir.
39  Function to compute the determinant from a square matrix, Det(A) of
40  dimension idim and order n. A working area ir is returned which is
41  needed by the Dfinv routine.
42 
43  @author T. Glebe
44 */
45 template <class Matrix, unsigned int n, unsigned int idim>
46 bool Dfactir(Matrix& rhs, typename Matrix::value_type& det, unsigned int* ir)
47  // int *n, float *a, int *idim, int *ir, int *ifail, float *det, int *jfail)
48 {
49 
50 #ifdef XXX
51  if (idim < n || n <= 0) {
52  return false;
53  }
54 #endif
55 
56 
57  /* Initialized data */
58  typename Matrix::value_type* a = rhs.Array();
59 
60  /* Local variables */
61  unsigned int nxch, i, j, k, l;
62  typename Matrix::value_type p, q, tf;
63 
64  /* Parameter adjustments */
65  a -= idim + 1;
66  --ir;
67 
68  /* Function Body */
69 
70  // fact.inc
71  nxch = 0;
72  det = 1.;
73  for (j = 1; j <= n; ++j) {
74  const unsigned int ji = j * idim;
75  const unsigned int jj = j + ji;
76 
77  k = j;
78  p = std::abs(a[jj]);
79 
80  if (j != n) {
81  for (i = j + 1; i <= n; ++i) {
82  q = std::abs(a[i + ji]);
83  if (q > p) {
84  k = i;
85  p = q;
86  }
87  } // for i
88 
89  if (k != j) {
90  for (l = 1; l <= n; ++l) {
91  const unsigned int li = l*idim;
92  const unsigned int jli = j + li;
93  const unsigned int kli = k + li;
94  tf = a[jli];
95  a[jli] = a[kli];
96  a[kli] = tf;
97  } // for l
98  ++nxch;
99  ir[nxch] = (j << 12) + k;
100  } // if k != j
101  } // if j!=n
102 
103  if (p <= 0.) {
104  det = 0;
105  return false;
106  }
107 
108  det *= a[jj];
109 #ifdef XXX
110  t = std::abs(det);
111  if (t < 1e-19 || t > 1e19) {
112  det = 0;
113  return false;
114  }
115 #endif
116 
117  a[jj] = 1. / a[jj];
118  if (j == n) {
119  continue;
120  }
121 
122  const unsigned int jm1 = j - 1;
123  const unsigned int jpi = (j + 1) * idim;
124  const unsigned int jjpi = j + jpi;
125 
126  for (k = j + 1; k <= n; ++k) {
127  const unsigned int ki = k * idim;
128  const unsigned int jki = j + ki;
129  const unsigned int kji = k + jpi;
130  if (j != 1) {
131  for (i = 1; i <= jm1; ++i) {
132  const unsigned int ii = i * idim;
133  a[jki] -= a[i + ki] * a[j + ii];
134  a[kji] -= a[i + jpi] * a[k + ii];
135  } // for i
136  }
137  a[jki] *= a[jj];
138  a[kji] -= a[jjpi] * a[k + ji];
139  } // for k
140  } // for j
141 
142  if (nxch % 2 != 0) {
143  det = -(det);
144  }
145  ir[n] = nxch;
146  return true;
147 } // end of Dfact
148 
149 
150  } // namespace Math
151 
152 } // namespace ROOT
153 
154 
155 
156 #endif /* ROOT_Math_Dfactir */
Namespace for new ROOT classes and functions.
Definition: TFoamSampler.h:19
bool Dfactir(Matrix &rhs, typename Matrix::value_type &det, unsigned int *ir)
Dfactir.
Definition: Dfactir.h:46
Namespace for new Math classes and functions.