45 const Int_t nrows = row_upb-row_lwb+1;
48 fU.
ResizeTo(row_lwb,row_lwb+nrows-1,row_lwb,row_lwb+nrows-1);
78 Error(
"TDecompChol(const TMatrixD &",
"matrix should be square");
111 Error(
"Decompose()",
"Matrix has not been set");
118 for (icol = 0; icol <
n; icol++) {
119 const Int_t rowOff = icol*
n;
123 for (irow = 0; irow < icol; irow++) {
124 const Int_t pos_ij = irow*n+icol;
125 ujj -= pU[pos_ij]*pU[pos_ij];
128 Error(
"Decompose()",
"matrix not positive definite");
132 pU[rowOff+icol] = ujj;
135 for (j = icol+1; j <
n; j++) {
136 for (i = 0; i < icol; i++) {
137 const Int_t rowOff2 = i*
n;
138 pU[rowOff+j] -= pU[rowOff2+j]*pU[rowOff2+icol];
141 for (j = icol+1; j <
n; j++)
146 for (irow = 0; irow <
n; irow++) {
147 const Int_t rowOff = irow*
n;
148 for (icol = 0; icol < irow; icol++)
149 pU[rowOff+icol] = 0.;
163 Error(
"GetMatrix()",
"Matrix is singular");
168 Error(
"GetMatrix()",
"Decomposition failed");
185 Error(
"SetMatrix(const TMatrixDSym &",
"matrix should be square");
207 Error(
"Solve()",
"Matrix is singular");
212 Error(
"Solve()",
"Decomposition failed");
218 Error(
"Solve(TVectorD &",
"vector and matrix incompatible");
229 for (i = 0; i <
n; i++) {
231 if (pU[off_i+i] <
fTol)
233 Error(
"Solve(TVectorD &b)",
"u[%d,%d]=%.4e < %.4e",i,i,pU[off_i+i],
fTol);
237 for (
Int_t j = 0; j < i; j++) {
239 r -= pU[off_j+i]*pb[j];
241 pb[i] = r/pU[off_i+i];
245 for (i = n-1; i >= 0; i--) {
248 for (
Int_t j = i+1; j <
n; j++)
249 r -= pU[off_i+j]*pb[j];
250 pb[i] = r/pU[off_i+i];
266 Error(
"Solve()",
"Matrix is singular");
271 Error(
"Solve()",
"Decomposition failed");
278 Error(
"Solve(TMatrixDColumn &cb",
"vector and matrix incompatible");
290 for (i = 0; i <
n; i++) {
292 const Int_t off_i2 = i*inc;
293 if (pU[off_i+i] <
fTol)
295 Error(
"Solve(TMatrixDColumn &cb)",
"u[%d,%d]=%.4e < %.4e",i,i,pU[off_i+i],
fTol);
299 for (
Int_t j = 0; j < i; j++) {
301 r -= pU[off_j+i]*pcb[j*inc];
303 pcb[off_i2] = r/pU[off_i+i];
307 for (i = n-1; i >= 0; i--) {
309 const Int_t off_i2 = i*inc;
311 for (
Int_t j = i+1; j <
n; j++)
312 r -= pU[off_i+j]*pcb[j*inc];
313 pcb[off_i2] = r/pU[off_i+i];
344 Error(
"Invert(TMatrixDSym &",
"Input matrix has wrong shape");
353 for (
Int_t icol = colLwb; icol <= colUpb && status; icol++) {
390 if (
this != &source) {
422 ::Error(
"NormalEqn",
"vectors b and std are not compatible");
430 mBw(irow) /= std(irow);
virtual TMatrixTBase< Element > & UnitMatrix()
Make a unit matrix (matrix need not be a square one).
virtual void Det(Double_t &d1, Double_t &d2)
Matrix determinant det = d1*TMath::Power(2.,d2) is square of diagProd of cholesky factor...
void Print(Option_t *opt="") const
Print class members .
virtual Bool_t MultiSolve(TMatrixD &B)
Solve set of equations with RHS in columns of B.
virtual const Element * GetMatrixArray() const
Decomposition Base class.
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
virtual TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t=-1)
Set size of the matrix to nrows x ncols New dynamic elements are created, the overlapping part of the...
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
const TMatrixDSym GetMatrix()
Reconstruct the original matrix using the decomposition parts.
TDecompBase & operator=(const TDecompBase &source)
Assignment operator.
virtual void SetMatrix(const TMatrixDSym &a)
Set the matrix to be decomposed, decomposition status is reset.
Bool_t AreCompatible(const TMatrixTBase< Element1 > &m1, const TMatrixTBase< Element2 > &m2, Int_t verbose=0)
Check that matrice sm1 and m2 areboth valid and have identical shapes .
TVectorT< Double_t > TVectorD
virtual Int_t GetNrows() const
TMatrixTRow< Double_t > TMatrixDRow
Element * GetMatrixArray()
TMatrixT< Double_t > TMatrixD
TVectorD NormalEqn(const TMatrixD &A, const TVectorD &b)
Solve min {(A .
Cholesky Decomposition class.
virtual Bool_t Solve(TVectorD &b)
Solve equations Ax=b assuming A has been factored by Cholesky.
void inv(rsa_NUMBER *, rsa_NUMBER *, rsa_NUMBER *)
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
void Print(Option_t *opt="") const
Print class members.
TMatrixTSym< Double_t > TMatrixDSym
virtual Bool_t Decompose()
Matrix A is decomposed in component U so that A = U^T * U If the decomposition succeeds, bit kDecomposed is set , otherwise kSingular.
virtual void Det(Double_t &d1, Double_t &d2)
Matrix determinant det = d1*TMath::Power(2.,d2)
void Print(Option_t *name="") const
Print the matrix as a table of elements.
const TMatrixTBase< Element > * GetMatrix() const
TDecompChol & operator=(const TDecompChol &source)
Assignment operator.
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Double_t Sqrt(Double_t x)