ROOT logo
// @(#)root/matrix:$Id: TDecompLU.cxx 23492 2008-04-23 20:42:49Z brun $
// Authors: Fons Rademakers, Eddy Offermann  Dec 2003

/*************************************************************************
 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

#include "TDecompLU.h"
#include "TMath.h"

ClassImp(TDecompLU)

///////////////////////////////////////////////////////////////////////////
//                                                                       //
// LU Decomposition class                                                //
//                                                                       //
// Decompose  a general n x n matrix A into P A = L U                    //
//                                                                       //
// where P is a permutation matrix, L is unit lower triangular and U     //
// is upper triangular.                                                  //
// L is stored in the strict lower triangular part of the matrix fLU.    //
// The diagonal elements of L are unity and are not stored.              //
// U is stored in the diagonal and upper triangular part of the matrix   //
// fU.                                                                   //
// P is stored in the index array fIndex : j = fIndex[i] indicates that  //
// row j and row i should be swapped .                                   //
//                                                                       //
// fSign gives the sign of the permutation, (-1)^n, where n is the       //
// number of interchanges in the permutation.                            //
//                                                                       //
// fLU has the same indexing range as matrix A .                         //
//                                                                       //
// The decomposition fails if a diagonal element of abs(fLU) is == 0,    //
// The matrix fUL is made invalid .                                      //
//                                                                       //
///////////////////////////////////////////////////////////////////////////

//______________________________________________________________________________
TDecompLU::TDecompLU()
{
// Default constructor

   fSign = 0.0;
   fNIndex = 0;
   fIndex = 0;
   fImplicitPivot = 0;
}

//______________________________________________________________________________
TDecompLU::TDecompLU(Int_t nrows)
{
// Constructor for (nrows x nrows) matrix

   fSign = 1.0;
   fNIndex = nrows;
   fIndex = new Int_t[fNIndex];
   memset(fIndex,0,fNIndex*sizeof(Int_t));
   fImplicitPivot = 0;
   fLU.ResizeTo(nrows,nrows);
}

//______________________________________________________________________________
TDecompLU::TDecompLU(Int_t row_lwb,Int_t row_upb)
{
// Constructor for ([row_lwb..row_upb] x [row_lwb..row_upb]) matrix

   const Int_t nrows = row_upb-row_lwb+1;
   fSign = 1.0;
   fNIndex = nrows;
   fIndex = new Int_t[fNIndex];
   memset(fIndex,0,fNIndex*sizeof(Int_t));
   fRowLwb = row_lwb;
   fColLwb = row_lwb;
   fImplicitPivot = 0;
   fLU.ResizeTo(row_lwb,row_lwb+nrows-1,row_lwb,row_lwb+nrows-1);
}

//______________________________________________________________________________
TDecompLU::TDecompLU(const TMatrixD &a,Double_t tol,Int_t implicit)
{
// Constructor for matrix a

   R__ASSERT(a.IsValid());

   if (a.GetNrows() != a.GetNcols() || a.GetRowLwb() != a.GetColLwb()) {
      Error("TDecompLU(const TMatrixD &","matrix should be square");
      return;
   }

   SetBit(kMatrixSet);
   fCondition = a.Norm1();
   fImplicitPivot = implicit;
   fTol = a.GetTol();
   if (tol > 0)
      fTol = tol;

   fSign = 1.0;
   fNIndex = a.GetNcols();
   fIndex = new Int_t[fNIndex];
   memset(fIndex,0,fNIndex*sizeof(Int_t));

   fRowLwb = a.GetRowLwb();
   fColLwb = a.GetColLwb();
   fLU.ResizeTo(a);
   fLU = a;
}

//______________________________________________________________________________
TDecompLU::TDecompLU(const TDecompLU &another) : TDecompBase(another)
{
// Copy constructor

   fNIndex = 0;
   fIndex  = 0;
   *this = another;
}

//______________________________________________________________________________
Bool_t TDecompLU::Decompose()
{
// Matrix A is decomposed in components U and L so that P * A = U * L
// If the decomposition succeeds, bit kDecomposed is set , otherwise kSingular

   if (TestBit(kDecomposed)) return kTRUE;

   if ( !TestBit(kMatrixSet) ) {
      Error("Decompose()","Matrix has not been set");
      return kFALSE;
   }

   Int_t nrZeros = 0;
   Bool_t ok;
   if (fImplicitPivot)
      ok = DecomposeLUCrout(fLU,fIndex,fSign,fTol,nrZeros);
   else
      ok = DecomposeLUGauss(fLU,fIndex,fSign,fTol,nrZeros);

   if (!ok) SetBit(kSingular);
   else     SetBit(kDecomposed);

   return ok;
}

//______________________________________________________________________________
const TMatrixD TDecompLU::GetMatrix()
{
// Reconstruct the original matrix using the decomposition parts

   if (TestBit(kSingular)) {
      Error("GetMatrix()","Matrix is singular");
      return TMatrixD();
   }
   if ( !TestBit(kDecomposed) ) {
      if (!Decompose()) {
         Error("GetMatrix()","Decomposition failed");
         return TMatrixD();
      }
   }

   TMatrixD mL = fLU;
   TMatrixD mU = fLU;
   Double_t * const pU = mU.GetMatrixArray();
   Double_t * const pL = mL.GetMatrixArray();
   const Int_t n = fLU.GetNcols();
   for (Int_t irow = 0; irow < n; irow++) {
      const Int_t off_row = irow*n;
         for (Int_t icol = 0; icol < n; icol++) {
            if (icol > irow)      pL[off_row+icol] = 0.;
            else if (icol < irow) pU[off_row+icol] = 0.;
            else                  pL[off_row+icol] = 1.;
         }
      }

   TMatrixD a = mL*mU;

   // swap rows

   Double_t * const pA = a.GetMatrixArray();
   for (Int_t i = n-1; i >= 0; i--) {
      const Int_t j = fIndex[i];
      if (j != i) {
         const Int_t off_j = j*n;
         const Int_t off_i = i*n;
         for (Int_t k = 0; k < n; k++) {
            const Double_t tmp = pA[off_j+k];
            pA[off_j+k] = pA[off_i+k];
            pA[off_i+k] = tmp;
         }
      }
   }

   return a;
}

//______________________________________________________________________________
void TDecompLU::SetMatrix(const TMatrixD &a)
{
// Set matrix to be decomposed

   R__ASSERT(a.IsValid());

   ResetStatus();
   if (a.GetNrows() != a.GetNcols() || a.GetRowLwb() != a.GetColLwb()) {
      Error("TDecompLU(const TMatrixD &","matrix should be square");
      return;
   }

   SetBit(kMatrixSet);
   fCondition = a.Norm1();

   fSign = 1.0;
   if (fNIndex != a.GetNcols()) {
      fNIndex = a.GetNcols();
      delete [] fIndex;
      fIndex = new Int_t[fNIndex];
      memset(fIndex,0,fNIndex*sizeof(Int_t));
   }

   fRowLwb = a.GetRowLwb();
   fColLwb = a.GetColLwb();
   fLU.ResizeTo(a);
   fLU = a;
}

//______________________________________________________________________________
Bool_t TDecompLU::Solve(TVectorD &b)
{
// Solve Ax=b assuming the LU form of A is stored in fLU, but assume b has *not*
// been transformed.  Solution returned in b.

   R__ASSERT(b.IsValid());
   if (TestBit(kSingular)) {
      Error("Solve()","Matrix is singular");
      return kFALSE;
   }
   if ( !TestBit(kDecomposed) ) {
      if (!Decompose()) {
         Error("Solve()","Decomposition failed");
         return kFALSE;
      }
   }

   if (fLU.GetNrows() != b.GetNrows() || fLU.GetRowLwb() != b.GetLwb()) {
      Error("Solve(TVectorD &","vector and matrix incompatible");
      return kFALSE;
   }

   const Int_t n = fLU.GetNrows();

   const Double_t *pLU = fLU.GetMatrixArray();
   Double_t *pb  = b.GetMatrixArray();

   Int_t i;

   // Check for zero diagonals
   for (i = 0; i < n ; i++) {
      const Int_t off_i = i*n;
      if (TMath::Abs(pLU[off_i+i]) < fTol) {
         Error("Solve(TVectorD &b)","LU[%d,%d]=%.4e < %.4e",i,i,pLU[off_i+i],fTol);
         return kFALSE;
      }
   }

   // Transform b allowing for leading zeros
   Int_t nonzero = -1;
   for (i = 0; i < n; i++) {
      const Int_t off_i = i*n;
      const Int_t iperm = fIndex[i];
      Double_t r = pb[iperm];
      pb[iperm] = pb[i];
      if (nonzero >= 0)
         for (Int_t j = nonzero; j < i; j++)
            r -= pLU[off_i+j]*pb[j];
      else if (r != 0.0)
         nonzero = i;
      pb[i] = r;
   }

   // Backward substitution
   for (i = n-1; i >= 0; i--) {
      const Int_t off_i = i*n;
      Double_t r = pb[i];
      for (Int_t j = i+1; j < n; j++)
         r -= pLU[off_i+j]*pb[j];
      pb[i] = r/pLU[off_i+i];
   }

   return kTRUE;
}

//______________________________________________________________________________
Bool_t TDecompLU::Solve(TMatrixDColumn &cb)
{
// Solve Ax=b assuming the LU form of A is stored in fLU, but assume b has *not*
// been transformed.  Solution returned in b.

   TMatrixDBase *b = const_cast<TMatrixDBase *>(cb.GetMatrix());
   R__ASSERT(b->IsValid());
   if (TestBit(kSingular)) {
      Error("Solve()","Matrix is singular");
      return kFALSE;
   }
   if ( !TestBit(kDecomposed) ) {
      if (!Decompose()) {
         Error("Solve()","Decomposition failed");
         return kFALSE;
      }
   }

   if (fLU.GetNrows() != b->GetNrows() || fLU.GetRowLwb() != b->GetRowLwb()) {
      Error("Solve(TMatrixDColumn &","vector and matrix incompatible");
      return kFALSE;
   }

   const Int_t     n   = fLU.GetNrows();
   const Double_t *pLU = fLU.GetMatrixArray();
 
   Int_t i;

   // Check for zero diagonals
   for (i = 0; i < n ; i++) {
      const Int_t off_i = i*n;
      if (TMath::Abs(pLU[off_i+i]) < fTol) {
         Error("Solve(TMatrixDColumn &cb)","LU[%d,%d]=%.4e < %.4e",i,i,pLU[off_i+i],fTol);
         return kFALSE;
      }
   }

   Double_t *pcb = cb.GetPtr();
   const Int_t     inc = cb.GetInc();

   // Transform b allowing for leading zeros
   Int_t nonzero = -1;
   for (i = 0; i < n; i++) {
      const Int_t off_i  = i*n;
      const Int_t off_i2 = i*inc;
      const Int_t iperm = fIndex[i];
      const Int_t off_iperm = iperm*inc;
      Double_t r = pcb[off_iperm];
      pcb[off_iperm] = pcb[off_i2];
      if (nonzero >=0)
         for (Int_t j = nonzero; j <= i-1; j++)
            r -= pLU[off_i+j]*pcb[j*inc];
      else if (r != 0.0)
         nonzero = i;
      pcb[off_i2] = r;
   }

   // Backward substitution
   for (i = n-1; i >= 0; i--) {
      const Int_t off_i  = i*n;
      const Int_t off_i2 = i*inc;
      Double_t r = pcb[off_i2];
      for (Int_t j = i+1; j < n; j++)
         r -= pLU[off_i+j]*pcb[j*inc];
      pcb[off_i2] = r/pLU[off_i+i];
   }

   return kTRUE;
}


//______________________________________________________________________________
Bool_t TDecompLU::TransSolve(TVectorD &b)
{
// Solve A^T x=b assuming the LU form of A^T is stored in fLU, but assume b has *not*
// been transformed.  Solution returned in b.

   R__ASSERT(b.IsValid());
   if (TestBit(kSingular)) {
      Error("TransSolve()","Matrix is singular");
      return kFALSE;
   }
   if ( !TestBit(kDecomposed) ) {
      if (!Decompose()) {
         Error("TransSolve()","Decomposition failed");
         return kFALSE;
      }
   }

   if (fLU.GetNrows() != b.GetNrows() || fLU.GetRowLwb() != b.GetLwb()) {
      Error("TransSolve(TVectorD &","vector and matrix incompatible");
      return kFALSE;
   }

   const Int_t n = fLU.GetNrows();

   const Double_t *pLU = fLU.GetMatrixArray();
   Double_t *pb  = b.GetMatrixArray();

   Int_t i;

   // Check for zero diagonals
   for (i = 0; i < n ; i++) {
      const Int_t off_i = i*n;
      if (TMath::Abs(pLU[off_i+i]) < fTol) {
         Error("TransSolve(TVectorD &b)","LU[%d,%d]=%.4e < %.4e",i,i,pLU[off_i+i],fTol);
         return kFALSE;
      }
   }

   // Forward Substitution
   for (i = 0; i < n; i++) {
      const Int_t off_i = i*n;
      Double_t r = pb[i];
      for (Int_t j = 0; j < i ; j++) {
         const Int_t off_j = j*n;
         r -= pLU[off_j+i]*pb[j];
      }
      pb[i] = r/pLU[off_i+i];
   }

   // Transform b allowing for leading zeros
   Int_t nonzero = -1;
   for (i = n-1 ; i >= 0; i--) {
      Double_t r = pb[i];
      if (nonzero >= 0) {
         for (Int_t j = i+1; j <= nonzero; j++) {
            const Int_t off_j = j*n;
            r -= pLU[off_j+i]*pb[j];
         }
      } else if (r != 0.0)
         nonzero = i;
      const Int_t iperm = fIndex[i];
      pb[i]     = pb[iperm];
      pb[iperm] = r;
   }

   return kTRUE;
}

//______________________________________________________________________________
Bool_t TDecompLU::TransSolve(TMatrixDColumn &cb)
{
// Solve A^T x=b assuming the LU form of A^T is stored in fLU, but assume b has *not*
// been transformed.  Solution returned in b.

   TMatrixDBase *b = const_cast<TMatrixDBase *>(cb.GetMatrix());
   R__ASSERT(b->IsValid());
   if (TestBit(kSingular)) {
      Error("TransSolve()","Matrix is singular");
      return kFALSE;
   }
   if ( !TestBit(kDecomposed) ) {
      if (!Decompose()) {
         Error("TransSolve()","Decomposition failed");
         return kFALSE;
      }
   }

   if (fLU.GetNrows() != b->GetNrows() || fLU.GetRowLwb() != b->GetRowLwb()) {
      Error("TransSolve(TMatrixDColumn &","vector and matrix incompatible");
      return kFALSE;
   }

   const Int_t n   = fLU.GetNrows();
   const Int_t lwb = fLU.GetRowLwb();
 
   const Double_t *pLU = fLU.GetMatrixArray();

   Int_t i;

   // Check for zero diagonals
   for (i = 0; i < n ; i++) {
      const Int_t off_i = i*n;
      if (TMath::Abs(pLU[off_i+i]) < fTol) {
         Error("TransSolve(TMatrixDColumn &cb)","LU[%d,%d]=%.4e < %.4e",i,i,pLU[off_i+i],fTol);
         return kFALSE;
      }
   }

   // Forward Substitution
   for (i = 0; i < n; i++) {
      const Int_t off_i = i*n;
      Double_t r = cb(i+lwb);
      for (Int_t j = 0; j < i ; j++) {
         const Int_t off_j = j*n;
         r -= pLU[off_j+i]*cb(j+lwb);
      }
      cb(i+lwb) = r/pLU[off_i+i];
   }

   // Transform b allowing for leading zeros
   Int_t nonzero = -1;
   for (i = n-1 ; i >= 0; i--) {
      Double_t r = cb(i+lwb);
      if (nonzero >= 0) {
         for (Int_t j = i+1; j <= nonzero; j++) {
            const Int_t off_j = j*n;
            r -= pLU[off_j+i]*cb(j+lwb);
         }
      } else if (r != 0.0)
         nonzero = i;
      const Int_t iperm = fIndex[i];
      cb(i+lwb)     = cb(iperm+lwb);
      cb(iperm+lwb) = r;
   }

   return kTRUE;
}

//______________________________________________________________________________
void TDecompLU::Det(Double_t &d1,Double_t &d2)
{
// Calculate determinant det = d1*TMath::Power(2.,d2) 

   if ( !TestBit(kDetermined) ) {
      if ( !TestBit(kDecomposed) )
         Decompose();
      TDecompBase::Det(d1,d2);
      fDet1 *= fSign;
      SetBit(kDetermined);
   }
   d1 = fDet1;
   d2 = fDet2;
}

//______________________________________________________________________________
Bool_t TDecompLU::Invert(TMatrixD &inv)
{
// For a matrix A(m,m), its inverse A_inv is defined as A * A_inv = A_inv * A = unit
// (m x m) Ainv is returned .

   if (inv.GetNrows()  != GetNrows()  || inv.GetNcols()  != GetNcols() ||
        inv.GetRowLwb() != GetRowLwb() || inv.GetColLwb() != GetColLwb()) {
      Error("Invert(TMatrixD &","Input matrix has wrong shape");
      return kFALSE;
   }

   inv.UnitMatrix();
   const Bool_t status = MultiSolve(inv);

   return status;
}

//______________________________________________________________________________
TMatrixD TDecompLU::Invert(Bool_t &status)
{
// For a matrix A(m,n), its inverse A_inv is defined as A * A_inv = A_inv * A = unit
// (n x m) Ainv is returned .

   const Int_t rowLwb = GetRowLwb();
   const Int_t rowUpb = rowLwb+GetNrows()-1;

   TMatrixD inv(rowLwb,rowUpb,rowLwb,rowUpb);
   inv.UnitMatrix();
   status = MultiSolve(inv);

   return inv;
}

//______________________________________________________________________________
void TDecompLU::Print(Option_t *opt) const
{
   //Print internals of this object
   TDecompBase::Print(opt);
   printf("fImplicitPivot = %d\n",fImplicitPivot);
   printf("fSign          = %f\n",fSign);
   printf("fIndex:\n");
   for (Int_t i = 0; i < fNIndex; i++)
      printf("[%d] = %d\n",i,fIndex[i]);
   fLU.Print("fLU");
}

//______________________________________________________________________________
TDecompLU &TDecompLU::operator=(const TDecompLU &source)
{
   //assignement operator
   if (this != &source) {
      TDecompBase::operator=(source);
      fLU.ResizeTo(source.fLU);
      fLU   = source.fLU;
      fSign = source.fSign;
      fImplicitPivot = source.fImplicitPivot;
      if (fNIndex != source.fNIndex) {
         if (fIndex)
            delete [] fIndex;
         fNIndex = source.fNIndex;
         fIndex = new Int_t[fNIndex];
      }
      memcpy(fIndex,source.fIndex,fNIndex*sizeof(Int_t));
   }
   return *this;
}

//______________________________________________________________________________
Bool_t TDecompLU::DecomposeLUCrout(TMatrixD &lu,Int_t *index,Double_t &sign,
                                   Double_t tol,Int_t &nrZeros)
{
// Crout/Doolittle algorithm of LU decomposing a square matrix, with implicit partial
// pivoting.  The decomposition is stored in fLU: U is explicit in the upper triag
// and L is in multiplier form in the subdiagionals .
// Row permutations are mapped out in fIndex. fSign, used for calculating the
// determinant, is +/- 1 for even/odd row permutations. .

   const Int_t     n     = lu.GetNcols();
   Double_t *pLU   = lu.GetMatrixArray();

   Double_t work[kWorkMax];
   Bool_t isAllocated = kFALSE;
   Double_t *scale = work;
   if (n > kWorkMax) {
      isAllocated = kTRUE;
      scale = new Double_t[n];
   }

   sign    = 1.0;
   nrZeros = 0;
   // Find implicit scaling factors for each row
   for (Int_t i = 0; i < n ; i++) {
      const Int_t off_i = i*n;
      Double_t max = 0.0;
      for (Int_t j = 0; j < n; j++) {
         const Double_t tmp = TMath::Abs(pLU[off_i+j]);
         if (tmp > max)
            max = tmp;
      }
      scale[i] = (max == 0.0 ? 0.0 : 1.0/max);
   }

   for (Int_t j = 0; j < n; j++) {
      const Int_t off_j = j*n;
      // Run down jth column from top to diag, to form the elements of U.
      for (Int_t i = 0; i < j; i++) {
         const Int_t off_i = i*n;
         Double_t r = pLU[off_i+j];
         for (Int_t k = 0; k < i; k++) {
            const Int_t off_k = k*n;
            r -= pLU[off_i+k]*pLU[off_k+j];
         }
         pLU[off_i+j] = r;
      }

      // Run down jth subdiag to form the residuals after the elimination of
      // the first j-1 subdiags.  These residuals divided by the appropriate
      // diagonal term will become the multipliers in the elimination of the jth.
      // subdiag. Find fIndex of largest scaled term in imax.

      Double_t max = 0.0;
      Int_t imax = 0;
      for (Int_t i = j; i < n; i++) {
         const Int_t off_i = i*n;
         Double_t r = pLU[off_i+j];
         for (Int_t k = 0; k < j; k++) {
            const Int_t off_k = k*n;
            r -= pLU[off_i+k]*pLU[off_k+j];
         }
         pLU[off_i+j] = r;
         const Double_t tmp = scale[i]*TMath::Abs(r);
         if (tmp >= max) {
            max = tmp;
            imax = i;
         }
      }

      // Permute current row with imax
      if (j != imax) {
         const Int_t off_imax = imax*n;
         for (Int_t k = 0; k < n; k++ ) {
            const Double_t tmp = pLU[off_imax+k];
            pLU[off_imax+k] = pLU[off_j+k];
            pLU[off_j+k]    = tmp;
         }
         sign = -sign;
         scale[imax] = scale[j];
      }
      index[j] = imax;

      // If diag term is not zero divide subdiag to form multipliers.
      if (pLU[off_j+j] != 0.0) {
         if (TMath::Abs(pLU[off_j+j]) < tol)
            nrZeros++;
         if (j != n-1) {
            const Double_t tmp = 1.0/pLU[off_j+j];
            for (Int_t i = j+1; i < n; i++) {
               const Int_t off_i = i*n;
               pLU[off_i+j] *= tmp;
            }
         }
      } else {
         ::Error("TDecompLU::DecomposeLUCrout","matrix is singular");
         return kFALSE;
      }
   }

   if (isAllocated)
      delete [] scale;

   return kTRUE;
}

//______________________________________________________________________________
Bool_t TDecompLU::DecomposeLUGauss(TMatrixD &lu,Int_t *index,Double_t &sign,
                                   Double_t tol,Int_t &nrZeros)
{
// LU decomposition using Gaussain Elimination with partial pivoting (See Golub &
// Van Loan, Matrix Computations, Algorithm 3.4.1) of a square matrix .
// The decomposition is stored in fLU: U is explicit in the upper triag and L is in
// multiplier form in the subdiagionals . Row permutations are mapped out in fIndex.
// fSign, used for calculating the determinant, is +/- 1 for even/odd row permutations.
// Since this algorithm uses partial pivoting without scaling like in Crout/Doolitle.
// it is somewhat faster but less precise .

   const Int_t     n   = lu.GetNcols();
   Double_t *pLU = lu.GetMatrixArray();

   sign    = 1.0;
   nrZeros = 0;

   index[n-1] = n-1;
   for (Int_t j = 0; j < n-1; j++) {
      const Int_t off_j = j*n;

      // Find maximum in the j-th column

      Double_t max = TMath::Abs(pLU[off_j+j]);
      Int_t i_pivot = j;

      for (Int_t i = j+1; i < n; i++) {
         const Int_t off_i = i*n;
         const Double_t mLUij = TMath::Abs(pLU[off_i+j]);

         if (mLUij > max) {
            max = mLUij;
            i_pivot = i;
         }
      }
  
      if (i_pivot != j) {
         const Int_t off_ipov = i_pivot*n;
         for (Int_t k = 0; k < n; k++ ) {
            const Double_t tmp = pLU[off_ipov+k];
            pLU[off_ipov+k] = pLU[off_j+k];
            pLU[off_j+k]    = tmp;
         }
         sign = -sign;
      }
      index[j] = i_pivot;

      const Double_t mLUjj = pLU[off_j+j];

      if (mLUjj != 0.0) {
         if (TMath::Abs(mLUjj) < tol)
            nrZeros++;
         for (Int_t i = j+1; i < n; i++) {
            const Int_t off_i = i*n;
            const Double_t mLUij = pLU[off_i+j]/mLUjj;
            pLU[off_i+j] = mLUij;

            for (Int_t k = j+1; k < n; k++) {
               const Double_t mLUik = pLU[off_i+k];
               const Double_t mLUjk = pLU[off_j+k];
               pLU[off_i+k] = mLUik-mLUij*mLUjk;
            }
         }
      } else {
         ::Error("TDecompLU::DecomposeLUGauss","matrix is singular");
         return kFALSE;
      }
   }

   return kTRUE;
}

//______________________________________________________________________________
Bool_t TDecompLU::InvertLU(TMatrixD &lu,Double_t tol,Double_t *det)
{
   // Calculate matrix inversion through in place forward/backward substitution

   if (det)
      *det = 0.0;

   if (lu.GetNrows() != lu.GetNcols() || lu.GetRowLwb() != lu.GetColLwb()) {
      ::Error("TDecompLU::InvertLU","matrix should be square");
      return kFALSE;
   }

   const Int_t     n   = lu.GetNcols();
   Double_t *pLU = lu.GetMatrixArray();

   Int_t worki[kWorkMax];
   Bool_t isAllocatedI = kFALSE;
   Int_t *index = worki;
   if (n > kWorkMax) {
      isAllocatedI = kTRUE;
      index = new Int_t[n];
   }

   Double_t sign = 1.0;
   Int_t nrZeros = 0;
   if (!DecomposeLUCrout(lu,index,sign,tol,nrZeros) || nrZeros > 0) {
      if (isAllocatedI)
         delete [] index;
      ::Error("TDecompLU::InvertLU","matrix is singular, %d diag elements < tolerance of %.4e",nrZeros,tol);
      return kFALSE;
   }

   if (det) {
      Double_t d1;
      Double_t d2;
      const TVectorD diagv = TMatrixDDiag_const(lu);
      DiagProd(diagv,tol,d1,d2);
      d1 *= sign;
      *det = d1*TMath::Power(2.0,d2);
   }

   //  Form inv(U).

   Int_t j;

   for (j = 0; j < n; j++) {
      const Int_t off_j = j*n;

      pLU[off_j+j] = 1./pLU[off_j+j];
      const Double_t mLU_jj = -pLU[off_j+j];

//    Compute elements 0:j-1 of j-th column.

      Double_t *pX = pLU+j;
      Int_t k;
      for (k = 0; k <= j-1; k++) {
         const Int_t off_k = k*n;
         if ( pX[off_k] != 0.0 ) {
            const Double_t tmp = pX[off_k];
            for (Int_t i = 0; i <= k-1; i++) {
               const Int_t off_i = i*n;
               pX[off_i] += tmp*pLU[off_i+k];
            }
            pX[off_k] *= pLU[off_k+k];
         }
      }
      for (k = 0; k <= j-1; k++) {
         const Int_t off_k = k*n;
         pX[off_k] *= mLU_jj;
      }
   }

   // Solve the equation inv(A)*L = inv(U) for inv(A).

   Double_t workd[kWorkMax];
   Bool_t isAllocatedD = kFALSE;
   Double_t *pWorkd = workd;
   if (n > kWorkMax) {
      isAllocatedD = kTRUE;
      pWorkd = new Double_t[n];
   }

   for (j = n-1; j >= 0; j--) {

      // Copy current column j of L to WORK and replace with zeros.
      for (Int_t i = j+1; i < n; i++) {
         const Int_t off_i = i*n;
         pWorkd[i] = pLU[off_i+j];
         pLU[off_i+j] = 0.0;
      }

      // Compute current column of inv(A).

      if (j < n-1) {
         const Double_t *mp = pLU+j+1;  // Matrix row ptr
         Double_t *tp = pLU+j;          // Target vector ptr

         for (Int_t irow = 0; irow < n; irow++) {
            Double_t sum = 0.;
            const Double_t *sp = pWorkd+j+1; // Source vector ptr
            for (Int_t icol = 0; icol < n-1-j ; icol++)
               sum += *mp++ * *sp++;
            *tp = -sum + *tp;
            mp += j+1;
            tp += n;
         }
      }
   }

   if (isAllocatedD)
      delete [] pWorkd;

   // Apply column interchanges.
   for (j = n-1; j >= 0; j--) {
      const Int_t jperm = index[j];
      if (jperm != j) {
         for (Int_t i = 0; i < n; i++) {
            const Int_t off_i = i*n;
            const Double_t tmp = pLU[off_i+jperm];
            pLU[off_i+jperm] = pLU[off_i+j];
            pLU[off_i+j]     = tmp;
         }
      }
   }

   if (isAllocatedI)
      delete [] index;

   return kTRUE;
}
 TDecompLU.cxx:1
 TDecompLU.cxx:2
 TDecompLU.cxx:3
 TDecompLU.cxx:4
 TDecompLU.cxx:5
 TDecompLU.cxx:6
 TDecompLU.cxx:7
 TDecompLU.cxx:8
 TDecompLU.cxx:9
 TDecompLU.cxx:10
 TDecompLU.cxx:11
 TDecompLU.cxx:12
 TDecompLU.cxx:13
 TDecompLU.cxx:14
 TDecompLU.cxx:15
 TDecompLU.cxx:16
 TDecompLU.cxx:17
 TDecompLU.cxx:18
 TDecompLU.cxx:19
 TDecompLU.cxx:20
 TDecompLU.cxx:21
 TDecompLU.cxx:22
 TDecompLU.cxx:23
 TDecompLU.cxx:24
 TDecompLU.cxx:25
 TDecompLU.cxx:26
 TDecompLU.cxx:27
 TDecompLU.cxx:28
 TDecompLU.cxx:29
 TDecompLU.cxx:30
 TDecompLU.cxx:31
 TDecompLU.cxx:32
 TDecompLU.cxx:33
 TDecompLU.cxx:34
 TDecompLU.cxx:35
 TDecompLU.cxx:36
 TDecompLU.cxx:37
 TDecompLU.cxx:38
 TDecompLU.cxx:39
 TDecompLU.cxx:40
 TDecompLU.cxx:41
 TDecompLU.cxx:42
 TDecompLU.cxx:43
 TDecompLU.cxx:44
 TDecompLU.cxx:45
 TDecompLU.cxx:46
 TDecompLU.cxx:47
 TDecompLU.cxx:48
 TDecompLU.cxx:49
 TDecompLU.cxx:50
 TDecompLU.cxx:51
 TDecompLU.cxx:52
 TDecompLU.cxx:53
 TDecompLU.cxx:54
 TDecompLU.cxx:55
 TDecompLU.cxx:56
 TDecompLU.cxx:57
 TDecompLU.cxx:58
 TDecompLU.cxx:59
 TDecompLU.cxx:60
 TDecompLU.cxx:61
 TDecompLU.cxx:62
 TDecompLU.cxx:63
 TDecompLU.cxx:64
 TDecompLU.cxx:65
 TDecompLU.cxx:66
 TDecompLU.cxx:67
 TDecompLU.cxx:68
 TDecompLU.cxx:69
 TDecompLU.cxx:70
 TDecompLU.cxx:71
 TDecompLU.cxx:72
 TDecompLU.cxx:73
 TDecompLU.cxx:74
 TDecompLU.cxx:75
 TDecompLU.cxx:76
 TDecompLU.cxx:77
 TDecompLU.cxx:78
 TDecompLU.cxx:79
 TDecompLU.cxx:80
 TDecompLU.cxx:81
 TDecompLU.cxx:82
 TDecompLU.cxx:83
 TDecompLU.cxx:84
 TDecompLU.cxx:85
 TDecompLU.cxx:86
 TDecompLU.cxx:87
 TDecompLU.cxx:88
 TDecompLU.cxx:89
 TDecompLU.cxx:90
 TDecompLU.cxx:91
 TDecompLU.cxx:92
 TDecompLU.cxx:93
 TDecompLU.cxx:94
 TDecompLU.cxx:95
 TDecompLU.cxx:96
 TDecompLU.cxx:97
 TDecompLU.cxx:98
 TDecompLU.cxx:99
 TDecompLU.cxx:100
 TDecompLU.cxx:101
 TDecompLU.cxx:102
 TDecompLU.cxx:103
 TDecompLU.cxx:104
 TDecompLU.cxx:105
 TDecompLU.cxx:106
 TDecompLU.cxx:107
 TDecompLU.cxx:108
 TDecompLU.cxx:109
 TDecompLU.cxx:110
 TDecompLU.cxx:111
 TDecompLU.cxx:112
 TDecompLU.cxx:113
 TDecompLU.cxx:114
 TDecompLU.cxx:115
 TDecompLU.cxx:116
 TDecompLU.cxx:117
 TDecompLU.cxx:118
 TDecompLU.cxx:119
 TDecompLU.cxx:120
 TDecompLU.cxx:121
 TDecompLU.cxx:122
 TDecompLU.cxx:123
 TDecompLU.cxx:124
 TDecompLU.cxx:125
 TDecompLU.cxx:126
 TDecompLU.cxx:127
 TDecompLU.cxx:128
 TDecompLU.cxx:129
 TDecompLU.cxx:130
 TDecompLU.cxx:131
 TDecompLU.cxx:132
 TDecompLU.cxx:133
 TDecompLU.cxx:134
 TDecompLU.cxx:135
 TDecompLU.cxx:136
 TDecompLU.cxx:137
 TDecompLU.cxx:138
 TDecompLU.cxx:139
 TDecompLU.cxx:140
 TDecompLU.cxx:141
 TDecompLU.cxx:142
 TDecompLU.cxx:143
 TDecompLU.cxx:144
 TDecompLU.cxx:145
 TDecompLU.cxx:146
 TDecompLU.cxx:147
 TDecompLU.cxx:148
 TDecompLU.cxx:149
 TDecompLU.cxx:150
 TDecompLU.cxx:151
 TDecompLU.cxx:152
 TDecompLU.cxx:153
 TDecompLU.cxx:154
 TDecompLU.cxx:155
 TDecompLU.cxx:156
 TDecompLU.cxx:157
 TDecompLU.cxx:158
 TDecompLU.cxx:159
 TDecompLU.cxx:160
 TDecompLU.cxx:161
 TDecompLU.cxx:162
 TDecompLU.cxx:163
 TDecompLU.cxx:164
 TDecompLU.cxx:165
 TDecompLU.cxx:166
 TDecompLU.cxx:167
 TDecompLU.cxx:168
 TDecompLU.cxx:169
 TDecompLU.cxx:170
 TDecompLU.cxx:171
 TDecompLU.cxx:172
 TDecompLU.cxx:173
 TDecompLU.cxx:174
 TDecompLU.cxx:175
 TDecompLU.cxx:176
 TDecompLU.cxx:177
 TDecompLU.cxx:178
 TDecompLU.cxx:179
 TDecompLU.cxx:180
 TDecompLU.cxx:181
 TDecompLU.cxx:182
 TDecompLU.cxx:183
 TDecompLU.cxx:184
 TDecompLU.cxx:185
 TDecompLU.cxx:186
 TDecompLU.cxx:187
 TDecompLU.cxx:188
 TDecompLU.cxx:189
 TDecompLU.cxx:190
 TDecompLU.cxx:191
 TDecompLU.cxx:192
 TDecompLU.cxx:193
 TDecompLU.cxx:194
 TDecompLU.cxx:195
 TDecompLU.cxx:196
 TDecompLU.cxx:197
 TDecompLU.cxx:198
 TDecompLU.cxx:199
 TDecompLU.cxx:200
 TDecompLU.cxx:201
 TDecompLU.cxx:202
 TDecompLU.cxx:203
 TDecompLU.cxx:204
 TDecompLU.cxx:205
 TDecompLU.cxx:206
 TDecompLU.cxx:207
 TDecompLU.cxx:208
 TDecompLU.cxx:209
 TDecompLU.cxx:210
 TDecompLU.cxx:211
 TDecompLU.cxx:212
 TDecompLU.cxx:213
 TDecompLU.cxx:214
 TDecompLU.cxx:215
 TDecompLU.cxx:216
 TDecompLU.cxx:217
 TDecompLU.cxx:218
 TDecompLU.cxx:219
 TDecompLU.cxx:220
 TDecompLU.cxx:221
 TDecompLU.cxx:222
 TDecompLU.cxx:223
 TDecompLU.cxx:224
 TDecompLU.cxx:225
 TDecompLU.cxx:226
 TDecompLU.cxx:227
 TDecompLU.cxx:228
 TDecompLU.cxx:229
 TDecompLU.cxx:230
 TDecompLU.cxx:231
 TDecompLU.cxx:232
 TDecompLU.cxx:233
 TDecompLU.cxx:234
 TDecompLU.cxx:235
 TDecompLU.cxx:236
 TDecompLU.cxx:237
 TDecompLU.cxx:238
 TDecompLU.cxx:239
 TDecompLU.cxx:240
 TDecompLU.cxx:241
 TDecompLU.cxx:242
 TDecompLU.cxx:243
 TDecompLU.cxx:244
 TDecompLU.cxx:245
 TDecompLU.cxx:246
 TDecompLU.cxx:247
 TDecompLU.cxx:248
 TDecompLU.cxx:249
 TDecompLU.cxx:250
 TDecompLU.cxx:251
 TDecompLU.cxx:252
 TDecompLU.cxx:253
 TDecompLU.cxx:254
 TDecompLU.cxx:255
 TDecompLU.cxx:256
 TDecompLU.cxx:257
 TDecompLU.cxx:258
 TDecompLU.cxx:259
 TDecompLU.cxx:260
 TDecompLU.cxx:261
 TDecompLU.cxx:262
 TDecompLU.cxx:263
 TDecompLU.cxx:264
 TDecompLU.cxx:265
 TDecompLU.cxx:266
 TDecompLU.cxx:267
 TDecompLU.cxx:268
 TDecompLU.cxx:269
 TDecompLU.cxx:270
 TDecompLU.cxx:271
 TDecompLU.cxx:272
 TDecompLU.cxx:273
 TDecompLU.cxx:274
 TDecompLU.cxx:275
 TDecompLU.cxx:276
 TDecompLU.cxx:277
 TDecompLU.cxx:278
 TDecompLU.cxx:279
 TDecompLU.cxx:280
 TDecompLU.cxx:281
 TDecompLU.cxx:282
 TDecompLU.cxx:283
 TDecompLU.cxx:284
 TDecompLU.cxx:285
 TDecompLU.cxx:286
 TDecompLU.cxx:287
 TDecompLU.cxx:288
 TDecompLU.cxx:289
 TDecompLU.cxx:290
 TDecompLU.cxx:291
 TDecompLU.cxx:292
 TDecompLU.cxx:293
 TDecompLU.cxx:294
 TDecompLU.cxx:295
 TDecompLU.cxx:296
 TDecompLU.cxx:297
 TDecompLU.cxx:298
 TDecompLU.cxx:299
 TDecompLU.cxx:300
 TDecompLU.cxx:301
 TDecompLU.cxx:302
 TDecompLU.cxx:303
 TDecompLU.cxx:304
 TDecompLU.cxx:305
 TDecompLU.cxx:306
 TDecompLU.cxx:307
 TDecompLU.cxx:308
 TDecompLU.cxx:309
 TDecompLU.cxx:310
 TDecompLU.cxx:311
 TDecompLU.cxx:312
 TDecompLU.cxx:313
 TDecompLU.cxx:314
 TDecompLU.cxx:315
 TDecompLU.cxx:316
 TDecompLU.cxx:317
 TDecompLU.cxx:318
 TDecompLU.cxx:319
 TDecompLU.cxx:320
 TDecompLU.cxx:321
 TDecompLU.cxx:322
 TDecompLU.cxx:323
 TDecompLU.cxx:324
 TDecompLU.cxx:325
 TDecompLU.cxx:326
 TDecompLU.cxx:327
 TDecompLU.cxx:328
 TDecompLU.cxx:329
 TDecompLU.cxx:330
 TDecompLU.cxx:331
 TDecompLU.cxx:332
 TDecompLU.cxx:333
 TDecompLU.cxx:334
 TDecompLU.cxx:335
 TDecompLU.cxx:336
 TDecompLU.cxx:337
 TDecompLU.cxx:338
 TDecompLU.cxx:339
 TDecompLU.cxx:340
 TDecompLU.cxx:341
 TDecompLU.cxx:342
 TDecompLU.cxx:343
 TDecompLU.cxx:344
 TDecompLU.cxx:345
 TDecompLU.cxx:346
 TDecompLU.cxx:347
 TDecompLU.cxx:348
 TDecompLU.cxx:349
 TDecompLU.cxx:350
 TDecompLU.cxx:351
 TDecompLU.cxx:352
 TDecompLU.cxx:353
 TDecompLU.cxx:354
 TDecompLU.cxx:355
 TDecompLU.cxx:356
 TDecompLU.cxx:357
 TDecompLU.cxx:358
 TDecompLU.cxx:359
 TDecompLU.cxx:360
 TDecompLU.cxx:361
 TDecompLU.cxx:362
 TDecompLU.cxx:363
 TDecompLU.cxx:364
 TDecompLU.cxx:365
 TDecompLU.cxx:366
 TDecompLU.cxx:367
 TDecompLU.cxx:368
 TDecompLU.cxx:369
 TDecompLU.cxx:370
 TDecompLU.cxx:371
 TDecompLU.cxx:372
 TDecompLU.cxx:373
 TDecompLU.cxx:374
 TDecompLU.cxx:375
 TDecompLU.cxx:376
 TDecompLU.cxx:377
 TDecompLU.cxx:378
 TDecompLU.cxx:379
 TDecompLU.cxx:380
 TDecompLU.cxx:381
 TDecompLU.cxx:382
 TDecompLU.cxx:383
 TDecompLU.cxx:384
 TDecompLU.cxx:385
 TDecompLU.cxx:386
 TDecompLU.cxx:387
 TDecompLU.cxx:388
 TDecompLU.cxx:389
 TDecompLU.cxx:390
 TDecompLU.cxx:391
 TDecompLU.cxx:392
 TDecompLU.cxx:393
 TDecompLU.cxx:394
 TDecompLU.cxx:395
 TDecompLU.cxx:396
 TDecompLU.cxx:397
 TDecompLU.cxx:398
 TDecompLU.cxx:399
 TDecompLU.cxx:400
 TDecompLU.cxx:401
 TDecompLU.cxx:402
 TDecompLU.cxx:403
 TDecompLU.cxx:404
 TDecompLU.cxx:405
 TDecompLU.cxx:406
 TDecompLU.cxx:407
 TDecompLU.cxx:408
 TDecompLU.cxx:409
 TDecompLU.cxx:410
 TDecompLU.cxx:411
 TDecompLU.cxx:412
 TDecompLU.cxx:413
 TDecompLU.cxx:414
 TDecompLU.cxx:415
 TDecompLU.cxx:416
 TDecompLU.cxx:417
 TDecompLU.cxx:418
 TDecompLU.cxx:419
 TDecompLU.cxx:420
 TDecompLU.cxx:421
 TDecompLU.cxx:422
 TDecompLU.cxx:423
 TDecompLU.cxx:424
 TDecompLU.cxx:425
 TDecompLU.cxx:426
 TDecompLU.cxx:427
 TDecompLU.cxx:428
 TDecompLU.cxx:429
 TDecompLU.cxx:430
 TDecompLU.cxx:431
 TDecompLU.cxx:432
 TDecompLU.cxx:433
 TDecompLU.cxx:434
 TDecompLU.cxx:435
 TDecompLU.cxx:436
 TDecompLU.cxx:437
 TDecompLU.cxx:438
 TDecompLU.cxx:439
 TDecompLU.cxx:440
 TDecompLU.cxx:441
 TDecompLU.cxx:442
 TDecompLU.cxx:443
 TDecompLU.cxx:444
 TDecompLU.cxx:445
 TDecompLU.cxx:446
 TDecompLU.cxx:447
 TDecompLU.cxx:448
 TDecompLU.cxx:449
 TDecompLU.cxx:450
 TDecompLU.cxx:451
 TDecompLU.cxx:452
 TDecompLU.cxx:453
 TDecompLU.cxx:454
 TDecompLU.cxx:455
 TDecompLU.cxx:456
 TDecompLU.cxx:457
 TDecompLU.cxx:458
 TDecompLU.cxx:459
 TDecompLU.cxx:460
 TDecompLU.cxx:461
 TDecompLU.cxx:462
 TDecompLU.cxx:463
 TDecompLU.cxx:464
 TDecompLU.cxx:465
 TDecompLU.cxx:466
 TDecompLU.cxx:467
 TDecompLU.cxx:468
 TDecompLU.cxx:469
 TDecompLU.cxx:470
 TDecompLU.cxx:471
 TDecompLU.cxx:472
 TDecompLU.cxx:473
 TDecompLU.cxx:474
 TDecompLU.cxx:475
 TDecompLU.cxx:476
 TDecompLU.cxx:477
 TDecompLU.cxx:478
 TDecompLU.cxx:479
 TDecompLU.cxx:480
 TDecompLU.cxx:481
 TDecompLU.cxx:482
 TDecompLU.cxx:483
 TDecompLU.cxx:484
 TDecompLU.cxx:485
 TDecompLU.cxx:486
 TDecompLU.cxx:487
 TDecompLU.cxx:488
 TDecompLU.cxx:489
 TDecompLU.cxx:490
 TDecompLU.cxx:491
 TDecompLU.cxx:492
 TDecompLU.cxx:493
 TDecompLU.cxx:494
 TDecompLU.cxx:495
 TDecompLU.cxx:496
 TDecompLU.cxx:497
 TDecompLU.cxx:498
 TDecompLU.cxx:499
 TDecompLU.cxx:500
 TDecompLU.cxx:501
 TDecompLU.cxx:502
 TDecompLU.cxx:503
 TDecompLU.cxx:504
 TDecompLU.cxx:505
 TDecompLU.cxx:506
 TDecompLU.cxx:507
 TDecompLU.cxx:508
 TDecompLU.cxx:509
 TDecompLU.cxx:510
 TDecompLU.cxx:511
 TDecompLU.cxx:512
 TDecompLU.cxx:513
 TDecompLU.cxx:514
 TDecompLU.cxx:515
 TDecompLU.cxx:516
 TDecompLU.cxx:517
 TDecompLU.cxx:518
 TDecompLU.cxx:519
 TDecompLU.cxx:520
 TDecompLU.cxx:521
 TDecompLU.cxx:522
 TDecompLU.cxx:523
 TDecompLU.cxx:524
 TDecompLU.cxx:525
 TDecompLU.cxx:526
 TDecompLU.cxx:527
 TDecompLU.cxx:528
 TDecompLU.cxx:529
 TDecompLU.cxx:530
 TDecompLU.cxx:531
 TDecompLU.cxx:532
 TDecompLU.cxx:533
 TDecompLU.cxx:534
 TDecompLU.cxx:535
 TDecompLU.cxx:536
 TDecompLU.cxx:537
 TDecompLU.cxx:538
 TDecompLU.cxx:539
 TDecompLU.cxx:540
 TDecompLU.cxx:541
 TDecompLU.cxx:542
 TDecompLU.cxx:543
 TDecompLU.cxx:544
 TDecompLU.cxx:545
 TDecompLU.cxx:546
 TDecompLU.cxx:547
 TDecompLU.cxx:548
 TDecompLU.cxx:549
 TDecompLU.cxx:550
 TDecompLU.cxx:551
 TDecompLU.cxx:552
 TDecompLU.cxx:553
 TDecompLU.cxx:554
 TDecompLU.cxx:555
 TDecompLU.cxx:556
 TDecompLU.cxx:557
 TDecompLU.cxx:558
 TDecompLU.cxx:559
 TDecompLU.cxx:560
 TDecompLU.cxx:561
 TDecompLU.cxx:562
 TDecompLU.cxx:563
 TDecompLU.cxx:564
 TDecompLU.cxx:565
 TDecompLU.cxx:566
 TDecompLU.cxx:567
 TDecompLU.cxx:568
 TDecompLU.cxx:569
 TDecompLU.cxx:570
 TDecompLU.cxx:571
 TDecompLU.cxx:572
 TDecompLU.cxx:573
 TDecompLU.cxx:574
 TDecompLU.cxx:575
 TDecompLU.cxx:576
 TDecompLU.cxx:577
 TDecompLU.cxx:578
 TDecompLU.cxx:579
 TDecompLU.cxx:580
 TDecompLU.cxx:581
 TDecompLU.cxx:582
 TDecompLU.cxx:583
 TDecompLU.cxx:584
 TDecompLU.cxx:585
 TDecompLU.cxx:586
 TDecompLU.cxx:587
 TDecompLU.cxx:588
 TDecompLU.cxx:589
 TDecompLU.cxx:590
 TDecompLU.cxx:591
 TDecompLU.cxx:592
 TDecompLU.cxx:593
 TDecompLU.cxx:594
 TDecompLU.cxx:595
 TDecompLU.cxx:596
 TDecompLU.cxx:597
 TDecompLU.cxx:598
 TDecompLU.cxx:599
 TDecompLU.cxx:600
 TDecompLU.cxx:601
 TDecompLU.cxx:602
 TDecompLU.cxx:603
 TDecompLU.cxx:604
 TDecompLU.cxx:605
 TDecompLU.cxx:606
 TDecompLU.cxx:607
 TDecompLU.cxx:608
 TDecompLU.cxx:609
 TDecompLU.cxx:610
 TDecompLU.cxx:611
 TDecompLU.cxx:612
 TDecompLU.cxx:613
 TDecompLU.cxx:614
 TDecompLU.cxx:615
 TDecompLU.cxx:616
 TDecompLU.cxx:617
 TDecompLU.cxx:618
 TDecompLU.cxx:619
 TDecompLU.cxx:620
 TDecompLU.cxx:621
 TDecompLU.cxx:622
 TDecompLU.cxx:623
 TDecompLU.cxx:624
 TDecompLU.cxx:625
 TDecompLU.cxx:626
 TDecompLU.cxx:627
 TDecompLU.cxx:628
 TDecompLU.cxx:629
 TDecompLU.cxx:630
 TDecompLU.cxx:631
 TDecompLU.cxx:632
 TDecompLU.cxx:633
 TDecompLU.cxx:634
 TDecompLU.cxx:635
 TDecompLU.cxx:636
 TDecompLU.cxx:637
 TDecompLU.cxx:638
 TDecompLU.cxx:639
 TDecompLU.cxx:640
 TDecompLU.cxx:641
 TDecompLU.cxx:642
 TDecompLU.cxx:643
 TDecompLU.cxx:644
 TDecompLU.cxx:645
 TDecompLU.cxx:646
 TDecompLU.cxx:647
 TDecompLU.cxx:648
 TDecompLU.cxx:649
 TDecompLU.cxx:650
 TDecompLU.cxx:651
 TDecompLU.cxx:652
 TDecompLU.cxx:653
 TDecompLU.cxx:654
 TDecompLU.cxx:655
 TDecompLU.cxx:656
 TDecompLU.cxx:657
 TDecompLU.cxx:658
 TDecompLU.cxx:659
 TDecompLU.cxx:660
 TDecompLU.cxx:661
 TDecompLU.cxx:662
 TDecompLU.cxx:663
 TDecompLU.cxx:664
 TDecompLU.cxx:665
 TDecompLU.cxx:666
 TDecompLU.cxx:667
 TDecompLU.cxx:668
 TDecompLU.cxx:669
 TDecompLU.cxx:670
 TDecompLU.cxx:671
 TDecompLU.cxx:672
 TDecompLU.cxx:673
 TDecompLU.cxx:674
 TDecompLU.cxx:675
 TDecompLU.cxx:676
 TDecompLU.cxx:677
 TDecompLU.cxx:678
 TDecompLU.cxx:679
 TDecompLU.cxx:680
 TDecompLU.cxx:681
 TDecompLU.cxx:682
 TDecompLU.cxx:683
 TDecompLU.cxx:684
 TDecompLU.cxx:685
 TDecompLU.cxx:686
 TDecompLU.cxx:687
 TDecompLU.cxx:688
 TDecompLU.cxx:689
 TDecompLU.cxx:690
 TDecompLU.cxx:691
 TDecompLU.cxx:692
 TDecompLU.cxx:693
 TDecompLU.cxx:694
 TDecompLU.cxx:695
 TDecompLU.cxx:696
 TDecompLU.cxx:697
 TDecompLU.cxx:698
 TDecompLU.cxx:699
 TDecompLU.cxx:700
 TDecompLU.cxx:701
 TDecompLU.cxx:702
 TDecompLU.cxx:703
 TDecompLU.cxx:704
 TDecompLU.cxx:705
 TDecompLU.cxx:706
 TDecompLU.cxx:707
 TDecompLU.cxx:708
 TDecompLU.cxx:709
 TDecompLU.cxx:710
 TDecompLU.cxx:711
 TDecompLU.cxx:712
 TDecompLU.cxx:713
 TDecompLU.cxx:714
 TDecompLU.cxx:715
 TDecompLU.cxx:716
 TDecompLU.cxx:717
 TDecompLU.cxx:718
 TDecompLU.cxx:719
 TDecompLU.cxx:720
 TDecompLU.cxx:721
 TDecompLU.cxx:722
 TDecompLU.cxx:723
 TDecompLU.cxx:724
 TDecompLU.cxx:725
 TDecompLU.cxx:726
 TDecompLU.cxx:727
 TDecompLU.cxx:728
 TDecompLU.cxx:729
 TDecompLU.cxx:730
 TDecompLU.cxx:731
 TDecompLU.cxx:732
 TDecompLU.cxx:733
 TDecompLU.cxx:734
 TDecompLU.cxx:735
 TDecompLU.cxx:736
 TDecompLU.cxx:737
 TDecompLU.cxx:738
 TDecompLU.cxx:739
 TDecompLU.cxx:740
 TDecompLU.cxx:741
 TDecompLU.cxx:742
 TDecompLU.cxx:743
 TDecompLU.cxx:744
 TDecompLU.cxx:745
 TDecompLU.cxx:746
 TDecompLU.cxx:747
 TDecompLU.cxx:748
 TDecompLU.cxx:749
 TDecompLU.cxx:750
 TDecompLU.cxx:751
 TDecompLU.cxx:752
 TDecompLU.cxx:753
 TDecompLU.cxx:754
 TDecompLU.cxx:755
 TDecompLU.cxx:756
 TDecompLU.cxx:757
 TDecompLU.cxx:758
 TDecompLU.cxx:759
 TDecompLU.cxx:760
 TDecompLU.cxx:761
 TDecompLU.cxx:762
 TDecompLU.cxx:763
 TDecompLU.cxx:764
 TDecompLU.cxx:765
 TDecompLU.cxx:766
 TDecompLU.cxx:767
 TDecompLU.cxx:768
 TDecompLU.cxx:769
 TDecompLU.cxx:770
 TDecompLU.cxx:771
 TDecompLU.cxx:772
 TDecompLU.cxx:773
 TDecompLU.cxx:774
 TDecompLU.cxx:775
 TDecompLU.cxx:776
 TDecompLU.cxx:777
 TDecompLU.cxx:778
 TDecompLU.cxx:779
 TDecompLU.cxx:780
 TDecompLU.cxx:781
 TDecompLU.cxx:782
 TDecompLU.cxx:783
 TDecompLU.cxx:784
 TDecompLU.cxx:785
 TDecompLU.cxx:786
 TDecompLU.cxx:787
 TDecompLU.cxx:788
 TDecompLU.cxx:789
 TDecompLU.cxx:790
 TDecompLU.cxx:791
 TDecompLU.cxx:792
 TDecompLU.cxx:793
 TDecompLU.cxx:794
 TDecompLU.cxx:795
 TDecompLU.cxx:796
 TDecompLU.cxx:797
 TDecompLU.cxx:798
 TDecompLU.cxx:799
 TDecompLU.cxx:800
 TDecompLU.cxx:801
 TDecompLU.cxx:802
 TDecompLU.cxx:803
 TDecompLU.cxx:804
 TDecompLU.cxx:805
 TDecompLU.cxx:806
 TDecompLU.cxx:807
 TDecompLU.cxx:808
 TDecompLU.cxx:809
 TDecompLU.cxx:810
 TDecompLU.cxx:811
 TDecompLU.cxx:812
 TDecompLU.cxx:813
 TDecompLU.cxx:814
 TDecompLU.cxx:815
 TDecompLU.cxx:816
 TDecompLU.cxx:817
 TDecompLU.cxx:818
 TDecompLU.cxx:819
 TDecompLU.cxx:820
 TDecompLU.cxx:821
 TDecompLU.cxx:822
 TDecompLU.cxx:823
 TDecompLU.cxx:824
 TDecompLU.cxx:825
 TDecompLU.cxx:826
 TDecompLU.cxx:827
 TDecompLU.cxx:828
 TDecompLU.cxx:829
 TDecompLU.cxx:830
 TDecompLU.cxx:831
 TDecompLU.cxx:832
 TDecompLU.cxx:833
 TDecompLU.cxx:834
 TDecompLU.cxx:835
 TDecompLU.cxx:836
 TDecompLU.cxx:837
 TDecompLU.cxx:838
 TDecompLU.cxx:839
 TDecompLU.cxx:840
 TDecompLU.cxx:841
 TDecompLU.cxx:842
 TDecompLU.cxx:843
 TDecompLU.cxx:844
 TDecompLU.cxx:845
 TDecompLU.cxx:846
 TDecompLU.cxx:847
 TDecompLU.cxx:848
 TDecompLU.cxx:849
 TDecompLU.cxx:850
 TDecompLU.cxx:851
 TDecompLU.cxx:852
 TDecompLU.cxx:853
 TDecompLU.cxx:854
 TDecompLU.cxx:855
 TDecompLU.cxx:856
 TDecompLU.cxx:857
 TDecompLU.cxx:858
 TDecompLU.cxx:859
 TDecompLU.cxx:860
 TDecompLU.cxx:861
 TDecompLU.cxx:862
 TDecompLU.cxx:863
 TDecompLU.cxx:864
 TDecompLU.cxx:865
 TDecompLU.cxx:866
 TDecompLU.cxx:867
 TDecompLU.cxx:868
 TDecompLU.cxx:869
 TDecompLU.cxx:870
 TDecompLU.cxx:871
 TDecompLU.cxx:872
 TDecompLU.cxx:873
 TDecompLU.cxx:874
 TDecompLU.cxx:875
 TDecompLU.cxx:876
 TDecompLU.cxx:877
 TDecompLU.cxx:878
 TDecompLU.cxx:879
 TDecompLU.cxx:880
 TDecompLU.cxx:881
 TDecompLU.cxx:882
 TDecompLU.cxx:883
 TDecompLU.cxx:884
 TDecompLU.cxx:885
 TDecompLU.cxx:886
 TDecompLU.cxx:887
 TDecompLU.cxx:888
 TDecompLU.cxx:889
 TDecompLU.cxx:890
 TDecompLU.cxx:891
 TDecompLU.cxx:892
 TDecompLU.cxx:893
 TDecompLU.cxx:894
 TDecompLU.cxx:895
 TDecompLU.cxx:896
 TDecompLU.cxx:897
 TDecompLU.cxx:898
 TDecompLU.cxx:899