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

/*************************************************************************
 * 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 "TDecompBK.h"
#include "TMath.h"
#include "TError.h"

ClassImp(TDecompBK)

///////////////////////////////////////////////////////////////////////////
//                                                                       //
// The Bunch-Kaufman diagonal pivoting method decomposes a real          //
// symmetric matrix A using                                              //
//                                                                       //
//     A = U*D*U^T                                                       //
//                                                                       //
//  where U is a product of permutation and unit upper triangular        //
//  matrices, U^T is the transpose of U, and D is symmetric and block    //
//  diagonal with 1-by-1 and 2-by-2 diagonal blocks.                     //
//                                                                       //
//     U = P(n-1)*U(n-1)* ... *P(k)U(k)* ...,                            //
//  i.e., U is a product of terms P(k)*U(k), where k decreases from n-1  //
//  to 0 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1//
//  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as    //
//  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such //
//  that if the diagonal block D(k) is of order s (s = 1 or 2), then     //
//                                                                       //
//             (   I    v    0   )   k-s                                 //
//     U(k) =  (   0    I    0   )   s                                   //
//             (   0    0    I   )   n-k                                 //
//                k-s   s   n-k                                          //
//                                                                       //
//  If s = 1, D(k) overwrites A(k,k), and v overwrites A(0:k-1,k).       //
//  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),//
//  and A(k,k), and v overwrites A(0:k-2,k-1:k).                         //
//                                                                       //
// fU contains on entry the symmetric matrix A of which only the upper   //
// triangular part is referenced . On exit fU contains the block diagonal//
// matrix D and the multipliers used to obtain the factor U, see above . //
//                                                                       //
// fIpiv if dimension n contains details of the interchanges and the     //
// the block structure of D . If (fIPiv(k) > 0, then rows and columns k  //
// and fIPiv(k) were interchanged and D(k,k) is a 1-by-1 diagonal block. //
// If IPiv(k) = fIPiv(k-1) < 0, rows and columns k-1 and -IPiv(k) were   //
// interchanged and D(k-1:k,k-1:k) is a 2-by-2 diagonal block.           //
//                                                                       //
///////////////////////////////////////////////////////////////////////////

//______________________________________________________________________________
TDecompBK::TDecompBK()
{
// Default constructor
   fNIpiv = 0;
   fIpiv  = 0;
}

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

   fNIpiv = nrows;
   fIpiv = new Int_t[fNIpiv];
   memset(fIpiv,0,fNIpiv*sizeof(Int_t));
   fU.ResizeTo(nrows,nrows);
}

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

   const Int_t nrows = row_upb-row_lwb+1;
   fNIpiv = nrows;
   fIpiv = new Int_t[fNIpiv];
   memset(fIpiv,0,fNIpiv*sizeof(Int_t));
   fColLwb = fRowLwb = row_lwb;
   fU.ResizeTo(nrows,nrows);
}

//______________________________________________________________________________
TDecompBK::TDecompBK(const TMatrixDSym &a,Double_t tol)
{
// Constructor for symmetric matrix A

   R__ASSERT(a.IsValid());

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

   fNIpiv = a.GetNcols();
   fIpiv = new Int_t[fNIpiv];
   memset(fIpiv,0,fNIpiv*sizeof(Int_t));

   const Int_t nRows = a.GetNrows();
   fColLwb = fRowLwb = a.GetRowLwb();
   fU.ResizeTo(nRows,nRows);
   memcpy(fU.GetMatrixArray(),a.GetMatrixArray(),nRows*nRows*sizeof(Double_t));
}

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

   fNIpiv = 0;
   fIpiv  = 0;
   *this = another;
}

//______________________________________________________________________________
Bool_t TDecompBK::Decompose()
{
// Matrix A is decomposed in components U and D so that A = U*D*U^T
// 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;
   }

   Bool_t ok = kTRUE;

// Initialize alpha for use in choosing pivot block size.
   const Double_t alpha = (1.+TMath::Sqrt(17.))/8.;

// Factorize a as u*d*u' using the upper triangle of a .
//  k is the main loop index, decreasing from n-1 to 0 in steps of 1 or 2

   const Int_t     n  = fU.GetNcols();
         Double_t *pU = fU.GetMatrixArray();
   TMatrixDDiag diag(fU);

   Int_t imax = 0;
   Int_t k = n-1;
   while (k >= 0) {
      Int_t kstep = 1;

      // determine rows and columns to be interchanged and whether
      // a 1-by-1 or 2-by-2 pivot block will be used

      const Double_t absakk = TMath::Abs(diag(k));

      // imax is the row-index of the largest off-diagonal element in
      // column k, and colmax is its absolute value

      Double_t colmax;
      if ( k > 0 ) {
         TVectorD vcol = TMatrixDColumn_const(fU,k);
         vcol.Abs();
         imax = TMath::LocMax(k,vcol.GetMatrixArray());
         colmax = vcol[imax];
      } else {
         colmax = 0.;
      }

      Int_t kp;
      if (TMath::Max(absakk,colmax) <= fTol) {
         // the block diagonal matrix will be singular
         kp = k;
         ok = kFALSE;
      } else {
         if (absakk >= alpha*colmax) {
           // no interchange, use 1-by-1 pivot block
            kp = k;
         } else {
            // jmax is the column-index of the largest off-diagonal
            // element in row imax, and rowmax is its absolute value
            TVectorD vrow = TMatrixDRow_const(fU,imax);
            vrow.Abs();
            Int_t jmax = imax+1+TMath::LocMax(k-imax,vrow.GetMatrixArray()+imax+1);
            Double_t rowmax = vrow[jmax];
            if (imax > 0) {
               TVectorD vcol = TMatrixDColumn_const(fU,imax);
               vcol.Abs();
               jmax = TMath::LocMax(imax,vcol.GetMatrixArray());
               rowmax = TMath::Max(rowmax,vcol[jmax]);
            }

            if (absakk >= alpha*colmax*(colmax/rowmax)) {
               // No interchange, use 1-by-1 pivot block
               kp = k;
            } else if( TMath::Abs(diag(imax)) >= alpha*rowmax) {
               // Interchange rows and columns k and imax, use 1-by-1 pivot block
               kp = imax;
            } else {
               // Interchange rows and columns k-1 and imax, use 2-by-2 pivot block
               kp = imax;
               kstep = 2;
            }
         }

         const Int_t kk = k-kstep+1;
         if (kp != kk) {
            // Interchange rows and columns kk and kp in the leading submatrix a(0:k,0:k)
            Double_t *c_kk = pU+kk;
            Double_t *c_kp = pU+kp;
            for (Int_t irow = 0; irow < kp; irow++) {
               const Double_t t = *c_kk;
               *c_kk = *c_kp;
               *c_kp = t;
               c_kk += n;
               c_kp += n;
            }

            c_kk = pU+(kp+1)*n+kk;
            Double_t *r_kp = pU+kp*n+kp+1;
            for (Int_t icol = 0; icol < kk-kp-1; icol++) {
               const Double_t t = *c_kk;
               *c_kk = *r_kp;
               *r_kp = t;
               c_kk += n;
               r_kp += 1;
            }

            Double_t t = diag(kk);
            diag(kk) = diag(kp);
            diag(kp) = t;
            if (kstep == 2) {
               t = pU[(k-1)*n+k];
               pU[(k-1)*n+k] = pU[kp*n+k];
               pU[kp*n+k]    = t;
            }
         }

         // Update the leading submatrix

         if (kstep == 1 && k > 0) {
            // 1-by-1 pivot block d(k): column k now holds w(k) = u(k)*d(k)
            // where u(k) is the k-th column of u

            // perform a rank-1 update of a(0:k-1,0:k-1) as
            // a := a - u(k)*d(k)*u(k)' = a - w(k)*1/d(k)*w(k)'

            const Double_t r1 = 1./diag(k);
            TMatrixDSub sub1(fU,0,k-1,0,k-1);
            sub1.Rank1Update(TMatrixDColumn_const(fU,k),-r1);

            // store u(k) in column k
            TMatrixDSub sub2(fU,0,k-1,k,k);
            sub2 *= r1;
         } else {
            // 2-by-2 pivot block d(k): columns k and k-1 now hold
            // ( w(k-1) w(k) ) = ( u(k-1) u(k) )*d(k)
            // where u(k) and u(k-1) are the k-th and (k-1)-th columns of u

            // perform a rank-2 update of a(0:k-2,0:k-2) as
            // a := a - ( u(k-1) u(k) )*d(k)*( u(k-1) u(k) )'
            //    = a - ( w(k-1) w(k) )*inv(d(k))*( w(k-1) w(k) )'

            if ( k > 1 ) {
                     Double_t *pU_k1 = pU+(k-1)*n;
                     Double_t d12    = pU_k1[k];
               const Double_t d22    = pU_k1[k-1]/d12;
               const Double_t d11    = diag(k)/d12;
               const Double_t t      = 1./(d11*d22-1.);
               d12 = t/d12;

               for (Int_t j = k-2; j >= 0; j--) {
                  Double_t *pU_j = pU+j*n;
                  const Double_t wkm1 = d12*(d11*pU_j[k-1]-pU_j[k]);
                  const Double_t wk   = d12*(d22*pU_j[k]-pU_j[k-1]);
                  for (Int_t i = j; i >= 0; i--) {
                     Double_t *pU_i = pU+i*n;
                     pU_i[j] -= (pU_i[k]*wk+pU_i[k-1]*wkm1);
                  }
                  pU_j[k]   = wk;
                  pU_j[k-1] = wkm1;
               }
            }
         }

         // Store details of the interchanges in fIpiv
         if (kstep == 1) {
            fIpiv[k] = (kp+1);
         } else {
            fIpiv[k]   = -(kp+1);
            fIpiv[k-1] = -(kp+1);
         }
      }

      k -= kstep;
   }

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

   fU.Shift(fRowLwb,fRowLwb);

   return ok;
}

//______________________________________________________________________________
void TDecompBK::SetMatrix(const TMatrixDSym &a)
{
// Set the matrix to be decomposed, decomposition status is reset.

   R__ASSERT(a.IsValid());

   ResetStatus();

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

   if (fNIpiv != a.GetNcols()) {
      fNIpiv = a.GetNcols();
      delete [] fIpiv;
      fIpiv = new Int_t[fNIpiv];
      memset(fIpiv,0,fNIpiv*sizeof(Int_t));
   }

   const Int_t nRows = a.GetNrows();
   fColLwb = fRowLwb = a.GetRowLwb();
   fU.ResizeTo(nRows,nRows);
   memcpy(fU.GetMatrixArray(),a.GetMatrixArray(),nRows*nRows*sizeof(Double_t));
}

//______________________________________________________________________________
Bool_t TDecompBK::Solve(TVectorD &b)
{
// Solve Ax=b assuming the BK form of A is stored in fU . 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 (fU.GetNrows() != b.GetNrows() || fU.GetRowLwb() != b.GetLwb()) {
      Error("Solve(TVectorD &","vector and matrix incompatible");
      return kFALSE;
   }

   const Int_t n = fU.GetNrows();

   TMatrixDDiag_const diag(fU);
   const Double_t *pU = fU.GetMatrixArray();
         Double_t *pb = b.GetMatrixArray();

   // solve a*x = b, where a = u*d*u'. First solve u*d*x = b, overwriting b with x.
   // k is the main loop index, decreasing from n-1 to 0 in steps of 1 or 2,
   // depending on the size of the diagonal blocks.

   Int_t k = n-1;
   while (k >= 0) {

      if (fIpiv[k] > 0) {

         //  1 x 1 diagonal block
         //  interchange rows k and ipiv(k).
         const Int_t kp = fIpiv[k]-1;
         if (kp != k) {
            const Double_t tmp = pb[k];
            pb[k]  = pb[kp];
            pb[kp] = tmp;
         }

         // multiply by inv(u(k)), where u(k) is the transformation
         // stored in column k of a.
         for (Int_t i = 0; i < k; i++)
            pb[i] -= pU[i*n+k]*pb[k];

         // multiply by the inverse of the diagonal block.
         pb[k] /= diag(k);
         k--;
      } else {

         // 2 x 2 diagonal block
         // interchange rows k-1 and -ipiv(k).
         const Int_t kp = -fIpiv[k]-1;
         if (kp != k-1) {
            const Double_t tmp = pb[k-1];
            pb[k-1]  = pb[kp];
            pb[kp]   = tmp;
         }

         // multiply by inv(u(k)), where u(k) is the transformation
         // stored in columns k-1 and k of a.
         Int_t i;
         for (i = 0; i < k-1; i++)
            pb[i] -= pU[i*n+k]*pb[k];
         for (i = 0; i < k-1; i++)
            pb[i] -= pU[i*n+k-1]*pb[k-1];

         // multiply by the inverse of the diagonal block.
         const Double_t *pU_k1 = pU+(k-1)*n;
         const Double_t ukm1k  = pU_k1[k];
         const Double_t ukm1   = pU_k1[k-1]/ukm1k;
         const Double_t uk     = diag(k)/ukm1k;
         const Double_t denom  = ukm1*uk-1.;
         const Double_t bkm1   = pb[k-1]/ukm1k;
         const Double_t bk     = pb[k]/ukm1k;
         pb[k-1] = (uk*bkm1-bk)/denom;
         pb[k]   = (ukm1*bk-bkm1)/denom;
         k -= 2;
      }
   }

   // Next solve u'*x = b, overwriting b with x.
   //
   //  k is the main loop index, increasing from 0 to n-1 in steps of
   //  1 or 2, depending on the size of the diagonal blocks.

   k = 0;
   while (k < n) {

      if (fIpiv[k] > 0) {
         // 1 x 1 diagonal block
         //  multiply by inv(u'(k)), where u(k) is the transformation
         //  stored in column k of a.
         for (Int_t i = 0; i < k; i++)
            pb[k] -= pU[i*n+k]*pb[i];

         // interchange elements k and ipiv(k).
         const Int_t kp = fIpiv[k]-1;
         if (kp != k) {
            const Double_t tmp = pb[k];
            pb[k]  = pb[kp];
            pb[kp] = tmp;
         }
         k++;
      } else {
         // 2 x 2 diagonal block
         // multiply by inv(u'(k+1)), where u(k+1) is the transformation
         // stored in columns k and k+1 of a.
         Int_t i ;
         for (i = 0; i < k; i++)
            pb[k] -= pU[i*n+k]*pb[i];
         for (i = 0; i < k; i++)
            pb[k+1] -= pU[i*n+k+1]*pb[i];

         // interchange elements k and -ipiv(k).
         const Int_t kp = -fIpiv[k]-1;
         if (kp != k) {
            const Double_t tmp = pb[k];
            pb[k]  = pb[kp];
            pb[kp] = tmp;
         }
         k += 2;
      }
   }

   return kTRUE;
}

//______________________________________________________________________________
Bool_t TDecompBK::Solve(TMatrixDColumn &cb)
{
// Solve Ax=b assuming the BK form of A is stored in fU . 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 (fU.GetNrows() != b->GetNrows() || fU.GetRowLwb() != b->GetRowLwb()) {
      Error("Solve(TMatrixDColumn &","vector and matrix incompatible");
      return kFALSE;
   }

   const Int_t n = fU.GetNrows();

   TMatrixDDiag_const diag(fU);
   const Double_t *pU  = fU.GetMatrixArray();
         Double_t *pcb = cb.GetPtr();
   const Int_t     inc = cb.GetInc();


  // solve a*x = b, where a = u*d*u'. First solve u*d*x = b, overwriting b with x.
  // k is the main loop index, decreasing from n-1 to 0 in steps of 1 or 2,
  // depending on the size of the diagonal blocks.

   Int_t k = n-1;
   while (k >= 0) {

      if (fIpiv[k] > 0) {

         //  1 x 1 diagonal block
         //  interchange rows k and ipiv(k).
         const Int_t kp = fIpiv[k]-1;
         if (kp != k) {
            const Double_t tmp = pcb[k*inc];
            pcb[k*inc]  = pcb[kp*inc];
            pcb[kp*inc] = tmp;
         }

         // multiply by inv(u(k)), where u(k) is the transformation
         // stored in column k of a.
         for (Int_t i = 0; i < k; i++)
            pcb[i*inc] -= pU[i*n+k]*pcb[k*inc];

         // multiply by the inverse of the diagonal block.
         pcb[k*inc] /= diag(k);
         k--;
      } else {

         // 2 x 2 diagonal block
         // interchange rows k-1 and -ipiv(k).
         const Int_t kp = -fIpiv[k]-1;
         if (kp != k-1) {
            const Double_t tmp = pcb[(k-1)*inc];
            pcb[(k-1)*inc] = pcb[kp*inc];
            pcb[kp*inc]    = tmp;
         }

         // multiply by inv(u(k)), where u(k) is the transformation
         // stored in columns k-1 and k of a.
         Int_t i;
         for (i = 0; i < k-1; i++)
            pcb[i*inc] -= pU[i*n+k]*pcb[k*inc];
         for (i = 0; i < k-1; i++)
            pcb[i*inc] -= pU[i*n+k-1]*pcb[(k-1)*inc];

         // multiply by the inverse of the diagonal block.
         const Double_t *pU_k1 = pU+(k-1)*n;
         const Double_t ukm1k  = pU_k1[k];
         const Double_t ukm1   = pU_k1[k-1]/ukm1k;
         const Double_t uk     = diag(k)/ukm1k;
         const Double_t denom  = ukm1*uk-1.;
         const Double_t bkm1   = pcb[(k-1)*inc]/ukm1k;
         const Double_t bk     = pcb[k*inc]/ukm1k;
         pcb[(k-1)*inc] = (uk*bkm1-bk)/denom;
         pcb[k*inc]     = (ukm1*bk-bkm1)/denom;
         k -= 2;
      }
   }

   // Next solve u'*x = b, overwriting b with x.
   //
   //  k is the main loop index, increasing from 0 to n-1 in steps of
   //  1 or 2, depending on the size of the diagonal blocks.

   k = 0;
   while (k < n) {

      if (fIpiv[k] > 0) {
         // 1 x 1 diagonal block
         //  multiply by inv(u'(k)), where u(k) is the transformation
         //  stored in column k of a.
         for (Int_t i = 0; i < k; i++)
            pcb[k*inc] -= pU[i*n+k]*pcb[i*inc];

         // interchange elements k and ipiv(k).
         const Int_t kp = fIpiv[k]-1;
         if (kp != k) {
            const Double_t tmp = pcb[k*inc];
            pcb[k*inc]  = pcb[kp*inc];
            pcb[kp*inc] = tmp;
         }
         k++;
      } else {
         // 2 x 2 diagonal block
         // multiply by inv(u'(k+1)), where u(k+1) is the transformation
         // stored in columns k and k+1 of a.
         Int_t i;
         for (i = 0; i < k; i++)
            pcb[k*inc] -= pU[i*n+k]*pcb[i*inc];
         for (i = 0; i < k; i++)
            pcb[(k+1)*inc] -= pU[i*n+k+1]*pcb[i*inc];

         // interchange elements k and -ipiv(k).
         const Int_t kp = -fIpiv[k]-1;
         if (kp != k) {
            const Double_t tmp = pcb[k*inc];
            pcb[k*inc]  = pcb[kp*inc];
            pcb[kp*inc] = tmp;
         }
         k += 2;
      }
   }

   return kTRUE;
}

//______________________________________________________________________________
Bool_t TDecompBK::Invert(TMatrixDSym &inv)
{
// For a symmetric matrix A(m,m), its inverse A_inv(m,m) is returned .

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

   inv.UnitMatrix();

   const Int_t colLwb = inv.GetColLwb();
   const Int_t colUpb = inv.GetColUpb();
   Bool_t status = kTRUE;
   for (Int_t icol = colLwb; icol <= colUpb && status; icol++) {
      TMatrixDColumn b(inv,icol);
      status &= Solve(b);
   }

   return status;
}

//______________________________________________________________________________
TMatrixDSym TDecompBK::Invert(Bool_t &status)
{
// For a symmetric matrix A(m,m), its inverse A_inv(m,m) is returned .

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

   TMatrixDSym inv(rowLwb,rowUpb);
   inv.UnitMatrix();
   status = Invert(inv);

   return inv;
}

//______________________________________________________________________________
void TDecompBK::Print(Option_t *opt) const
{
// Print the class members

   TDecompBase::Print(opt);
   printf("fIpiv:\n");
   for (Int_t i = 0; i < fNIpiv; i++)
      printf("[%d] = %d\n",i,fIpiv[i]);
   fU.Print("fU");
}

//______________________________________________________________________________
TDecompBK &TDecompBK::operator=(const TDecompBK &source)
{
// Assigment operator

   if (this != &source) {
      TDecompBase::operator=(source);
      fU.ResizeTo(source.fU);
      fU   = source.fU;
      if (fNIpiv != source.fNIpiv) {
         if (fIpiv)
            delete [] fIpiv;
         fNIpiv = source.fNIpiv;
         fIpiv = new Int_t[fNIpiv];
      }
      memcpy(fIpiv,source.fIpiv,fNIpiv*sizeof(Int_t));
   }
   return *this;
}
 TDecompBK.cxx:1
 TDecompBK.cxx:2
 TDecompBK.cxx:3
 TDecompBK.cxx:4
 TDecompBK.cxx:5
 TDecompBK.cxx:6
 TDecompBK.cxx:7
 TDecompBK.cxx:8
 TDecompBK.cxx:9
 TDecompBK.cxx:10
 TDecompBK.cxx:11
 TDecompBK.cxx:12
 TDecompBK.cxx:13
 TDecompBK.cxx:14
 TDecompBK.cxx:15
 TDecompBK.cxx:16
 TDecompBK.cxx:17
 TDecompBK.cxx:18
 TDecompBK.cxx:19
 TDecompBK.cxx:20
 TDecompBK.cxx:21
 TDecompBK.cxx:22
 TDecompBK.cxx:23
 TDecompBK.cxx:24
 TDecompBK.cxx:25
 TDecompBK.cxx:26
 TDecompBK.cxx:27
 TDecompBK.cxx:28
 TDecompBK.cxx:29
 TDecompBK.cxx:30
 TDecompBK.cxx:31
 TDecompBK.cxx:32
 TDecompBK.cxx:33
 TDecompBK.cxx:34
 TDecompBK.cxx:35
 TDecompBK.cxx:36
 TDecompBK.cxx:37
 TDecompBK.cxx:38
 TDecompBK.cxx:39
 TDecompBK.cxx:40
 TDecompBK.cxx:41
 TDecompBK.cxx:42
 TDecompBK.cxx:43
 TDecompBK.cxx:44
 TDecompBK.cxx:45
 TDecompBK.cxx:46
 TDecompBK.cxx:47
 TDecompBK.cxx:48
 TDecompBK.cxx:49
 TDecompBK.cxx:50
 TDecompBK.cxx:51
 TDecompBK.cxx:52
 TDecompBK.cxx:53
 TDecompBK.cxx:54
 TDecompBK.cxx:55
 TDecompBK.cxx:56
 TDecompBK.cxx:57
 TDecompBK.cxx:58
 TDecompBK.cxx:59
 TDecompBK.cxx:60
 TDecompBK.cxx:61
 TDecompBK.cxx:62
 TDecompBK.cxx:63
 TDecompBK.cxx:64
 TDecompBK.cxx:65
 TDecompBK.cxx:66
 TDecompBK.cxx:67
 TDecompBK.cxx:68
 TDecompBK.cxx:69
 TDecompBK.cxx:70
 TDecompBK.cxx:71
 TDecompBK.cxx:72
 TDecompBK.cxx:73
 TDecompBK.cxx:74
 TDecompBK.cxx:75
 TDecompBK.cxx:76
 TDecompBK.cxx:77
 TDecompBK.cxx:78
 TDecompBK.cxx:79
 TDecompBK.cxx:80
 TDecompBK.cxx:81
 TDecompBK.cxx:82
 TDecompBK.cxx:83
 TDecompBK.cxx:84
 TDecompBK.cxx:85
 TDecompBK.cxx:86
 TDecompBK.cxx:87
 TDecompBK.cxx:88
 TDecompBK.cxx:89
 TDecompBK.cxx:90
 TDecompBK.cxx:91
 TDecompBK.cxx:92
 TDecompBK.cxx:93
 TDecompBK.cxx:94
 TDecompBK.cxx:95
 TDecompBK.cxx:96
 TDecompBK.cxx:97
 TDecompBK.cxx:98
 TDecompBK.cxx:99
 TDecompBK.cxx:100
 TDecompBK.cxx:101
 TDecompBK.cxx:102
 TDecompBK.cxx:103
 TDecompBK.cxx:104
 TDecompBK.cxx:105
 TDecompBK.cxx:106
 TDecompBK.cxx:107
 TDecompBK.cxx:108
 TDecompBK.cxx:109
 TDecompBK.cxx:110
 TDecompBK.cxx:111
 TDecompBK.cxx:112
 TDecompBK.cxx:113
 TDecompBK.cxx:114
 TDecompBK.cxx:115
 TDecompBK.cxx:116
 TDecompBK.cxx:117
 TDecompBK.cxx:118
 TDecompBK.cxx:119
 TDecompBK.cxx:120
 TDecompBK.cxx:121
 TDecompBK.cxx:122
 TDecompBK.cxx:123
 TDecompBK.cxx:124
 TDecompBK.cxx:125
 TDecompBK.cxx:126
 TDecompBK.cxx:127
 TDecompBK.cxx:128
 TDecompBK.cxx:129
 TDecompBK.cxx:130
 TDecompBK.cxx:131
 TDecompBK.cxx:132
 TDecompBK.cxx:133
 TDecompBK.cxx:134
 TDecompBK.cxx:135
 TDecompBK.cxx:136
 TDecompBK.cxx:137
 TDecompBK.cxx:138
 TDecompBK.cxx:139
 TDecompBK.cxx:140
 TDecompBK.cxx:141
 TDecompBK.cxx:142
 TDecompBK.cxx:143
 TDecompBK.cxx:144
 TDecompBK.cxx:145
 TDecompBK.cxx:146
 TDecompBK.cxx:147
 TDecompBK.cxx:148
 TDecompBK.cxx:149
 TDecompBK.cxx:150
 TDecompBK.cxx:151
 TDecompBK.cxx:152
 TDecompBK.cxx:153
 TDecompBK.cxx:154
 TDecompBK.cxx:155
 TDecompBK.cxx:156
 TDecompBK.cxx:157
 TDecompBK.cxx:158
 TDecompBK.cxx:159
 TDecompBK.cxx:160
 TDecompBK.cxx:161
 TDecompBK.cxx:162
 TDecompBK.cxx:163
 TDecompBK.cxx:164
 TDecompBK.cxx:165
 TDecompBK.cxx:166
 TDecompBK.cxx:167
 TDecompBK.cxx:168
 TDecompBK.cxx:169
 TDecompBK.cxx:170
 TDecompBK.cxx:171
 TDecompBK.cxx:172
 TDecompBK.cxx:173
 TDecompBK.cxx:174
 TDecompBK.cxx:175
 TDecompBK.cxx:176
 TDecompBK.cxx:177
 TDecompBK.cxx:178
 TDecompBK.cxx:179
 TDecompBK.cxx:180
 TDecompBK.cxx:181
 TDecompBK.cxx:182
 TDecompBK.cxx:183
 TDecompBK.cxx:184
 TDecompBK.cxx:185
 TDecompBK.cxx:186
 TDecompBK.cxx:187
 TDecompBK.cxx:188
 TDecompBK.cxx:189
 TDecompBK.cxx:190
 TDecompBK.cxx:191
 TDecompBK.cxx:192
 TDecompBK.cxx:193
 TDecompBK.cxx:194
 TDecompBK.cxx:195
 TDecompBK.cxx:196
 TDecompBK.cxx:197
 TDecompBK.cxx:198
 TDecompBK.cxx:199
 TDecompBK.cxx:200
 TDecompBK.cxx:201
 TDecompBK.cxx:202
 TDecompBK.cxx:203
 TDecompBK.cxx:204
 TDecompBK.cxx:205
 TDecompBK.cxx:206
 TDecompBK.cxx:207
 TDecompBK.cxx:208
 TDecompBK.cxx:209
 TDecompBK.cxx:210
 TDecompBK.cxx:211
 TDecompBK.cxx:212
 TDecompBK.cxx:213
 TDecompBK.cxx:214
 TDecompBK.cxx:215
 TDecompBK.cxx:216
 TDecompBK.cxx:217
 TDecompBK.cxx:218
 TDecompBK.cxx:219
 TDecompBK.cxx:220
 TDecompBK.cxx:221
 TDecompBK.cxx:222
 TDecompBK.cxx:223
 TDecompBK.cxx:224
 TDecompBK.cxx:225
 TDecompBK.cxx:226
 TDecompBK.cxx:227
 TDecompBK.cxx:228
 TDecompBK.cxx:229
 TDecompBK.cxx:230
 TDecompBK.cxx:231
 TDecompBK.cxx:232
 TDecompBK.cxx:233
 TDecompBK.cxx:234
 TDecompBK.cxx:235
 TDecompBK.cxx:236
 TDecompBK.cxx:237
 TDecompBK.cxx:238
 TDecompBK.cxx:239
 TDecompBK.cxx:240
 TDecompBK.cxx:241
 TDecompBK.cxx:242
 TDecompBK.cxx:243
 TDecompBK.cxx:244
 TDecompBK.cxx:245
 TDecompBK.cxx:246
 TDecompBK.cxx:247
 TDecompBK.cxx:248
 TDecompBK.cxx:249
 TDecompBK.cxx:250
 TDecompBK.cxx:251
 TDecompBK.cxx:252
 TDecompBK.cxx:253
 TDecompBK.cxx:254
 TDecompBK.cxx:255
 TDecompBK.cxx:256
 TDecompBK.cxx:257
 TDecompBK.cxx:258
 TDecompBK.cxx:259
 TDecompBK.cxx:260
 TDecompBK.cxx:261
 TDecompBK.cxx:262
 TDecompBK.cxx:263
 TDecompBK.cxx:264
 TDecompBK.cxx:265
 TDecompBK.cxx:266
 TDecompBK.cxx:267
 TDecompBK.cxx:268
 TDecompBK.cxx:269
 TDecompBK.cxx:270
 TDecompBK.cxx:271
 TDecompBK.cxx:272
 TDecompBK.cxx:273
 TDecompBK.cxx:274
 TDecompBK.cxx:275
 TDecompBK.cxx:276
 TDecompBK.cxx:277
 TDecompBK.cxx:278
 TDecompBK.cxx:279
 TDecompBK.cxx:280
 TDecompBK.cxx:281
 TDecompBK.cxx:282
 TDecompBK.cxx:283
 TDecompBK.cxx:284
 TDecompBK.cxx:285
 TDecompBK.cxx:286
 TDecompBK.cxx:287
 TDecompBK.cxx:288
 TDecompBK.cxx:289
 TDecompBK.cxx:290
 TDecompBK.cxx:291
 TDecompBK.cxx:292
 TDecompBK.cxx:293
 TDecompBK.cxx:294
 TDecompBK.cxx:295
 TDecompBK.cxx:296
 TDecompBK.cxx:297
 TDecompBK.cxx:298
 TDecompBK.cxx:299
 TDecompBK.cxx:300
 TDecompBK.cxx:301
 TDecompBK.cxx:302
 TDecompBK.cxx:303
 TDecompBK.cxx:304
 TDecompBK.cxx:305
 TDecompBK.cxx:306
 TDecompBK.cxx:307
 TDecompBK.cxx:308
 TDecompBK.cxx:309
 TDecompBK.cxx:310
 TDecompBK.cxx:311
 TDecompBK.cxx:312
 TDecompBK.cxx:313
 TDecompBK.cxx:314
 TDecompBK.cxx:315
 TDecompBK.cxx:316
 TDecompBK.cxx:317
 TDecompBK.cxx:318
 TDecompBK.cxx:319
 TDecompBK.cxx:320
 TDecompBK.cxx:321
 TDecompBK.cxx:322
 TDecompBK.cxx:323
 TDecompBK.cxx:324
 TDecompBK.cxx:325
 TDecompBK.cxx:326
 TDecompBK.cxx:327
 TDecompBK.cxx:328
 TDecompBK.cxx:329
 TDecompBK.cxx:330
 TDecompBK.cxx:331
 TDecompBK.cxx:332
 TDecompBK.cxx:333
 TDecompBK.cxx:334
 TDecompBK.cxx:335
 TDecompBK.cxx:336
 TDecompBK.cxx:337
 TDecompBK.cxx:338
 TDecompBK.cxx:339
 TDecompBK.cxx:340
 TDecompBK.cxx:341
 TDecompBK.cxx:342
 TDecompBK.cxx:343
 TDecompBK.cxx:344
 TDecompBK.cxx:345
 TDecompBK.cxx:346
 TDecompBK.cxx:347
 TDecompBK.cxx:348
 TDecompBK.cxx:349
 TDecompBK.cxx:350
 TDecompBK.cxx:351
 TDecompBK.cxx:352
 TDecompBK.cxx:353
 TDecompBK.cxx:354
 TDecompBK.cxx:355
 TDecompBK.cxx:356
 TDecompBK.cxx:357
 TDecompBK.cxx:358
 TDecompBK.cxx:359
 TDecompBK.cxx:360
 TDecompBK.cxx:361
 TDecompBK.cxx:362
 TDecompBK.cxx:363
 TDecompBK.cxx:364
 TDecompBK.cxx:365
 TDecompBK.cxx:366
 TDecompBK.cxx:367
 TDecompBK.cxx:368
 TDecompBK.cxx:369
 TDecompBK.cxx:370
 TDecompBK.cxx:371
 TDecompBK.cxx:372
 TDecompBK.cxx:373
 TDecompBK.cxx:374
 TDecompBK.cxx:375
 TDecompBK.cxx:376
 TDecompBK.cxx:377
 TDecompBK.cxx:378
 TDecompBK.cxx:379
 TDecompBK.cxx:380
 TDecompBK.cxx:381
 TDecompBK.cxx:382
 TDecompBK.cxx:383
 TDecompBK.cxx:384
 TDecompBK.cxx:385
 TDecompBK.cxx:386
 TDecompBK.cxx:387
 TDecompBK.cxx:388
 TDecompBK.cxx:389
 TDecompBK.cxx:390
 TDecompBK.cxx:391
 TDecompBK.cxx:392
 TDecompBK.cxx:393
 TDecompBK.cxx:394
 TDecompBK.cxx:395
 TDecompBK.cxx:396
 TDecompBK.cxx:397
 TDecompBK.cxx:398
 TDecompBK.cxx:399
 TDecompBK.cxx:400
 TDecompBK.cxx:401
 TDecompBK.cxx:402
 TDecompBK.cxx:403
 TDecompBK.cxx:404
 TDecompBK.cxx:405
 TDecompBK.cxx:406
 TDecompBK.cxx:407
 TDecompBK.cxx:408
 TDecompBK.cxx:409
 TDecompBK.cxx:410
 TDecompBK.cxx:411
 TDecompBK.cxx:412
 TDecompBK.cxx:413
 TDecompBK.cxx:414
 TDecompBK.cxx:415
 TDecompBK.cxx:416
 TDecompBK.cxx:417
 TDecompBK.cxx:418
 TDecompBK.cxx:419
 TDecompBK.cxx:420
 TDecompBK.cxx:421
 TDecompBK.cxx:422
 TDecompBK.cxx:423
 TDecompBK.cxx:424
 TDecompBK.cxx:425
 TDecompBK.cxx:426
 TDecompBK.cxx:427
 TDecompBK.cxx:428
 TDecompBK.cxx:429
 TDecompBK.cxx:430
 TDecompBK.cxx:431
 TDecompBK.cxx:432
 TDecompBK.cxx:433
 TDecompBK.cxx:434
 TDecompBK.cxx:435
 TDecompBK.cxx:436
 TDecompBK.cxx:437
 TDecompBK.cxx:438
 TDecompBK.cxx:439
 TDecompBK.cxx:440
 TDecompBK.cxx:441
 TDecompBK.cxx:442
 TDecompBK.cxx:443
 TDecompBK.cxx:444
 TDecompBK.cxx:445
 TDecompBK.cxx:446
 TDecompBK.cxx:447
 TDecompBK.cxx:448
 TDecompBK.cxx:449
 TDecompBK.cxx:450
 TDecompBK.cxx:451
 TDecompBK.cxx:452
 TDecompBK.cxx:453
 TDecompBK.cxx:454
 TDecompBK.cxx:455
 TDecompBK.cxx:456
 TDecompBK.cxx:457
 TDecompBK.cxx:458
 TDecompBK.cxx:459
 TDecompBK.cxx:460
 TDecompBK.cxx:461
 TDecompBK.cxx:462
 TDecompBK.cxx:463
 TDecompBK.cxx:464
 TDecompBK.cxx:465
 TDecompBK.cxx:466
 TDecompBK.cxx:467
 TDecompBK.cxx:468
 TDecompBK.cxx:469
 TDecompBK.cxx:470
 TDecompBK.cxx:471
 TDecompBK.cxx:472
 TDecompBK.cxx:473
 TDecompBK.cxx:474
 TDecompBK.cxx:475
 TDecompBK.cxx:476
 TDecompBK.cxx:477
 TDecompBK.cxx:478
 TDecompBK.cxx:479
 TDecompBK.cxx:480
 TDecompBK.cxx:481
 TDecompBK.cxx:482
 TDecompBK.cxx:483
 TDecompBK.cxx:484
 TDecompBK.cxx:485
 TDecompBK.cxx:486
 TDecompBK.cxx:487
 TDecompBK.cxx:488
 TDecompBK.cxx:489
 TDecompBK.cxx:490
 TDecompBK.cxx:491
 TDecompBK.cxx:492
 TDecompBK.cxx:493
 TDecompBK.cxx:494
 TDecompBK.cxx:495
 TDecompBK.cxx:496
 TDecompBK.cxx:497
 TDecompBK.cxx:498
 TDecompBK.cxx:499
 TDecompBK.cxx:500
 TDecompBK.cxx:501
 TDecompBK.cxx:502
 TDecompBK.cxx:503
 TDecompBK.cxx:504
 TDecompBK.cxx:505
 TDecompBK.cxx:506
 TDecompBK.cxx:507
 TDecompBK.cxx:508
 TDecompBK.cxx:509
 TDecompBK.cxx:510
 TDecompBK.cxx:511
 TDecompBK.cxx:512
 TDecompBK.cxx:513
 TDecompBK.cxx:514
 TDecompBK.cxx:515
 TDecompBK.cxx:516
 TDecompBK.cxx:517
 TDecompBK.cxx:518
 TDecompBK.cxx:519
 TDecompBK.cxx:520
 TDecompBK.cxx:521
 TDecompBK.cxx:522
 TDecompBK.cxx:523
 TDecompBK.cxx:524
 TDecompBK.cxx:525
 TDecompBK.cxx:526
 TDecompBK.cxx:527
 TDecompBK.cxx:528
 TDecompBK.cxx:529
 TDecompBK.cxx:530
 TDecompBK.cxx:531
 TDecompBK.cxx:532
 TDecompBK.cxx:533
 TDecompBK.cxx:534
 TDecompBK.cxx:535
 TDecompBK.cxx:536
 TDecompBK.cxx:537
 TDecompBK.cxx:538
 TDecompBK.cxx:539
 TDecompBK.cxx:540
 TDecompBK.cxx:541
 TDecompBK.cxx:542
 TDecompBK.cxx:543
 TDecompBK.cxx:544
 TDecompBK.cxx:545
 TDecompBK.cxx:546
 TDecompBK.cxx:547
 TDecompBK.cxx:548
 TDecompBK.cxx:549
 TDecompBK.cxx:550
 TDecompBK.cxx:551
 TDecompBK.cxx:552
 TDecompBK.cxx:553
 TDecompBK.cxx:554
 TDecompBK.cxx:555
 TDecompBK.cxx:556
 TDecompBK.cxx:557
 TDecompBK.cxx:558
 TDecompBK.cxx:559
 TDecompBK.cxx:560
 TDecompBK.cxx:561
 TDecompBK.cxx:562
 TDecompBK.cxx:563
 TDecompBK.cxx:564
 TDecompBK.cxx:565
 TDecompBK.cxx:566
 TDecompBK.cxx:567
 TDecompBK.cxx:568
 TDecompBK.cxx:569
 TDecompBK.cxx:570
 TDecompBK.cxx:571
 TDecompBK.cxx:572
 TDecompBK.cxx:573
 TDecompBK.cxx:574
 TDecompBK.cxx:575
 TDecompBK.cxx:576
 TDecompBK.cxx:577
 TDecompBK.cxx:578
 TDecompBK.cxx:579
 TDecompBK.cxx:580
 TDecompBK.cxx:581
 TDecompBK.cxx:582
 TDecompBK.cxx:583
 TDecompBK.cxx:584
 TDecompBK.cxx:585
 TDecompBK.cxx:586
 TDecompBK.cxx:587
 TDecompBK.cxx:588
 TDecompBK.cxx:589
 TDecompBK.cxx:590
 TDecompBK.cxx:591
 TDecompBK.cxx:592
 TDecompBK.cxx:593
 TDecompBK.cxx:594
 TDecompBK.cxx:595
 TDecompBK.cxx:596
 TDecompBK.cxx:597
 TDecompBK.cxx:598
 TDecompBK.cxx:599
 TDecompBK.cxx:600
 TDecompBK.cxx:601
 TDecompBK.cxx:602
 TDecompBK.cxx:603
 TDecompBK.cxx:604
 TDecompBK.cxx:605
 TDecompBK.cxx:606
 TDecompBK.cxx:607
 TDecompBK.cxx:608
 TDecompBK.cxx:609
 TDecompBK.cxx:610
 TDecompBK.cxx:611
 TDecompBK.cxx:612
 TDecompBK.cxx:613
 TDecompBK.cxx:614
 TDecompBK.cxx:615
 TDecompBK.cxx:616
 TDecompBK.cxx:617
 TDecompBK.cxx:618
 TDecompBK.cxx:619
 TDecompBK.cxx:620
 TDecompBK.cxx:621
 TDecompBK.cxx:622
 TDecompBK.cxx:623
 TDecompBK.cxx:624
 TDecompBK.cxx:625
 TDecompBK.cxx:626
 TDecompBK.cxx:627
 TDecompBK.cxx:628
 TDecompBK.cxx:629
 TDecompBK.cxx:630
 TDecompBK.cxx:631
 TDecompBK.cxx:632
 TDecompBK.cxx:633
 TDecompBK.cxx:634
 TDecompBK.cxx:635
 TDecompBK.cxx:636
 TDecompBK.cxx:637
 TDecompBK.cxx:638
 TDecompBK.cxx:639
 TDecompBK.cxx:640
 TDecompBK.cxx:641
 TDecompBK.cxx:642
 TDecompBK.cxx:643
 TDecompBK.cxx:644
 TDecompBK.cxx:645
 TDecompBK.cxx:646
 TDecompBK.cxx:647
 TDecompBK.cxx:648
 TDecompBK.cxx:649
 TDecompBK.cxx:650
 TDecompBK.cxx:651
 TDecompBK.cxx:652
 TDecompBK.cxx:653
 TDecompBK.cxx:654
 TDecompBK.cxx:655
 TDecompBK.cxx:656
 TDecompBK.cxx:657
 TDecompBK.cxx:658
 TDecompBK.cxx:659
 TDecompBK.cxx:660
 TDecompBK.cxx:661
 TDecompBK.cxx:662
 TDecompBK.cxx:663
 TDecompBK.cxx:664
 TDecompBK.cxx:665
 TDecompBK.cxx:666
 TDecompBK.cxx:667
 TDecompBK.cxx:668
 TDecompBK.cxx:669
 TDecompBK.cxx:670