/usr/include/root/Math/Dfactir.h is in libroot-math-smatrix-dev 5.34.30-0ubuntu8.
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 | // @(#)root/smatrix:$Id$
// Authors: T. Glebe, L. Moneta 2005
#ifndef ROOT_Math_Dfactir
#define ROOT_Math_Dfactir
// ********************************************************************
//
// 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, needed for Dfinv()
// Code was taken from CERNLIB::kernlib dfact function, translated
// from FORTRAN to C++ and optimized.
//
// changes:
// 02 Apr 2001 (TG) creation
//
// ********************************************************************
#include <cmath>
namespace ROOT {
namespace Math {
/** Dfactir.
Function to compute the determinant from a square matrix, Det(A) of
dimension idim and order n. A working area ir is returned which is
needed by the Dfinv routine.
@author T. Glebe
*/
template <class Matrix, unsigned int n, unsigned int idim>
bool Dfactir(Matrix& rhs, typename Matrix::value_type& det, unsigned int* ir)
// int *n, float *a, int *idim, int *ir, int *ifail, float *det, int *jfail)
{
#ifdef XXX
if (idim < n || n <= 0) {
return false;
}
#endif
/* Initialized data */
typename Matrix::value_type* a = rhs.Array();
/* Local variables */
unsigned int nxch, i, j, k, l;
typename Matrix::value_type p, q, tf;
/* Parameter adjustments */
a -= idim + 1;
--ir;
/* 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(a[jj]);
if (j != n) {
for (i = j + 1; i <= n; ++i) {
q = std::abs(a[i + ji]);
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 = a[jli];
a[jli] = a[kli];
a[kli] = tf;
} // for l
++nxch;
ir[nxch] = (j << 12) + k;
} // if k != j
} // if j!=n
if (p <= 0.) {
det = 0;
return false;
}
det *= a[jj];
#ifdef XXX
t = std::abs(det);
if (t < 1e-19 || t > 1e19) {
det = 0;
return false;
}
#endif
a[jj] = 1. / a[jj];
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;
a[jki] -= a[i + ki] * a[j + ii];
a[kji] -= a[i + jpi] * a[k + ii];
} // for i
}
a[jki] *= a[jj];
a[kji] -= a[jjpi] * a[k + ji];
} // for k
} // for j
if (nxch % 2 != 0) {
det = -(det);
}
ir[n] = nxch;
return true;
} // end of Dfact
} // namespace Math
} // namespace ROOT
#endif /* ROOT_Math_Dfactir */
|