Logo ROOT   6.12/07
Reference Guide
Dsfact.h
Go to the documentation of this file.
1 // @(#)root/smatrix:$Id$
2 // Authors: T. Glebe, L. Moneta 2005
3 
4 #ifndef ROOT_Math_Dsfact
5 #define ROOT_Math_Dsfact
6 // ********************************************************************
7 //
8 // source:
9 //
10 // type: source code
11 //
12 // created: 22. Mar 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 symmetric, positive definite matrix.
23 // Code was taken from CERNLIB::kernlib dsfact function, translated
24 // from FORTRAN to C++ and optimized.
25 //
26 // changes:
27 // 22 Mar 2001 (TG) creation
28 // 18 Apr 2001 (TG) removed internal copying of array.
29 //
30 // ********************************************************************
31 
33 
34 namespace ROOT {
35 
36  namespace Math {
37 
38 
39 
40 
41 /** Dsfact.
42  Compute determinant of a symmetric, positive definite matrix of dimension
43  $idim$ and order $n$.
44 
45  @author T. Glebe
46 */
47 
48 template <unsigned int n, unsigned int idim =n>
49 class SDeterminant {
50 
51 public:
52 template <class T>
53 static bool Dsfact(MatRepStd<T,n,idim>& rhs, T& det) {
54 
55 #ifdef XXX
56  /* Function Body */
57  if (idim < n || n <= 0) {
58  return false;
59  }
60 #endif
61 
62 #ifdef OLD_IMPL
63  typename MatrixRep::value_type* a = rhs.Array();
64 #endif
65 
66 #ifdef XXX
67  const typename MatrixRep::value_type* A = rhs.Array();
68  typename MatrixRep::value_type array[MatrixRep::kSize];
69  typename MatrixRep::value_type* a = array;
70 
71  // copy contents of matrix to working place
72  for(unsigned int i=0; i<MatrixRep::kSize; ++i) {
73  array[i] = A[i];
74  }
75 #endif
76 
77  /* Local variables */
78  unsigned int i, j, l;
79 
80  /* Parameter adjustments */
81  // a -= idim + 1;
82  const int arrayOffset = -(idim+1);
83  /* sfactd.inc */
84  det = 1.;
85  for (j = 1; j <= n; ++j) {
86  const unsigned int ji = j * idim;
87  const unsigned int jj = j + ji;
88 
89  if (rhs[jj + arrayOffset] <= 0.) {
90  det = 0.;
91  return false;
92  }
93 
94  const unsigned int jp1 = j + 1;
95  const unsigned int jpi = jp1 * idim;
96 
97  det *= rhs[jj + arrayOffset];
98  rhs[jj + arrayOffset] = 1. / rhs[jj + arrayOffset];
99 
100  for (l = jp1; l <= n; ++l) {
101  rhs[j + l * idim + arrayOffset] = rhs[jj + arrayOffset] * rhs[l + ji + arrayOffset];
102 
103  const unsigned int lj = l + jpi;
104 
105  for (i = 1; i <= j; ++i) {
106  rhs[lj + arrayOffset] -= rhs[l + i * idim + arrayOffset] * rhs[i + jpi + arrayOffset];
107  } // for i
108  } // for l
109  } // for j
110 
111  return true;
112 } // end of Dsfact
113 
114 
115  // t.b.d re-implement methods for symmetric
116  // symmetric function (copy in a general one)
117  template <class T>
118  static bool Dsfact(MatRepSym<T,n> & rhs, T & det) {
119  // not very efficient but need to re-do Dsinv for new storage of
120  // symmetric matrices
121  MatRepStd<T,n> tmp;
122  for (unsigned int i = 0; i< n*n; ++i)
123  tmp[i] = rhs[i];
124  if (! SDeterminant<n>::Dsfact(tmp,det) ) return false;
125 // // recopy the data
126 // for (int i = 0; i< idim*n; ++i)
127 // rhs[i] = tmp[i];
128 
129  return true;
130  }
131 
132 
133 }; // end of clas Sdeterminant
134 
135  } // namespace Math
136 
137 } // namespace ROOT
138 
139 #endif /* ROOT_Math_Dsfact */
140 
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 double A[]
Expression wrapper class for Matrix objects.
Definition: Expression.h:134
auto * a
Definition: textangle.C:12
static bool Dsfact(MatRepStd< T, n, idim > &rhs, T &det)
Definition: Dsfact.h:53
Namespace for new Math classes and functions.
auto * l
Definition: textangle.C:4
const Int_t n
Definition: legend1.C:16
static bool Dsfact(MatRepSym< T, n > &rhs, T &det)
Definition: Dsfact.h:118