ROOT logo
// @(#)root/matrix:$Id: TDecompSVD.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.             *
 *************************************************************************/

///////////////////////////////////////////////////////////////////////////
//                                                                       //
// Single Value Decomposition class                                      //
//                                                                       //
// For an (m x n) matrix A with m >= n, the singular value decomposition //
// is                                                                    //
// an (m x m) orthogonal matrix fU, an (m x n) diagonal matrix fS, and   //
// an (n x n) orthogonal matrix fV so that A = U*S*V'.                   //
//                                                                       //
// If the row/column index of A starts at (rowLwb,colLwb) then the       //
// decomposed matrices/vectors start at :                                //
//  fU   : (rowLwb,colLwb)                                               //
//  fV   : (colLwb,colLwb)                                               //
//  fSig : (colLwb)                                                      //
//                                                                       //
// The diagonal matrix fS is stored in the singular values vector fSig . //
// The singular values, fSig[k] = S[k][k], are ordered so that           //
// fSig[0] >= fSig[1] >= ... >= fSig[n-1].                               //
//                                                                       //
// The singular value decompostion always exists, so the decomposition   //
// will (as long as m >=n) never fail. If m < n, the user should add     //
// sufficient zero rows to A , so that m == n                            //
//                                                                       //
// Here fTol is used to set the threshold on the minimum allowed value   //
// of the singular values:                                               //
//  min_singular = fTol*max(fSig[i])                                     //
//                                                                       //
///////////////////////////////////////////////////////////////////////////

#include "TDecompSVD.h"
#include "TMath.h"
#include "TArrayD.h"

ClassImp(TDecompSVD)

//______________________________________________________________________________
TDecompSVD::TDecompSVD(Int_t nrows,Int_t ncols)
{
// Constructor for (nrows x ncols) matrix

   if (nrows < ncols) {
      Error("TDecompSVD(Int_t,Int_t","matrix rows should be >= columns");
      return;
   }
   fU.ResizeTo(nrows,nrows);
   fSig.ResizeTo(ncols);
   fV.ResizeTo(nrows,ncols); // In the end we only need the nColxnCol part
}

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

   const Int_t nrows = row_upb-row_lwb+1;
   const Int_t ncols = col_upb-col_lwb+1;

   if (nrows < ncols) {
      Error("TDecompSVD(Int_t,Int_t,Int_t,Int_t","matrix rows should be >= columns");
      return;
   }
   fRowLwb = row_lwb;
   fColLwb = col_lwb;
   fU.ResizeTo(nrows,nrows);
   fSig.ResizeTo(ncols);
   fV.ResizeTo(nrows,ncols); // In the end we only need the nColxnCol part
}

//______________________________________________________________________________
TDecompSVD::TDecompSVD(const TMatrixD &a,Double_t tol)
{
// Constructor for general matrix A .

   R__ASSERT(a.IsValid());
   if (a.GetNrows() < a.GetNcols()) {
      Error("TDecompSVD(const TMatrixD &","matrix rows should be >= columns");
      return;
   }

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

   fRowLwb = a.GetRowLwb();
   fColLwb = a.GetColLwb();
   const Int_t nRow = a.GetNrows();
   const Int_t nCol = a.GetNcols();

   fU.ResizeTo(nRow,nRow);
   fSig.ResizeTo(nCol);
   fV.ResizeTo(nRow,nCol); // In the end we only need the nColxnCol part

   fU.UnitMatrix();
   memcpy(fV.GetMatrixArray(),a.GetMatrixArray(),nRow*nCol*sizeof(Double_t));
}

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

   *this = another;
}

//______________________________________________________________________________
Bool_t TDecompSVD::Decompose()
{
// SVD decomposition of matrix
// 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;
   }

   const Int_t nCol   = this->GetNcols();
   const Int_t rowLwb = this->GetRowLwb();
   const Int_t colLwb = this->GetColLwb();

   TVectorD offDiag;
   Double_t work[kWorkMax];
   if (nCol > kWorkMax) offDiag.ResizeTo(nCol);
   else                 offDiag.Use(nCol,work);

   // step 1: bidiagonalization of A
   if (!Bidiagonalize(fV,fU,fSig,offDiag))
      return kFALSE;

   // step 2: diagonalization of bidiagonal matrix
   if (!Diagonalize(fV,fU,fSig,offDiag))
      return kFALSE;

   // step 3: order singular values and perform permutations
   SortSingular(fV,fU,fSig);
   fV.ResizeTo(nCol,nCol); fV.Shift(colLwb,colLwb);
   fSig.Shift(colLwb);
   fU.Transpose(fU);       fU.Shift(rowLwb,colLwb);
   SetBit(kDecomposed);

   return kTRUE;
}

