86 const Int_t nrows = row_upb-row_lwb+1;
136 Error(
"Decompose()",
"Matrix has not been set");
181 if (absakk >= alpha*colmax) {
198 if (absakk >= alpha*colmax*(colmax/rowmax)) {
201 }
else if(
TMath::Abs(diag(imax)) >= alpha*rowmax) {
211 const Int_t kk = k-kstep+1;
216 for (
Int_t irow = 0; irow < kp; irow++) {
224 c_kk = pU+(kp+1)*n+kk;
226 for (
Int_t icol = 0; icol < kk-kp-1; icol++) {
239 pU[(k-1)*n+k] = pU[kp*n+k];
246 if (kstep == 1 && k > 0) {
272 const Double_t d22 = pU_k1[k-1]/d12;
277 for (
Int_t j = k-2; j >= 0; j--) {
279 const Double_t wkm1 = d12*(d11*pU_j[k-1]-pU_j[k]);
280 const Double_t wk = d12*(d22*pU_j[k]-pU_j[k-1]);
281 for (
Int_t i = j; i >= 0; i--) {
283 pU_i[j] -= (pU_i[k]*wk+pU_i[k-1]*wkm1);
296 fIpiv[k-1] = -(kp+1);
343 Error(
"Solve()",
"Matrix is singular");
348 Error(
"Solve()",
"Decomposition failed");
354 Error(
"Solve(TVectorD &",
"vector and matrix incompatible");
384 for (
Int_t i = 0; i < k; i++)
385 pb[i] -= pU[i*n+k]*pb[k];
404 for (i = 0; i < k-1; i++)
405 pb[i] -= pU[i*n+k]*pb[k];
406 for (i = 0; i < k-1; i++)
407 pb[i] -= pU[i*n+k-1]*pb[k-1];
412 const Double_t ukm1 = pU_k1[k-1]/ukm1k;
415 const Double_t bkm1 = pb[k-1]/ukm1k;
417 pb[k-1] = (uk*bkm1-bk)/denom;
418 pb[k] = (ukm1*bk-bkm1)/denom;
435 for (
Int_t i = 0; i < k; i++)
436 pb[k] -= pU[i*n+k]*pb[i];
451 for (i = 0; i < k; i++)
452 pb[k] -= pU[i*n+k]*pb[i];
453 for (i = 0; i < k; i++)
454 pb[k+1] -= pU[i*n+k+1]*pb[i];
478 Error(
"Solve()",
"Matrix is singular");
483 Error(
"Solve()",
"Decomposition failed");
489 Error(
"Solve(TMatrixDColumn &",
"vector and matrix incompatible");
515 pcb[k*inc] = pcb[kp*inc];
521 for (
Int_t i = 0; i < k; i++)
522 pcb[i*inc] -= pU[i*n+k]*pcb[k*inc];
525 pcb[k*inc] /= diag(k);
533 const Double_t tmp = pcb[(k-1)*inc];
534 pcb[(k-1)*inc] = pcb[kp*inc];
541 for (i = 0; i < k-1; i++)
542 pcb[i*inc] -= pU[i*n+k]*pcb[k*inc];
543 for (i = 0; i < k-1; i++)
544 pcb[i*inc] -= pU[i*n+k-1]*pcb[(k-1)*inc];
549 const Double_t ukm1 = pU_k1[k-1]/ukm1k;
552 const Double_t bkm1 = pcb[(k-1)*inc]/ukm1k;
553 const Double_t bk = pcb[k*inc]/ukm1k;
554 pcb[(k-1)*inc] = (uk*bkm1-bk)/denom;
555 pcb[k*inc] = (ukm1*bk-bkm1)/denom;
572 for (
Int_t i = 0; i < k; i++)
573 pcb[k*inc] -= pU[i*n+k]*pcb[i*inc];
579 pcb[k*inc] = pcb[kp*inc];
588 for (i = 0; i < k; i++)
589 pcb[k*inc] -= pU[i*n+k]*pcb[i*inc];
590 for (i = 0; i < k; i++)
591 pcb[(k+1)*inc] -= pU[i*n+k+1]*pcb[i*inc];
597 pcb[k*inc] = pcb[kp*inc];
613 Error(
"Invert(TMatrixDSym &",
"Input matrix has wrong shape");
622 for (
Int_t icol = colLwb; icol <= colUpb && status; icol++) {
653 printf(
"[%d] = %d\n",i,
fIpiv[i]);
662 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)
Return index of array with the maximum element.
virtual const Element * GetMatrixArray() const
virtual Bool_t Solve(TVectorD &b)
Solve Ax=b assuming the BK form of A is stored in fU . Solution returned in b.
Decomposition Base class.
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
TMatrixTColumn_const< Double_t > TMatrixDColumn_const
virtual const Element * GetMatrixArray() 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.
The Bunch-Kaufman diagonal pivoting method decomposes a real symmetric matrix A using.
TDecompBase & operator=(const TDecompBase &source)
Assignment operator.
TMatrixTRow_const< Double_t > TMatrixDRow_const
Element * GetMatrixArray()
TDecompBK & operator=(const TDecompBK &source)
Assignment operator.
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.
TDecompBK()
Default constructor.
virtual void SetMatrix(const TMatrixDSym &a)
Set the matrix to be decomposed, decomposition status is reset.
virtual Int_t GetNrows() const
void Print(Option_t *name="") const
Print the matrix as a table of elements.
const TMatrixTBase< Element > * GetMatrix() const
Short_t Max(Short_t a, Short_t b)
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
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...
void Print(Option_t *opt="") const
Print the class members.