ROOT   Reference Guide
Dfactir.h
Go to the documentation of this file.
1// @(#)root/smatrix:$Id$
2// Authors: T. Glebe, L. Moneta 2005
3
4#ifndef ROOT_Math_Dfactir
5#define ROOT_Math_Dfactir
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, needed for Dfinv()
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
33namespace ROOT {
34
35 namespace Math {
36
37
38/** Dfactir.
39 Function to compute the determinant from a square matrix, Det(A) of
40 dimension idim and order n. A working area ir is returned which is
41 needed by the Dfinv routine.
42
43 @author T. Glebe
44*/
45template <class Matrix, unsigned int n, unsigned int idim>
46bool Dfactir(Matrix& rhs, typename Matrix::value_type& det, unsigned int* ir)
47 // int *n, float *a, int *idim, int *ir, int *ifail, float *det, int *jfail)
48{
49
50#ifdef XXX
51 if (idim < n || n <= 0) {
52 return false;
53 }
54#endif
55
56
57 /* Initialized data */
58 typename Matrix::value_type* a = rhs.Array();
59
60 /* Local variables */
61 unsigned int nxch, i, j, k, l;
62 typename Matrix::value_type p, q, tf;
63
65 a -= idim + 1;
66 --ir;
67
68 /* Function Body */
69
70 // fact.inc
71 nxch = 0;
72 det = 1.;
73 for (j = 1; j <= n; ++j) {
74 const unsigned int ji = j * idim;
75 const unsigned int jj = j + ji;
76
77 k = j;
78 p = std::abs(a[jj]);
79
80 if (j != n) {
81 for (i = j + 1; i <= n; ++i) {
82 q = std::abs(a[i + ji]);
83 if (q > p) {
84 k = i;
85 p = q;
86 }
87 } // for i
88
89 if (k != j) {
90 for (l = 1; l <= n; ++l) {
91 const unsigned int li = l*idim;
92 const unsigned int jli = j + li;
93 const unsigned int kli = k + li;
94 tf = a[jli];
95 a[jli] = a[kli];
96 a[kli] = tf;
97 } // for l
98 ++nxch;
99 ir[nxch] = (j << 12) + k;
100 } // if k != j
101 } // if j!=n
102
103 if (p <= 0.) {
104 det = 0;
105 return false;
106 }
107
108 det *= a[jj];
109#ifdef XXX
110 t = std::abs(det);
111 if (t < 1e-19 || t > 1e19) {
112 det = 0;
113 return false;
114 }
115#endif
116
117 a[jj] = 1. / a[jj];
118 if (j == n) {
119 continue;
120 }
121
122 const unsigned int jm1 = j - 1;
123 const unsigned int jpi = (j + 1) * idim;
124 const unsigned int jjpi = j + jpi;
125
126 for (k = j + 1; k <= n; ++k) {
127 const unsigned int ki = k * idim;
128 const unsigned int jki = j + ki;
129 const unsigned int kji = k + jpi;
130 if (j != 1) {
131 for (i = 1; i <= jm1; ++i) {
132 const unsigned int ii = i * idim;
133 a[jki] -= a[i + ki] * a[j + ii];
134 a[kji] -= a[i + jpi] * a[k + ii];
135 } // for i
136 }
137 a[jki] *= a[jj];
138 a[kji] -= a[jjpi] * a[k + ji];
139 } // for k
140 } // for j
141
142 if (nxch % 2 != 0) {
143 det = -(det);
144 }
145 ir[n] = nxch;
146 return true;
147} // end of Dfact
148
149
150 } // namespace Math
151
152} // namespace ROOT
153
154
155
156#endif /* ROOT_Math_Dfactir */
float * q
Definition: THbookFile.cxx:87
const Int_t n
Definition: legend1.C:16
Namespace for new Math classes and functions.
bool Dfactir(Matrix &rhs, typename Matrix::value_type &det, unsigned int *ir)
Dfactir.
Definition: Dfactir.h:46
VSD Structures.
Definition: StringConv.hxx:21
auto * l
Definition: textangle.C:4
auto * a
Definition: textangle.C:12