//______________________________________________________________________________
Bool_t TDecompSVD::Bidiagonalize(TMatrixD &v,TMatrixD &u,TVectorD &sDiag,TVectorD &oDiag)
{
// Bidiagonalize the (m x n) - matrix a (stored in v) through a series of Householder
// transformations applied to the left (Q^T) and to the right (H) of a ,
// so that A = Q . C . H^T with matrix C bidiagonal. Q and H are orthogonal matrices .
//
// Output:
//   v     - (n x n) - matrix H in the (n x n) part of v
//   u     - (m x m) - matrix Q^T
//   sDiag - diagonal of the (m x n) C
//   oDiag - off-diagonal elements of matrix C
//
//  Test code for the output:
//    const Int_t nRow = v.GetNrows();
//    const Int_t nCol = v.GetNcols();
//    TMatrixD H(v); H.ResizeTo(nCol,nCol);
//    TMatrixD E1(nCol,nCol); E1.UnitMatrix();
//    TMatrixD Ht(TMatrixDBase::kTransposed,H);
//    Bool_t ok = kTRUE;
//    ok &= VerifyMatrixIdentity(Ht * H,E1,kTRUE,1.0e-13);
//    ok &= VerifyMatrixIdentity(H * Ht,E1,kTRUE,1.0e-13);
//    TMatrixD E2(nRow,nRow); E2.UnitMatrix();
//    TMatrixD Qt(u);
//    TMatrixD Q(TMatrixDBase::kTransposed,Qt);
//    ok &= VerifyMatrixIdentity(Q * Qt,E2,kTRUE,1.0e-13);
//    TMatrixD C(nRow,nCol);
//    TMatrixDDiag(C) = sDiag;
//    for (Int_t i = 0; i < nCol-1; i++)
//      C(i,i+1) = oDiag(i+1);
//    TMatrixD A = Q*C*Ht;
//    ok &= VerifyMatrixIdentity(A,a,kTRUE,1.0e-13);

   const Int_t nRow_v = v.GetNrows();
   const Int_t nCol_v = v.GetNcols();
   const Int_t nCol_u = u.GetNcols();

   TArrayD ups(nCol_v);
   TArrayD betas(nCol_v);

   for (Int_t i = 0; i < nCol_v; i++) {
      // Set up Householder Transformation q(i)

      Double_t up,beta;
      if (i < nCol_v-1 || nRow_v > nCol_v) {
         const TVectorD vc_i = TMatrixDColumn_const(v,i);
         //if (!DefHouseHolder(vc_i,i,i+1,up,beta))
         //  return kFALSE;
         DefHouseHolder(vc_i,i,i+1,up,beta);

         // Apply q(i) to v
         for (Int_t j = i; j < nCol_v; j++) {
            TMatrixDColumn vc_j = TMatrixDColumn(v,j);
            ApplyHouseHolder(vc_i,up,beta,i,i+1,vc_j);
         }

         // Apply q(i) to u
         for (Int_t j = 0; j < nCol_u; j++)
         {
            TMatrixDColumn uc_j = TMatrixDColumn(u,j);
            ApplyHouseHolder(vc_i,up,beta,i,i+1,uc_j);
         }
      }
      if (i < nCol_v-2) {
         // set up Householder Transformation h(i)
         const TVectorD vr_i = TMatrixDRow_const(v,i);

         //if (!DefHouseHolder(vr_i,i+1,i+2,up,beta))
         //  return kFALSE;
         DefHouseHolder(vr_i,i+1,i+2,up,beta);

         // save h(i)
         ups[i]   = up;
         betas[i] = beta;

         // apply h(i) to v
         for (Int_t j = i; j < nRow_v; j++) {
            TMatrixDRow vr_j = TMatrixDRow(v,j);
            ApplyHouseHolder(vr_i,up,beta,i+1,i+2,vr_j);

            // save elements i+2,...in row j of matrix v
            if (j == i) {
               for (Int_t k = i+2; k < nCol_v; k++)
                  vr_j(k) = vr_i(k);
            }
         }
      }
   }

   // copy diagonal of transformed matrix v to sDiag and upper parallel v to oDiag
   if (nCol_v > 1) {
      sDiag = TMatrixDDiag(v);
      for (Int_t i = 1; i < nCol_v; i++)
         oDiag(i) = v(i-1,i);
   }
   oDiag(0) = 0.;

   // construct product matrix h = h(1)*h(2)*...*h(nCol_v-1), h(nCol_v-1) = I

   TVectorD vr_i(nCol_v);
   for (Int_t i = nCol_v-1; i >= 0; i--) {
      if (i < nCol_v-1)
         vr_i = TMatrixDRow_const(v,i);
      TMatrixDRow(v,i) = 0.0;
      v(i,i) = 1.;

      if (i < nCol_v-2) {
         for (Int_t k = i; k < nCol_v; k++) {
            // householder transformation on k-th column
            TMatrixDColumn vc_k = TMatrixDColumn(v,k);
            ApplyHouseHolder(vr_i,ups[i],betas[i],i+1,i+2,vc_k);
         }
      }
   }

   return kTRUE;
}

