82 const Int_t nrows = row_upb-row_lwb+1;
132 Error(
"Decompose()",
"Matrix has not been set");
177 if (absakk >= alpha*colmax) {
194 if (absakk >= alpha*colmax*(colmax/rowmax)) {
197 }
else if(
TMath::Abs(diag(imax)) >= alpha*rowmax) {
207 const Int_t kk = k-kstep+1;
212 for (
Int_t irow = 0; irow < kp; irow++) {
220 c_kk = pU+(kp+1)*n+kk;
222 for (
Int_t icol = 0; icol < kk-kp-1; icol++) {
235 pU[(k-1)*n+k] = pU[kp*n+k];
242 if (kstep == 1 && k > 0) {
268 const Double_t d22 = pU_k1[k-1]/d12;
273 for (
Int_t j = k-2; j >= 0; j--) {
275 const Double_t wkm1 = d12*(d11*pU_j[k-1]-pU_j[k]);
276 const Double_t wk = d12*(d22*pU_j[k]-pU_j[k-1]);
277 for (
Int_t i = j; i >= 0; i--) {
279 pU_i[j] -= (pU_i[k]*wk+pU_i[k-1]*wkm1);
292 fIpiv[k-1] = -(kp+1);
339 Error(
"Solve()",
"Matrix is singular");
344 Error(
"Solve()",
"Decomposition failed");
350 Error(
"Solve(TVectorD &",
"vector and matrix incompatible");
380 for (
Int_t i = 0; i < k; i++)
381 pb[i] -= pU[i*n+k]*pb[k];
400 for (i = 0; i < k-1; i++)
401 pb[i] -= pU[i*n+k]*pb[k];
402 for (i = 0; i < k-1; i++)
403 pb[i] -= pU[i*n+k-1]*pb[k-1];
408 const Double_t ukm1 = pU_k1[k-1]/ukm1k;
411 const Double_t bkm1 = pb[k-1]/ukm1k;
413 pb[k-1] = (uk*bkm1-bk)/denom;
414 pb[k] = (ukm1*bk-bkm1)/denom;
431 for (
Int_t i = 0; i < k; i++)
432 pb[k] -= pU[i*n+k]*pb[i];
447 for (i = 0; i < k; i++)
448 pb[k] -= pU[i*n+k]*pb[i];
449 for (i = 0; i < k; i++)
450 pb[k+1] -= pU[i*n+k+1]*pb[i];
474 Error(
"Solve()",
"Matrix is singular");
479 Error(
"Solve()",
"Decomposition failed");
485 Error(
"Solve(TMatrixDColumn &",
"vector and matrix incompatible");
511 pcb[k*inc] = pcb[kp*inc];
517 for (
Int_t i = 0; i < k; i++)
518 pcb[i*inc] -= pU[i*n+k]*pcb[k*inc];
521 pcb[k*inc] /= diag(k);
529 const Double_t tmp = pcb[(k-1)*inc];
530 pcb[(k-1)*inc] = pcb[kp*inc];
537 for (i = 0; i < k-1; i++)
538 pcb[i*inc] -= pU[i*n+k]*pcb[k*inc];
539 for (i = 0; i < k-1; i++)
540 pcb[i*inc] -= pU[i*n+k-1]*pcb[(k-1)*inc];
545 const Double_t ukm1 = pU_k1[k-1]/ukm1k;
548 const Double_t bkm1 = pcb[(k-1)*inc]/ukm1k;
549 const Double_t bk = pcb[k*inc]/ukm1k;
550 pcb[(k-1)*inc] = (uk*bkm1-bk)/denom;
551 pcb[k*inc] = (ukm1*bk-bkm1)/denom;
568 for (
Int_t i = 0; i < k; i++)
569 pcb[k*inc] -= pU[i*n+k]*pcb[i*inc];
575 pcb[k*inc] = pcb[kp*inc];
584 for (i = 0; i < k; i++)
585 pcb[k*inc] -= pU[i*n+k]*pcb[i*inc];
586 for (i = 0; i < k; i++)
587 pcb[(k+1)*inc] -= pU[i*n+k+1]*pcb[i*inc];
593 pcb[k*inc] = pcb[kp*inc];
609 Error(
"Invert(TMatrixDSym &",
"Input matrix has wrong shape");
618 for (
Int_t icol = colLwb; icol <= colUpb &&
status; icol++) {
658 if (
this != &source) {
void Rank1Update(const TVectorT< Element > &vec, Element alpha=1.0)
Perform a rank 1 operation on the matrix: A += alpha * v * v^T.
virtual TMatrixTBase< Element > & UnitMatrix()
Make a unit matrix (matrix need not be a square one).
virtual TMatrixTBase< Element > & Shift(Int_t row_shift, Int_t col_shift)
Shift the row index by adding row_shift and the column index by adding col_shift, respectively...
Long64_t LocMax(Long64_t n, const T *a)
virtual Bool_t Solve(TVectorD &b)
Solve Ax=b assuming the BK form of A is stored in fU . Solution returned in b.
void Print(Option_t *opt="") const
Print class members.
double inv(double x)
For comparisons.
TMatrixTColumn_const< Double_t > TMatrixDColumn_const
ClassImp(TDecompBK) TDecompBK
Default constructor.
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.
TDecompBase & operator=(const TDecompBase &source)
Assignment operator.
TMatrixTRow_const< Double_t > TMatrixDRow_const
void Print(Option_t *opt="") const
Print the class members.
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
virtual const Element * GetMatrixArray() const
Element * GetMatrixArray()
void Print(Option_t *name="") const
Print the matrix as a table of elements.
TDecompBK & operator=(const TDecompBK &source)
Assigment operator.
unsigned int r1[N_CITIES]
virtual Int_t GetNrows() const
Bool_t TestBit(UInt_t f) const
virtual const Element * GetMatrixArray() const
const TMatrixTBase< Element > * GetMatrix() const
virtual void SetMatrix(const TMatrixDSym &a)
Set the matrix to be decomposed, decomposition status is reset.
ClassImp(TMCParticle) void TMCParticle printf(": p=(%7.3f,%7.3f,%9.3f) ;", fPx, fPy, fPz)
Short_t Max(Short_t a, Short_t b)
TVectorT< Element > & Abs()
Take an absolute value of a vector, i.e. apply Abs() to each element.
Double_t Sqrt(Double_t x)
virtual Bool_t Decompose()
Matrix A is decomposed in components U and D so that A = U*D*U^T If the decomposition succeeds...