Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
TMatrixTSym.cxx
Go to the documentation of this file.
1// @(#)root/matrix:$Id$
2// Authors: Fons Rademakers, Eddy Offermann Nov 2003
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/** \class TMatrixTSym
13 \ingroup Matrix
14
15TMatrixTSym
16
17Template class of a symmetric matrix in the linear algebra package.
18
19See the \ref Matrix page for the documentation of the linear algebra package
20
21Note that in this implementation both matrix element m[i][j] and
22m[j][i] are updated and stored in memory. However, when making the
23object persistent only the upper right triangle is stored.
24
25*/
26
27#include "TMatrixTSym.h"
28#include "TBuffer.h"
29#include "TMatrixTLazy.h"
31#include "TDecompLU.h"
32#include "TMatrixDSymEigen.h"
33#include "TMath.h"
34
35
36
37////////////////////////////////////////////////////////////////////////////////
38
39template<class Element>
45////////////////////////////////////////////////////////////////////////////////
46
47template<class Element>
53
54////////////////////////////////////////////////////////////////////////////////
55/// option=
56/// - "F": array elements contains the matrix stored column-wise
57/// like in Fortran, so a[i,j] = elements[i+no_rows*j],
58/// - else it is supposed that array elements are stored row-wise
59/// a[i,j] = elements[i*no_cols+j]
60///
61/// array elements are copied
62
63template<class Element>
65{
66 Allocate(no_rows,no_rows);
67 SetMatrixArray(elements,option);
68}
70////////////////////////////////////////////////////////////////////////////////
71/// array elements are copied
73template<class Element>
75{
76 const Int_t no_rows = row_upb-row_lwb+1;
78 SetMatrixArray(elements,option);
79}
80
81////////////////////////////////////////////////////////////////////////////////
83template<class Element>
85{
86 R__ASSERT(another.IsValid());
87 Allocate(another.GetNrows(),another.GetNcols(),another.GetRowLwb(),another.GetColLwb());
88 *this = another;
89}
90
91////////////////////////////////////////////////////////////////////////////////
92/// Create a matrix applying a specific operation to the prototype.
93/// Example: TMatrixTSym<Element> a(10,12); ...; TMatrixTSym<Element> b(TMatrixT::kTransposed, a);
94/// Supported operations are: kZero, kUnit, kTransposed, kInverted and kAtA.
95
96template<class Element>
98{
99 R__ASSERT(prototype.IsValid());
100
101 switch(op) {
102 case kZero:
103 Allocate(prototype.GetNrows(),prototype.GetNcols(),
104 prototype.GetRowLwb(),prototype.GetColLwb(),1);
105 break;
106
107 case kUnit:
108 Allocate(prototype.GetNrows(),prototype.GetNcols(),
109 prototype.GetRowLwb(),prototype.GetColLwb(),1);
110 this->UnitMatrix();
111 break;
112
113 case kTransposed:
114 Allocate(prototype.GetNcols(), prototype.GetNrows(),
115 prototype.GetColLwb(),prototype.GetRowLwb());
116 Transpose(prototype);
117 break;
118
119 case kInverted:
121 Allocate(prototype.GetNrows(),prototype.GetNcols(),
122 prototype.GetRowLwb(),prototype.GetColLwb(),1);
123 *this = prototype;
124 // Since the user can not control the tolerance of this newly created matrix
125 // we put it to the smallest possible number
126 const Element oldTol = this->SetTol(std::numeric_limits<Element>::min());
127 this->Invert();
128 this->SetTol(oldTol);
129 break;
131
132 case kAtA:
133 Allocate(prototype.GetNcols(),prototype.GetNcols(),prototype.GetColLwb(),prototype.GetColLwb(),1);
134 TMult(prototype);
135 break;
137 default:
138 Error("TMatrixTSym(EMatrixCreatorOp1,const TMatrixTSym)",
139 "operation %d not yet implemented", op);
140 }
141}
142
143////////////////////////////////////////////////////////////////////////////////
144
145template<class Element>
148 R__ASSERT(prototype.IsValid());
149
150 switch(op) {
151 case kAtA:
152 Allocate(prototype.GetNcols(),prototype.GetNcols(),prototype.GetColLwb(),prototype.GetColLwb(),1);
153 TMult(prototype);
154 break;
155
156 default:
157 Error("TMatrixTSym(EMatrixCreatorOp1,const TMatrixT)",
158 "operation %d not yet implemented", op);
159 }
160}
161
162////////////////////////////////////////////////////////////////////////////////
163
164template<class Element>
167 R__ASSERT(a.IsValid());
168 R__ASSERT(b.IsValid());
170 switch(op) {
171 case kPlus:
173 Allocate(a.GetNcols(),a.GetNcols(),a.GetColLwb(),a.GetColLwb(),1);
174 Plus(a,b);
175 break;
177
178 case kMinus:
179 {
180 Allocate(a.GetNcols(),a.GetNcols(),a.GetColLwb(),a.GetColLwb(),1);
181 Minus(a,b);
182 break;
183 }
184
185 default:
186 Error("TMatrixTSym(EMatrixCreatorOp2)", "operation %d not yet implemented", op);
187 }
188}
189
190////////////////////////////////////////////////////////////////////////////////
191
192template<class Element>
194{
195 Allocate(lazy_constructor.GetRowUpb()-lazy_constructor.GetRowLwb()+1,
196 lazy_constructor.GetRowUpb()-lazy_constructor.GetRowLwb()+1,
197 lazy_constructor.GetRowLwb(),lazy_constructor.GetRowLwb(),1);
198 lazy_constructor.FillIn(*this);
200 Error("TMatrixTSym(TMatrixTSymLazy)","matrix not symmetric");
201 }
202}
203
204////////////////////////////////////////////////////////////////////////////////
205/// delete data pointer m, if it was assigned on the heap
206
207template<class Element>
209{
210 if (m) {
211 if (size > this->kSizeMax)
212 delete [] m;
213 m = nullptr;
214 }
215}
216
217////////////////////////////////////////////////////////////////////////////////
218/// return data pointer . if requested size <= kSizeMax, assign pointer
219/// to the stack space
220
221template<class Element>
223{
224 if (size == 0) return nullptr;
225 else {
226 if ( size <= this->kSizeMax )
227 return fDataStack;
228 else {
229 Element *heap = new Element[size];
230 return heap;
231 }
232 }
233}
234
235////////////////////////////////////////////////////////////////////////////////
236/// copy copySize doubles from *oldp to *newp . However take care of the
237/// situation where both pointers are assigned to the same stack space
238
239template<class Element>
242{
243 if (copySize == 0 || oldp == newp)
244 return 0;
245 else {
246 if ( newSize <= this->kSizeMax && oldSize <= this->kSizeMax ) {
247 // both pointers are inside fDataStack, be careful with copy direction !
248 if (newp > oldp) {
249 for (Int_t i = copySize-1; i >= 0; i--)
250 newp[i] = oldp[i];
251 } else {
252 for (Int_t i = 0; i < copySize; i++)
253 newp[i] = oldp[i];
254 }
255 }
256 else
257 memcpy(newp,oldp,copySize*sizeof(Element));
258 }
259 return 0;
260}
261
262////////////////////////////////////////////////////////////////////////////////
263/// Allocate new matrix. Arguments are number of rows, columns, row
264/// lowerbound (0 default) and column lowerbound (0 default).
265
266template<class Element>
268 Int_t init,Int_t /*nr_nonzeros*/)
269{
270 this->fIsOwner = kTRUE;
271 this->fTol = std::numeric_limits<Element>::epsilon();
272 this->fNrows = 0;
273 this->fNcols = 0;
274 this->fRowLwb = 0;
275 this->fColLwb = 0;
276 this->fNelems = 0;
277 fElements = nullptr;
278
279 if (no_rows < 0 || no_cols < 0)
280 {
281 Error("Allocate","no_rows=%d no_cols=%d",no_rows,no_cols);
282 this->Invalidate();
283 return;
284 }
285
286 this->MakeValid();
287 this->fNrows = no_rows;
288 this->fNcols = no_cols;
289 this->fRowLwb = row_lwb;
290 this->fColLwb = col_lwb;
291 this->fNelems = this->fNrows*this->fNcols;
292
293 if (this->fNelems > 0) {
294 fElements = New_m(this->fNelems);
295 if (init)
296 memset(fElements,0,this->fNelems*sizeof(Element));
297 } else
298 fElements = nullptr;
299}
300
301////////////////////////////////////////////////////////////////////////////////
302/// Symmetric matrix summation. Replace this matrix with C such that C = A + B.
303
304template<class Element>
306{
307 if (gMatrixCheck) {
308 if (!AreCompatible(a,b)) {
309 Error("Plus","matrices not compatible");
310 return;
311 }
312
313 if (this->GetMatrixArray() == a.GetMatrixArray()) {
314 Error("Plus","this->GetMatrixArray() == a.GetMatrixArray()");
315 return;
316 }
317
318 if (this->GetMatrixArray() == b.GetMatrixArray()) {
319 Error("Plus","this->GetMatrixArray() == b.GetMatrixArray()");
320 return;
321 }
322 }
323
324 const Element * ap = a.GetMatrixArray();
325 const Element * bp = b.GetMatrixArray();
326 Element * cp = this->GetMatrixArray();
327 const Element * const cp_last = cp+this->fNelems;
328
329 while (cp < cp_last) {
330 *cp = *ap++ + *bp++;
331 cp++;
332 }
333}
334
335////////////////////////////////////////////////////////////////////////////////
336/// Symmetric matrix subtraction. Replace this matrix with C such that C = A - B.
337
338template<class Element>
340{
341 if (gMatrixCheck) {
342 if (!AreCompatible(a,b)) {
343 Error("Minus","matrices not compatible");
344 return;
345 }
346
347 if (this->GetMatrixArray() == a.GetMatrixArray()) {
348 Error("Minus","this->GetMatrixArray() == a.GetMatrixArray()");
349 return;
350 }
351
352 if (this->GetMatrixArray() == b.GetMatrixArray()) {
353 Error("Minus","this->GetMatrixArray() == b.GetMatrixArray()");
354 return;
355 }
356 }
357
358 const Element * ap = a.GetMatrixArray();
359 const Element * bp = b.GetMatrixArray();
360 Element * cp = this->GetMatrixArray();
361 const Element * const cp_last = cp+this->fNelems;
362
363 while (cp < cp_last) {
364 *cp = *ap++ - *bp++;
365 cp++;
366 }
367}
368
369////////////////////////////////////////////////////////////////////////////////
370/// Replace this matrix with C such that C = A' * A. In other words,
371/// c[i,j] = SUM{ a[k,i] * a[k,j] }.
372
373template<class Element>
375{
376 R__ASSERT(a.IsValid());
377
378#ifdef CBLAS
379 const Element *ap = a.GetMatrixArray();
380 Element *cp = this->GetMatrixArray();
381 if (typeid(Element) == typeid(Double_t))
382 cblas_dgemm (CblasRowMajor,CblasTrans,CblasNoTrans,this->fNrows,this->fNcols,a.GetNrows(),
383 1.0,ap,a.GetNcols(),ap,a.GetNcols(),1.0,cp,this->fNcols);
384 else if (typeid(Element) != typeid(Float_t))
385 cblas_sgemm (CblasRowMajor,CblasTrans,CblasNoTrans,fNrows,fNcols,a.GetNrows(),
386 1.0,ap,a.GetNcols(),ap,a.GetNcols(),1.0,cp,fNcols);
387 else
388 Error("TMult","type %s not implemented in BLAS library",typeid(Element));
389#else
390 const Int_t nb = a.GetNoElements();
391 const Int_t ncolsa = a.GetNcols();
392 const Int_t ncolsb = ncolsa;
393 const Element * const ap = a.GetMatrixArray();
394 const Element * const bp = ap;
395 Element * cp = this->GetMatrixArray();
396
397 const Element *acp0 = ap; // Pointer to A[i,0];
398 while (acp0 < ap+a.GetNcols()) {
399 for (const Element *bcp = bp; bcp < bp+ncolsb; ) { // Pointer to the j-th column of A, Start bcp = A[0,0]
400 const Element *acp = acp0; // Pointer to the i-th column of A, reset to A[0,i]
401 Element cij = 0;
402 while (bcp < bp+nb) { // Scan the i-th column of A and
403 cij += *acp * *bcp; // the j-th col of A
404 acp += ncolsa;
405 bcp += ncolsb;
406 }
407 *cp++ = cij;
408 bcp -= nb-1; // Set bcp to the (j+1)-th col
409 }
410 acp0++; // Set acp0 to the (i+1)-th col
411 }
412
413 R__ASSERT(cp == this->GetMatrixArray()+this->fNelems && acp0 == ap+ncolsa);
414#endif
415}
416
417////////////////////////////////////////////////////////////////////////////////
418/// Matrix multiplication, with A symmetric
419/// Replace this matrix with C such that C = A' * A = A * A = A * A'
420
421template<class Element>
423{
424 R__ASSERT(a.IsValid());
425
426#ifdef CBLAS
427 const Element *ap = a.GetMatrixArray();
428 Element *cp = this->GetMatrixArray();
429 if (typeid(Element) == typeid(Double_t))
430 cblas_dsymm (CblasRowMajor,CblasLeft,CblasUpper,this->fNrows,this->fNcols,1.0,
431 ap,a.GetNcols(),ap,a.GetNcols(),0.0,cp,this->fNcols);
432 else if (typeid(Element) != typeid(Float_t))
434 ap1,a.GetNcols(),ap,a.GetNcols(),0.0,cp,fNcols);
435 else
436 Error("TMult","type %s not implemented in BLAS library",typeid(Element));
437#else
438 const Int_t nb = a.GetNoElements();
439 const Int_t ncolsa = a.GetNcols();
440 const Int_t ncolsb = ncolsa;
441 const Element * const ap = a.GetMatrixArray();
442 const Element * const bp = ap;
443 Element * cp = this->GetMatrixArray();
444
445 const Element *acp0 = ap; // Pointer to A[i,0];
446 while (acp0 < ap+a.GetNcols()) {
447 for (const Element *bcp = bp; bcp < bp+ncolsb; ) { // Pointer to the j-th column of A, Start bcp = A[0,0]
448 const Element *acp = acp0; // Pointer to the i-th column of A, reset to A[0,i]
449 Element cij = 0;
450 while (bcp < bp+nb) { // Scan the i-th column of A and
451 cij += *acp * *bcp; // the j-th col of A
452 acp += ncolsa;
453 bcp += ncolsb;
454 }
455 *cp++ = cij;
456 bcp -= nb-1; // Set bcp to the (j+1)-th col
457 }
458 acp0++; // Set acp0 to the (i+1)-th col
459 }
460
461 R__ASSERT(cp == this->GetMatrixArray()+this->fNelems && acp0 == ap+ncolsa);
462#endif
463}
464
465////////////////////////////////////////////////////////////////////////////////
466
467template<class Element>
469{
471 {
472 Error("Use","row_upb=%d < row_lwb=%d",row_upb,row_lwb);
473 return *this;
474 }
475
476 this->Clear();
477 this->fNrows = row_upb-row_lwb+1;
478 this->fNcols = this->fNrows;
479 this->fRowLwb = row_lwb;
480 this->fColLwb = row_lwb;
481 this->fNelems = this->fNrows*this->fNcols;
482 fElements = data;
483 this->fIsOwner = kFALSE;
484
485 return *this;
486}
487
488////////////////////////////////////////////////////////////////////////////////
489/// Get submatrix [row_lwb..row_upb][row_lwb..row_upb]; The indexing range of the
490/// returned matrix depends on the argument option:
491///
492/// option == "S" : return [0..row_upb-row_lwb+1][0..row_upb-row_lwb+1] (default)
493/// else : return [row_lwb..row_upb][row_lwb..row_upb]
494
495template<class Element>
497 TMatrixTSym<Element> &target,Option_t *option) const
498{
499 if (gMatrixCheck) {
500 R__ASSERT(this->IsValid());
501 if (row_lwb < this->fRowLwb || row_lwb > this->fRowLwb+this->fNrows-1) {
502 Error("GetSub","row_lwb out of bounds");
503 return target;
504 }
505 if (row_upb < this->fRowLwb || row_upb > this->fRowLwb+this->fNrows-1) {
506 Error("GetSub","row_upb out of bounds");
507 return target;
508 }
509 if (row_upb < row_lwb) {
510 Error("GetSub","row_upb < row_lwb");
511 return target;
512 }
513 }
514
515 TString opt(option);
516 opt.ToUpper();
517 const Int_t shift = (opt.Contains("S")) ? 1 : 0;
518
521 if (shift) {
522 row_lwb_sub = 0;
524 } else {
527 }
528
531
532 if (target.GetRowIndexArray() && target.GetColIndexArray()) {
533 for (Int_t irow = 0; irow < nrows_sub; irow++) {
534 for (Int_t icol = 0; icol < nrows_sub; icol++) {
536 }
537 }
538 } else {
539 const Element *ap = this->GetMatrixArray()+(row_lwb-this->fRowLwb)*this->fNrows+(row_lwb-this->fRowLwb);
540 Element *bp = target.GetMatrixArray();
541
542 for (Int_t irow = 0; irow < nrows_sub; irow++) {
543 const Element *ap_sub = ap;
544 for (Int_t icol = 0; icol < nrows_sub; icol++) {
545 *bp++ = *ap_sub++;
546 }
547 ap += this->fNrows;
548 }
549 }
550
551 return target;
552}
553
554////////////////////////////////////////////////////////////////////////////////
555/// Get submatrix [row_lwb..row_upb][col_lwb..col_upb]; The indexing range of the
556/// returned matrix depends on the argument option:
557///
558/// option == "S" : return [0..row_upb-row_lwb+1][0..col_upb-col_lwb+1] (default)
559/// else : return [row_lwb..row_upb][col_lwb..col_upb]
560
561template<class Element>
564{
565 if (gMatrixCheck) {
566 R__ASSERT(this->IsValid());
567 if (row_lwb < this->fRowLwb || row_lwb > this->fRowLwb+this->fNrows-1) {
568 Error("GetSub","row_lwb out of bounds");
569 return target;
570 }
571 if (col_lwb < this->fColLwb || col_lwb > this->fColLwb+this->fNcols-1) {
572 Error("GetSub","col_lwb out of bounds");
573 return target;
574 }
575 if (row_upb < this->fRowLwb || row_upb > this->fRowLwb+this->fNrows-1) {
576 Error("GetSub","row_upb out of bounds");
577 return target;
578 }
579 if (col_upb < this->fColLwb || col_upb > this->fColLwb+this->fNcols-1) {
580 Error("GetSub","col_upb out of bounds");
581 return target;
582 }
583 if (row_upb < row_lwb || col_upb < col_lwb) {
584 Error("GetSub","row_upb < row_lwb || col_upb < col_lwb");
585 return target;
586 }
587 }
588
589 TString opt(option);
590 opt.ToUpper();
591 const Int_t shift = (opt.Contains("S")) ? 1 : 0;
592
593 const Int_t row_lwb_sub = (shift) ? 0 : row_lwb;
594 const Int_t row_upb_sub = (shift) ? row_upb-row_lwb : row_upb;
595 const Int_t col_lwb_sub = (shift) ? 0 : col_lwb;
596 const Int_t col_upb_sub = (shift) ? col_upb-col_lwb : col_upb;
597
601
602 if (target.GetRowIndexArray() && target.GetColIndexArray()) {
603 for (Int_t irow = 0; irow < nrows_sub; irow++) {
604 for (Int_t icol = 0; icol < ncols_sub; icol++) {
606 }
607 }
608 } else {
609 const Element *ap = this->GetMatrixArray()+(row_lwb-this->fRowLwb)*this->fNcols+(col_lwb-this->fColLwb);
610 Element *bp = target.GetMatrixArray();
611
612 for (Int_t irow = 0; irow < nrows_sub; irow++) {
613 const Element *ap_sub = ap;
614 for (Int_t icol = 0; icol < ncols_sub; icol++) {
615 *bp++ = *ap_sub++;
616 }
617 ap += this->fNcols;
618 }
619 }
620
621 return target;
622}
623
624////////////////////////////////////////////////////////////////////////////////
625/// Insert matrix source starting at [row_lwb][row_lwb], thereby overwriting the part
626/// [row_lwb..row_lwb+nrows_source][row_lwb..row_lwb+nrows_source];
627
628template<class Element>
630{
631 if (gMatrixCheck) {
632 R__ASSERT(this->IsValid());
633 R__ASSERT(source.IsValid());
634
635 if (!source.IsSymmetric()) {
636 Error("SetSub","source matrix is not symmetric");
637 return *this;
638 }
639 if (row_lwb < this->fRowLwb || row_lwb > this->fRowLwb+this->fNrows-1) {
640 Error("SetSub","row_lwb outof bounds");
641 return *this;
642 }
643 if (row_lwb+source.GetNrows() > this->fRowLwb+this->fNrows) {
644 Error("SetSub","source matrix too large");
645 return *this;
646 }
647 }
648
649 const Int_t nRows_source = source.GetNrows();
650
651 if (source.GetRowIndexArray() && source.GetColIndexArray()) {
652 const Int_t rowlwb_s = source.GetRowLwb();
653 for (Int_t irow = 0; irow < nRows_source; irow++) {
654 for (Int_t icol = 0; icol < nRows_source; icol++) {
656 }
657 }
658 } else {
659 const Element *bp = source.GetMatrixArray();
660 Element *ap = this->GetMatrixArray()+(row_lwb-this->fRowLwb)*this->fNrows+(row_lwb-this->fRowLwb);
661
662 for (Int_t irow = 0; irow < nRows_source; irow++) {
663 Element *ap_sub = ap;
664 for (Int_t icol = 0; icol < nRows_source; icol++) {
665 *ap_sub++ = *bp++;
666 }
667 ap += this->fNrows;
668 }
669 }
670
671 return *this;
672}
673
674////////////////////////////////////////////////////////////////////////////////
675/// Insert matrix source starting at [row_lwb][col_lwb] in a symmetric fashion, thereby overwriting the part
676/// [row_lwb..row_lwb+nrows_source][row_lwb..row_lwb+nrows_source];
677
678template<class Element>
680{
681 if (gMatrixCheck) {
682 R__ASSERT(this->IsValid());
683 R__ASSERT(source.IsValid());
684
685 if (row_lwb < this->fRowLwb || row_lwb > this->fRowLwb+this->fNrows-1) {
686 Error("SetSub","row_lwb out of bounds");
687 return *this;
688 }
689 if (col_lwb < this->fColLwb || col_lwb > this->fColLwb+this->fNcols-1) {
690 Error("SetSub","col_lwb out of bounds");
691 return *this;
692 }
693
694 if (row_lwb+source.GetNrows() > this->fRowLwb+this->fNrows || col_lwb+source.GetNcols() > this->fRowLwb+this->fNrows) {
695 Error("SetSub","source matrix too large");
696 return *this;
697 }
698 if (col_lwb+source.GetNcols() > this->fRowLwb+this->fNrows || row_lwb+source.GetNrows() > this->fRowLwb+this->fNrows) {
699 Error("SetSub","source matrix too large");
700 return *this;
701 }
702 }
703
704 const Int_t nRows_source = source.GetNrows();
705 const Int_t nCols_source = source.GetNcols();
706
707 const Int_t rowlwb_s = source.GetRowLwb();
708 const Int_t collwb_s = source.GetColLwb();
709 if (row_lwb >= col_lwb) {
710 // lower triangle
711 Int_t irow;
712 for (irow = 0; irow < nRows_source; irow++) {
713 for (Int_t icol = 0; col_lwb+icol <= row_lwb+irow &&
714 icol < nCols_source; icol++) {
716 }
717 }
718
719 // upper triangle
720 for (irow = 0; irow < nCols_source; irow++) {
722 icol >= 0; icol--) {
724 }
725 }
726 } else {
727
728 }
729
730 return *this;
731}
732
733////////////////////////////////////////////////////////////////////////////////
734
735template<class Element>
737{
740 Error("SetMatrixArray","Matrix is not symmetric after Set");
741 }
742
743 return *this;
744}
745
746////////////////////////////////////////////////////////////////////////////////
747
748template<class Element>
750{
751 if (row_shift != col_shift) {
752 Error("Shift","row_shift != col_shift");
753 return *this;
754 }
756}
757
758////////////////////////////////////////////////////////////////////////////////
759/// Set size of the matrix to nrows x ncols
760/// New dynamic elements are created, the overlapping part of the old ones are
761/// copied to the new structures, then the old elements are deleted.
762
763template<class Element>
765{
766 R__ASSERT(this->IsValid());
767 if (!this->fIsOwner) {
768 Error("ResizeTo(Int_t,Int_t)","Not owner of data array,cannot resize");
769 return *this;
770 }
771
772 if (nrows != ncols) {
773 Error("ResizeTo(Int_t,Int_t)","nrows != ncols");
774 return *this;
775 }
776
777 if (this->fNelems > 0) {
778 if (this->fNrows == nrows && this->fNcols == ncols)
779 return *this;
780 else if (nrows == 0 || ncols == 0) {
781 this->fNrows = nrows; this->fNcols = ncols;
782 this->Clear();
783 return *this;
784 }
785
786 Element *elements_old = GetMatrixArray();
787 const Int_t nelems_old = this->fNelems;
788 const Int_t nrows_old = this->fNrows;
789 const Int_t ncols_old = this->fNcols;
790
791 Allocate(nrows,ncols);
792 R__ASSERT(this->IsValid());
793
794 Element *elements_new = GetMatrixArray();
795 // new memory should be initialized but be careful not to wipe out the stack
796 // storage. Initialize all when old or new storage was on the heap
797 if (this->fNelems > this->kSizeMax || nelems_old > this->kSizeMax)
798 memset(elements_new,0,this->fNelems*sizeof(Element));
799 else if (this->fNelems > nelems_old)
800 memset(elements_new+nelems_old,0,(this->fNelems-nelems_old)*sizeof(Element));
801
802 // Copy overlap
803 const Int_t ncols_copy = TMath::Min(this->fNcols,ncols_old);
804 const Int_t nrows_copy = TMath::Min(this->fNrows,nrows_old);
805
806 const Int_t nelems_new = this->fNelems;
807 if (ncols_old < this->fNcols) {
808 for (Int_t i = nrows_copy-1; i >= 0; i--) {
809 Memcpy_m(elements_new+i*this->fNcols,elements_old+i*ncols_old,ncols_copy,
811 if (this->fNelems <= this->kSizeMax && nelems_old <= this->kSizeMax)
812 memset(elements_new+i*this->fNcols+ncols_copy,0,(this->fNcols-ncols_copy)*sizeof(Element));
813 }
814 } else {
815 for (Int_t i = 0; i < nrows_copy; i++)
816 Memcpy_m(elements_new+i*this->fNcols,elements_old+i*ncols_old,ncols_copy,
818 }
819
820 Delete_m(nelems_old,elements_old);
821 } else {
822 Allocate(nrows,ncols,0,0,1);
823 }
824
825 return *this;
826}
827
828////////////////////////////////////////////////////////////////////////////////
829/// Set size of the matrix to [row_lwb:row_upb] x [col_lwb:col_upb]
830/// New dynamic elements are created, the overlapping part of the old ones are
831/// copied to the new structures, then the old elements are deleted.
832
833template<class Element>
835 Int_t /*nr_nonzeros*/)
836{
837 R__ASSERT(this->IsValid());
838 if (!this->fIsOwner) {
839 Error("ResizeTo(Int_t,Int_t,Int_t,Int_t)","Not owner of data array,cannot resize");
840 return *this;
841 }
842
843 if (row_lwb != col_lwb) {
844 Error("ResizeTo(Int_t,Int_t,Int_t,Int_t)","row_lwb != col_lwb");
845 return *this;
846 }
847 if (row_upb != col_upb) {
848 Error("ResizeTo(Int_t,Int_t,Int_t,Int_t)","row_upb != col_upb");
849 return *this;
850 }
851
852 const Int_t new_nrows = row_upb-row_lwb+1;
853 const Int_t new_ncols = col_upb-col_lwb+1;
854
855 if (this->fNelems > 0) {
856
857 if (this->fNrows == new_nrows && this->fNcols == new_ncols &&
858 this->fRowLwb == row_lwb && this->fColLwb == col_lwb)
859 return *this;
860 else if (new_nrows == 0 || new_ncols == 0) {
861 this->fNrows = new_nrows; this->fNcols = new_ncols;
862 this->fRowLwb = row_lwb; this->fColLwb = col_lwb;
863 this->Clear();
864 return *this;
865 }
866
867 Element *elements_old = GetMatrixArray();
868 const Int_t nelems_old = this->fNelems;
869 const Int_t nrows_old = this->fNrows;
870 const Int_t ncols_old = this->fNcols;
871 const Int_t rowLwb_old = this->fRowLwb;
872 const Int_t colLwb_old = this->fColLwb;
873
875 R__ASSERT(this->IsValid());
876
877 Element *elements_new = GetMatrixArray();
878 // new memory should be initialized but be careful not to wipe out the stack
879 // storage. Initialize all when old or new storage was on the heap
880 if (this->fNelems > this->kSizeMax || nelems_old > this->kSizeMax)
881 memset(elements_new,0,this->fNelems*sizeof(Element));
882 else if (this->fNelems > nelems_old)
883 memset(elements_new+nelems_old,0,(this->fNelems-nelems_old)*sizeof(Element));
884
885 // Copy overlap
886 const Int_t rowLwb_copy = TMath::Max(this->fRowLwb,rowLwb_old);
887 const Int_t colLwb_copy = TMath::Max(this->fColLwb,colLwb_old);
888 const Int_t rowUpb_copy = TMath::Min(this->fRowLwb+this->fNrows-1,rowLwb_old+nrows_old-1);
889 const Int_t colUpb_copy = TMath::Min(this->fColLwb+this->fNcols-1,colLwb_old+ncols_old-1);
890
893
894 if (nrows_copy > 0 && ncols_copy > 0) {
896 const Int_t colNewOff = colLwb_copy-this->fColLwb;
897 if (ncols_old < this->fNcols) {
898 for (Int_t i = nrows_copy-1; i >= 0; i--) {
900 const Int_t iRowNew = rowLwb_copy+i-this->fRowLwb;
901 Memcpy_m(elements_new+iRowNew*this->fNcols+colNewOff,
903 if (this->fNelems <= this->kSizeMax && nelems_old <= this->kSizeMax)
905 (this->fNcols-ncols_copy)*sizeof(Element));
906 }
907 } else {
908 for (Int_t i = 0; i < nrows_copy; i++) {
910 const Int_t iRowNew = rowLwb_copy+i-this->fRowLwb;
911 Memcpy_m(elements_new+iRowNew*this->fNcols+colNewOff,
913 }
914 }
915 }
916
917 Delete_m(nelems_old,elements_old);
918 } else {
920 }
921
922 return *this;
923}
924
925////////////////////////////////////////////////////////////////////////////////
926
927template<class Element>
929{
930 const TMatrixT<Element> &tmp = *this;
931 TDecompLU lu(tmp,this->fTol);
932 Double_t d1,d2;
933 lu.Det(d1,d2);
934 return d1*TMath::Power(2.0,d2);
935}
936
937////////////////////////////////////////////////////////////////////////////////
938
939template<class Element>
941{
942 const TMatrixT<Element> &tmp = *this;
943 TDecompLU lu(tmp,this->fTol);
944 lu.Det(d1,d2);
945}
946
947////////////////////////////////////////////////////////////////////////////////
948/// Invert the matrix and calculate its determinant
949/// Notice that the LU decomposition is used instead of Bunch-Kaufman
950/// Bunch-Kaufman guarantees a symmetric inverted matrix but is slower than LU .
951/// The user can access Bunch-Kaufman through the TDecompBK class .
952
953template<class Element>
955{
956 R__ASSERT(this->IsValid());
957 TMatrixD tmp(*this);
958 if (TDecompLU::InvertLU(tmp,Double_t(this->fTol),det)) {
959 const Double_t *p1 = tmp.GetMatrixArray();
960 Element *p2 = this->GetMatrixArray();
961 for (Int_t i = 0; i < this->GetNoElements(); i++)
962 p2[i] = p1[i];
963 }
964
965 return *this;
966}
967
968////////////////////////////////////////////////////////////////////////////////
969/// Invert the matrix and calculate its determinant
970
971template<class Element>
973{
974 R__ASSERT(this->IsValid());
975
976 const Char_t nRows = Char_t(this->GetNrows());
977 switch (nRows) {
978 case 1:
979 {
980 Element *pM = this->GetMatrixArray();
981 if (*pM == 0.) {
982 Error("InvertFast","matrix is singular");
983 *det = 0;
984 } else {
985 *det = *pM;
986 *pM = 1.0/(*pM);
987 }
988 return *this;
989 }
990 case 2:
991 {
992 TMatrixTSymCramerInv::Inv2x2<Element>(*this,det);
993 return *this;
994 }
995 case 3:
996 {
997 TMatrixTSymCramerInv::Inv3x3<Element>(*this,det);
998 return *this;
999 }
1000 case 4:
1001 {
1002 TMatrixTSymCramerInv::Inv4x4<Element>(*this,det);
1003 return *this;
1004 }
1005 case 5:
1006 {
1007 TMatrixTSymCramerInv::Inv5x5<Element>(*this,det);
1008 return *this;
1009 }
1010 case 6:
1011 {
1012 TMatrixTSymCramerInv::Inv6x6<Element>(*this,det);
1013 return *this;
1014 }
1015
1016 default:
1017 {
1018 TMatrixD tmp(*this);
1019 if (TDecompLU::InvertLU(tmp,Double_t(this->fTol),det)) {
1020 const Double_t *p1 = tmp.GetMatrixArray();
1021 Element *p2 = this->GetMatrixArray();
1022 for (Int_t i = 0; i < this->GetNoElements(); i++)
1023 p2[i] = p1[i];
1024 }
1025 return *this;
1026 }
1027 }
1028}
1029
1030////////////////////////////////////////////////////////////////////////////////
1031/// Transpose a matrix.
1032
1033template<class Element>
1035{
1036 if (gMatrixCheck) {
1037 R__ASSERT(this->IsValid());
1038 R__ASSERT(source.IsValid());
1039
1040 if (this->fNrows != source.GetNcols() || this->fRowLwb != source.GetColLwb())
1041 {
1042 Error("Transpose","matrix has wrong shape");
1043 return *this;
1044 }
1045 }
1046
1047 *this = source;
1048 return *this;
1049}
1050
1051////////////////////////////////////////////////////////////////////////////////
1052/// Perform a rank 1 operation on the matrix:
1053/// A += alpha * v * v^T
1054
1055template<class Element>
1057{
1058 if (gMatrixCheck) {
1059 R__ASSERT(this->IsValid());
1060 R__ASSERT(v.IsValid());
1061 if (v.GetNoElements() < this->fNrows) {
1062 Error("Rank1Update","vector too short");
1063 return *this;
1064 }
1065 }
1066
1067 const Element * const pv = v.GetMatrixArray();
1068 Element *trp = this->GetMatrixArray(); // pointer to UR part and diagonal, traverse row-wise
1069 Element *tcp = trp; // pointer to LL part, traverse col-wise
1070 for (Int_t i = 0; i < this->fNrows; i++) {
1071 trp += i; // point to [i,i]
1072 tcp += i*this->fNcols; // point to [i,i]
1073 const Element tmp = alpha*pv[i];
1074 for (Int_t j = i; j < this->fNcols; j++) {
1075 if (j > i) *tcp += tmp*pv[j];
1076 *trp++ += tmp*pv[j];
1077 tcp += this->fNcols;
1078 }
1079 tcp -= this->fNelems-1; // point to [0,i]
1080 }
1081
1082 return *this;
1083}
1084
1085////////////////////////////////////////////////////////////////////////////////
1086/// Calculate B * (*this) * B^T , final matrix will be (nrowsb x nrowsb)
1087/// This is a similarity transform when B is orthogonal . It is more
1088/// efficient than applying the actual multiplication because this
1089/// routine realizes that the final matrix is symmetric .
1090
1091template<class Element>
1093{
1094 if (gMatrixCheck) {
1095 R__ASSERT(this->IsValid());
1096 R__ASSERT(b.IsValid());
1097 if (this->fNcols != b.GetNcols() || this->fColLwb != b.GetColLwb()) {
1098 Error("Similarity(const TMatrixT &)","matrices incompatible");
1099 return *this;
1100 }
1101 }
1102
1103 const Int_t ncolsa = this->fNcols;
1104 const Int_t nb = b.GetNoElements();
1105 const Int_t nrowsb = b.GetNrows();
1106 const Int_t ncolsb = b.GetNcols();
1107
1108 const Element * const bp = b.GetMatrixArray();
1109
1110 Element work[kWorkMax];
1112 Element *bap = work;
1113 if (nrowsb*ncolsa > kWorkMax) {
1115 bap = new Element[nrowsb*ncolsa];
1116 }
1117
1118 AMultB(bp,nb,ncolsb,this->fElements,this->fNelems,this->fNcols,bap);
1119
1120 if (nrowsb != this->fNrows)
1121 this->ResizeTo(nrowsb,nrowsb);
1122
1123#ifdef CBLAS
1124 Element *cp = this->GetMatrixArray();
1125 if (typeid(Element) == typeid(Double_t))
1126 cblas_dgemm (CblasRowMajor,CblasNoTrans,CblasTrans,this->fNrows,this->fNcols,ba.GetNcols(),
1127 1.0,bap,ba.GetNcols(),bp,b.GetNcols(),1.0,cp,this->fNcols);
1128 else if (typeid(Element) != typeid(Float_t))
1129 cblas_sgemm (CblasRowMajor,CblasNoTrans,CblasTrans,this->fNrows,this->fNcols,ba.GetNcols(),
1130 1.0,bap,ba.GetNcols(),bp,b.GetNcols(),1.0,cp,this->fNcols);
1131 else
1132 Error("Similarity","type %s not implemented in BLAS library",typeid(Element));
1133#else
1134 const Int_t nba = nrowsb*ncolsa;
1135 const Int_t ncolsba = ncolsa;
1136 const Element * bi1p = bp;
1137 Element * cp = this->GetMatrixArray();
1138 Element * const cp0 = cp;
1139
1140 Int_t ishift = 0;
1141 const Element *barp0 = bap;
1142 while (barp0 < bap+nba) {
1143 const Element *brp0 = bi1p;
1144 while (brp0 < bp+nb) {
1145 const Element *barp = barp0;
1146 const Element *brp = brp0;
1147 Element cij = 0;
1148 while (brp < brp0+ncolsb)
1149 cij += *barp++ * *brp++;
1150 *cp++ = cij;
1151 brp0 += ncolsb;
1152 }
1153 barp0 += ncolsba;
1154 bi1p += ncolsb;
1155 cp += ++ishift;
1156 }
1157
1158 R__ASSERT(cp == cp0+this->fNelems+ishift && barp0 == bap+nba);
1159
1160 cp = cp0;
1161 for (Int_t irow = 0; irow < this->fNrows; irow++) {
1162 const Int_t rowOff1 = irow*this->fNrows;
1163 for (Int_t icol = 0; icol < irow; icol++) {
1164 const Int_t rowOff2 = icol*this->fNrows;
1166 }
1167 }
1168#endif
1169
1170 if (isAllocated)
1171 delete [] bap;
1172
1173 return *this;
1174}
1175
1176////////////////////////////////////////////////////////////////////////////////
1177/// Calculate B * (*this) * B^T , final matrix will be (nrowsb x nrowsb)
1178/// This is a similarity transform when B is orthogonal . It is more
1179/// efficient than applying the actual multiplication because this
1180/// routine realizes that the final matrix is symmetric .
1181
1182template<class Element>
1184{
1185 if (gMatrixCheck) {
1186 R__ASSERT(this->IsValid());
1187 R__ASSERT(b.IsValid());
1188 if (this->fNcols != b.GetNcols() || this->fColLwb != b.GetColLwb()) {
1189 Error("Similarity(const TMatrixTSym &)","matrices incompatible");
1190 return *this;
1191 }
1192 }
1193
1194#ifdef CBLAS
1195 const Int_t nrowsb = b.GetNrows();
1196 const Int_t ncolsa = this->GetNcols();
1197
1198 Element work[kWorkMax];
1200 Element *abtp = work;
1201 if (this->fNcols > kWorkMax) {
1203 abtp = new Element[this->fNcols];
1204 }
1205
1207
1208 const Element *bp = b.GetMatrixArray();
1209 Element *cp = this->GetMatrixArray();
1210 if (typeid(Element) == typeid(Double_t))
1211 cblas_dsymm (CblasRowMajor,CblasLeft,CblasUpper,this->fNrows,this->fNcols,1.0,
1212 bp,b.GetNcols(),abtp,abt.GetNcols(),0.0,cp,this->fNcols);
1213 else if (typeid(Element) != typeid(Float_t))
1214 cblas_ssymm (CblasRowMajor,CblasLeft,CblasUpper,this->fNrows,this->fNcols,1.0,
1215 bp,b.GetNcols(),abtp,abt.GetNcols(),0.0,cp,this->fNcols);
1216 else
1217 Error("Similarity","type %s not implemented in BLAS library",typeid(Element));
1218
1219 if (isAllocated)
1220 delete [] abtp;
1221#else
1222 const Int_t ncolsa = this->GetNcols();
1223 const Int_t nb = b.GetNoElements();
1224 const Int_t nrowsb = b.GetNrows();
1225 const Int_t ncolsb = b.GetNcols();
1226
1227 const Element * const bp = b.GetMatrixArray();
1228
1229 Element work[kWorkMax];
1231 Element *bap = work;
1232 if (nrowsb*ncolsa > kWorkMax) {
1234 bap = new Element[nrowsb*ncolsa];
1235 }
1236
1237 AMultB(bp,nb,ncolsb,this->fElements,this->fNelems,this->fNcols,bap);
1238
1239 const Int_t nba = nrowsb*ncolsa;
1240 const Int_t ncolsba = ncolsa;
1241 const Element * bi1p = bp;
1242 Element * cp = this->GetMatrixArray();
1243 Element * const cp0 = cp;
1244
1245 Int_t ishift = 0;
1246 const Element *barp0 = bap;
1247 while (barp0 < bap+nba) {
1248 const Element *brp0 = bi1p;
1249 while (brp0 < bp+nb) {
1250 const Element *barp = barp0;
1251 const Element *brp = brp0;
1252 Element cij = 0;
1253 while (brp < brp0+ncolsb)
1254 cij += *barp++ * *brp++;
1255 *cp++ = cij;
1256 brp0 += ncolsb;
1257 }
1258 barp0 += ncolsba;
1259 bi1p += ncolsb;
1260 cp += ++ishift;
1261 }
1262
1263 R__ASSERT(cp == cp0+this->fNelems+ishift && barp0 == bap+nba);
1264
1265 cp = cp0;
1266 for (Int_t irow = 0; irow < this->fNrows; irow++) {
1267 const Int_t rowOff1 = irow*this->fNrows;
1268 for (Int_t icol = 0; icol < irow; icol++) {
1269 const Int_t rowOff2 = icol*this->fNrows;
1271 }
1272 }
1273
1274 if (isAllocated)
1275 delete [] bap;
1276#endif
1277
1278 return *this;
1279}
1280
1281////////////////////////////////////////////////////////////////////////////////
1282/// Calculate scalar v * (*this) * v^T
1283
1284template<class Element>
1286{
1287 if (gMatrixCheck) {
1288 R__ASSERT(this->IsValid());
1289 R__ASSERT(v.IsValid());
1290 if (this->fNcols != v.GetNrows() || this->fColLwb != v.GetLwb()) {
1291 Error("Similarity(const TVectorT &)","vector and matrix incompatible");
1292 return -1.;
1293 }
1294 }
1295
1296 const Element *mp = this->GetMatrixArray(); // Matrix row ptr
1297 const Element *vp = v.GetMatrixArray(); // vector ptr
1298
1299 Element sum1 = 0;
1300 const Element * const vp_first = vp;
1301 const Element * const vp_last = vp+v.GetNrows();
1302 while (vp < vp_last) {
1303 Element sum2 = 0;
1304 for (const Element *sp = vp_first; sp < vp_last; )
1305 sum2 += *mp++ * *sp++;
1306 sum1 += sum2 * *vp++;
1307 }
1308
1309 R__ASSERT(mp == this->GetMatrixArray()+this->GetNoElements());
1310
1311 return sum1;
1312}
1313
1314////////////////////////////////////////////////////////////////////////////////
1315/// Calculate B^T * (*this) * B , final matrix will be (ncolsb x ncolsb)
1316/// It is more efficient than applying the actual multiplication because this
1317/// routine realizes that the final matrix is symmetric .
1318
1319template<class Element>
1321{
1322 if (gMatrixCheck) {
1323 R__ASSERT(this->IsValid());
1324 R__ASSERT(b.IsValid());
1325 if (this->fNrows != b.GetNrows() || this->fRowLwb != b.GetRowLwb()) {
1326 Error("SimilarityT(const TMatrixT &)","matrices incompatible");
1327 return *this;
1328 }
1329 }
1330
1331 const Int_t ncolsb = b.GetNcols();
1332 const Int_t ncolsa = this->GetNcols();
1333
1334 Element work[kWorkMax];
1336 Element *btap = work;
1337 if (ncolsb*ncolsa > kWorkMax) {
1339 btap = new Element[ncolsb*ncolsa];
1340 }
1341
1343 bta.TMult(b,*this);
1344
1345 if (ncolsb != this->fNcols)
1346 this->ResizeTo(ncolsb,ncolsb);
1347
1348#ifdef CBLAS
1349 const Element *bp = b.GetMatrixArray();
1350 Element *cp = this->GetMatrixArray();
1351 if (typeid(Element) == typeid(Double_t))
1352 cblas_dgemm (CblasRowMajor,CblasNoTrans,CblasNoTrans,this->fNrows,this->fNcols,bta.GetNcols(),
1353 1.0,btap,bta.GetNcols(),bp,b.GetNcols(),1.0,cp,this->fNcols);
1354 else if (typeid(Element) != typeid(Float_t))
1355 cblas_sgemm (CblasRowMajor,CblasNoTrans,CblasNoTrans,this->fNrows,this->fNcols,bta.GetNcols(),
1356 1.0,btap,bta.GetNcols(),bp,b.GetNcols(),1.0,cp,this->fNcols);
1357 else
1358 Error("similarityT","type %s not implemented in BLAS library",typeid(Element));
1359#else
1360 const Int_t nbta = bta.GetNoElements();
1361 const Int_t nb = b.GetNoElements();
1362 const Int_t ncolsbta = bta.GetNcols();
1363 const Element * const bp = b.GetMatrixArray();
1364 Element * cp = this->GetMatrixArray();
1365 Element * const cp0 = cp;
1366
1367 Int_t ishift = 0;
1368 const Element *btarp0 = btap; // Pointer to A[i,0];
1369 const Element *bcp0 = bp;
1370 while (btarp0 < btap+nbta) {
1371 for (const Element *bcp = bcp0; bcp < bp+ncolsb; ) { // Pointer to the j-th column of B, Start bcp = B[0,0]
1372 const Element *btarp = btarp0; // Pointer to the i-th row of A, reset to A[i,0]
1373 Element cij = 0;
1374 while (bcp < bp+nb) { // Scan the i-th row of A and
1375 cij += *btarp++ * *bcp; // the j-th col of B
1376 bcp += ncolsb;
1377 }
1378 *cp++ = cij;
1379 bcp -= nb-1; // Set bcp to the (j+1)-th col
1380 }
1381 btarp0 += ncolsbta; // Set ap to the (i+1)-th row
1382 bcp0++;
1383 cp += ++ishift;
1384 }
1385
1386 R__ASSERT(cp == cp0+this->fNelems+ishift && btarp0 == btap+nbta);
1387
1388 cp = cp0;
1389 for (Int_t irow = 0; irow < this->fNrows; irow++) {
1390 const Int_t rowOff1 = irow*this->fNrows;
1391 for (Int_t icol = 0; icol < irow; icol++) {
1392 const Int_t rowOff2 = icol*this->fNrows;
1394 }
1395 }
1396#endif
1397
1398 if (isAllocated)
1399 delete [] btap;
1400
1401 return *this;
1402}
1403
1404////////////////////////////////////////////////////////////////////////////////
1405
1406template<class Element>
1408{
1409 if (gMatrixCheck && !AreCompatible(*this,source)) {
1410 Error("operator=","matrices not compatible");
1411 return *this;
1412 }
1413
1414 if (this->GetMatrixArray() != source.GetMatrixArray()) {
1416 memcpy(this->GetMatrixArray(),source.fElements,this->fNelems*sizeof(Element));
1417 }
1418 return *this;
1419}
1420
1421////////////////////////////////////////////////////////////////////////////////
1422
1423template<class Element>
1425{
1426 R__ASSERT(this->IsValid());
1427
1428 if (lazy_constructor.fRowUpb != this->GetRowUpb() ||
1429 lazy_constructor.fRowLwb != this->GetRowLwb()) {
1430 Error("operator=(const TMatrixTSymLazy&)", "matrix is incompatible with "
1431 "the assigned Lazy matrix");
1432 return *this;
1433 }
1434
1435 lazy_constructor.FillIn(*this);
1436 return *this;
1437}
1438
1439////////////////////////////////////////////////////////////////////////////////
1440/// Assign val to every element of the matrix.
1441
1442template<class Element>
1444{
1445 R__ASSERT(this->IsValid());
1446
1447 Element *ep = fElements;
1448 const Element * const ep_last = ep+this->fNelems;
1449 while (ep < ep_last)
1450 *ep++ = val;
1451
1452 return *this;
1453}
1454
1455////////////////////////////////////////////////////////////////////////////////
1456/// Add val to every element of the matrix.
1457
1458template<class Element>
1460{
1461 R__ASSERT(this->IsValid());
1462
1463 Element *ep = fElements;
1464 const Element * const ep_last = ep+this->fNelems;
1465 while (ep < ep_last)
1466 *ep++ += val;
1467
1468 return *this;
1469}
1470
1471////////////////////////////////////////////////////////////////////////////////
1472/// Subtract val from every element of the matrix.
1473
1474template<class Element>
1476{
1477 R__ASSERT(this->IsValid());
1478
1479 Element *ep = fElements;
1480 const Element * const ep_last = ep+this->fNelems;
1481 while (ep < ep_last)
1482 *ep++ -= val;
1483
1484 return *this;
1485}
1486
1487////////////////////////////////////////////////////////////////////////////////
1488/// Multiply every element of the matrix with val.
1489
1490template<class Element>
1492{
1493 R__ASSERT(this->IsValid());
1494
1495 Element *ep = fElements;
1496 const Element * const ep_last = ep+this->fNelems;
1497 while (ep < ep_last)
1498 *ep++ *= val;
1499
1500 return *this;
1501}
1502
1503////////////////////////////////////////////////////////////////////////////////
1504/// Add the source matrix.
1505
1506template<class Element>
1508{
1509 if (gMatrixCheck && !AreCompatible(*this,source)) {
1510 Error("operator+=","matrices not compatible");
1511 return *this;
1512 }
1513
1514 const Element *sp = source.GetMatrixArray();
1515 Element *tp = this->GetMatrixArray();
1516 const Element * const tp_last = tp+this->fNelems;
1517 while (tp < tp_last)
1518 *tp++ += *sp++;
1519
1520 return *this;
1521}
1522
1523////////////////////////////////////////////////////////////////////////////////
1524/// Subtract the source matrix.
1525
1526template<class Element>
1528{
1529 if (gMatrixCheck && !AreCompatible(*this,source)) {
1530 Error("operator-=","matrices not compatible");
1531 return *this;
1532 }
1533
1534 const Element *sp = source.GetMatrixArray();
1535 Element *tp = this->GetMatrixArray();
1536 const Element * const tp_last = tp+this->fNelems;
1537 while (tp < tp_last)
1538 *tp++ -= *sp++;
1539
1540 return *this;
1541}
1542
1543////////////////////////////////////////////////////////////////////////////////
1544
1545template<class Element>
1547{
1548 R__ASSERT(this->IsValid());
1549
1550 Element val = 0;
1551 Element *trp = this->GetMatrixArray(); // pointer to UR part and diagonal, traverse row-wise
1552 Element *tcp = trp; // pointer to LL part, traverse col-wise
1553 for (Int_t i = 0; i < this->fNrows; i++) {
1554 trp += i; // point to [i,i]
1555 tcp += i*this->fNcols; // point to [i,i]
1556 for (Int_t j = i; j < this->fNcols; j++) {
1557 action.Operation(val);
1558 if (j > i) *tcp = val;
1559 *trp++ = val;
1560 tcp += this->fNcols;
1561 }
1562 tcp -= this->fNelems-1; // point to [0,i]
1563 }
1564
1565 return *this;
1566}
1567
1568////////////////////////////////////////////////////////////////////////////////
1569/// Apply action to each element of the matrix. To action the location
1570/// of the current element is passed.
1571
1572template<class Element>
1574{
1575 R__ASSERT(this->IsValid());
1576
1577 Element val = 0;
1578 Element *trp = this->GetMatrixArray(); // pointer to UR part and diagonal, traverse row-wise
1579 Element *tcp = trp; // pointer to LL part, traverse col-wise
1580 for (Int_t i = 0; i < this->fNrows; i++) {
1581 action.fI = i+this->fRowLwb;
1582 trp += i; // point to [i,i]
1583 tcp += i*this->fNcols; // point to [i,i]
1584 for (Int_t j = i; j < this->fNcols; j++) {
1585 action.fJ = j+this->fColLwb;
1586 action.Operation(val);
1587 if (j > i) *tcp = val;
1588 *trp++ = val;
1589 tcp += this->fNcols;
1590 }
1591 tcp -= this->fNelems-1; // point to [0,i]
1592 }
1593
1594 return *this;
1595}
1596
1597////////////////////////////////////////////////////////////////////////////////
1598/// randomize matrix element values but keep matrix symmetric
1599
1600template<class Element>
1602{
1603 if (gMatrixCheck) {
1604 R__ASSERT(this->IsValid());
1605 if (this->fNrows != this->fNcols || this->fRowLwb != this->fColLwb) {
1606 Error("Randomize(Element,Element,Element &","matrix should be square");
1607 return *this;
1608 }
1609 }
1610
1611 const Element scale = beta-alpha;
1612 const Element shift = alpha/scale;
1613
1614 Element *ep = GetMatrixArray();
1615 for (Int_t i = 0; i < this->fNrows; i++) {
1616 const Int_t off = i*this->fNcols;
1617 for (Int_t j = 0; j <= i; j++) {
1618 ep[off+j] = scale*(Drand(seed)+shift);
1619 if (i != j) {
1620 ep[j*this->fNcols+i] = ep[off+j];
1621 }
1622 }
1623 }
1624
1625 return *this;
1626}
1627
1628////////////////////////////////////////////////////////////////////////////////
1629/// randomize matrix element values but keep matrix symmetric positive definite
1630
1631template<class Element>
1633{
1634 if (gMatrixCheck) {
1635 R__ASSERT(this->IsValid());
1636 if (this->fNrows != this->fNcols || this->fRowLwb != this->fColLwb) {
1637 Error("RandomizeSym(Element,Element,Element &","matrix should be square");
1638 return *this;
1639 }
1640 }
1641
1642 const Element scale = beta-alpha;
1643 const Element shift = alpha/scale;
1644
1645 Element *ep = GetMatrixArray();
1646 Int_t i;
1647 for (i = 0; i < this->fNrows; i++) {
1648 const Int_t off = i*this->fNcols;
1649 for (Int_t j = 0; j <= i; j++)
1650 ep[off+j] = scale*(Drand(seed)+shift);
1651 }
1652
1653 for (i = this->fNrows-1; i >= 0; i--) {
1654 const Int_t off1 = i*this->fNcols;
1655 for (Int_t j = i; j >= 0; j--) {
1656 const Int_t off2 = j*this->fNcols;
1657 ep[off1+j] *= ep[off2+j];
1658 for (Int_t k = j-1; k >= 0; k--) {
1659 ep[off1+j] += ep[off1+k]*ep[off2+k];
1660 }
1661 if (i != j)
1662 ep[off2+i] = ep[off1+j];
1663 }
1664 }
1665
1666 return *this;
1667}
1668
1669////////////////////////////////////////////////////////////////////////////////
1670/// Return a matrix containing the eigen-vectors ordered by descending eigen-values.
1671/// For full functionality use TMatrixDSymEigen .
1672
1673template<class Element>
1675{
1676 TMatrixDSym tmp = *this;
1678 eigenValues.ResizeTo(this->fNrows);
1679 eigenValues = eigen.GetEigenValues();
1680 return eigen.GetEigenVectors();
1681}
1682
1683////////////////////////////////////////////////////////////////////////////////
1684/// Check to see if two matrices are identical.
1685
1686template<class Element>
1688{
1689 if (!AreCompatible(m1,m2)) return kFALSE;
1690 return (memcmp(m1.GetMatrixArray(),m2.GetMatrixArray(),
1691 m1.GetNoElements()*sizeof(Element)) == 0);
1692}
1693
1694////////////////////////////////////////////////////////////////////////////////
1695
1696template<class Element>
1703
1704////////////////////////////////////////////////////////////////////////////////
1705
1706template<class Element>
1708{
1710 target += val;
1711 return target;
1712}
1713
1714////////////////////////////////////////////////////////////////////////////////
1715
1716template<class Element>
1718{
1719 return operator+(source1,val);
1720}
1721
1722////////////////////////////////////////////////////////////////////////////////
1723
1724template<class Element>
1731
1732////////////////////////////////////////////////////////////////////////////////
1733
1734template<class Element>
1736{
1738 target -= val;
1739 return target;
1740}
1741
1742////////////////////////////////////////////////////////////////////////////////
1743
1744template<class Element>
1746{
1747 return Element(-1.0)*operator-(source1,val);
1748}
1749
1750////////////////////////////////////////////////////////////////////////////////
1751
1752template<class Element>
1754{
1756 target *= val;
1757 return target;
1758}
1759
1760////////////////////////////////////////////////////////////////////////////////
1761
1762template<class Element>
1764{
1765 return operator*(source1,val);
1766}
1767
1768////////////////////////////////////////////////////////////////////////////////
1769/// Logical AND
1770
1771template<class Element>
1773{
1775
1777 Error("operator&&(const TMatrixTSym&,const TMatrixTSym&)","matrices not compatible");
1778 return target;
1779 }
1780
1781 target.ResizeTo(source1);
1782
1783 const Element *sp1 = source1.GetMatrixArray();
1784 const Element *sp2 = source2.GetMatrixArray();
1785 Element *tp = target.GetMatrixArray();
1786 const Element * const tp_last = tp+target.GetNoElements();
1787 while (tp < tp_last)
1788 *tp++ = (*sp1++ != 0.0 && *sp2++ != 0.0);
1789
1790 return target;
1791}
1792
1793////////////////////////////////////////////////////////////////////////////////
1794/// Logical Or
1795
1796template<class Element>
1798{
1800
1802 Error("operator||(const TMatrixTSym&,const TMatrixTSym&)","matrices not compatible");
1803 return target;
1804 }
1805
1806 target.ResizeTo(source1);
1807
1808 const Element *sp1 = source1.GetMatrixArray();
1809 const Element *sp2 = source2.GetMatrixArray();
1810 Element *tp = target.GetMatrixArray();
1811 const Element * const tp_last = tp+target.GetNoElements();
1812 while (tp < tp_last)
1813 *tp++ = (*sp1++ != 0.0 || *sp2++ != 0.0);
1814
1815 return target;
1816}
1817
1818////////////////////////////////////////////////////////////////////////////////
1819/// source1 > source2
1820
1821template<class Element>
1823{
1825
1827 Error("operator>(const TMatrixTSym&,const TMatrixTSym&)","matrices not compatible");
1828 return target;
1829 }
1830
1831 target.ResizeTo(source1);
1832
1833 const Element *sp1 = source1.GetMatrixArray();
1834 const Element *sp2 = source2.GetMatrixArray();
1835 Element *tp = target.GetMatrixArray();
1836 const Element * const tp_last = tp+target.GetNoElements();
1837 while (tp < tp_last) {
1838 *tp++ = (*sp1) > (*sp2); sp1++; sp2++;
1839 }
1840
1841 return target;
1842}
1843
1844////////////////////////////////////////////////////////////////////////////////
1845/// source1 >= source2
1846
1847template<class Element>
1849{
1851
1853 Error("operator>=(const TMatrixTSym&,const TMatrixTSym&)","matrices not compatible");
1854 return target;
1855 }
1856
1857 target.ResizeTo(source1);
1858
1859 const Element *sp1 = source1.GetMatrixArray();
1860 const Element *sp2 = source2.GetMatrixArray();
1861 Element *tp = target.GetMatrixArray();
1862 const Element * const tp_last = tp+target.GetNoElements();
1863 while (tp < tp_last) {
1864 *tp++ = (*sp1) >= (*sp2); sp1++; sp2++;
1865 }
1866
1867 return target;
1868}
1869
1870////////////////////////////////////////////////////////////////////////////////
1871/// source1 <= source2
1872
1873template<class Element>
1875{
1877
1879 Error("operator<=(const TMatrixTSym&,const TMatrixTSym&)","matrices not compatible");
1880 return target;
1881 }
1882
1883 target.ResizeTo(source1);
1884
1885 const Element *sp1 = source1.GetMatrixArray();
1886 const Element *sp2 = source2.GetMatrixArray();
1887 Element *tp = target.GetMatrixArray();
1888 const Element * const tp_last = tp+target.GetNoElements();
1889 while (tp < tp_last) {
1890 *tp++ = (*sp1) <= (*sp2); sp1++; sp2++;
1891 }
1892
1893 return target;
1894}
1895
1896////////////////////////////////////////////////////////////////////////////////
1897/// source1 < source2
1898
1899template<class Element>
1901{
1903
1905 Error("operator<(const TMatrixTSym&,const TMatrixTSym&)","matrices not compatible");
1906 return target;
1907 }
1908
1909 target.ResizeTo(source1);
1910
1911 const Element *sp1 = source1.GetMatrixArray();
1912 const Element *sp2 = source2.GetMatrixArray();
1913 Element *tp = target.GetMatrixArray();
1914 const Element * const tp_last = tp+target.GetNoElements();
1915 while (tp < tp_last) {
1916 *tp++ = (*sp1) < (*sp2); sp1++; sp2++;
1917 }
1918
1919 return target;
1920}
1921
1922////////////////////////////////////////////////////////////////////////////////
1923/// Modify addition: target += scalar * source.
1924
1925template<class Element>
1927{
1929 ::Error("Add","matrices not compatible");
1930 return target;
1931 }
1932
1933 const Int_t nrows = target.GetNrows();
1934 const Int_t ncols = target.GetNcols();
1935 const Int_t nelems = target.GetNoElements();
1936 const Element *sp = source.GetMatrixArray();
1937 Element *trp = target.GetMatrixArray(); // pointer to UR part and diagonal, traverse row-wise
1938 Element *tcp = target.GetMatrixArray(); // pointer to LL part, traverse col-wise
1939 for (Int_t i = 0; i < nrows; i++) {
1940 sp += i;
1941 trp += i; // point to [i,i]
1942 tcp += i*ncols; // point to [i,i]
1943 for (Int_t j = i; j < ncols; j++) {
1944 const Element tmp = scalar * *sp++;
1945 if (j > i) *tcp += tmp;
1946 *trp++ += tmp;
1947 tcp += ncols;
1948 }
1949 tcp -= nelems-1; // point to [0,i]
1950 }
1951
1952 return target;
1953}
1954
1955////////////////////////////////////////////////////////////////////////////////
1956/// Multiply target by the source, element-by-element.
1957
1958template<class Element>
1960{
1962 ::Error("ElementMult","matrices not compatible");
1963 return target;
1964 }
1965
1966 const Int_t nrows = target.GetNrows();
1967 const Int_t ncols = target.GetNcols();
1968 const Int_t nelems = target.GetNoElements();
1969 const Element *sp = source.GetMatrixArray();
1970 Element *trp = target.GetMatrixArray(); // pointer to UR part and diagonal, traverse row-wise
1971 Element *tcp = target.GetMatrixArray(); // pointer to LL part, traverse col-wise
1972 for (Int_t i = 0; i < nrows; i++) {
1973 sp += i;
1974 trp += i; // point to [i,i]
1975 tcp += i*ncols; // point to [i,i]
1976 for (Int_t j = i; j < ncols; j++) {
1977 if (j > i) *tcp *= *sp;
1978 *trp++ *= *sp++;
1979 tcp += ncols;
1980 }
1981 tcp -= nelems-1; // point to [0,i]
1982 }
1983
1984 return target;
1985}
1986
1987////////////////////////////////////////////////////////////////////////////////
1988/// Multiply target by the source, element-by-element.
1989
1990template<class Element>
1992{
1994 ::Error("ElementDiv","matrices not compatible");
1995 return target;
1996 }
1997
1998 const Int_t nrows = target.GetNrows();
1999 const Int_t ncols = target.GetNcols();
2000 const Int_t nelems = target.GetNoElements();
2001 const Element *sp = source.GetMatrixArray();
2002 Element *trp = target.GetMatrixArray(); // pointer to UR part and diagonal, traverse row-wise
2003 Element *tcp = target.GetMatrixArray(); // pointer to LL part, traverse col-wise
2004 for (Int_t i = 0; i < nrows; i++) {
2005 sp += i;
2006 trp += i; // point to [i,i]
2007 tcp += i*ncols; // point to [i,i]
2008 for (Int_t j = i; j < ncols; j++) {
2009 if (*sp != 0.0) {
2010 if (j > i) *tcp /= *sp;
2011 *trp++ /= *sp++;
2012 } else {
2013 const Int_t irow = (sp-source.GetMatrixArray())/source.GetNcols();
2014 const Int_t icol = (sp-source.GetMatrixArray())%source.GetNcols();
2015 Error("ElementDiv","source (%d,%d) is zero",irow,icol);
2016 trp++;
2017 }
2018 tcp += ncols;
2019 }
2020 tcp -= nelems-1; // point to [0,i]
2021 }
2022
2023 return target;
2024}
2025
2026////////////////////////////////////////////////////////////////////////////////
2027/// Stream an object of class TMatrixTSym.
2028
2029template<class Element>
2031{
2032 if (R__b.IsReading()) {
2033 UInt_t R__s, R__c;
2034 Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
2035 Clear();
2036 R__b.ReadClassBuffer(TMatrixTBase<Element>::Class(),this,R__v,R__s,R__c);
2037 fElements = new Element[this->fNelems];
2038 Int_t i;
2039 for (i = 0; i < this->fNrows; i++) {
2040 R__b.ReadFastArray(fElements+i*this->fNcols+i,this->fNcols-i);
2041 }
2042 // copy to Lower left triangle
2043 for (i = 0; i < this->fNrows; i++) {
2044 for (Int_t j = 0; j < i; j++) {
2045 fElements[i*this->fNcols+j] = fElements[j*this->fNrows+i];
2046 }
2047 }
2048 if (this->fNelems <= this->kSizeMax) {
2049 memcpy(fDataStack,fElements,this->fNelems*sizeof(Element));
2050 delete [] fElements;
2051 fElements = fDataStack;
2052 }
2053 } else {
2054 R__b.WriteClassBuffer(TMatrixTBase<Element>::Class(),this);
2055 // Only write the Upper right triangle
2056 for (Int_t i = 0; i < this->fNrows; i++) {
2057 R__b.WriteFastArray(fElements+i*this->fNcols+i,this->fNcols-i);
2058 }
2059 }
2060}
2061
2062#include "TMatrixFSymfwd.h"
2063
2064template class TMatrixTSym<Float_t>;
2065
2066template Bool_t operator== <Float_t>(const TMatrixFSym &source1,const TMatrixFSym &source2);
2067template TMatrixFSym operator+ <Float_t>(const TMatrixFSym &source1,const TMatrixFSym &source2);
2068template TMatrixFSym operator+ <Float_t>(const TMatrixFSym &source1, Float_t val);
2069template TMatrixFSym operator+ <Float_t>( Float_t val ,const TMatrixFSym &source2);
2070template TMatrixFSym operator- <Float_t>(const TMatrixFSym &source1,const TMatrixFSym &source2);
2071template TMatrixFSym operator- <Float_t>(const TMatrixFSym &source1, Float_t val);
2072template TMatrixFSym operator- <Float_t>( Float_t val ,const TMatrixFSym &source2);
2073template TMatrixFSym operator* <Float_t>(const TMatrixFSym &source, Float_t val );
2074template TMatrixFSym operator* <Float_t>( Float_t val, const TMatrixFSym &source );
2075template TMatrixFSym operator&& <Float_t>(const TMatrixFSym &source1,const TMatrixFSym &source2);
2076template TMatrixFSym operator|| <Float_t>(const TMatrixFSym &source1,const TMatrixFSym &source2);
2077template TMatrixFSym operator> <Float_t>(const TMatrixFSym &source1,const TMatrixFSym &source2);
2078template TMatrixFSym operator>= <Float_t>(const TMatrixFSym &source1,const TMatrixFSym &source2);
2081
2085
2086#include "TMatrixDSymfwd.h"
2087
2088template class TMatrixTSym<Double_t>;
2089
2090template Bool_t operator== <Double_t>(const TMatrixDSym &source1,const TMatrixDSym &source2);
2091template TMatrixDSym operator+ <Double_t>(const TMatrixDSym &source1,const TMatrixDSym &source2);
2092template TMatrixDSym operator+ <Double_t>(const TMatrixDSym &source1, Double_t val);
2093template TMatrixDSym operator+ <Double_t>( Double_t val ,const TMatrixDSym &source2);
2094template TMatrixDSym operator- <Double_t>(const TMatrixDSym &source1,const TMatrixDSym &source2);
2095template TMatrixDSym operator- <Double_t>(const TMatrixDSym &source1, Double_t val);
2096template TMatrixDSym operator- <Double_t>( Double_t val ,const TMatrixDSym &source2);
2097template TMatrixDSym operator* <Double_t>(const TMatrixDSym &source, Double_t val );
2098template TMatrixDSym operator* <Double_t>( Double_t val, const TMatrixDSym &source );
2099template TMatrixDSym operator&& <Double_t>(const TMatrixDSym &source1,const TMatrixDSym &source2);
2100template TMatrixDSym operator|| <Double_t>(const TMatrixDSym &source1,const TMatrixDSym &source2);
2101template TMatrixDSym operator> <Double_t>(const TMatrixDSym &source1,const TMatrixDSym &source2);
2102template TMatrixDSym operator>= <Double_t>(const TMatrixDSym &source1,const TMatrixDSym &source2);
2105
#define b(i)
Definition RSha256.hxx:100
#define a(i)
Definition RSha256.hxx:99
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
bool Bool_t
Boolean (0=false, 1=true) (bool)
Definition RtypesCore.h:77
int Int_t
Signed integer 4 bytes (int)
Definition RtypesCore.h:59
short Version_t
Class version identifier (short)
Definition RtypesCore.h:79
char Char_t
Character 1 byte (char)
Definition RtypesCore.h:51
float Float_t
Float 4 bytes (float)
Definition RtypesCore.h:71
constexpr Bool_t kFALSE
Definition RtypesCore.h:108
double Double_t
Double 8 bytes.
Definition RtypesCore.h:73
constexpr Bool_t kTRUE
Definition RtypesCore.h:107
const char Option_t
Option string (const char)
Definition RtypesCore.h:80
@ kPlus
Definition TAttMarker.h:56
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
Bool_t operator<=(const TDatime &d1, const TDatime &d2)
Definition TDatime.h:108
Bool_t operator>(const TDatime &d1, const TDatime &d2)
Definition TDatime.h:110
Bool_t operator>=(const TDatime &d1, const TDatime &d2)
Definition TDatime.h:112
Bool_t operator<(const TDatime &d1, const TDatime &d2)
Definition TDatime.h:106
#define R__ASSERT(e)
Checks condition e and reports a fatal error if it's false.
Definition TError.h:125
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition TError.cxx:208
Option_t Option_t option
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t target
R__EXTERN Int_t gMatrixCheck
template TMatrixFSparse & ElementDiv< Float_t >(TMatrixFSparse &target, const TMatrixFSparse &source)
template TMatrixDSparse & Add< Double_t >(TMatrixDSparse &target, Double_t scalar, const TMatrixDSparse &source)
template TMatrixDSparse & ElementDiv< Double_t >(TMatrixDSparse &target, const TMatrixDSparse &source)
template TMatrixFSparse & ElementMult< Float_t >(TMatrixFSparse &target, const TMatrixFSparse &source)
template TMatrixDSparse & ElementMult< Double_t >(TMatrixDSparse &target, const TMatrixDSparse &source)
template TMatrixFSparse & Add< Float_t >(TMatrixFSparse &target, Float_t scalar, const TMatrixFSparse &source)
template TMatrixFSym operator<=< Float_t >(const TMatrixFSym &source1, const TMatrixFSym &source2)
template TMatrixFSym operator<<Float_t >(const TMatrixFSym &source1, const TMatrixFSym &source2)
template TMatrixDSym operator<<Double_t >(const TMatrixDSym &source1, const TMatrixDSym &source2)
template TMatrixDSym operator<=< Double_t >(const TMatrixDSym &source1, const TMatrixDSym &source2)
Double_t Drand(Double_t &ix)
Random number generator [0....1] with seed ix.
TTime operator*(const TTime &t1, const TTime &t2)
Definition TTime.h:85
TTime operator-(const TTime &t1, const TTime &t2)
Definition TTime.h:83
Buffer base class used for serializing objects.
Definition TBuffer.h:43
LU Decomposition class.
Definition TDecompLU.h:24
static Bool_t InvertLU(TMatrixD &a, Double_t tol, Double_t *det=nullptr)
Calculate matrix inversion through in place forward/backward substitution.
TMatrixDSymEigen.
TMatrixTBase.
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.
virtual TMatrixTBase< Element > & SetMatrixArray(const Element *data, Option_t *option="")
Copy array data to matrix .
TMatrixTSym< Element > & operator*=(Element val)
Multiply every element of the matrix with val.
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...
TMatrixTSym< Element > & operator+=(Element val)
Add val to every element of the matrix.
virtual TMatrixTSym< Element > & RandomizePD(Element alpha, Element beta, Double_t &seed)
randomize matrix element values but keep matrix symmetric positive definite
TMatrixTBase< Element > & Randomize(Element alpha, Element beta, Double_t &seed) override
randomize matrix element values but keep matrix symmetric
TMatrixTSym< Element > & Transpose(const TMatrixTSym< Element > &source)
Transpose a matrix.
Element * New_m(Int_t size)
return data pointer .
TMatrixTSym< Element > & Use(Int_t row_lwb, Int_t row_upb, Element *data)
TMatrixTSym< Element > & Rank1Update(const TVectorT< Element > &v, Element alpha=1.0)
Perform a rank 1 operation on the matrix: A += alpha * v * v^T.
void Minus(const TMatrixTSym< Element > &a, const TMatrixTSym< Element > &b)
Symmetric matrix subtraction. Replace this matrix with C such that C = A - B.
void TMult(const TMatrixT< Element > &a)
Replace this matrix with C such that C = A' * A.
void Delete_m(Int_t size, Element *&)
delete data pointer m, if it was assigned on the heap
void Allocate(Int_t nrows, Int_t ncols, Int_t row_lwb=0, Int_t col_lwb=0, Int_t init=0, Int_t=-1)
Allocate new matrix.
TMatrixTSym< Element > & GetSub(Int_t row_lwb, Int_t row_upb, TMatrixTSym< Element > &target, Option_t *option="S") const
Get submatrix [row_lwb..row_upb][row_lwb..row_upb]; The indexing range of the returned matrix depends...
Int_t Memcpy_m(Element *newp, const Element *oldp, Int_t copySize, Int_t newSize, Int_t oldSize)
copy copySize doubles from *oldp to *newp .
const TMatrixT< Element > EigenVectors(TVectorT< Element > &eigenValues) const
Return a matrix containing the eigen-vectors ordered by descending eigen-values.
TMatrixTSym< Element > & operator-=(Element val)
Subtract val from every element of the matrix.
void Plus(const TMatrixTSym< Element > &a, const TMatrixTSym< Element > &b)
Symmetric matrix summation. Replace this matrix with C such that C = A + B.
TMatrixTSym< Element > & InvertFast(Double_t *det=nullptr)
Invert the matrix and calculate its determinant.
TMatrixTSym< Element > & SimilarityT(const TMatrixT< Element > &n)
Calculate B^T * (*this) * B , final matrix will be (ncolsb x ncolsb) It is more efficient than applyi...
TMatrixTSym< Element > & Similarity(const TMatrixT< Element > &n)
Calculate B * (*this) * B^T , final matrix will be (nrowsb x nrowsb) This is a similarity transform w...
TMatrixTSym< Element > & operator=(const TMatrixTSym< Element > &source)
TMatrixTSym< Element > & SetSub(Int_t row_lwb, const TMatrixTBase< Element > &source)
Insert matrix source starting at [row_lwb][row_lwb], thereby overwriting the part [row_lwb....
TMatrixTBase< Element > & SetMatrixArray(const Element *data, Option_t *option="") override
Copy array data to matrix .
TMatrixTBase< Element > & Shift(Int_t row_shift, Int_t col_shift) override
Shift the row index by adding row_shift and the column index by adding col_shift, respectively.
TMatrixTBase< Element > & Apply(const TElementActionT< Element > &action) override
Apply action to each matrix element.
Double_t Determinant() const override
TMatrixTSym< Element > & Invert(Double_t *det=nullptr)
Invert the matrix and calculate its determinant Notice that the LU decomposition is used instead of B...
void Streamer(TBuffer &) override
Stream an object of class TMatrixTSym.
TObject & operator=(const TObject &rhs) noexcept
TObject assignment operator.
Definition TObject.h:299
Basic string class.
Definition TString.h:138
void ToUpper()
Change string to upper case.
Definition TString.cxx:1202
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition TString.h:640
bool operator==(const RConcurrentHashColl::HashValue &lhs, const RConcurrentHashColl::HashValue &rhs)
ELogLevel operator+(ELogLevel severity, int offset)
Definition RLogger.hxx:42
Short_t Max(Short_t a, Short_t b)
Returns the largest of a and b.
Definition TMathBase.h:251
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition TMath.h:732
Short_t Min(Short_t a, Short_t b)
Returns the smallest of a and b.
Definition TMathBase.h:199
void AMultB(const Element *const ap, Int_t na, Int_t ncolsa, const Element *const bp, Int_t nb, Int_t ncolsb, Element *cp)
Elementary routine to calculate matrix multiplication A*B.
TMatrixT< Element > & ElementMult(TMatrixT< Element > &target, const TMatrixT< Element > &source)
Multiply target by the source, element-by-element.
TMatrixT< Element > operator||(const TMatrixT< Element > &source1, const TMatrixT< Element > &source2)
Logical OR.
TMatrixT< Element > & Add(TMatrixT< Element > &target, Element scalar, const TMatrixT< Element > &source)
Modify addition: target += scalar * source.
TMatrixT< Element > & ElementDiv(TMatrixT< Element > &target, const TMatrixT< Element > &source)
Divide target by the source, element-by-element.
TMatrixT< Element > operator&&(const TMatrixT< Element > &source1, const TMatrixT< Element > &source2)
Logical AND.
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 .
TMarker m
Definition textangle.C:8