//______________________________________________________________________________
Bool_t TDecompSVD::Diagonalize(TMatrixD &v,TMatrixD &u,TVectorD &sDiag,TVectorD &oDiag)
{
// Diagonalizes in an iterative fashion the bidiagonal matrix C as described through
// sDiag and oDiag, so that S' = U'^T . C . V' is diagonal. U' and V' are orthogonal
// matrices .
//
// Output:
//   v     - (n x n) - matrix H . V' in the (n x n) part of v
//   u     - (m x m) - matrix U'^T . Q^T
//   sDiag - diagonal of the (m x n) S'
//
//   return convergence flag:  0 -> no convergence
//                             1 -> convergence
//
//  Test code for the output:
//    const Int_t nRow = v.GetNrows();
//    const Int_t nCol = v.GetNcols();
//    TMatrixD tmp = v; tmp.ResizeTo(nCol,nCol);
//    TMatrixD Vprime  = Ht*tmp;
//    TMatrixD Vprimet(TMatrixDBase::kTransposed,Vprime);
//    TMatrixD Uprimet = u*Q;
//    TMatrixD Uprime(TMatrixDBase::kTransposed,Uprimet);
//    TMatrixD Sprime(nRow,nCol);
//    TMatrixDDiag(Sprime) = sDiag;
//    ok &= VerifyMatrixIdentity(Uprimet * C * Vprime,Sprime,kTRUE,1.0e-13);
//    ok &= VerifyMatrixIdentity(Q*Uprime * Sprime * Vprimet * Ht,a,kTRUE,1.0e-13);

   Bool_t ok    = kTRUE;
   Int_t niter  = 0;
   Double_t bmx = sDiag(0);

   const Int_t nCol_v = v.GetNcols();

   if (nCol_v > 1) {
      for (Int_t i = 1; i < nCol_v; i++)
         bmx = TMath::Max(TMath::Abs(sDiag(i))+TMath::Abs(oDiag(i)),bmx);
   }

   const Double_t eps = std::numeric_limits<double>::epsilon();

   const Int_t niterm = 10*nCol_v;
   for (Int_t k = nCol_v-1; k >= 0; k--) {
      loop:
         if (k != 0) {
            // since sDiag(k) == 0 perform Givens transform with result oDiag[k] = 0
            if (TMath::Abs(sDiag(k)) < eps*bmx)
               Diag_1(v,sDiag,oDiag,k);

            // find l (1 <= l <=k) so that either oDiag(l) = 0 or sDiag(l-1) = 0.
            // In the latter case transform oDiag(l) to zero. In both cases the matrix
            // splits and the bottom right minor begins with row l.
            // If no such l is found set l = 1.

            Int_t elzero = 0;
            Int_t l = 0;
            for (Int_t ll = k; ll >= 0 ; ll--) {
               l = ll;
               if (l == 0) {
                  elzero = 0;
                  break;
               } else if (TMath::Abs(oDiag(l)) < eps*bmx) {
                  elzero = 1;
                  break;
               } else if (TMath::Abs(sDiag(l-1)) < eps*bmx)
                  elzero = 0;
            }
            if (l > 0 && !elzero)
               Diag_2(sDiag,oDiag,k,l);
            if (l != k) {
               // one more QR pass with order k
               Diag_3(v,u,sDiag,oDiag,k,l);
               niter++;
               if (niter <= niterm) goto loop;
               ::Error("TDecompSVD::Diagonalize","no convergence after %d steps",niter);
               ok = kFALSE;
            }
         }

         if (sDiag(k) < 0.) {
            // for negative singular values perform change of sign
            sDiag(k) = -sDiag(k);
            TMatrixDColumn(v,k) *= -1.0;
         }
      // order is decreased by one in next pass
   }

   return ok;
}

