root/Math/MatrixInversion.icc

0001 // @(#)root/smatrix:$Id$
0002 // Authors: CLHEP authors, L. Moneta    2006
0003
0004 #ifndef ROOT_Math_MatrixInversion_icc
0005 #define ROOT_Math_MatrixInversion_icc
0006
0007 #ifndef ROOT_Math_Dinv
0008 #error "Do not use MatrixInversion.icc directly. #include \"Math/Dinv.h\" instead."
0009 #endif // ROOT_Math_Dinv
0010
0011
0012 #include "Math/SVector.h"
0013 #include "Math/Math.h"
0014 #include <limits>
0015
0016
0017 // inversion algorithms for matrices
0018 // taken  from CLHEP (L. Moneta May 2006)
0019
0020 namespace ROOT {
0021
0022   namespace Math {
0023
0024
0025   /** General Inversion for a symmetric matrix
0026       Bunch-Kaufman diagonal pivoting method
0027       It is described in J.R. Bunch, L. Kaufman (1977).
0028       "Some Stable Methods for Calculating Inertia and Solving Symmetric
0029       Linear Systems", Math. Comp. 31, p. 162-179. or in Gene H. Golub,
0030       /Charles F. van Loan, "Matrix Computations" (the second edition
0031       has a bug.) and implemented in "lapack"
0032       Mario Stanke, 09/97
0033
0034   */
0035
0036
0037
0038 template <unsigned int idim, unsigned int N>
0039 template<class T>
0040 void Inverter<idim,N>::InvertBunchKaufman(MatRepSym<T,idim> & rhs, int &ifail) {
0041
0042    typedef T value_type;
0043    //typedef double value_type;  // for making inversions in double's
0044
0045
0046    int i, j, k, s;
0047    int pivrow;
0048
0049    const int nrow = MatRepSym<T,idim>::kRows;
0050
0051    // Establish the two working-space arrays needed:  x and piv are
0052    // used as pointers to arrays of doubles and ints respectively, each
0053    // of length nrow.  We do not want to reallocate each time through
0054    // unless the size needs to grow.  We do not want to leak memory, even
0055    // by having a new without a delete that is only done once.
0056
0057
0058    SVector<T, MatRepSym<T,idim>::kRows> xvec;
0059    SVector<int, MatRepSym<T,idim>::kRows> pivv;
0060
0061    typedef int* pivIter;
0062    typedef T* mIter;
0063
0064
0065    // Note - resize should do nothing if the size is already larger than nrow,
0066    //        but on VC++ there are indications that it does so we check.
0067    // Note - the data elements in a vector are guaranteed to be contiguous,
0068    //        so x[i] and piv[i] are optimally fast.
0069    mIter   x   = xvec.begin();
0070    // x[i] is used as helper storage, needs to have at least size nrow.
0071    pivIter piv = pivv.begin();
0072    // piv[i] is used to store details of exchanges
0073
0074    value_type temp1, temp2;
0075    mIter ip, mjj, iq;
0076    value_type lambda, sigma;
0077    const value_type alpha = .6404; // = (1+sqrt(17))/8
0078    // LM (04/2009) remove this useless check (it is not in LAPACK) which fails inversion of
0079    // a matrix with  values < epsilon in the diagonal
0080    //
0081    //const double epsilon = 32*std::numeric_limits<T>::epsilon();
0082    // whenever a sum of two doubles is below or equal to epsilon
0083    // it is set to zero.
0084    // this constant could be set to zero but then the algorithm
0085    // doesn't necessarily detect that a matrix is singular
0086
0087    for (i = 0; i < nrow; i++)
0088       piv[i] = i+1;
0089
0090    ifail = 0;
0091
0092    // compute the factorization P*A*P^T = L * D * L^T
0093    // L is unit lower triangular, D is direct sum of 1x1 and 2x2 matrices
0094    // L and D^-1 are stored in A = *this, P is stored in piv[]
0095
0096    for (j=1; j < nrow; j+=s)  // main loop over columns
0097    {
0098       mjj = rhs.Array() + j*(j-1)/2 + j-1;
0099       lambda = 0;           // compute lambda = max of A(j+1:n,j)
0100       pivrow = j+1;
0101       ip = rhs.Array() + (j+1)*j/2 + j-1;
0102       for (i=j+1; i <= nrow ; ip += i++)
0103          if (std::abs(*ip) > lambda)
0104          {
0105             lambda = std::abs(*ip);
0106             pivrow = i;
0107          }
0108
0109       if (lambda == 0 )
0110       {
0111          if (*mjj == 0)
0112          {
0113             ifail = 1;
0114             return;
0115          }
0116          s=1;
0117          *mjj = 1.0f / *mjj;
0118       }
0119       else
0120       {
0121          if (std::abs(*mjj) >= lambda*alpha)
0122          {
0123             s=1;
0124             pivrow=j;
0125          }
0126          else
0127          {
0128             sigma = 0;  // compute sigma = max A(pivrow, j:pivrow-1)
0129             ip = rhs.Array() + pivrow*(pivrow-1)/2+j-1;
0130             for (k=j; k < pivrow; k++)
0131             {
0132                if (std::abs(*ip) > sigma)
0133                   sigma = std::abs(*ip);
0134                ip++;
0135             }
0136             // sigma cannot be zero because it is at least lambda which is not zero
0137             if ( std::abs(*mjj) >= alpha * lambda * (lambda/ sigma) )
0138             {
0139                s=1;
0140                pivrow = j;
0141             }
0142             else if (std::abs(*(rhs.Array()+pivrow*(pivrow-1)/2+pivrow-1))
0143                      >= alpha * sigma)
0144                s=1;
0145             else
0146                s=2;
0147          }
0148          if (pivrow == j)  // no permutation necessary
0149          {
0150             piv[j-1] = pivrow;
0151             if (*mjj == 0)
0152             {
0153                ifail=1;
0154                return;
0155             }
0156             temp2 = *mjj = 1.0f/ *mjj; // invert D(j,j)
0157
0158             // update A(j+1:n, j+1,n)
0159             for (i=j+1; i <= nrow; i++)
0160             {
0161                temp1 = *(rhs.Array() + i*(i-1)/2 + j-1) * temp2;
0162                ip = rhs.Array()+i*(i-1)/2+j;
0163                for (k=j+1; k<=i; k++)
0164                {
0165                   *ip -= static_cast<T> ( temp1 * *(rhs.Array() + k*(k-1)/2 + j-1) );
0166 //                   if (std::abs(*ip) <= epsilon)
0167 //                      *ip=0;
0168                   ip++;
0169                }
0170             }
0171             // update L
0172             ip = rhs.Array() + (j+1)*j/2 + j-1;
0173             for (i=j+1; i <= nrow; ip += i++)
0174                *ip *= static_cast<T> ( temp2 );
0175          }
0176          else if (s==1) // 1x1 pivot
0177          {
0178             piv[j-1] = pivrow;
0179
0180             // interchange rows and columns j and pivrow in
0181             // submatrix (j:n,j:n)
0182             ip = rhs.Array() + pivrow*(pivrow-1)/2 + j;
0183             for (i=j+1; i < pivrow; i++, ip++)
0184             {
0185                temp1 = *(rhs.Array() + i*(i-1)/2 + j-1);
0186                *(rhs.Array() + i*(i-1)/2 + j-1)= *ip;
0187                *ip = static_cast<T> ( temp1 );
0188             }
0189             temp1 = *mjj;
0190             *mjj = *(rhs.Array()+pivrow*(pivrow-1)/2+pivrow-1);
0191             *(rhs.Array()+pivrow*(pivrow-1)/2+pivrow-1) = static_cast<T> (temp1 );
0192             ip = rhs.Array() + (pivrow+1)*pivrow/2 + j-1;
0193             iq = ip + pivrow-j;
0194             for (i = pivrow+1; i <= nrow; ip += i, iq += i++)
0195             {
0196                temp1 = *iq;
0197                *iq = *ip;
0198                *ip = static_cast<T>( temp1 );
0199             }
0200
0201             if (*mjj == 0)
0202             {
0203                ifail = 1;
0204                return;
0205             }
0206             temp2 = *mjj = 1.0f / *mjj; // invert D(j,j)
0207
0208             // update A(j+1:n, j+1:n)
0209             for (i = j+1; i <= nrow; i++)
0210             {
0211                temp1 = *(rhs.Array() + i*(i-1)/2 + j-1) * temp2;
0212                ip = rhs.Array()+i*(i-1)/2+j;
0213                for (k=j+1; k<=i; k++)
0214                {
0215                   *ip -= static_cast<T> (temp1 * *(rhs.Array() + k*(k-1)/2 + j-1) );
0216 //                   if (std::abs(*ip) <= epsilon)
0217 //                      *ip=0;
0218                   ip++;
0219                }
0220             }
0221             // update L
0222             ip = rhs.Array() + (j+1)*j/2 + j-1;
0223             for (i=j+1; i<=nrow; ip += i++)
0224                *ip *= static_cast<T>( temp2 );
0225          }
0226          else // s=2, ie use a 2x2 pivot
0227          {
0228             piv[j-1] = -pivrow;
0229             piv[j] = 0; // that means this is the second row of a 2x2 pivot
0230
0231             if (j+1 != pivrow)
0232             {
0233                // interchange rows and columns j+1 and pivrow in
0234                // submatrix (j:n,j:n)
0235                ip = rhs.Array() + pivrow*(pivrow-1)/2 + j+1;
0236                for (i=j+2; i < pivrow; i++, ip++)
0237                {
0238                   temp1 = *(rhs.Array() + i*(i-1)/2 + j);
0239                   *(rhs.Array() + i*(i-1)/2 + j) = *ip;
0240                   *ip = static_cast<T>( temp1 );
0241                }
0242                temp1 = *(mjj + j + 1);
0243                *(mjj + j + 1) =
0244                   *(rhs.Array() + pivrow*(pivrow-1)/2 + pivrow-1);
0245                *(rhs.Array() + pivrow*(pivrow-1)/2 + pivrow-1) = static_cast<T>( temp1 );
0246                temp1 = *(mjj + j);
0247                *(mjj + j) = *(rhs.Array() + pivrow*(pivrow-1)/2 + j-1);
0248                *(rhs.Array() + pivrow*(pivrow-1)/2 + j-1) = static_cast<T>( temp1 );
0249                ip = rhs.Array() + (pivrow+1)*pivrow/2 + j;
0250                iq = ip + pivrow-(j+1);
0251                for (i = pivrow+1; i <= nrow; ip += i, iq += i++)
0252                {
0253                   temp1 = *iq;
0254                   *iq = *ip;
0255                   *ip = static_cast<T>( temp1 );
0256                }
0257             }
0258             // invert D(j:j+1,j:j+1)
0259             temp2 = *mjj * *(mjj + j + 1) - *(mjj + j) * *(mjj + j);
0260             if (temp2 == 0)
0261                std::cerr
0262                   << "SymMatrix::bunch_invert: error in pivot choice"
0263                   << std::endl;
0264             temp2 = 1. / temp2;
0265             // this quotient is guaranteed to exist by the choice
0266             // of the pivot
0267             temp1 = *mjj;
0268             *mjj = static_cast<T>( *(mjj + j + 1) * temp2 );
0269             *(mjj + j + 1) = static_cast<T>( temp1 * temp2 );
0270             *(mjj + j) = static_cast<T>( - *(mjj + j) * temp2 );
0271
0272             if (j < nrow-1) // otherwise do nothing
0273             {
0274                // update A(j+2:n, j+2:n)
0275                for (i=j+2; i <= nrow ; i++)
0276                {
0277                   ip = rhs.Array() + i*(i-1)/2 + j-1;
0278                   temp1 = *ip * *mjj + *(ip + 1) * *(mjj + j);
0279 //                   if (std::abs(temp1 ) <= epsilon)
0280 //                      temp1 = 0;
0281                   temp2 = *ip * *(mjj + j) + *(ip + 1) * *(mjj + j + 1);
0282 //                   if (std::abs(temp2 ) <= epsilon)
0283 //                      temp2 = 0;
0284                   for (k = j+2; k <= i ; k++)
0285                   {
0286                      ip = rhs.Array() + i*(i-1)/2 + k-1;
0287                      iq = rhs.Array() + k*(k-1)/2 + j-1;
0288                      *ip -= static_cast<T>( temp1 * *iq + temp2 * *(iq+1) );
0289 //                      if (std::abs(*ip) <= epsilon)
0290 //                         *ip = 0;
0291                   }
0292                }
0293                // update L
0294                for (i=j+2; i <= nrow ; i++)
0295                {
0296                   ip = rhs.Array() + i*(i-1)/2 + j-1;
0297                   temp1 = *ip * *mjj + *(ip+1) * *(mjj + j);
0298 //                   if (std::abs(temp1) <= epsilon)
0299 //                      temp1 = 0;
0300                   *(ip+1) = *ip * *(mjj + j)
0301                      + *(ip+1) * *(mjj + j + 1);
0302 //                   if (std::abs(*(ip+1)) <= epsilon)
0303 //                      *(ip+1) = 0;
0304                   *ip = static_cast<T>( temp1 );
0305                }
0306             }
0307          }
0308       }
0309    } // end of main loop over columns
0310
0311    if (j == nrow) // the last pivot is 1x1
0312    {
0313       mjj = rhs.Array() + j*(j-1)/2 + j-1;
0314       if (*mjj == 0)
0315       {
0316          ifail = 1;
0317          return;
0318       }
0319       else
0320          *mjj = 1.0f  / *mjj;
0321    } // end of last pivot code
0322
0323    // computing the inverse from the factorization
0324
0325    for (j = nrow ; j >= 1 ; j -= s) // loop over columns
0326    {
0327       mjj = rhs.Array() + j*(j-1)/2 + j-1;
0328       if (piv[j-1] > 0) // 1x1 pivot, compute column j of inverse
0329       {
0330          s = 1;
0331          if (j < nrow)
0332          {
0333             ip = rhs.Array() + (j+1)*j/2 + j-1;
0334             for (i=0; i < nrow-j; ip += 1+j+i++)
0335                x[i] = *ip;
0336             for (i=j+1; i<=nrow ; i++)
0337             {
0338                temp2=0;
0339                ip = rhs.Array() + i*(i-1)/2 + j;
0340                for (k=0; k <= i-j-1; k++)
0341                   temp2 += *ip++ * x[k];
0342                for (ip += i-1; k < nrow-j; ip += 1+j+k++)
0343                   temp2 += *ip * x[k];
0344                *(rhs.Array()+ i*(i-1)/2 + j-1) = static_cast<T>( -temp2 );
0345             }
0346             temp2 = 0;
0347             ip = rhs.Array() + (j+1)*j/2 + j-1;
0348             for (k=0; k < nrow-j; ip += 1+j+k++)
0349                temp2 += x[k] * *ip;
0350             *mjj -= static_cast<T>( temp2 );
0351          }
0352       }
0353       else //2x2 pivot, compute columns j and j-1 of the inverse
0354       {
0355          if (piv[j-1] != 0)
0356             std::cerr << "error in piv" << piv[j-1] << std::endl;
0357          s=2;
0358          if (j < nrow)
0359          {
0360             ip = rhs.Array() + (j+1)*j/2 + j-1;
0361             for (i=0; i < nrow-j; ip += 1+j+i++)
0362                x[i] = *ip;
0363             for (i=j+1; i<=nrow ; i++)
0364             {
0365                temp2 = 0;
0366                ip = rhs.Array() + i*(i-1)/2 + j;
0367                for (k=0; k <= i-j-1; k++)
0368                   temp2 += *ip++ * x[k];
0369                for (ip += i-1; k < nrow-j; ip += 1+j+k++)
0370                   temp2 += *ip * x[k];
0371                *(rhs.Array()+ i*(i-1)/2 + j-1) = static_cast<T>( -temp2 );
0372             }
0373             temp2 = 0;
0374             ip = rhs.Array() + (j+1)*j/2 + j-1;
0375             for (k=0; k < nrow-j; ip += 1+j+k++)
0376                temp2 += x[k] * *ip;
0377             *mjj -= static_cast<T>( temp2 );
0378             temp2 = 0;
0379             ip = rhs.Array() + (j+1)*j/2 + j-2;
0380             for (i=j+1; i <= nrow; ip += i++)
0381                temp2 += *ip * *(ip+1);
0382             *(mjj-1) -= static_cast<T>( temp2 );
0383             ip = rhs.Array() + (j+1)*j/2 + j-2;
0384             for (i=0; i < nrow-j; ip += 1+j+i++)
0385                x[i] = *ip;
0386             for (i=j+1; i <= nrow ; i++)
0387             {
0388                temp2 = 0;
0389                ip = rhs.Array() + i*(i-1)/2 + j;
0390                for (k=0; k <= i-j-1; k++)
0391                   temp2 += *ip++ * x[k];
0392                for (ip += i-1; k < nrow-j; ip += 1+j+k++)
0393                   temp2 += *ip * x[k];
0394                *(rhs.Array()+ i*(i-1)/2 + j-2)= static_cast<T>( -temp2 );
0395             }
0396             temp2 = 0;
0397             ip = rhs.Array() + (j+1)*j/2 + j-2;
0398             for (k=0; k < nrow-j; ip += 1+j+k++)
0399                temp2 += x[k] * *ip;
0400             *(mjj-j) -= static_cast<T>( temp2 );
0401          }
0402       }
0403
0404       // interchange rows and columns j and piv[j-1]
0405       // or rows and columns j and -piv[j-2]
0406
0407       pivrow = (piv[j-1]==0)? -piv[j-2] : piv[j-1];
0408       ip = rhs.Array() + pivrow*(pivrow-1)/2 + j;
0409       for (i=j+1;i < pivrow; i++, ip++)
0410       {
0411          temp1 = *(rhs.Array() + i*(i-1)/2 + j-1);
0412          *(rhs.Array() + i*(i-1)/2 + j-1) = *ip;
0413          *ip = static_cast<T>( temp1 );
0414       }
0415       temp1 = *mjj;
0416       *mjj = *(rhs.Array() + pivrow*(pivrow-1)/2 + pivrow-1);
0417       *(rhs.Array() + pivrow*(pivrow-1)/2 + pivrow-1) = static_cast<T>( temp1 );
0418       if (s==2)
0419       {
0420          temp1 = *(mjj-1);
0421          *(mjj-1) = *( rhs.Array() + pivrow*(pivrow-1)/2 + j-2);
0422          *( rhs.Array() + pivrow*(pivrow-1)/2 + j-2) = static_cast<T>( temp1 );
0423       }
0424
0425       ip = rhs.Array() + (pivrow+1)*pivrow/2 + j-1;  // &A(i,j)
0426       iq = ip + pivrow-j;
0427       for (i = pivrow+1; i <= nrow; ip += i, iq += i++)
0428       {
0429          temp1 = *iq;
0430          *iq = *ip;
0431          *ip = static_cast<T>(temp1);
0432       }
0433    } // end of loop over columns (in computing inverse from factorization)
0434
0435    return; // inversion successful
0436
0437 }
0438
0439
0440
0441 /**
0442    LU factorization : code originally from CERNLIB dfact routine and ported in C++ for CLHEP
0443 */
0444
0445 template <unsigned int idim, unsigned int n>
0446 template<class T>
0447 int Inverter<idim,n>::DfactMatrix(MatRepStd<T,idim,n> & rhs, T &det, unsigned int *ir) {
0448
0449    if (idim != n) return   -1;
0450
0451    int ifail, jfail;
0452
0453    typedef T* mIter;
0454
0455    typedef T value_type;
0456    //typedef double value_type;  // for making inversions in double's
0457
0458
0459    value_type tf;
0460    value_type g1 = 1.0e-19, g2 = 1.0e19;
0461
0462    value_type p, q, t;
0463    value_type s11, s12;
0464
0465    // LM (04.09) : remove useless check on epsilon and set it to zero
0466    const value_type epsilon = 0.0;
0467    //double epsilon = 8*std::numeric_limits<T>::epsilon();
0468    // could be set to zero (like it was before)
0469    // but then the algorithm often doesn't detect
0470    // that a matrix is singular
0471
0472    int normal = 0, imposs = -1;
0473    int jrange = 0, jover = 1, junder = -1;
0474    ifail = normal;
0475    jfail = jrange;
0476    int nxch = 0;
0477    det = 1.0;
0478    mIter mj = rhs.Array();
0479    mIter mjj = mj;
0480    for (unsigned int j=1;j<=n;j++) {
0481       unsigned int k = j;
0482       p = (std::abs(*mjj));
0483       if (j!=n) {
0484          mIter mij = mj + n + j - 1;
0485          for (unsigned int i=j+1;i<=n;i++) {
0486             q = (std::abs(*(mij)));
0487             if (q > p) {
0488                k = i;
0489                p = q;
0490             }
0491             mij += n;
0492          }
0493          if (k==j) {
0494             if (p <= epsilon) {
0495                det = 0;
0496                ifail = imposs;
0497                jfail = jrange;
0498                return ifail;
0499             }
0500             det = -det; // in this case the sign of the determinant
0501             // must not change. So I change it twice.
0502          }
0503          mIter mjl = mj;
0504          mIter mkl = rhs.Array() + (k-1)*n;
0505          for (unsigned int l=1;l<=n;l++) {
0506             tf = *mjl;
0507             *(mjl++) = *mkl;
0508             *(mkl++) = static_cast<T>(tf);
0509          }
0510          nxch = nxch + 1;  // this makes the determinant change its sign
0511          ir[nxch] = (((j)<<12)+(k));
0512       } else {
0513          if (p <= epsilon) {
0514             det = 0.0;
0515             ifail = imposs;
0516             jfail = jrange;
0517             return ifail;
0518          }
0519       }
0520       det *= *mjj;
0521       *mjj = 1.0f / *mjj;
0522       t = (std::abs(det));
0523       if (t < g1) {
0524          det = 0.0;
0525          if (jfail == jrange) jfail = junder;
0526       } else if (t > g2) {
0527          det = 1.0;
0528          if (jfail==jrange) jfail = jover;
0529       }
0530       if (j!=n) {
0531          mIter mk = mj + n;
0532          mIter mkjp = mk + j;
0533          mIter mjk = mj + j;
0534          for (k=j+1;k<=n;k++) {
0535             s11 = - (*mjk);
0536             s12 = - (*mkjp);
0537             if (j!=1) {
0538                mIter mik = rhs.Array() + k - 1;
0539                mIter mijp = rhs.Array() + j;
0540                mIter mki = mk;
0541                mIter mji = mj;
0542                for (unsigned int i=1;i<j;i++) {
0543                   s11 += (*mik) * (*(mji++));
0544                   s12 += (*mijp) * (*(mki++));
0545                   mik += n;
0546                   mijp += n;
0547                }
0548             }
0549             // cast to avoid warnings from double to float conversions
0550             *(mjk++) = static_cast<T>( - s11 * (*mjj) );
0551             *(mkjp) = static_cast<T> ( -(((*(mjj+1)))*((*(mkjp-1)))+(s12)) );
0552             mk += n;
0553             mkjp += n;
0554          }
0555       }
0556       mj += n;
0557       mjj += (n+1);
0558    }
0559    if (nxch%2==1) det = -det;
0560    if (jfail !=jrange) det = 0.0;
0561    ir[n] = nxch;
0562    return 0;
0563 }
0564
0565
0566
0567     /**
0568     Inversion for General square matrices.
0569     Code from  dfinv routine from CERNLIB
0570     Assumed first the LU decomposition via DfactMatrix function
0571
0572     taken from CLHEP : L. Moneta May 2006
0573     */
0574
0575 template <unsigned int idim, unsigned int n>
0576 template<class T>
0577 int Inverter<idim,n>::DfinvMatrix(MatRepStd<T,idim,n> & rhs,unsigned int * ir) {
0578
0579
0580    typedef T* mIter;
0581
0582    typedef T value_type;
0583    //typedef double value_type;  // for making inversions in double's
0584
0585
0586    if (idim != n) return   -1;
0587
0588
0589    value_type s31, s32;
0590    value_type s33, s34;
0591
0592    mIter m11 = rhs.Array();
0593    mIter m12 = m11 + 1;
0594    mIter m21 = m11 + n;
0595    mIter m22 = m12 + n;
0596    *m21 = -(*m22) * (*m11) * (*m21);
0597    *m12 = -(*m12);
0598    if (n>2) {
0599       mIter mi = rhs.Array() + 2 * n;
0600       mIter mii= rhs.Array() + 2 * n + 2;
0601       mIter mimim = rhs.Array() + n + 1;
0602       for (unsigned int i=3;i<=n;i++) {
0603          unsigned int im2 = i - 2;
0604          mIter mj = rhs.Array();
0605          mIter mji = mj + i - 1;
0606          mIter mij = mi;
0607          for (unsigned int j=1;j<=im2;j++) {
0608             s31 = 0.0;
0609             s32 = *mji;
0610             mIter mkj = mj + j - 1;
0611             mIter mik = mi + j - 1;
0612             mIter mjkp = mj + j;
0613             mIter mkpi = mj + n + i - 1;
0614             for (unsigned int k=j;k<=im2;k++) {
0615                s31 += (*mkj) * (*(mik++));
0616                s32 += (*(mjkp++)) * (*mkpi);
0617                mkj += n;
0618                mkpi += n;
0619             }
0620             *mij = static_cast<T>( -(*mii) * (((*(mij-n)))*( (*(mii-1)))+(s31)) );
0621             *mji = static_cast<T> ( -s32 );
0622             mj += n;
0623             mji += n;
0624             mij++;
0625          }
0626          *(mii-1) = -(*mii) * (*mimim) * (*(mii-1));
0627          *(mimim+1) = -(*(mimim+1));
0628          mi += n;
0629          mimim += (n+1);
0630          mii += (n+1);
0631       }
0632    }
0633    mIter mi = rhs.Array();
0634    mIter mii = rhs.Array();
0635    for (unsigned  int i=1;i<n;i++) {
0636       unsigned int ni = n - i;
0637       mIter mij = mi;
0638       //int j;
0639       for (unsigned j=1; j<=i;j++) {
0640          s33 = *mij;
0641          mIter mikj = mi + n + j - 1;
0642          mIter miik = mii + 1;
0643          mIter min_end = mi + n;
0644          for (;miik<min_end;) {
0645             s33 += (*mikj) * (*(miik++));
0646             mikj += n;
0647          }
0648          *(mij++) = static_cast<T> ( s33 );
0649       }
0650       for (unsigned j=1;j<=ni;j++) {
0651          s34 = 0.0;
0652          mIter miik = mii + j;
0653          mIter mikij = mii + j * n + j;
0654          for (unsigned int k=j;k<=ni;k++) {
0655             s34 += *mikij * (*(miik++));
0656             mikij += n;
0657          }
0658          *(mii+j) = s34;
0659       }
0660       mi += n;
0661       mii += (n+1);
0662    }
0663    unsigned int nxch = ir[n];
0664    if (nxch==0) return 0;
0665    for (unsigned int mm=1;mm<=nxch;mm++) {
0666       unsigned int k = nxch - mm + 1;
0667       int ij = ir[k];
0668       int i = ij >> 12;
0669       int j = ij%4096;
0670       mIter mki = rhs.Array() + i - 1;
0671       mIter mkj = rhs.Array() + j - 1;
0672       for (k=1; k<=n;k++) {
0673          // 2/24/05 David Sachs fix of improper swap bug that was present
0674          // for many years:
0675          T ti = *mki; // 2/24/05
0676          *mki = *mkj;
0677          *mkj = ti; // 2/24/05
0678          mki += n;
0679          mkj += n;
0680       }
0681    }
0682    return 0;
0683 }
0684
0685
0686   }  // end namespace Math
0687 }    // end namespace ROOT
0688
0689 #endif