Logo ROOT   6.14/05
Reference Guide
TDecompBK.cxx
Go to the documentation of this file.
1 // @(#)root/matrix:$Id$
2 // Authors: Fons Rademakers, Eddy Offermann Sep 2004
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 #include "TDecompBK.h"
13 #include "TMath.h"
14 #include "TError.h"
15 
17 
18 /** \class TDecompBK
19  \ingroup Matrix
20 
21  The Bunch-Kaufman diagonal pivoting method decomposes a real
22  symmetric matrix A using
23 
24 ~~~
25  A = U*D*U^T
26 ~~~
27 
28  where U is a product of permutation and unit upper triangular
29  matrices, U^T is the transpose of U, and D is symmetric and block
30  diagonal with 1-by-1 and 2-by-2 diagonal blocks.
31 
32  U = P(n-1)*U(n-1)* ... *P(k)U(k)* ...,
33  i.e., U is a product of terms P(k)*U(k), where k decreases from n-1
34  to 0 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
35  and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
36  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
37  that if the diagonal block D(k) is of order s (s = 1 or 2), then
38 
39 ~~~
40  ( I v 0 ) k-s
41  U(k) = ( 0 I 0 ) s
42  ( 0 0 I ) n-k
43  k-s s n-k
44 ~~~
45 
46  If s = 1, D(k) overwrites A(k,k), and v overwrites A(0:k-1,k).
47  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
48  and A(k,k), and v overwrites A(0:k-2,k-1:k).
49 
50  fU contains on entry the symmetric matrix A of which only the upper
51  triangular part is referenced . On exit fU contains the block diagonal
52  matrix D and the multipliers used to obtain the factor U, see above .
53 
54  fIpiv if dimension n contains details of the interchanges and the
55  the block structure of D . If (fIPiv(k) > 0, then rows and columns k
56  and fIPiv(k) were interchanged and D(k,k) is a 1-by-1 diagonal block.
57  If IPiv(k) = fIPiv(k-1) < 0, rows and columns k-1 and -IPiv(k) were
58  interchanged and D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
59 */
60 
61 ////////////////////////////////////////////////////////////////////////////////
62 /// Default constructor
63 
65 {
66  fNIpiv = 0;
67  fIpiv = 0;
68 }
69 
70 ////////////////////////////////////////////////////////////////////////////////
71 /// Constructor for (nrows x nrows) symmetric matrix
72 
74 {
75  fNIpiv = nrows;
76  fIpiv = new Int_t[fNIpiv];
77  memset(fIpiv,0,fNIpiv*sizeof(Int_t));
78  fU.ResizeTo(nrows,nrows);
79 }
80 
81 ////////////////////////////////////////////////////////////////////////////////
82 /// Constructor for ([row_lwb..row_upb] x [row_lwb..row_upb]) symmetric matrix
83 
85 {
86  const Int_t nrows = row_upb-row_lwb+1;
87  fNIpiv = nrows;
88  fIpiv = new Int_t[fNIpiv];
89  memset(fIpiv,0,fNIpiv*sizeof(Int_t));
90  fColLwb = fRowLwb = row_lwb;
91  fU.ResizeTo(nrows,nrows);
92 }
93 
94 ////////////////////////////////////////////////////////////////////////////////
95 /// Constructor for symmetric matrix A
96 
98 {
99  R__ASSERT(a.IsValid());
100 
102  fCondition = a.Norm1();
103  fTol = a.GetTol();
104  if (tol > 0)
105  fTol = tol;
106 
107  fNIpiv = a.GetNcols();
108  fIpiv = new Int_t[fNIpiv];
109  memset(fIpiv,0,fNIpiv*sizeof(Int_t));
110 
111  const Int_t nRows = a.GetNrows();
112  fColLwb = fRowLwb = a.GetRowLwb();
113  fU.ResizeTo(nRows,nRows);
114  memcpy(fU.GetMatrixArray(),a.GetMatrixArray(),nRows*nRows*sizeof(Double_t));
115 }
116 
117 ////////////////////////////////////////////////////////////////////////////////
118 /// Copy constructor
119 
120 TDecompBK::TDecompBK(const TDecompBK &another) : TDecompBase(another)
121 {
122  fNIpiv = 0;
123  fIpiv = 0;
124  *this = another;
125 }
126 
127 ////////////////////////////////////////////////////////////////////////////////
128 /// Matrix A is decomposed in components U and D so that A = U*D*U^T
129 /// If the decomposition succeeds, bit kDecomposed is set , otherwise kSingular
130 
132 {
133  if (TestBit(kDecomposed)) return kTRUE;
134 
135  if ( !TestBit(kMatrixSet) ) {
136  Error("Decompose()","Matrix has not been set");
137  return kFALSE;
138  }
139 
140  Bool_t ok = kTRUE;
141 
142 // Initialize alpha for use in choosing pivot block size.
143  const Double_t alpha = (1.+TMath::Sqrt(17.))/8.;
144 
145 // Factorize a as u*d*u' using the upper triangle of a .
146 // k is the main loop index, decreasing from n-1 to 0 in steps of 1 or 2
147 
148  const Int_t n = fU.GetNcols();
149  Double_t *pU = fU.GetMatrixArray();
150  TMatrixDDiag diag(fU);
151 
152  Int_t imax = 0;
153  Int_t k = n-1;
154  while (k >= 0) {
155  Int_t kstep = 1;
156 
157  // determine rows and columns to be interchanged and whether
158  // a 1-by-1 or 2-by-2 pivot block will be used
159 
160  const Double_t absakk = TMath::Abs(diag(k));
161 
162  // imax is the row-index of the largest off-diagonal element in
163  // column k, and colmax is its absolute value
164 
165  Double_t colmax;
166  if ( k > 0 ) {
167  TVectorD vcol = TMatrixDColumn_const(fU,k);
168  vcol.Abs();
169  imax = TMath::LocMax(k,vcol.GetMatrixArray());
170  colmax = vcol[imax];
171  } else {
172  colmax = 0.;
173  }
174 
175  Int_t kp;
176  if (TMath::Max(absakk,colmax) <= fTol) {
177  // the block diagonal matrix will be singular
178  kp = k;
179  ok = kFALSE;
180  } else {
181  if (absakk >= alpha*colmax) {
182  // no interchange, use 1-by-1 pivot block
183  kp = k;
184  } else {
185  // jmax is the column-index of the largest off-diagonal
186  // element in row imax, and rowmax is its absolute value
187  TVectorD vrow = TMatrixDRow_const(fU,imax);
188  vrow.Abs();
189  Int_t jmax = imax+1+TMath::LocMax(k-imax,vrow.GetMatrixArray()+imax+1);
190  Double_t rowmax = vrow[jmax];
191  if (imax > 0) {
192  TVectorD vcol = TMatrixDColumn_const(fU,imax);
193  vcol.Abs();
194  jmax = TMath::LocMax(imax,vcol.GetMatrixArray());
195  rowmax = TMath::Max(rowmax,vcol[jmax]);
196  }
197 
198  if (absakk >= alpha*colmax*(colmax/rowmax)) {
199  // No interchange, use 1-by-1 pivot block
200  kp = k;
201  } else if( TMath::Abs(diag(imax)) >= alpha*rowmax) {
202  // Interchange rows and columns k and imax, use 1-by-1 pivot block
203  kp = imax;
204  } else {
205  // Interchange rows and columns k-1 and imax, use 2-by-2 pivot block
206  kp = imax;
207  kstep = 2;
208  }
209  }
210 
211  const Int_t kk = k-kstep+1;
212  if (kp != kk) {
213  // Interchange rows and columns kk and kp in the leading submatrix a(0:k,0:k)
214  Double_t *c_kk = pU+kk;
215  Double_t *c_kp = pU+kp;
216  for (Int_t irow = 0; irow < kp; irow++) {
217  const Double_t t = *c_kk;
218  *c_kk = *c_kp;
219  *c_kp = t;
220  c_kk += n;
221  c_kp += n;
222  }
223 
224  c_kk = pU+(kp+1)*n+kk;
225  Double_t *r_kp = pU+kp*n+kp+1;
226  for (Int_t icol = 0; icol < kk-kp-1; icol++) {
227  const Double_t t = *c_kk;
228  *c_kk = *r_kp;
229  *r_kp = t;
230  c_kk += n;
231  r_kp += 1;
232  }
233 
234  Double_t t = diag(kk);
235  diag(kk) = diag(kp);
236  diag(kp) = t;
237  if (kstep == 2) {
238  t = pU[(k-1)*n+k];
239  pU[(k-1)*n+k] = pU[kp*n+k];
240  pU[kp*n+k] = t;
241  }
242  }
243 
244  // Update the leading submatrix
245 
246  if (kstep == 1 && k > 0) {
247  // 1-by-1 pivot block d(k): column k now holds w(k) = u(k)*d(k)
248  // where u(k) is the k-th column of u
249 
250  // perform a rank-1 update of a(0:k-1,0:k-1) as
251  // a := a - u(k)*d(k)*u(k)' = a - w(k)*1/d(k)*w(k)'
252 
253  const Double_t r1 = 1./diag(k);
254  TMatrixDSub sub1(fU,0,k-1,0,k-1);
255  sub1.Rank1Update(TMatrixDColumn_const(fU,k),-r1);
256 
257  // store u(k) in column k
258  TMatrixDSub sub2(fU,0,k-1,k,k);
259  sub2 *= r1;
260  } else {
261  // 2-by-2 pivot block d(k): columns k and k-1 now hold
262  // ( w(k-1) w(k) ) = ( u(k-1) u(k) )*d(k)
263  // where u(k) and u(k-1) are the k-th and (k-1)-th columns of u
264 
265  // perform a rank-2 update of a(0:k-2,0:k-2) as
266  // a := a - ( u(k-1) u(k) )*d(k)*( u(k-1) u(k) )'
267  // = a - ( w(k-1) w(k) )*inv(d(k))*( w(k-1) w(k) )'
268 
269  if ( k > 1 ) {
270  Double_t *pU_k1 = pU+(k-1)*n;
271  Double_t d12 = pU_k1[k];
272  const Double_t d22 = pU_k1[k-1]/d12;
273  const Double_t d11 = diag(k)/d12;
274  const Double_t t = 1./(d11*d22-1.);
275  d12 = t/d12;
276 
277  for (Int_t j = k-2; j >= 0; j--) {
278  Double_t *pU_j = pU+j*n;
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--) {
282  Double_t *pU_i = pU+i*n;
283  pU_i[j] -= (pU_i[k]*wk+pU_i[k-1]*wkm1);
284  }
285  pU_j[k] = wk;
286  pU_j[k-1] = wkm1;
287  }
288  }
289  }
290 
291  // Store details of the interchanges in fIpiv
292  if (kstep == 1) {
293  fIpiv[k] = (kp+1);
294  } else {
295  fIpiv[k] = -(kp+1);
296  fIpiv[k-1] = -(kp+1);
297  }
298  }
299 
300  k -= kstep;
301  }
302 
303  if (!ok) SetBit(kSingular);
304  else SetBit(kDecomposed);
305 
307 
308  return ok;
309 }
310 
311 ////////////////////////////////////////////////////////////////////////////////
312 /// Set the matrix to be decomposed, decomposition status is reset.
313 
315 {
316  R__ASSERT(a.IsValid());
317 
318  ResetStatus();
319 
321  fCondition = a.Norm1();
322 
323  if (fNIpiv != a.GetNcols()) {
324  fNIpiv = a.GetNcols();
325  delete [] fIpiv;
326  fIpiv = new Int_t[fNIpiv];
327  memset(fIpiv,0,fNIpiv*sizeof(Int_t));
328  }
329 
330  const Int_t nRows = a.GetNrows();
331  fColLwb = fRowLwb = a.GetRowLwb();
332  fU.ResizeTo(nRows,nRows);
333  memcpy(fU.GetMatrixArray(),a.GetMatrixArray(),nRows*nRows*sizeof(Double_t));
334 }
335 
336 ////////////////////////////////////////////////////////////////////////////////
337 /// Solve Ax=b assuming the BK form of A is stored in fU . Solution returned in b.
338 
340 {
341  R__ASSERT(b.IsValid());
342  if (TestBit(kSingular)) {
343  Error("Solve()","Matrix is singular");
344  return kFALSE;
345  }
346  if ( !TestBit(kDecomposed) ) {
347  if (!Decompose()) {
348  Error("Solve()","Decomposition failed");
349  return kFALSE;
350  }
351  }
352 
353  if (fU.GetNrows() != b.GetNrows() || fU.GetRowLwb() != b.GetLwb()) {
354  Error("Solve(TVectorD &","vector and matrix incompatible");
355  return kFALSE;
356  }
357 
358  const Int_t n = fU.GetNrows();
359 
360  TMatrixDDiag_const diag(fU);
361  const Double_t *pU = fU.GetMatrixArray();
362  Double_t *pb = b.GetMatrixArray();
363 
364  // solve a*x = b, where a = u*d*u'. First solve u*d*x = b, overwriting b with x.
365  // k is the main loop index, decreasing from n-1 to 0 in steps of 1 or 2,
366  // depending on the size of the diagonal blocks.
367 
368  Int_t k = n-1;
369  while (k >= 0) {
370 
371  if (fIpiv[k] > 0) {
372 
373  // 1 x 1 diagonal block
374  // interchange rows k and ipiv(k).
375  const Int_t kp = fIpiv[k]-1;
376  if (kp != k) {
377  const Double_t tmp = pb[k];
378  pb[k] = pb[kp];
379  pb[kp] = tmp;
380  }
381 
382  // multiply by inv(u(k)), where u(k) is the transformation
383  // stored in column k of a.
384  for (Int_t i = 0; i < k; i++)
385  pb[i] -= pU[i*n+k]*pb[k];
386 
387  // multiply by the inverse of the diagonal block.
388  pb[k] /= diag(k);
389  k--;
390  } else {
391 
392  // 2 x 2 diagonal block
393  // interchange rows k-1 and -ipiv(k).
394  const Int_t kp = -fIpiv[k]-1;
395  if (kp != k-1) {
396  const Double_t tmp = pb[k-1];
397  pb[k-1] = pb[kp];
398  pb[kp] = tmp;
399  }
400 
401  // multiply by inv(u(k)), where u(k) is the transformation
402  // stored in columns k-1 and k of a.
403  Int_t i;
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];
408 
409  // multiply by the inverse of the diagonal block.
410  const Double_t *pU_k1 = pU+(k-1)*n;
411  const Double_t ukm1k = pU_k1[k];
412  const Double_t ukm1 = pU_k1[k-1]/ukm1k;
413  const Double_t uk = diag(k)/ukm1k;
414  const Double_t denom = ukm1*uk-1.;
415  const Double_t bkm1 = pb[k-1]/ukm1k;
416  const Double_t bk = pb[k]/ukm1k;
417  pb[k-1] = (uk*bkm1-bk)/denom;
418  pb[k] = (ukm1*bk-bkm1)/denom;
419  k -= 2;
420  }
421  }
422 
423  // Next solve u'*x = b, overwriting b with x.
424  //
425  // k is the main loop index, increasing from 0 to n-1 in steps of
426  // 1 or 2, depending on the size of the diagonal blocks.
427 
428  k = 0;
429  while (k < n) {
430 
431  if (fIpiv[k] > 0) {
432  // 1 x 1 diagonal block
433  // multiply by inv(u'(k)), where u(k) is the transformation
434  // stored in column k of a.
435  for (Int_t i = 0; i < k; i++)
436  pb[k] -= pU[i*n+k]*pb[i];
437 
438  // interchange elements k and ipiv(k).
439  const Int_t kp = fIpiv[k]-1;
440  if (kp != k) {
441  const Double_t tmp = pb[k];
442  pb[k] = pb[kp];
443  pb[kp] = tmp;
444  }
445  k++;
446  } else {
447  // 2 x 2 diagonal block
448  // multiply by inv(u'(k+1)), where u(k+1) is the transformation
449  // stored in columns k and k+1 of a.
450  Int_t 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];
455 
456  // interchange elements k and -ipiv(k).
457  const Int_t kp = -fIpiv[k]-1;
458  if (kp != k) {
459  const Double_t tmp = pb[k];
460  pb[k] = pb[kp];
461  pb[kp] = tmp;
462  }
463  k += 2;
464  }
465  }
466 
467  return kTRUE;
468 }
469 
470 ////////////////////////////////////////////////////////////////////////////////
471 /// Solve Ax=b assuming the BK form of A is stored in fU . Solution returned in b.
472 
474 {
475  TMatrixDBase *b = const_cast<TMatrixDBase *>(cb.GetMatrix());
476  R__ASSERT(b->IsValid());
477  if (TestBit(kSingular)) {
478  Error("Solve()","Matrix is singular");
479  return kFALSE;
480  }
481  if ( !TestBit(kDecomposed) ) {
482  if (!Decompose()) {
483  Error("Solve()","Decomposition failed");
484  return kFALSE;
485  }
486  }
487 
488  if (fU.GetNrows() != b->GetNrows() || fU.GetRowLwb() != b->GetRowLwb()) {
489  Error("Solve(TMatrixDColumn &","vector and matrix incompatible");
490  return kFALSE;
491  }
492 
493  const Int_t n = fU.GetNrows();
494 
495  TMatrixDDiag_const diag(fU);
496  const Double_t *pU = fU.GetMatrixArray();
497  Double_t *pcb = cb.GetPtr();
498  const Int_t inc = cb.GetInc();
499 
500 
501  // solve a*x = b, where a = u*d*u'. First solve u*d*x = b, overwriting b with x.
502  // k is the main loop index, decreasing from n-1 to 0 in steps of 1 or 2,
503  // depending on the size of the diagonal blocks.
504 
505  Int_t k = n-1;
506  while (k >= 0) {
507 
508  if (fIpiv[k] > 0) {
509 
510  // 1 x 1 diagonal block
511  // interchange rows k and ipiv(k).
512  const Int_t kp = fIpiv[k]-1;
513  if (kp != k) {
514  const Double_t tmp = pcb[k*inc];
515  pcb[k*inc] = pcb[kp*inc];
516  pcb[kp*inc] = tmp;
517  }
518 
519  // multiply by inv(u(k)), where u(k) is the transformation
520  // stored in column k of a.
521  for (Int_t i = 0; i < k; i++)
522  pcb[i*inc] -= pU[i*n+k]*pcb[k*inc];
523 
524  // multiply by the inverse of the diagonal block.
525  pcb[k*inc] /= diag(k);
526  k--;
527  } else {
528 
529  // 2 x 2 diagonal block
530  // interchange rows k-1 and -ipiv(k).
531  const Int_t kp = -fIpiv[k]-1;
532  if (kp != k-1) {
533  const Double_t tmp = pcb[(k-1)*inc];
534  pcb[(k-1)*inc] = pcb[kp*inc];
535  pcb[kp*inc] = tmp;
536  }
537 
538  // multiply by inv(u(k)), where u(k) is the transformation
539  // stored in columns k-1 and k of a.
540  Int_t i;
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];
545 
546  // multiply by the inverse of the diagonal block.
547  const Double_t *pU_k1 = pU+(k-1)*n;
548  const Double_t ukm1k = pU_k1[k];
549  const Double_t ukm1 = pU_k1[k-1]/ukm1k;
550  const Double_t uk = diag(k)/ukm1k;
551  const Double_t denom = ukm1*uk-1.;
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;
556  k -= 2;
557  }
558  }
559 
560  // Next solve u'*x = b, overwriting b with x.
561  //
562  // k is the main loop index, increasing from 0 to n-1 in steps of
563  // 1 or 2, depending on the size of the diagonal blocks.
564 
565  k = 0;
566  while (k < n) {
567 
568  if (fIpiv[k] > 0) {
569  // 1 x 1 diagonal block
570  // multiply by inv(u'(k)), where u(k) is the transformation
571  // stored in column k of a.
572  for (Int_t i = 0; i < k; i++)
573  pcb[k*inc] -= pU[i*n+k]*pcb[i*inc];
574 
575  // interchange elements k and ipiv(k).
576  const Int_t kp = fIpiv[k]-1;
577  if (kp != k) {
578  const Double_t tmp = pcb[k*inc];
579  pcb[k*inc] = pcb[kp*inc];
580  pcb[kp*inc] = tmp;
581  }
582  k++;
583  } else {
584  // 2 x 2 diagonal block
585  // multiply by inv(u'(k+1)), where u(k+1) is the transformation
586  // stored in columns k and k+1 of a.
587  Int_t i;
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];
592 
593  // interchange elements k and -ipiv(k).
594  const Int_t kp = -fIpiv[k]-1;
595  if (kp != k) {
596  const Double_t tmp = pcb[k*inc];
597  pcb[k*inc] = pcb[kp*inc];
598  pcb[kp*inc] = tmp;
599  }
600  k += 2;
601  }
602  }
603 
604  return kTRUE;
605 }
606 
607 ////////////////////////////////////////////////////////////////////////////////
608 /// For a symmetric matrix A(m,m), its inverse A_inv(m,m) is returned .
609 
611 {
612  if (inv.GetNrows() != GetNrows() || inv.GetRowLwb() != GetRowLwb()) {
613  Error("Invert(TMatrixDSym &","Input matrix has wrong shape");
614  return kFALSE;
615  }
616 
617  inv.UnitMatrix();
618 
619  const Int_t colLwb = inv.GetColLwb();
620  const Int_t colUpb = inv.GetColUpb();
621  Bool_t status = kTRUE;
622  for (Int_t icol = colLwb; icol <= colUpb && status; icol++) {
623  TMatrixDColumn b(inv,icol);
624  status &= Solve(b);
625  }
626 
627  return status;
628 }
629 
630 ////////////////////////////////////////////////////////////////////////////////
631 /// For a symmetric matrix A(m,m), its inverse A_inv(m,m) is returned .
632 
634 {
635  const Int_t rowLwb = GetRowLwb();
636  const Int_t rowUpb = rowLwb+GetNrows()-1;
637 
638  TMatrixDSym inv(rowLwb,rowUpb);
639  inv.UnitMatrix();
640  status = Invert(inv);
641 
642  return inv;
643 }
644 
645 ////////////////////////////////////////////////////////////////////////////////
646 /// Print the class members
647 
648 void TDecompBK::Print(Option_t *opt) const
649 {
650  TDecompBase::Print(opt);
651  printf("fIpiv:\n");
652  for (Int_t i = 0; i < fNIpiv; i++)
653  printf("[%d] = %d\n",i,fIpiv[i]);
654  fU.Print("fU");
655 }
656 
657 ////////////////////////////////////////////////////////////////////////////////
658 /// Assignment operator
659 
661 {
662  if (this != &source) {
663  TDecompBase::operator=(source);
664  fU.ResizeTo(source.fU);
665  fU = source.fU;
666  if (fNIpiv != source.fNIpiv) {
667  if (fIpiv)
668  delete [] fIpiv;
669  fNIpiv = source.fNIpiv;
670  fIpiv = new Int_t[fNIpiv];
671  }
672  if (fIpiv) memcpy(fIpiv,source.fIpiv,fNIpiv*sizeof(Int_t));
673  }
674  return *this;
675 }
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...
Int_t GetRowLwb() const
Definition: TMatrixTBase.h:120
Long64_t LocMax(Long64_t n, const T *a)
Return index of array with the maximum element.
Definition: TMath.h:999
Int_t GetInc() const
Double_t fCondition
Definition: TDecompBase.h:39
Int_t GetRowLwb() const
Definition: TDecompBase.h:73
Int_t fRowLwb
Definition: TDecompBase.h:40
const char Option_t
Definition: RtypesCore.h:62
virtual const Element * GetMatrixArray() const
Definition: TMatrixT.h:222
Int_t fNIpiv
Definition: TDecompBK.h:30
virtual Bool_t Solve(TVectorD &b)
Solve Ax=b assuming the BK form of A is stored in fU . Solution returned in b.
Definition: TDecompBK.cxx:339
Int_t GetNcols() const
Definition: TMatrixTBase.h:125
Int_t GetLwb() const
Definition: TVectorT.h:73
Decomposition Base class.
Definition: TDecompBase.h:33
#define R__ASSERT(e)
Definition: TError.h:96
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
Definition: TObject.h:172
Int_t GetNrows() const
Definition: TVectorT.h:75
TMatrixTColumn_const< Double_t > TMatrixDColumn_const
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
virtual const Element * GetMatrixArray() const
Definition: TMatrixTSym.h:189
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...
Definition: TMatrixT.cxx:1210
Short_t Abs(Short_t d)
Definition: TMathBase.h:108
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition: TObject.cxx:694
Double_t fTol
Definition: TDecompBase.h:36
The Bunch-Kaufman diagonal pivoting method decomposes a real symmetric matrix A using.
Definition: TDecompBK.h:26
Int_t GetColLwb() const
Definition: TMatrixTBase.h:123
Element * GetPtr() const
TDecompBase & operator=(const TDecompBase &source)
Assignment operator.
TMatrixTRow_const< Double_t > TMatrixDRow_const
Int_t fColLwb
Definition: TDecompBase.h:41
Element Norm1() const
Definition: TMatrixTBase.h:176
Element * GetMatrixArray()
Definition: TVectorT.h:78
Element GetTol() const
Definition: TMatrixTBase.h:127
TDecompBK & operator=(const TDecompBK &source)
Assignment operator.
Definition: TDecompBK.cxx:660
void inv(rsa_NUMBER *, rsa_NUMBER *, rsa_NUMBER *)
Definition: rsaaux.cxx:949
auto * a
Definition: textangle.C:12
TMatrixDSym Invert()
Definition: TDecompBK.h:64
Bool_t IsValid() const
Definition: TVectorT.h:83
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:880
void Print(Option_t *opt="") const
Print class members.
void ResetStatus()
Definition: TDecompBase.h:43
Int_t GetColUpb() const
Definition: TMatrixTBase.h:124
Linear Algebra Package.
Int_t GetNrows() const
Definition: TMatrixTBase.h:122
const Bool_t kFALSE
Definition: RtypesCore.h:88
TDecompBK()
Default constructor.
Definition: TDecompBK.cxx:64
TMatrixD fU
Definition: TDecompBK.h:32
#define ClassImp(name)
Definition: Rtypes.h:359
double Double_t
Definition: RtypesCore.h:55
virtual void SetMatrix(const TMatrixDSym &a)
Set the matrix to be decomposed, decomposition status is reset.
Definition: TDecompBK.cxx:314
virtual Int_t GetNrows() const
Definition: TDecompBK.h:45
Int_t * fIpiv
Definition: TDecompBK.h:31
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)
Definition: TMathBase.h:200
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
Definition: TRolke.cxx:630
TVectorT< Element > & Abs()
Take an absolute value of a vector, i.e. apply Abs() to each element.
Definition: TVectorT.cxx:462
Double_t Sqrt(Double_t x)
Definition: TMath.h:690
virtual Bool_t Decompose()
Matrix A is decomposed in components U and D so that A = U*D*U^T If the decomposition succeeds...
Definition: TDecompBK.cxx:131
const Bool_t kTRUE
Definition: RtypesCore.h:87
const Int_t n
Definition: legend1.C:16
void Print(Option_t *opt="") const
Print the class members.
Definition: TDecompBK.cxx:648
Bool_t IsValid() const
Definition: TMatrixTBase.h:145