ROOT   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
34namespace ROOT {
35
36 namespace Math {
37
38
39
40
41/** Dsfact.
42 Compute determinant of a symmetric, positive definite matrix of dimension
43 \f$idim\f$ and order \f$n\f$.
44
45 @author T. Glebe
46*/
47
48template <unsigned int n, unsigned int idim =n>
50
51public:
52template <class T>
53static 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 class Sdeterminant
134
135 } // namespace Math
136
137} // namespace ROOT
138
139#endif /* ROOT_Math_Dsfact */
140
@ kSize
Definition: TStructNode.h:26
Expression wrapper class for Matrix objects.
MatRepSym Matrix storage representation for a symmetric matrix of dimension NxN This class is a templ...
static bool Dsfact(MatRepStd< T, n, idim > &rhs, T &det)
Definition: Dsfact.h:53
static bool Dsfact(MatRepSym< T, n > &rhs, T &det)
Definition: Dsfact.h:118
const Int_t n
Definition: legend1.C:16
Namespace for new Math classes and functions.
static double A[]
double T(double x)
Definition: ChebyshevPol.h:34
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.
TLine l
Definition: textangle.C:4
TArc a
Definition: textangle.C:12