root/Math/Dfinv.h

0001 // @(#)root/smatrix:$Id$
0002 // Authors: T. Glebe, L. Moneta    2005
0003
0004 #ifndef ROOT_Math_Dfinv
0005 #define ROOT_Math_Dfinv
0006 // ********************************************************************
0007 //
0008 // source:
0009 //
0010 // type:      source code
0011 //
0012 // created:   03. 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: Matrix inversion
0023 //              Code was taken from CERNLIB::kernlib dfinv function, translated
0024 //              from FORTRAN to C++ and optimized.
0025 //
0026 // changes:
0027 // 03 Apr 2001 (TG) creation
0028 //
0029 // ********************************************************************
0030
0031
0032 namespace ROOT {
0033
0034   namespace Math {
0035
0036
0037
0038
0039 /** Dfinv.
0040     Function to compute the inverse of a square matrix (\f$A^{-1}\f$) of
0041     dimension \f$idim\f$ and order \f$n\f$. The routine Dfactir must be called
0042     before Dfinv!
0043
0044     @author T. Glebe
0045 */
0046 template <class Matrix, unsigned int n, unsigned int idim>
0047 bool Dfinv(Matrix& rhs, unsigned int* ir) {
0048 #ifdef XXX
0049   if (idim < n || n <= 0 || n==1) {
0050     return false;
0051   }
0052 #endif
0053
0054   typename Matrix::value_type* a = rhs.Array();
0055
0056   /* Local variables */
0057   unsigned int nxch, i, j, k, m, ij;
0058   unsigned int im2, nm1, nmi;
0059   typename Matrix::value_type s31, s34, ti;
0060
0061   /* Parameter adjustments */
0062   a -= idim + 1;
0063   --ir;
0064
0065   /* Function Body */
0066
0067   /* finv.inc */
0068
0069   a[idim + 2] = -a[(idim << 1) + 2] * (a[idim + 1] * a[idim + 2]);
0070   a[(idim << 1) + 1] = -a[(idim << 1) + 1];
0071
0072    if (n != 2) {
0073       for (i = 3; i <= n; ++i) {
0074          const unsigned int ii   = i * idim;
0075          const unsigned int iii  = i + ii;
0076          const unsigned int imi  = ii - idim;
0077          const unsigned int iimi = i + imi;
0078          im2 = i - 2;
0079          for (j = 1; j <= im2; ++j) {
0080             const unsigned int ji  = j * idim;
0081             const unsigned int jii = j + ii;
0082             s31 = 0.;
0083             for (k = j; k <= im2; ++k) {
0084                s31 += a[k + ji] * a[i + k * idim];
0085                a[jii] += a[j + (k + 1) * idim] * a[k + 1 + ii];
0086             } // for k
0087             a[i + ji] = -a[iii] * (a[i - 1 + ji] * a[iimi] + s31);
0088             a[jii] *= -1;
0089          } // for j
0090          a[iimi] = -a[iii] * (a[i - 1 + imi] * a[iimi]);
0091          a[i - 1 + ii] *= -1;
0092       } // for i
0093    } // if n!=2
0094
0095    nm1 = n - 1;
0096    for (i = 1; i <= nm1; ++i) {
0097       const unsigned int ii = i * idim;
0098       nmi = n - i;
0099       for (j = 1; j <= i; ++j) {
0100          const unsigned int ji  = j * idim;
0101          const unsigned int iji = i + ji;
0102          for (k = 1; k <= nmi; ++k) {
0103             a[iji] += a[i + k + ji] * a[i + (i + k) * idim];
0104          } // for k
0105       } // for j
0106
0107       for (j = 1; j <= nmi; ++j) {
0108          const unsigned int ji = j * idim;
0109          s34 = 0.;
0110          for (k = j; k <= nmi; ++k) {
0111             s34 += a[i + k + ii + ji] * a[i + (i + k) * idim];
0112          } // for k
0113          a[i + ii + ji] = s34;
0114       } // for j
0115    } // for i
0116
0117    nxch = ir[n];
0118    if (nxch == 0) {
0119       return true;
0120    }
0121
0122    for (m = 1; m <= nxch; ++m) {
0123       k = nxch - m + 1;
0124       ij = ir[k];
0125       i = ij / 4096;
0126       j = ij % 4096;
0127       const unsigned int ii = i * idim;
0128       const unsigned int ji = j * idim;
0129       for (k = 1; k <= n; ++k) {
0130          ti = a[k + ii];
0131          a[k + ii] = a[k + ji];
0132          a[k + ji] = ti;
0133       } // for k
0134    } // for m
0135
0136    return true;
0137 } // Dfinv
0138
0139
0140   }  // namespace Math
0141
0142 }  // namespace ROOT
0143
0144
0145
0146 #endif  /* ROOT_Math_Dfinv */