/usr/include/root/Math/Dfact.h is in libroot-math-smatrix-dev 5.34.19+dfsg-1.2.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 | // @(#)root/smatrix:$Id$
// Authors: T. Glebe, L. Moneta 2005
#ifndef ROOT_Math_Dfact
#define ROOT_Math_Dfact
// ********************************************************************
//
// source:
//
// type: source code
//
// created: 02. Apr 2001
//
// author: Thorsten Glebe
// HERA-B Collaboration
// Max-Planck-Institut fuer Kernphysik
// Saupfercheckweg 1
// 69117 Heidelberg
// Germany
// E-mail: T.Glebe@mpi-hd.mpg.de
//
// Description: Determinant of a square matrix
// Code was taken from CERNLIB::kernlib dfact function, translated
// from FORTRAN to C++ and optimized.
//
// changes:
// 02 Apr 2001 (TG) creation
//
// ********************************************************************
#include <cmath>
#ifndef ROOT_Math_MatrixRepresentationsStatic
#include "Math/MatrixRepresentationsStatic.h"
#endif
namespace ROOT {
namespace Math {
/**
Detrminant for a general squared matrix
Function to compute the determinant from a square matrix (\f$ \det(A)\f$) of
dimension idim and order n.
@author T. Glebe
*/
template <unsigned int n, unsigned int idim = n>
class Determinant {
public:
template <class T>
static bool Dfact(MatRepStd<T,n,idim>& rhs, T& det) {
#ifdef XXX
if (idim < n || n <= 0) {
return false;
}
#endif
/* Initialized data */
// const typename MatrixRep::value_type* A = rhs.Array();
//typename MatrixRep::value_type* a = rhs.Array();
/* Local variables */
unsigned int nxch, i, j, k, l;
//static typename MatrixRep::value_type p, q, tf;
T p, q, tf;
/* Parameter adjustments */
// a -= idim + 1;
const int arrayOffset = - int(idim+1);
/* Function Body */
// fact.inc
nxch = 0;
det = 1.;
for (j = 1; j <= n; ++j) {
const unsigned int ji = j * idim;
const unsigned int jj = j + ji;
k = j;
p = std::abs(rhs[jj + arrayOffset]);
if (j != n) {
for (i = j + 1; i <= n; ++i) {
q = std::abs(rhs[i + ji + arrayOffset]);
if (q > p) {
k = i;
p = q;
}
} // for i
if (k != j) {
for (l = 1; l <= n; ++l) {
const unsigned int li = l*idim;
const unsigned int jli = j + li;
const unsigned int kli = k + li;
tf = rhs[jli + arrayOffset];
rhs[jli + arrayOffset] = rhs[kli + arrayOffset];
rhs[kli + arrayOffset] = tf;
} // for l
++nxch;
} // if k != j
} // if j!=n
if (p <= 0.) {
det = 0;
return false;
}
det *= rhs[jj + arrayOffset];
#ifdef XXX
t = std::abs(det);
if (t < 1e-19 || t > 1e19) {
det = 0;
return false;
}
#endif
// using 1.0f removes a warning on Windows (1.0f is still the same as 1.0)
rhs[jj + arrayOffset] = 1.0f / rhs[jj + arrayOffset];
if (j == n) {
continue;
}
const unsigned int jm1 = j - 1;
const unsigned int jpi = (j + 1) * idim;
const unsigned int jjpi = j + jpi;
for (k = j + 1; k <= n; ++k) {
const unsigned int ki = k * idim;
const unsigned int jki = j + ki;
const unsigned int kji = k + jpi;
if (j != 1) {
for (i = 1; i <= jm1; ++i) {
const unsigned int ii = i * idim;
rhs[jki + arrayOffset] -= rhs[i + ki + arrayOffset] * rhs[j + ii + arrayOffset];
rhs[kji + arrayOffset] -= rhs[i + jpi + arrayOffset] * rhs[k + ii + arrayOffset];
} // for i
}
rhs[jki + arrayOffset] *= rhs[jj + arrayOffset];
rhs[kji + arrayOffset] -= rhs[jjpi + arrayOffset] * rhs[k + ji + arrayOffset];
} // for k
} // for j
if (nxch % 2 != 0) {
det = -(det);
}
return true;
} // end of Dfact
// t.b.d re-implement methods for symmetric
// symmetric function (copy in a general one)
template <class T>
static bool Dfact(MatRepSym<T,n> & rhs, T & det) {
// not very efficient but need to re-do Dsinv for new storage of
// symmetric matrices
MatRepStd<T,n> tmp;
for (unsigned int i = 0; i< n*n; ++i)
tmp[i] = rhs[i];
if (! Determinant<n>::Dfact(tmp,det) ) return false;
// // recopy the data
// for (int i = 0; i< idim*n; ++i)
// rhs[i] = tmp[i];
return true;
}
};
} // namespace Math
} // namespace ROOT
#endif /* ROOT_Math_Dfact */
|