ROOT  6.06/09
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 
32 #ifndef ROOT_Math_MatrixRepresentationsStatic
34 #endif
35 
36 namespace ROOT {
37 
38  namespace Math {
39 
40 
41 
42 
43 /** Dsfact.
44  Compute determinant of a symmetric, positive definite matrix of dimension
45  $idim$ and order $n$.
46 
47  @author T. Glebe
48 */
49 
50 template <unsigned int n, unsigned int idim =n>
51 class SDeterminant {
52 
53 public:
54 template <class T>
55 static bool Dsfact(MatRepStd<T,n,idim>& rhs, T& det) {
56 
57 #ifdef XXX
58  /* Function Body */
59  if (idim < n || n <= 0) {
60  return false;
61  }
62 #endif
63 
64 #ifdef OLD_IMPL
65  typename MatrixRep::value_type* a = rhs.Array();
66 #endif
67 
68 #ifdef XXX
69  const typename MatrixRep::value_type* A = rhs.Array();
70  typename MatrixRep::value_type array[MatrixRep::kSize];
71  typename MatrixRep::value_type* a = array;
72 
73  // copy contents of matrix to working place
74  for(unsigned int i=0; i<MatrixRep::kSize; ++i) {
75  array[i] = A[i];
76  }
77 #endif
78 
79  /* Local variables */
80  unsigned int i, j, l;
81 
82  /* Parameter adjustments */
83  // a -= idim + 1;
84  const int arrayOffset = -(idim+1);
85  /* sfactd.inc */
86  det = 1.;
87  for (j = 1; j <= n; ++j) {
88  const unsigned int ji = j * idim;
89  const unsigned int jj = j + ji;
90 
91  if (rhs[jj + arrayOffset] <= 0.) {
92  det = 0.;
93  return false;
94  }
95 
96  const unsigned int jp1 = j + 1;
97  const unsigned int jpi = jp1 * idim;
98 
99  det *= rhs[jj + arrayOffset];
100  rhs[jj + arrayOffset] = 1. / rhs[jj + arrayOffset];
101 
102  for (l = jp1; l <= n; ++l) {
103  rhs[j + l * idim + arrayOffset] = rhs[jj + arrayOffset] * rhs[l + ji + arrayOffset];
104 
105  const unsigned int lj = l + jpi;
106 
107  for (i = 1; i <= j; ++i) {
108  rhs[lj + arrayOffset] -= rhs[l + i * idim + arrayOffset] * rhs[i + jpi + arrayOffset];
109  } // for i
110  } // for l
111  } // for j
112 
113  return true;
114 } // end of Dsfact
115 
116 
117  // t.b.d re-implement methods for symmetric
118  // symmetric function (copy in a general one)
119  template <class T>
120  static bool Dsfact(MatRepSym<T,n> & rhs, T & det) {
121  // not very efficient but need to re-do Dsinv for new storage of
122  // symmetric matrices
123  MatRepStd<T,n> tmp;
124  for (unsigned int i = 0; i< n*n; ++i)
125  tmp[i] = rhs[i];
126  if (! SDeterminant<n>::Dsfact(tmp,det) ) return false;
127 // // recopy the data
128 // for (int i = 0; i< idim*n; ++i)
129 // rhs[i] = tmp[i];
130 
131  return true;
132  }
133 
134 
135 }; // end of clas Sdeterminant
136 
137  } // namespace Math
138 
139 } // namespace ROOT
140 
141 #endif /* ROOT_Math_Dsfact */
142 
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
TArc * a
Definition: textangle.C:12
static double A[]
Expression wrapper class for Matrix objects.
Definition: Expression.h:134
TLine * l
Definition: textangle.C:4
static bool Dsfact(MatRepStd< T, n, idim > &rhs, T &det)
Definition: Dsfact.h:55
Namespace for new Math classes and functions.
const Int_t n
Definition: legend1.C:16
static bool Dsfact(MatRepSym< T, n > &rhs, T &det)
Definition: Dsfact.h:120