root/Math/Dfactir.h

0001 // @(#)root/smatrix:$Id$
0002 // Authors: T. Glebe, L. Moneta    2005
0003
0004 #ifndef ROOT_Math_Dfactir
0005 #define ROOT_Math_Dfactir
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, needed for Dfinv()
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 namespace ROOT {
0034
0035   namespace Math {
0036
0037
0038 /** Dfactir.
0039     Function to compute the determinant from a square matrix, Det(A) of
0040     dimension idim and order n. A working area ir is returned which is
0041     needed by the Dfinv routine.
0042
0043     @author T. Glebe
0044 */
0045 template <class Matrix, unsigned int n, unsigned int idim>
0046 bool Dfactir(Matrix& rhs, typename Matrix::value_type& det, unsigned int* ir)
0047   // int *n, float *a, int *idim, int *ir, int *ifail, float *det, int *jfail)
0048 {
0049
0050 #ifdef XXX
0051   if (idim < n || n <= 0) {
0052     return false;
0053   }
0054 #endif
0055
0056
0057   /* Initialized data */
0058    typename Matrix::value_type* a = rhs.Array();
0059
0060    /* Local variables */
0061    unsigned int nxch, i, j, k, l;
0062    typename Matrix::value_type p, q, tf;
0063
0064    /* Parameter adjustments */
0065    a -= idim + 1;
0066    --ir;
0067
0068    /* Function Body */
0069
0070    // fact.inc
0071    nxch = 0;
0072    det = 1.;
0073    for (j = 1; j <= n; ++j) {
0074       const unsigned int ji = j * idim;
0075       const unsigned int jj = j + ji;
0076
0077       k = j;
0078       p = std::abs(a[jj]);
0079
0080       if (j != n) {
0081          for (i = j + 1; i <= n; ++i) {
0082             q = std::abs(a[i + ji]);
0083             if (q > p) {
0084                k = i;
0085                p = q;
0086             }
0087          } // for i
0088
0089          if (k != j) {
0090             for (l = 1; l <= n; ++l) {
0091                const unsigned int li = l*idim;
0092                const unsigned int jli = j + li;
0093                const unsigned int kli = k + li;
0094                tf = a[jli];
0095                a[jli] = a[kli];
0096                a[kli] = tf;
0097             } // for l
0098             ++nxch;
0099             ir[nxch] = (j << 12) + k;
0100          } // if k != j
0101       } // if j!=n
0102
0103       if (p <= 0.) {
0104          det = 0;
0105          return false;
0106       }
0107
0108       det *= a[jj];
0109 #ifdef XXX
0110       t = std::abs(det);
0111       if (t < 1e-19 || t > 1e19) {
0112          det = 0;
0113          return false;
0114       }
0115 #endif
0116
0117       a[jj] = 1. / a[jj];
0118       if (j == n) {
0119          continue;
0120       }
0121
0122       const unsigned int jm1 = j - 1;
0123       const unsigned int jpi = (j + 1) * idim;
0124       const unsigned int jjpi = j + jpi;
0125
0126       for (k = j + 1; k <= n; ++k) {
0127          const unsigned int ki  = k * idim;
0128          const unsigned int jki = j + ki;
0129          const unsigned int kji = k + jpi;
0130          if (j != 1) {
0131             for (i = 1; i <= jm1; ++i) {
0132                const unsigned int ii = i * idim;
0133                a[jki] -= a[i + ki] * a[j + ii];
0134                a[kji] -= a[i + jpi] * a[k + ii];
0135             } // for i
0136          }
0137          a[jki] *= a[jj];
0138          a[kji] -= a[jjpi] * a[k + ji];
0139       } // for k
0140    } // for j
0141
0142    if (nxch % 2 != 0) {
0143       det = -(det);
0144    }
0145    ir[n] = nxch;
0146    return true;
0147 } // end of Dfact
0148
0149
0150   }  // namespace Math
0151
0152 }  // namespace ROOT
0153
0154
0155
0156 #endif /* ROOT_Math_Dfactir */