//______________________________________________________________________________
void TDecompSVD::Diag_1(TMatrixD &v,TVectorD &sDiag,TVectorD &oDiag,Int_t k)
{
// Step 1 in the matrix diagonalization

   const Int_t nCol_v = v.GetNcols();

   TMatrixDColumn vc_k = TMatrixDColumn(v,k);
   for (Int_t i = k-1; i >= 0; i--) {
      TMatrixDColumn vc_i = TMatrixDColumn(v,i);
      Double_t h,cs,sn;
      if (i == k-1)
         DefAplGivens(sDiag[i],oDiag[i+1],cs,sn);
      else
         DefAplGivens(sDiag[i],h,cs,sn);
      if (i > 0) {
         h = 0.;
         ApplyGivens(oDiag[i],h,cs,sn);
      }
      for (Int_t j = 0; j < nCol_v; j++)
         ApplyGivens(vc_i(j),vc_k(j),cs,sn);
   }
}

//______________________________________________________________________________
void TDecompSVD::Diag_2(TVectorD &sDiag,TVectorD &oDiag,Int_t k,Int_t l)
{
// Step 2 in the matrix diagonalization

   for (Int_t i = l; i <= k; i++) {
      Double_t h,cs,sn;
      if (i == l)
         DefAplGivens(sDiag(i),oDiag(i),cs,sn);
      else
         DefAplGivens(sDiag(i),h,cs,sn);
      if (i < k) {
         h = 0.;
         ApplyGivens(oDiag(i+1),h,cs,sn);
      }
   }
}

//______________________________________________________________________________
void TDecompSVD::Diag_3(TMatrixD &v,TMatrixD &u,TVectorD &sDiag,TVectorD &oDiag,Int_t k,Int_t l)
{
// Step 3 in the matrix diagonalization

   Double_t *pS = sDiag.GetMatrixArray();
   Double_t *pO = oDiag.GetMatrixArray();

   // determine shift parameter

   const Double_t psk1 = pS[k-1];
   const Double_t psk  = pS[k];
   const Double_t pok1 = pO[k-1];
   const Double_t pok  = pO[k];
   const Double_t psl  = pS[l];

   Double_t f;
   if (psl == 0.0 || pok == 0.0 || psk1 == 0.0) {
      const Double_t b = ((psk1-psk)*(psk1+psk)+pok1*pok1)/2.0;
      const Double_t c = (psk*pok1)*(psk*pok1);

      Double_t shift = 0.0;
      if ((b != 0.0) | (c != 0.0)) {
         shift = TMath::Sqrt(b*b+c);
         if (b < 0.0)
            shift = -shift;
         shift = c/(b+shift);
      }

      f = (psl+psk)*(psl-psk)+shift;
   } else {
      f = ((psk1-psk)*(psk1+psk)+(pok1-pok)*(pok1+pok))/(2.*pok*psk1);
      const Double_t g = TMath::Hypot(1.,f);
      const Double_t t = (f >= 0.) ? f+g : f-g;

      f = ((psl-psk)*(psl+psk)+pok*(psk1/t-pok))/psl;
   }

   const Int_t nCol_v = v.GetNcols();
   const Int_t nCol_u = u.GetNcols();

   Double_t h,cs,sn;
   Int_t j;
   for (Int_t i = l; i <= k-1; i++) {
      if (i == l)
         // define r[l]
         DefGivens(f,pO[i+1],cs,sn);
      else
         // define r[i]
         DefAplGivens(pO[i],h,cs,sn);

      ApplyGivens(pS[i],pO[i+1],cs,sn);
      h = 0.;
      ApplyGivens(h,pS[i+1],cs,sn);

      TMatrixDColumn vc_i  = TMatrixDColumn(v,i);
      TMatrixDColumn vc_i1 = TMatrixDColumn(v,i+1);
      for (j = 0; j < nCol_v; j++)
         ApplyGivens(vc_i(j),vc_i1(j),cs,sn);
      // define t[i]
      DefAplGivens(pS[i],h,cs,sn);
      ApplyGivens(pO[i+1],pS[i+1],cs,sn);
      if (i < k-1) {
         h = 0.;
         ApplyGivens(h,pO[i+2],cs,sn);
      }

      TMatrixDRow ur_i  = TMatrixDRow(u,i);
      TMatrixDRow ur_i1 = TMatrixDRow(u,i+1);
      for (j = 0; j < nCol_u; j++)
         ApplyGivens(ur_i(j),ur_i1(j),cs,sn);
   }
}

