/usr/include/root/Math/Dsinv.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 | // @(#)root/smatrix:$Id$
// Authors: T. Glebe, L. Moneta 2005
#ifndef ROOT_Math_Dsinv
#define ROOT_Math_Dsinv
// ********************************************************************
//
// source:
//
// type: source code
//
// created: 22. Mar 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: Inversion of a symmetric, positive definite matrix.
// Code was taken from CERNLIB::kernlib dsinv function, translated
// from FORTRAN to C++ and optimized.
//
// changes:
// 22 Mar 2001 (TG) creation
//
// ********************************************************************
namespace ROOT {
namespace Math {
/** Dsinv.
Compute inverse of a symmetric, positive definite matrix of dimension
$idim$ and order $n$.
@author T. Glebe
*/
template <class T, int n, int idim>
class SInverter
{
public:
template <class MatrixRep>
inline static bool Dsinv(MatrixRep& rhs) {
/* Local variables */
int i, j, k, l;
T s31, s32;
int jm1, jp1;
/* Parameter adjustments */
const int arrayOffset = -1*(idim + 1);
/* Function Body */
if (idim < n || n <= 1) {
return false;
}
/* sfact.inc */
for (j = 1; j <= n; ++j) {
const int ja = j * idim;
const int jj = j + ja;
const int ja1 = ja + idim;
if (rhs[jj + arrayOffset] <= 0.) { return false; }
rhs[jj + arrayOffset] = 1. / rhs[jj + arrayOffset];
if (j == n) { break; }
for (l = j + 1; l <= n; ++l) {
rhs[j + (l * idim) + arrayOffset] = rhs[jj + arrayOffset] * rhs[l + ja + arrayOffset];
const int lj = l + ja1;
for (i = 1; i <= j; ++i) {
rhs[lj + arrayOffset] -= rhs[l + (i * idim) + arrayOffset] * rhs[i + ja1 + arrayOffset];
}
}
}
/* sfinv.inc */
// idim << 1 is equal to idim * 2
// compiler will compute the arguments!
rhs[((idim << 1) + 1) + arrayOffset] = -rhs[((idim << 1) + 1) + arrayOffset];
rhs[idim + 2 + arrayOffset] = rhs[((idim << 1)) + 1 + arrayOffset] * rhs[((idim << 1)) + 2 + arrayOffset];
if(n > 2) {
for (j = 3; j <= n; ++j) {
const int jm2 = j - 2;
const int ja = j * idim;
const int jj = j + ja;
const int j1 = j - 1 + ja;
for (k = 1; k <= jm2; ++k) {
s31 = rhs[k + ja + arrayOffset];
for (i = k; i <= jm2; ++i) {
s31 += rhs[k + ((i + 1) * idim) + arrayOffset] * rhs[i + 1 + ja + arrayOffset];
} // for i
rhs[k + ja + arrayOffset] = -s31;
rhs[j + (k * idim) + arrayOffset] = -s31 * rhs[jj + arrayOffset];
} // for k
rhs[j1 + arrayOffset] *= -1;
// rhs[j1] = -rhs[j1];
rhs[jj - idim + arrayOffset] = rhs[j1 + arrayOffset] * rhs[jj + arrayOffset];
} // for j
} // if (n>2)
j = 1;
do {
const int jad = j * idim;
const int jj = j + jad;
jp1 = j + 1;
for (i = jp1; i <= n; ++i) {
rhs[jj + arrayOffset] += rhs[j + (i * idim) + arrayOffset] * rhs[i + jad + arrayOffset];
} // for i
jm1 = j;
j = jp1;
const int ja = j * idim;
for (k = 1; k <= jm1; ++k) {
s32 = 0.;
for (i = j; i <= n; ++i) {
s32 += rhs[k + (i * idim) + arrayOffset] * rhs[i + ja + arrayOffset];
} // for i
//rhs[k + ja + arrayOffset] = rhs[j + (k * idim) + arrayOffset] = s32;
rhs[k + ja + arrayOffset] = s32;
} // for k
} while(j < n);
return true;
}
// for symmetric matrices
static bool Dsinv(MatRepSym<T,n> & rhs) {
// not very efficient but need to re-do Dsinv for new storage of
// symmetric matrices
MatRepStd<T,n,n> tmp;
for (int i = 0; i< n*n; ++i)
tmp[i] = rhs[i];
// call dsinv
if (! SInverter<T,n,n>::Dsinv(tmp) ) return false;
//if (! Inverter<n>::Dinv(tmp) ) return false;
// recopy the data
for (int i = 0; i< n*n; ++i)
rhs[i] = tmp[i];
return true;
}
}; // end of Dsinv
} // namespace Math
} // namespace ROOT
#endif /* ROOT_Math_Dsinv */
|