root/Math/Dfact.h

0001 // @(#)root/smatrix:$Id$
0002 // Authors: T. Glebe, L. Moneta    2005
0003
0004 #ifndef ROOT_Math_Dfact
0005 #define ROOT_Math_Dfact
0006 // ********************************************************************
0007 //
0008 // source:
0009 //
0010 // type:      source code
0011 //
0012 // created:   02. Apr 2001
0013 //
0014 // author:    Thorsten Glebe
0015 //            HERA-B Collaboration
0016 //            Max-Planck-Institut fuer Kernphysik
0017 //            Saupfercheckweg 1
0018 //            69117 Heidelberg
0019 //            Germany
0020 //            E-mail: T.Glebe@mpi-hd.mpg.de
0021 //
0022 // Description: Determinant of a square matrix
0023 //              Code was taken from CERNLIB::kernlib dfact function, translated
0024 //              from FORTRAN to C++ and optimized.
0025 //
0026 // changes:
0027 // 02 Apr 2001 (TG) creation
0028 //
0029 // ********************************************************************
0030
0031 #include <cmath>
0032
0033 #include "Math/MatrixRepresentationsStatic.h"
0034
0035 namespace ROOT {
0036
0037   namespace Math {
0038
0039
0040
0041 /**
0042     Detrminant for a general squared matrix
0043     Function to compute the determinant from a square matrix (\f$ \det(A)\f$) of
0044     dimension idim and order n.
0045
0046     @author T. Glebe
0047 */
0048 template <unsigned int n, unsigned int idim = n>
0049 class Determinant {
0050 public:
0051
0052 template <class T>
0053 static bool Dfact(MatRepStd<T,n,idim>& rhs, T& det) {
0054
0055 #ifdef XXX
0056   if (idim < n || n <= 0) {
0057     return false;
0058   }
0059 #endif
0060
0061
0062   /* Initialized data */
0063   //  const typename MatrixRep::value_type* A = rhs.Array();
0064   //typename MatrixRep::value_type* a = rhs.Array();
0065
0066   /* Local variables */
0067   unsigned int nxch, i, j, k, l;
0068   //static typename MatrixRep::value_type p, q, tf;
0069   T p, q, tf;
0070
0071   /* Parameter adjustments */
0072   //  a -= idim + 1;
0073   const int arrayOffset = - int(idim+1);
0074   /* Function Body */
0075
0076   // fact.inc
0077
0078    nxch = 0;
0079    det = 1.;
0080    for (j = 1; j <= n; ++j) {
0081       const unsigned int ji = j * idim;
0082       const unsigned int jj = j + ji;
0083
0084       k = j;
0085       p = std::abs(rhs[jj + arrayOffset]);
0086
0087       if (j != n) {
0088          for (i = j + 1; i <= n; ++i) {
0089             q = std::abs(rhs[i + ji + arrayOffset]);
0090             if (q > p) {
0091                k = i;
0092                p = q;
0093             }
0094          } // for i
0095          if (k != j) {
0096             for (l = 1; l <= n; ++l) {
0097                const unsigned int li = l*idim;
0098                const unsigned int jli = j + li;
0099                const unsigned int kli = k + li;
0100                tf = rhs[jli + arrayOffset];
0101                rhs[jli + arrayOffset] = rhs[kli + arrayOffset];
0102                rhs[kli + arrayOffset] = tf;
0103             } // for l
0104             ++nxch;
0105          } // if k != j
0106       } // if j!=n
0107
0108       if (p <= 0.) {
0109          det = 0;
0110          return false;
0111       }
0112
0113       det *= rhs[jj + arrayOffset];
0114 #ifdef XXX
0115       t = std::abs(det);
0116       if (t < 1e-19 || t > 1e19) {
0117          det = 0;
0118          return false;
0119       }
0120 #endif
0121       // using 1.0f removes a warning on Windows (1.0f is still the same  as 1.0)
0122       rhs[jj + arrayOffset] = 1.0f / rhs[jj + arrayOffset];
0123       if (j == n) {
0124          continue;
0125       }
0126
0127       const unsigned int jm1 = j - 1;
0128       const unsigned int jpi = (j + 1) * idim;
0129       const unsigned int jjpi = j + jpi;
0130
0131       for (k = j + 1; k <= n; ++k) {
0132          const unsigned int ki  = k * idim;
0133          const unsigned int jki = j + ki;
0134          const unsigned int kji = k + jpi;
0135          if (j != 1) {
0136             for (i = 1; i <= jm1; ++i) {
0137                const unsigned int ii = i * idim;
0138                rhs[jki + arrayOffset] -= rhs[i + ki + arrayOffset] * rhs[j + ii + arrayOffset];
0139                rhs[kji + arrayOffset] -= rhs[i + jpi + arrayOffset] * rhs[k + ii + arrayOffset];
0140             } // for i
0141          }
0142          rhs[jki + arrayOffset] *= rhs[jj + arrayOffset];
0143          rhs[kji + arrayOffset] -= rhs[jjpi + arrayOffset] * rhs[k + ji + arrayOffset];
0144       } // for k
0145    } // for j
0146
0147    if (nxch % 2 != 0) {
0148       det = -(det);
0149   }
0150   return true;
0151 } // end of Dfact
0152
0153
0154    // t.b.d re-implement methods for symmetric
0155   // symmetric function (copy in a general  one)
0156   template <class T>
0157   static bool Dfact(MatRepSym<T,n> & rhs, T & det) {
0158     // not very efficient but need to re-do Dsinv for new storage of
0159     // symmetric matrices
0160     MatRepStd<T,n> tmp;
0161     for (unsigned int i = 0; i< n*n; ++i)
0162       tmp[i] = rhs[i];
0163     if (! Determinant<n>::Dfact(tmp,det) ) return false;
0164 //     // recopy the data
0165 //     for (int i = 0; i< idim*n; ++i)
0166 //       rhs[i] = tmp[i];
0167
0168     return true;
0169   }
0170
0171 };
0172
0173
0174   }  // namespace Math
0175
0176 }  // namespace ROOT
0177
0178
0179
0180 #endif /* ROOT_Math_Dfact */