//______________________________________________________________________________
void TDecompSVD::SortSingular(TMatrixD &v,TMatrixD &u,TVectorD &sDiag)
{
// Perform a permutation transformation on the diagonal matrix S', so that
// matrix S'' = U''^T . S' . V''  has diagonal elements ordered such that they
// do not increase.
//
// Output:
//   v     - (n x n) - matrix H . V' . V'' in the (n x n) part of v
//   u     - (m x m) - matrix U''^T . U'^T . Q^T
//   sDiag - diagonal of the (m x n) S''

   const Int_t nCol_v = v.GetNcols();
   const Int_t nCol_u = u.GetNcols();

   Double_t *pS = sDiag.GetMatrixArray();
   Double_t *pV = v.GetMatrixArray();
   Double_t *pU = u.GetMatrixArray();

   // order singular values

   Int_t i,j;
   if (nCol_v > 1) {
      while (1) {
         Bool_t found = kFALSE;
         i = 1;
         while (!found && i < nCol_v) {
            if (pS[i] > pS[i-1])
               found = kTRUE;
            else
               i++;
         }
         if (!found) break;
         for (i = 1; i < nCol_v; i++) {
            Double_t t = pS[i-1];
            Int_t k = i-1;
            for (j = i; j < nCol_v; j++) {
               if (t < pS[j]) {
                  t = pS[j];
                  k = j;
               }
            }
            if (k != i-1) {
               // perform permutation on singular values
               pS[k]   = pS[i-1];
               pS[i-1] = t;
               // perform permutation on matrix v
               for (j = 0; j < nCol_v; j++) {
                  const Int_t off_j = j*nCol_v;
                  t             = pV[off_j+k];
                  pV[off_j+k]   = pV[off_j+i-1];
                  pV[off_j+i-1] = t;
               }
               // perform permutation on vector u
               for (j = 0; j < nCol_u; j++) {
                  const Int_t off_k  = k*nCol_u;
                  const Int_t off_i1 = (i-1)*nCol_u;
                  t            = pU[off_k+j];
                  pU[off_k+j]  = pU[off_i1+j];
                  pU[off_i1+j] = t;
               }
            }
         }
      }
   }
}

//______________________________________________________________________________
const TMatrixD TDecompSVD::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();
      }
   }

   const Int_t nRows  = fU.GetNrows();
   const Int_t nCols  = fV.GetNcols();
   const Int_t colLwb = this->GetColLwb();
   TMatrixD s(nRows,nCols); s.Shift(colLwb,colLwb);
   TMatrixDDiag diag(s); diag = fSig;
   const TMatrixD vt(TMatrixD::kTransposed,fV);
   return fU * s * vt;
}

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

   R__ASSERT(a.IsValid());

   ResetStatus();
   if (a.GetNrows() < a.GetNcols()) {
      Error("TDecompSVD(const TMatrixD &","matrix rows should be >= columns");
      return;
   }

   SetBit(kMatrixSet);
   fCondition = -1.0;

   fRowLwb = a.GetRowLwb();
   fColLwb = a.GetColLwb();
   const Int_t nRow = a.GetNrows();
   const Int_t nCol = a.GetNcols();

   fU.ResizeTo(nRow,nRow);
   fSig.ResizeTo(nCol);
   fV.ResizeTo(nRow,nCol); // In the end we only need the nColxnCol part

   fU.UnitMatrix();
   memcpy(fV.GetMatrixArray(),a.GetMatrixArray(),nRow*nCol*sizeof(Double_t));
}

