```// @(#)root/matrix:\$Id\$
// Authors: Fons Rademakers, Eddy Offermann  Dec 2003

/*************************************************************************
*                                                                       *
* 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) {
for (Int_t i = 1; i < nCol_v; i++)
oDiag(i) = v(i-1,i);
}
oDiag(0) = 0.;
sDiag = TMatrixDDiag(v);

// 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;
}

//______________________________________________________________________________
Int_t  TDecompSVD::GetNrows  () const
{
return fU.GetNrows();
}

Int_t TDecompSVD::GetNcols  () const
{
return fV.GetNcols();
}

//______________________________________________________________________________
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  = fU.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+fU.GetNrows()-1;
TMatrixD inv(rowLwb,rowUpb,colLwb,colLwb+fU.GetNrows()-1);
inv.UnitMatrix();
status = MultiSolve(inv);
inv.ResizeTo(rowLwb,rowLwb+fV.GetNcols()-1,colLwb,colLwb+fU.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
TDecompSVD.cxx:925
TDecompSVD.cxx:926
TDecompSVD.cxx:927
TDecompSVD.cxx:928
TDecompSVD.cxx:929
TDecompSVD.cxx:930
TDecompSVD.cxx:931
TDecompSVD.cxx:932
TDecompSVD.cxx:933
TDecompSVD.cxx:934
TDecompSVD.cxx:935