/usr/include/root/Math/Dsfact.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 | // @(#)root/smatrix:$Id$
// Authors: T. Glebe, L. Moneta 2005
#ifndef ROOT_Math_Dsfact
#define ROOT_Math_Dsfact
// ********************************************************************
//
// 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: Determinant of a symmetric, positive definite matrix.
// Code was taken from CERNLIB::kernlib dsfact function, translated
// from FORTRAN to C++ and optimized.
//
// changes:
// 22 Mar 2001 (TG) creation
// 18 Apr 2001 (TG) removed internal copying of array.
//
// ********************************************************************
namespace ROOT {
namespace Math {
/** Dsfact.
Compute determinant of a symmetric, positive definite matrix of dimension
$idim$ and order $n$.
@author T. Glebe
*/
template <unsigned int n, unsigned int idim =n>
class SDeterminant {
public:
template <class T>
static bool Dsfact(MatRepStd<T,n,idim>& rhs, T& det) {
#ifdef XXX
/* Function Body */
if (idim < n || n <= 0) {
return false;
}
#endif
#ifdef OLD_IMPL
typename MatrixRep::value_type* a = rhs.Array();
#endif
#ifdef XXX
const typename MatrixRep::value_type* A = rhs.Array();
typename MatrixRep::value_type array[MatrixRep::kSize];
typename MatrixRep::value_type* a = array;
// copy contents of matrix to working place
for(unsigned int i=0; i<MatrixRep::kSize; ++i) {
array[i] = A[i];
}
#endif
/* Local variables */
unsigned int i, j, l;
/* Parameter adjustments */
// a -= idim + 1;
const int arrayOffset = -(idim+1);
/* sfactd.inc */
det = 1.;
for (j = 1; j <= n; ++j) {
const unsigned int ji = j * idim;
const unsigned int jj = j + ji;
if (rhs[jj + arrayOffset] <= 0.) {
det = 0.;
return false;
}
const unsigned int jp1 = j + 1;
const unsigned int jpi = jp1 * idim;
det *= rhs[jj + arrayOffset];
rhs[jj + arrayOffset] = 1. / rhs[jj + arrayOffset];
for (l = jp1; l <= n; ++l) {
rhs[j + l * idim + arrayOffset] = rhs[jj + arrayOffset] * rhs[l + ji + arrayOffset];
const unsigned int lj = l + jpi;
for (i = 1; i <= j; ++i) {
rhs[lj + arrayOffset] -= rhs[l + i * idim + arrayOffset] * rhs[i + jpi + arrayOffset];
} // for i
} // for l
} // for j
return true;
} // end of Dsfact
// t.b.d re-implement methods for symmetric
// symmetric function (copy in a general one)
template <class T>
static bool Dsfact(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 (! SDeterminant<n>::Dsfact(tmp,det) ) return false;
// // recopy the data
// for (int i = 0; i< idim*n; ++i)
// rhs[i] = tmp[i];
return true;
}
}; // end of clas Sdeterminant
} // namespace Math
} // namespace ROOT
#endif /* ROOT_Math_Dsfact */
|