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