//______________________________________________________________________________
Bool_t TDecompSVD::Solve(TVectorD &b)
{
// Solve Ax=b assuming the SVD form of A is stored . Solution returned in b.
// If A is of size (m x n), input vector b should be of size (m), however,
// the solution, returned in b, will be in the first (n) elements .
//
// For m > n , x  is the least-squares solution of min(A . x - b)

   R__ASSERT(b.IsValid());
   if (TestBit(kSingular)) {
      return kFALSE;
   }
   if ( !TestBit(kDecomposed) ) {
      if (!Decompose()) {
         return kFALSE;
      }
   }

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

   // We start with fU fSig fV^T x = b, and turn it into  fV^T x = fSig^-1 fU^T b
   // Form tmp = fSig^-1 fU^T b but ignore diagonal elements with
   // fSig(i) < fTol * max(fSig)

   const Int_t    lwb       = fU.GetColLwb();
   const Int_t    upb       = lwb+fV.GetNcols()-1;
   const Double_t threshold = fSig(lwb)*fTol;

   TVectorD tmp(lwb,upb);
   for (Int_t irow = lwb; irow <= upb; irow++) {
      Double_t r = 0.0;
      if (fSig(irow) > threshold) {
         const TVectorD uc_i = TMatrixDColumn(fU,irow);
         r = uc_i * b;
         r /= fSig(irow);
      }
      tmp(irow) = r;
   }

   if (b.GetNrows() > fV.GetNrows()) {
      TVectorD tmp2;
      tmp2.Use(lwb,upb,b.GetMatrixArray());
      tmp2 = fV*tmp;
   } else
      b = fV*tmp;

   return kTRUE;
}

//______________________________________________________________________________
Bool_t TDecompSVD::Solve(TMatrixDColumn &cb)
{
// Solve Ax=b assuming the SVD form of A is stored . Solution returned in the
// matrix column cb b.
// If A is of size (m x n), input vector b should be of size (m), however,
// the solution, returned in b, will be in the first (n) elements .
//
// For m > n , x  is the least-squares solution of min(A . x - b)

   TMatrixDBase *b = const_cast<TMatrixDBase *>(cb.GetMatrix());
   R__ASSERT(b->IsValid());
   if (TestBit(kSingular)) {
      return kFALSE;
   }
   if ( !TestBit(kDecomposed) ) {
      if (!Decompose()) {
         return kFALSE;
      }
   }

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

   // We start with fU fSig fV^T x = b, and turn it into  fV^T x = fSig^-1 fU^T b
   // Form tmp = fSig^-1 fU^T b but ignore diagonal elements in
   // fSig(i) < fTol * max(fSig)

   const Int_t    lwb       = fU.GetColLwb();
   const Int_t    upb       = lwb+fV.GetNcols()-1;
   const Double_t threshold = fSig(lwb)*fTol;

   TVectorD tmp(lwb,upb);
   const TVectorD vb = cb;
   for (Int_t irow = lwb; irow <= upb; irow++) {
      Double_t r = 0.0;
      if (fSig(irow) > threshold) {
         const TVectorD uc_i = TMatrixDColumn(fU,irow);
         r = uc_i * vb;
         r /= fSig(irow);
      }
      tmp(irow) = r;
   }

   if (b->GetNrows() > fV.GetNrows()) {
      const TVectorD tmp2 = fV*tmp;
      TVectorD tmp3(cb);
      tmp3.SetSub(tmp2.GetLwb(),tmp2);
      cb = tmp3;
   } else
      cb = fV*tmp;

   return kTRUE;
}

//______________________________________________________________________________
Bool_t TDecompSVD::TransSolve(TVectorD &b)
{
// Solve A^T x=b assuming the SVD form of A is stored . Solution returned in b.

   R__ASSERT(b.IsValid());
   if (TestBit(kSingular)) {
      return kFALSE;
   }
   if ( !TestBit(kDecomposed) ) {
      if (!Decompose()) {
         return kFALSE;
      }
   }

   if (fU.GetNrows() != fV.GetNrows() || fU.GetRowLwb() != fV.GetRowLwb()) {
      Error("TransSolve(TVectorD &","matrix should be square");
      return kFALSE;
   }

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

   // We start with fV fSig fU^T x = b, and turn it into  fU^T x = fSig^-1 fV^T b
   // Form tmp = fSig^-1 fV^T b but ignore diagonal elements in
   // fSig(i) < fTol * max(fSig)

   const Int_t    lwb       = fU.GetColLwb();
   const Int_t    upb       = lwb+fV.GetNcols()-1;
   const Double_t threshold = fSig(lwb)*fTol;

   TVectorD tmp(lwb,upb);
   for (Int_t i = lwb; i <= upb; i++) {
      Double_t r = 0.0;
      if (fSig(i) > threshold) {
         const TVectorD vc = TMatrixDColumn(fV,i);
         r = vc * b;
         r /= fSig(i);
      }
      tmp(i) = r;
   }
   b = fU*tmp;

   return kTRUE;
}

