Logo ROOT  
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
35namespace 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*/
48template <unsigned int n, unsigned int idim = n>
50public:
51
52template <class T>
53static 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 */
winID h TVirtualViewer3D TVirtualGLPainter p
float * q
Definition: THbookFile.cxx:89
Detrminant for a general squared matrix Function to compute the determinant from a square matrix ( ) ...
Definition: Dfact.h:49
static bool Dfact(MatRepStd< T, n, idim > &rhs, T &det)
Definition: Dfact.h:53
static bool Dfact(MatRepSym< T, n > &rhs, T &det)
Definition: Dfact.h:157
Expression wrapper class for Matrix objects.
MatRepSym Matrix storage representation for a symmetric matrix of dimension NxN This class is a templ...
RVec< PromoteType< T > > abs(const RVec< T > &v)
Definition: RVec.hxx:1739
const Int_t n
Definition: legend1.C:16
Namespace for new Math classes and functions.
double T(double x)
Definition: ChebyshevPol.h:34
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.
auto * l
Definition: textangle.C:4