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