Logo ROOT   6.10/09
Reference Guide
Dfact.h
Go to the documentation of this file.
1 // @(#)root/smatrix:$Id$
2 // Authors: T. Glebe, L. Moneta 2005
3 
4 #ifndef ROOT_Math_Dfact
5 #define ROOT_Math_Dfact
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
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 
34 
35 namespace ROOT {
36 
37  namespace Math {
38 
39 
40 
41 /**
42  Detrminant for a general squared matrix
43  Function to compute the determinant from a square matrix (\f$ \det(A)\f$) of
44  dimension idim and order n.
45 
46  @author T. Glebe
47 */
48 template <unsigned int n, unsigned int idim = n>
49 class Determinant {
50 public:
51 
52 template <class T>
53 static bool Dfact(MatRepStd<T,n,idim>& rhs, T& det) {
54 
55 #ifdef XXX
56  if (idim < n || n <= 0) {
57  return false;
58  }
59 #endif
60 
61 
62  /* Initialized data */
63  // const typename MatrixRep::value_type* A = rhs.Array();
64  //typename MatrixRep::value_type* a = rhs.Array();
65 
66  /* Local variables */
67  unsigned int nxch, i, j, k, l;
68  //static typename MatrixRep::value_type p, q, tf;
69  T p, q, tf;
70 
71  /* Parameter adjustments */
72  // a -= idim + 1;
73  const int arrayOffset = - int(idim+1);
74  /* Function Body */
75 
76  // fact.inc
77 
78  nxch = 0;
79  det = 1.;
80  for (j = 1; j <= n; ++j) {
81  const unsigned int ji = j * idim;
82  const unsigned int jj = j + ji;
83 
84  k = j;
85  p = std::abs(rhs[jj + arrayOffset]);
86 
87  if (j != n) {
88  for (i = j + 1; i <= n; ++i) {
89  q = std::abs(rhs[i + ji + arrayOffset]);
90  if (q > p) {
91  k = i;
92  p = q;
93  }
94  } // for i
95  if (k != j) {
96  for (l = 1; l <= n; ++l) {
97  const unsigned int li = l*idim;
98  const unsigned int jli = j + li;
99  const unsigned int kli = k + li;
100  tf = rhs[jli + arrayOffset];
101  rhs[jli + arrayOffset] = rhs[kli + arrayOffset];
102  rhs[kli + arrayOffset] = tf;
103  } // for l
104  ++nxch;
105  } // if k != j
106  } // if j!=n
107 
108  if (p <= 0.) {
109  det = 0;
110  return false;
111  }
112 
113  det *= rhs[jj + arrayOffset];
114 #ifdef XXX
115  t = std::abs(det);
116  if (t < 1e-19 || t > 1e19) {
117  det = 0;
118  return false;
119  }
120 #endif
121  // using 1.0f removes a warning on Windows (1.0f is still the same as 1.0)
122  rhs[jj + arrayOffset] = 1.0f / rhs[jj + arrayOffset];
123  if (j == n) {
124  continue;
125  }
126 
127  const unsigned int jm1 = j - 1;
128  const unsigned int jpi = (j + 1) * idim;
129  const unsigned int jjpi = j + jpi;
130 
131  for (k = j + 1; k <= n; ++k) {
132  const unsigned int ki = k * idim;
133  const unsigned int jki = j + ki;
134  const unsigned int kji = k + jpi;
135  if (j != 1) {
136  for (i = 1; i <= jm1; ++i) {
137  const unsigned int ii = i * idim;
138  rhs[jki + arrayOffset] -= rhs[i + ki + arrayOffset] * rhs[j + ii + arrayOffset];
139  rhs[kji + arrayOffset] -= rhs[i + jpi + arrayOffset] * rhs[k + ii + arrayOffset];
140  } // for i
141  }
142  rhs[jki + arrayOffset] *= rhs[jj + arrayOffset];
143  rhs[kji + arrayOffset] -= rhs[jjpi + arrayOffset] * rhs[k + ji + arrayOffset];
144  } // for k
145  } // for j
146 
147  if (nxch % 2 != 0) {
148  det = -(det);
149  }
150  return true;
151 } // end of Dfact
152 
153 
154  // t.b.d re-implement methods for symmetric
155  // symmetric function (copy in a general one)
156  template <class T>
157  static bool Dfact(MatRepSym<T,n> & rhs, T & det) {
158  // not very efficient but need to re-do Dsinv for new storage of
159  // symmetric matrices
160  MatRepStd<T,n> tmp;
161  for (unsigned int i = 0; i< n*n; ++i)
162  tmp[i] = rhs[i];
163  if (! Determinant<n>::Dfact(tmp,det) ) return false;
164 // // recopy the data
165 // for (int i = 0; i< idim*n; ++i)
166 // rhs[i] = tmp[i];
167 
168  return true;
169  }
170 
171 };
172 
173 
174  } // namespace Math
175 
176 } // namespace ROOT
177 
178 
179 
180 #endif /* ROOT_Math_Dfact */
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
double T(double x)
Definition: ChebyshevPol.h:34
MatRepSym Matrix storage representation for a symmetric matrix of dimension NxN This class is a templ...
Definition: HelperOps.h:35
static bool Dfact(MatRepSym< T, n > &rhs, T &det)
Definition: Dfact.h:157
Expression wrapper class for Matrix objects.
Definition: Expression.h:134
TLine * l
Definition: textangle.C:4
Detrminant for a general squared matrix Function to compute the determinant from a square matrix ( ) ...
Definition: Dfact.h:49
Namespace for new Math classes and functions.
static bool Dfact(MatRepStd< T, n, idim > &rhs, T &det)
Definition: Dfact.h:53
float * q
Definition: THbookFile.cxx:87
const Int_t n
Definition: legend1.C:16