Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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 = nullptr;
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
121{
122 fNIpiv = 0;
123 fIpiv = nullptr;
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 ) {
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);
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
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
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
649{
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) {
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}
#define b(i)
Definition RSha256.hxx:100
#define a(i)
Definition RSha256.hxx:99
constexpr Bool_t kFALSE
Definition RtypesCore.h:101
constexpr Bool_t kTRUE
Definition RtypesCore.h:100
const char Option_t
Definition RtypesCore.h:66
#define ClassImp(name)
Definition Rtypes.h:377
#define R__ASSERT(e)
Definition TError.h:118
TMatrixTColumn_const< Double_t > TMatrixDColumn_const
TMatrixTRow_const< Double_t > TMatrixDRow_const
The Bunch-Kaufman diagonal pivoting method decomposes a real symmetric matrix A using.
Definition TDecompBK.h:27
Int_t GetNrows() const override
Definition TDecompBK.h:45
TDecompBK & operator=(const TDecompBK &source)
Assignment operator.
TDecompBK()
Default constructor.
Definition TDecompBK.cxx:64
TMatrixDSym Invert()
Definition TDecompBK.h:64
Bool_t Decompose() override
Matrix A is decomposed in components U and D so that A = U*D*U^T If the decomposition succeeds,...
Int_t fNIpiv
Definition TDecompBK.h:30
Bool_t Solve(TVectorD &b) override
Solve Ax=b assuming the BK form of A is stored in fU . Solution returned in b.
TMatrixD fU
Definition TDecompBK.h:32
virtual void SetMatrix(const TMatrixDSym &a)
Set the matrix to be decomposed, decomposition status is reset.
void Print(Option_t *opt="") const override
Print the class members.
Int_t * fIpiv
Definition TDecompBK.h:31
Decomposition Base class.
Definition TDecompBase.h:34
Int_t GetRowLwb() const
Definition TDecompBase.h:73
void ResetStatus()
Definition TDecompBase.h:43
Int_t fRowLwb
Definition TDecompBase.h:40
Double_t fTol
Definition TDecompBase.h:36
TDecompBase & operator=(const TDecompBase &source)
Assignment operator.
Double_t fCondition
Definition TDecompBase.h:39
Int_t fColLwb
Definition TDecompBase.h:41
void Print(Option_t *opt="") const override
Print class members.
Int_t GetNrows() const
void Print(Option_t *name="") const override
Print the matrix as a table of elements.
Int_t GetRowLwb() const
Int_t GetNcols() const
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.
const TMatrixTBase< Element > * GetMatrix() const
Element * GetPtr() const
void Rank1Update(const TVectorT< Element > &vec, Element alpha=1.0)
Perform a rank 1 operation on the matrix: A += alpha * v * v^T.
const Element * GetMatrixArray() const override
Definition TMatrixT.h:225
TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t=-1) override
Set size of the matrix to nrows x ncols New dynamic elements are created, the overlapping part of the...
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
Definition TObject.h:201
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition TObject.cxx:780
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:976
TVectorT< Element > & Abs()
Take an absolute value of a vector, i.e. apply Abs() to each element.
Definition TVectorT.cxx:464
Element * GetMatrixArray()
Definition TVectorT.h:76
const Int_t n
Definition legend1.C:16
Short_t Max(Short_t a, Short_t b)
Returns the largest of a and b.
Definition TMathBase.h:250
Long64_t LocMax(Long64_t n, const T *a)
Returns index of array with the maximum element.
Definition TMath.h:1012
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:662
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:123
void inv(rsa_NUMBER *, rsa_NUMBER *, rsa_NUMBER *)
Definition rsaaux.cxx:949