//______________________________________________________________________________
Bool_t TDecompSVD::TransSolve(TMatrixDColumn &cb)
{
// Solve A^T x=b assuming the SVD form of A is stored . Solution returned in b.

   TMatrixDBase *b = const_cast<TMatrixDBase *>(cb.GetMatrix());
   R__ASSERT(b->IsValid());
   if (TestBit(kSingular)) {
      return kFALSE;
   }
   if ( !TestBit(kDecomposed) ) {
      if (!Decompose()) {
         return kFALSE;
      }
   }

   if (fU.GetNrows() != fV.GetNrows() || fU.GetRowLwb() != fV.GetRowLwb()) {
      Error("TransSolve(TMatrixDColumn &","matrix should be square");
      return kFALSE;
   }

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

   // We start with fV fSig fU^T x = b, and turn it into  fU^T x = fSig^-1 fV^T b
   // Form tmp = fSig^-1 fV^T b but ignore diagonal elements in
   // fSig(i) < fTol * max(fSig)

   const Int_t    lwb       = fU.GetColLwb();
   const Int_t    upb       = lwb+fV.GetNcols()-1;
   const Double_t threshold = fSig(lwb)*fTol;

   const TVectorD vb = cb;
   TVectorD tmp(lwb,upb);
   for (Int_t i = lwb; i <= upb; i++) {
      Double_t r = 0.0;
      if (fSig(i) > threshold) {
         const TVectorD vc = TMatrixDColumn(fV,i);
         r = vc * vb;
         r /= fSig(i);
      }
      tmp(i) = r;
   }
   cb = fU*tmp;

   return kTRUE;
}

//______________________________________________________________________________
Double_t TDecompSVD::Condition()
{
// Matrix condition number

   if ( !TestBit(kCondition) ) {
      fCondition = -1;
      if (TestBit(kSingular))
         return fCondition;
      if ( !TestBit(kDecomposed) ) {
         if (!Decompose())
            return fCondition;
      }
      const Int_t colLwb = GetColLwb();
      const Int_t nCols  = GetNcols();
      const Double_t max = fSig(colLwb);
      const Double_t min = fSig(colLwb+nCols-1);
      fCondition = (min > 0.0) ? max/min : -1.0;
      SetBit(kCondition);
   }
   return fCondition;
}

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

   if ( !TestBit(kDetermined) ) {
      if ( !TestBit(kDecomposed) )
         Decompose();
      if (TestBit(kSingular)) {
         fDet1 = 0.0;
         fDet2 = 0.0;
      } else {
         DiagProd(fSig,fTol,fDet1,fDet2);
      }
      SetBit(kDetermined);
   }
   d1 = fDet1;
   d2 = fDet2;
}

//______________________________________________________________________________
Bool_t TDecompSVD::Invert(TMatrixD &inv)
{
// For a matrix A(m,n), its inverse A_inv is defined as A * A_inv = A_inv * A = unit
// The user should always supply a matrix of size (m x m) !
// If m > n , only the (n x m) part of the returned (pseudo inverse) matrix
// should be used .

   const Int_t rowLwb = GetRowLwb();
   const Int_t colLwb = GetColLwb();
   const Int_t nRows  = GetNrows();

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

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

   return status;
}

//______________________________________________________________________________
TMatrixD TDecompSVD::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 colLwb = GetColLwb();
   const Int_t rowUpb = rowLwb+GetNrows()-1;
   TMatrixD inv(rowLwb,rowUpb,colLwb,colLwb+GetNrows()-1);
   inv.UnitMatrix();
   status = MultiSolve(inv);
   inv.ResizeTo(rowLwb,rowLwb+GetNcols()-1,colLwb,colLwb+GetNrows()-1);

   return inv;
}

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

   TDecompBase::Print(opt);
   fU.Print("fU");
   fV.Print("fV");
   fSig.Print("fSig");
}

//______________________________________________________________________________
TDecompSVD &TDecompSVD::operator=(const TDecompSVD &source)
{
// Assignment operator

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