Logo ROOT  

Reference Guide
 
Loading...
Searching...
No Matches
TMatrixT.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 TMatrixT
13 \ingroup Matrix
14
15TMatrixT
16
17Template class of a general matrix in the linear algebra package
18
19See the \ref Matrix page for the documentation of the linear algebra package
20
21*/
22
23#include <typeinfo>
24
25#include "TMatrixT.h"
26#include "TBuffer.h"
27#include "TMatrixTSym.h"
28#include "TMatrixTLazy.h"
29#include "TMatrixTCramerInv.h"
30#include "TDecompLU.h"
31#include "TMatrixDEigen.h"
32#include "TMath.h"
33
34////////////////////////////////////////////////////////////////////////////////
35/// Constructor for (nrows x ncols) matrix
36
37template <class Element>
38TMatrixT<Element>::TMatrixT(Int_t nrows, Int_t ncols)
39{
40 Allocate(nrows, ncols, 0, 0, 1);
41}
42
43////////////////////////////////////////////////////////////////////////////////
44/// Constructor for ([row_lwb..row_upb] x [col_lwb..col_upb]) matrix
45
46template <class Element>
47TMatrixT<Element>::TMatrixT(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb)
48{
49 Allocate(row_upb - row_lwb + 1, col_upb - col_lwb + 1, row_lwb, col_lwb, 1);
50}
51
52////////////////////////////////////////////////////////////////////////////////
53/// option="F": array elements contains the matrix stored column-wise
54/// like in Fortran, so a[i,j] = elements[i+no_rows*j],
55/// else it is supposed that array elements are stored row-wise
56/// a[i,j] = elements[i*no_cols+j]
57///
58/// array elements are copied
59
60template <class Element>
61TMatrixT<Element>::TMatrixT(Int_t no_rows, Int_t no_cols, const Element *elements, Option_t *option)
62{
63 Allocate(no_rows, no_cols);
64 TMatrixTBase<Element>::SetMatrixArray(elements, option);
65}
66
67////////////////////////////////////////////////////////////////////////////////
68/// array elements are copied
69
70template <class Element>
77
78////////////////////////////////////////////////////////////////////////////////
79/// Copy constructor
80
81template <class Element>
82TMatrixT<Element>::TMatrixT(const TMatrixT<Element> &another) : TMatrixTBase<Element>(another)
83{
84 R__ASSERT(another.IsValid());
85 Allocate(another.GetNrows(), another.GetNcols(), another.GetRowLwb(), another.GetColLwb());
86 *this = another;
87}
88
89////////////////////////////////////////////////////////////////////////////////
90/// Copy constructor of a symmetric matrix
91
92template <class Element>
93TMatrixT<Element>::TMatrixT(const TMatrixTSym<Element> &another)
94{
95 R__ASSERT(another.IsValid());
96 Allocate(another.GetNrows(), another.GetNcols(), another.GetRowLwb(), another.GetColLwb());
97 *this = another;
98}
99
100////////////////////////////////////////////////////////////////////////////////
101/// Copy constructor of a sparse matrix
102
103template <class Element>
104TMatrixT<Element>::TMatrixT(const TMatrixTSparse<Element> &another)
105{
106 R__ASSERT(another.IsValid());
107 Allocate(another.GetNrows(), another.GetNcols(), another.GetRowLwb(), another.GetColLwb());
108 *this = another;
109}
110
111////////////////////////////////////////////////////////////////////////////////
112/// Constructor of matrix applying a specific operation to the prototype.
113/// Example: TMatrixT<Element> a(10,12); ...; TMatrixT<Element> b(TMatrixT::kTransposed, a);
114/// Supported operations are: kZero, kUnit, kTransposed, kInverted and kAtA.
115
116template <class Element>
117TMatrixT<Element>::TMatrixT(EMatrixCreatorsOp1 op, const TMatrixT<Element> &prototype)
118{
119 R__ASSERT(prototype.IsValid());
120
121 switch (op) {
122 case kZero:
123 Allocate(prototype.GetNrows(), prototype.GetNcols(), prototype.GetRowLwb(), prototype.GetColLwb(), 1);
124 break;
125
126 case kUnit:
127 Allocate(prototype.GetNrows(), prototype.GetNcols(), prototype.GetRowLwb(), prototype.GetColLwb(), 1);
128 this->UnitMatrix();
129 break;
130
131 case kTransposed:
132 Allocate(prototype.GetNcols(), prototype.GetNrows(), prototype.GetColLwb(), prototype.GetRowLwb());
133 Transpose(prototype);
134 break;
135
136 case kInverted: {
137 Allocate(prototype.GetNrows(), prototype.GetNcols(), prototype.GetRowLwb(), prototype.GetColLwb(), 1);
138 *this = prototype;
139 // Since the user can not control the tolerance of this newly created matrix
140 // we put it to the smallest possible number
141 const Element oldTol = this->SetTol(std::numeric_limits<Element>::min());
142 this->Invert();
143 this->SetTol(oldTol);
144 break;
145 }
146
147 case kAtA:
148 Allocate(prototype.GetNcols(), prototype.GetNcols(), prototype.GetColLwb(), prototype.GetColLwb(), 1);
149 TMult(prototype, prototype);
150 break;
151
152 default: Error("TMatrixT(EMatrixCreatorOp1)", "operation %d not yet implemented", op);
153 }
154}
155
156////////////////////////////////////////////////////////////////////////////////
157/// Constructor of matrix applying a specific operation to two prototypes.
158/// Example: TMatrixT<Element> a(10,12), b(12,5); ...; TMatrixT<Element> c(a, TMatrixT::kMult, b);
159/// Supported operations are: kMult (a*b), kTransposeMult (a'*b), kInvMult (a^(-1)*b);
160/// Whenever kInvMult is invoked and b is not squared, additional memory is allocated for a^(-1)
161
162template <class Element>
163TMatrixT<Element>::TMatrixT(const TMatrixT<Element> &a, EMatrixCreatorsOp2 op, const TMatrixT<Element> &b)
164{
165 R__ASSERT(a.IsValid());
166 R__ASSERT(b.IsValid());
167
168 switch (op) {
169 case kMult:
170 Allocate(a.GetNrows(), b.GetNcols(), a.GetRowLwb(), b.GetColLwb(), 1);
171 Mult(a, b);
172 break;
173
174 case kTransposeMult:
175 Allocate(a.GetNcols(), b.GetNcols(), a.GetColLwb(), b.GetColLwb(), 1);
176 TMult(a, b);
177 break;
178
179 case kMultTranspose:
180 Allocate(a.GetNrows(), b.GetNrows(), a.GetRowLwb(), b.GetRowLwb(), 1);
181 MultT(a, b);
182 break;
183
184 case kInvMult: {
185 Allocate(a.GetNrows(), b.GetNcols(), a.GetRowLwb(), b.GetColLwb(), 1);
186 // if size(a) == size(b), perform in place computation
187 if (a.GetNrows() == b.GetNcols()) {
188 *this = a;
189 const Element oldTol = this->SetTol(std::numeric_limits<Element>::min());
190 this->Invert();
191 this->SetTol(oldTol);
192 *this *= b;
193 } else {
194 TMatrixT<Element> ainv = a;
195 const Element oldTol = ainv.SetTol(std::numeric_limits<Element>::min());
196 ainv.Invert();
197 ainv.SetTol(oldTol);
198 Mult(ainv, b);
199 }
200 break;
201 }
202
203 case kPlus: {
204 Allocate(a.GetNrows(), a.GetNcols(), a.GetRowLwb(), a.GetColLwb(), 1);
205 Plus(a, b);
206 break;
207 }
208
209 case kMinus: {
210 Allocate(a.GetNrows(), a.GetNcols(), a.GetRowLwb(), a.GetColLwb(), 1);
211 Minus(a, b);
212 break;
213 }
214
215 default: Error("TMatrixT(EMatrixCreatorOp2)", "operation %d not yet implemented", op);
216 }
217}
218
219////////////////////////////////////////////////////////////////////////////////
220/// Constructor of matrix applying a specific operation to two prototypes.
221/// Example: TMatrixT<Element> a(10,12), b(12,5); ...; TMatrixT<Element> c(a, TMatrixT::kMult, b);
222/// Supported operations are: kMult (a*b), kTransposeMult (a'*b), kInvMult (a^(-1)*b)
223
224template <class Element>
225TMatrixT<Element>::TMatrixT(const TMatrixT<Element> &a, EMatrixCreatorsOp2 op, const TMatrixTSym<Element> &b)
226{
227 R__ASSERT(a.IsValid());
228 R__ASSERT(b.IsValid());
229
230 switch (op) {
231 case kMult:
232 Allocate(a.GetNrows(), b.GetNcols(), a.GetRowLwb(), b.GetColLwb(), 1);
233 Mult(a, b);
234 break;
235
236 case kTransposeMult:
237 Allocate(a.GetNcols(), b.GetNcols(), a.GetColLwb(), b.GetColLwb(), 1);
238 TMult(a, b);
239 break;
240
241 case kMultTranspose:
242 Allocate(a.GetNrows(), b.GetNrows(), a.GetRowLwb(), b.GetRowLwb(), 1);
243 MultT(a, b);
244 break;
245
246 case kInvMult: {
247 Allocate(a.GetNrows(), a.GetNcols(), a.GetRowLwb(), a.GetColLwb(), 1);
248 *this = a;
249 const Element oldTol = this->SetTol(std::numeric_limits<Element>::min());
250 this->Invert();
251 this->SetTol(oldTol);
252 *this *= b;
253 break;
254 }
255
256 case kPlus: {
257 Allocate(a.GetNrows(), a.GetNcols(), a.GetRowLwb(), a.GetColLwb(), 1);
258 Plus(a, b);
259 break;
260 }
261
262 case kMinus: {
263 Allocate(a.GetNrows(), a.GetNcols(), a.GetRowLwb(), a.GetColLwb(), 1);
264 Minus(a, b);
265 break;
266 }
267
268 default: Error("TMatrixT(EMatrixCreatorOp2)", "operation %d not yet implemented", op);
269 }
270}
271
272////////////////////////////////////////////////////////////////////////////////
273/// Constructor of matrix applying a specific operation to two prototypes.
274/// Example: TMatrixT<Element> a(10,12), b(12,5); ...; TMatrixT<Element> c(a, TMatrixT::kMult, b);
275/// Supported operations are: kMult (a*b), kTransposeMult (a'*b), kInvMult (a^(-1)*b)
276/// Whenever kInvMult is invoked and b is not squared, additional memory is allocated for a^(-1)
277
278template <class Element>
279TMatrixT<Element>::TMatrixT(const TMatrixTSym<Element> &a, EMatrixCreatorsOp2 op, const TMatrixT<Element> &b)
280{
281 R__ASSERT(a.IsValid());
282 R__ASSERT(b.IsValid());
283
284 switch (op) {
285 case kMult:
286 Allocate(a.GetNrows(), b.GetNcols(), a.GetRowLwb(), b.GetColLwb(), 1);
287 Mult(a, b);
288 break;
289
290 case kTransposeMult:
291 Allocate(a.GetNcols(), b.GetNcols(), a.GetColLwb(), b.GetColLwb(), 1);
292 TMult(a, b);
293 break;
294
295 case kMultTranspose:
296 Allocate(a.GetNrows(), b.GetNrows(), a.GetRowLwb(), b.GetRowLwb(), 1);
297 MultT(a, b);
298 break;
299
300 case kInvMult: {
301 Allocate(a.GetNrows(), b.GetNcols(), a.GetRowLwb(), b.GetColLwb(), 1);
302 // if size(a) == size(b), perform in place computation
303 if (a.GetNrows() == b.GetNcols()) {
304 *this = a;
305 const Element oldTol = this->SetTol(std::numeric_limits<Element>::min());
306 this->Invert();
307 this->SetTol(oldTol);
308 *this *= b;
309 } else {
310 TMatrixTSym<Element> ainv = a;
311 const Element oldTol = ainv.SetTol(std::numeric_limits<Element>::min());
312 ainv.Invert();
313 ainv.SetTol(oldTol);
314 Mult(ainv, b);
315 }
316 break;
317 }
318
319 case kPlus: {
320 Allocate(a.GetNrows(), a.GetNcols(), a.GetRowLwb(), a.GetColLwb(), 1);
321 Plus(a, b);
322 break;
323 }
324
325 case kMinus: {
326 Allocate(a.GetNrows(), a.GetNcols(), a.GetRowLwb(), a.GetColLwb(), 1);
327 Minus(a, b);
328 break;
329 }
330
331 default: Error("TMatrixT(EMatrixCreatorOp2)", "operation %d not yet implemented", op);
332 }
333}
334
335////////////////////////////////////////////////////////////////////////////////
336/// Constructor of matrix applying a specific operation to two prototypes.
337/// Example: TMatrixT<Element> a(10,12), b(12,5); ...; TMatrixT<Element> c(a, TMatrixT::kMult, b);
338/// Supported operations are: kMult (a*b), kTransposeMult (a'*b), kInvMult (a^(-1)*b)
339
340template <class Element>
341TMatrixT<Element>::TMatrixT(const TMatrixTSym<Element> &a, EMatrixCreatorsOp2 op, const TMatrixTSym<Element> &b)
342{
343 R__ASSERT(a.IsValid());
344 R__ASSERT(b.IsValid());
345
346 switch (op) {
347 case kMult:
348 Allocate(a.GetNrows(), b.GetNcols(), a.GetRowLwb(), b.GetColLwb(), 1);
349 Mult(a, b);
350 break;
351
352 case kTransposeMult:
353 Allocate(a.GetNcols(), b.GetNcols(), a.GetColLwb(), b.GetColLwb(), 1);
354 TMult(a, b);
355 break;
356
357 case kMultTranspose:
358 Allocate(a.GetNrows(), b.GetNrows(), a.GetRowLwb(), b.GetRowLwb(), 1);
359 MultT(a, b);
360 break;
361
362 case kInvMult: {
363 Allocate(a.GetNrows(), a.GetNcols(), a.GetRowLwb(), a.GetColLwb(), 1);
364 *this = a;
365 const Element oldTol = this->SetTol(std::numeric_limits<Element>::min());
366 this->Invert();
367 this->SetTol(oldTol);
368 *this *= b;
369 break;
370 }
371
372 case kPlus: {
373 Allocate(a.GetNrows(), a.GetNcols(), a.GetRowLwb(), a.GetColLwb(), 1);
374 Plus(*dynamic_cast<const TMatrixT<Element> *>(&a), b);
375 break;
376 }
377
378 case kMinus: {
379 Allocate(a.GetNrows(), a.GetNcols(), a.GetRowLwb(), a.GetColLwb(), 1);
380 Minus(*dynamic_cast<const TMatrixT<Element> *>(&a), b);
381 break;
382 }
383
384 default: Error("TMatrixT(EMatrixCreatorOp2)", "operation %d not yet implemented", op);
385 }
386}
387
388////////////////////////////////////////////////////////////////////////////////
389/// Constructor using the TMatrixTLazy class
390
391template <class Element>
392TMatrixT<Element>::TMatrixT(const TMatrixTLazy<Element> &lazy_constructor)
393{
394 Allocate(lazy_constructor.GetRowUpb() - lazy_constructor.GetRowLwb() + 1,
395 lazy_constructor.GetColUpb() - lazy_constructor.GetColLwb() + 1, lazy_constructor.GetRowLwb(),
396 lazy_constructor.GetColLwb(), 1);
397 lazy_constructor.FillIn(*this);
398}
399
400////////////////////////////////////////////////////////////////////////////////
401/// Delete data pointer m, if it was assigned on the heap
402
403template <class Element>
404void TMatrixT<Element>::Delete_m(Int_t size, Element *&m)
405{
406 if (m) {
407 if (size > this->kSizeMax)
408 delete[] m;
409 m = nullptr;
410 }
411}
412
413////////////////////////////////////////////////////////////////////////////////
414/// Return data pointer . if requested size <= kSizeMax, assign pointer
415/// to the stack space
416
417template <class Element>
418Element *TMatrixT<Element>::New_m(Int_t size)
419{
420 if (size == 0)
421 return nullptr;
422 else {
423 if (size <= this->kSizeMax)
424 return fDataStack;
425 else {
426 Element *heap = new Element[size];
427 return heap;
428 }
429 }
430}
431
432////////////////////////////////////////////////////////////////////////////////
433/// Copy copySize doubles from *oldp to *newp . However take care of the
434/// situation where both pointers are assigned to the same stack space
435
436template <class Element>
437Int_t TMatrixT<Element>::Memcpy_m(Element *newp, const Element *oldp, Int_t copySize, Int_t newSize, Int_t oldSize)
438{
439 if (copySize == 0 || oldp == newp)
440 return 0;
441 else {
442 if (newSize <= this->kSizeMax && oldSize <= this->kSizeMax) {
443 // both pointers are inside fDataStack, be careful with copy direction !
444 if (newp > oldp) {
445 for (Int_t i = copySize - 1; i >= 0; i--)
446 newp[i] = oldp[i];
447 } else {
448 for (Int_t i = 0; i < copySize; i++)
449 newp[i] = oldp[i];
450 }
451 } else
452 memcpy(newp, oldp, copySize * sizeof(Element));
453 }
454 return 0;
455}
456
457////////////////////////////////////////////////////////////////////////////////
458/// Allocate new matrix. Arguments are number of rows, columns, row
459/// lowerbound (0 default) and column lowerbound (0 default).
460
461template <class Element>
462void TMatrixT<Element>::Allocate(Int_t no_rows, Int_t no_cols, Int_t row_lwb, Int_t col_lwb, Int_t init,
463 Int_t /*nr_nonzeros*/)
464{
465 this->fIsOwner = kTRUE;
466 this->fTol = std::numeric_limits<Element>::epsilon();
467 fElements = nullptr;
468 this->fNrows = 0;
469 this->fNcols = 0;
470 this->fRowLwb = 0;
471 this->fColLwb = 0;
472 this->fNelems = 0;
473
474 if (no_rows < 0 || no_cols < 0) {
475 Error("Allocate", "no_rows=%d no_cols=%d", no_rows, no_cols);
476 this->Invalidate();
477 return;
478 }
479
480 this->MakeValid();
481 this->fNrows = no_rows;
482 this->fNcols = no_cols;
483 this->fRowLwb = row_lwb;
484 this->fColLwb = col_lwb;
485 this->fNelems = this->fNrows * this->fNcols;
486
487 // Check if fNelems does not have an overflow.
488 if (((Long64_t)this->fNrows) * this->fNcols != this->fNelems) {
489 Error("Allocate", "too large: no_rows=%d no_cols=%d", no_rows, no_cols);
490 this->Invalidate();
491 return;
492 }
493
494 if (this->fNelems > 0) {
495 fElements = New_m(this->fNelems);
496 if (init)
497 memset(fElements, 0, this->fNelems * sizeof(Element));
498 } else
499 fElements = nullptr;
500}
501
502////////////////////////////////////////////////////////////////////////////////
503/// General matrix summation. Replace this matrix with C such that C = A + B.
504
505template <class Element>
506void TMatrixT<Element>::Plus(const TMatrixT<Element> &a, const TMatrixT<Element> &b)
507{
508 if (gMatrixCheck) {
509 if (!AreCompatible(a, b)) {
510 Error("Plus", "matrices not compatible");
511 return;
512 }
513
514 if (this->GetMatrixArray() == a.GetMatrixArray()) {
515 Error("Plus", "this->GetMatrixArray() == a.GetMatrixArray()");
516 return;
517 }
518
519 if (this->GetMatrixArray() == b.GetMatrixArray()) {
520 Error("Plus", "this->GetMatrixArray() == b.GetMatrixArray()");
521 return;
522 }
523 }
524
525 const Element *ap = a.GetMatrixArray();
526 const Element *bp = b.GetMatrixArray();
527 Element *cp = this->GetMatrixArray();
528 const Element *const cp_last = cp + this->fNelems;
529
530 while (cp < cp_last) {
531 *cp = *ap++ + *bp++;
532 cp++;
533 }
534}
535
536////////////////////////////////////////////////////////////////////////////////
537/// General matrix summation. Replace this matrix with C such that C = A + B.
538
539template <class Element>
540void TMatrixT<Element>::Plus(const TMatrixT<Element> &a, const TMatrixTSym<Element> &b)
541{
542 if (gMatrixCheck) {
543 if (!AreCompatible(a, b)) {
544 Error("Plus", "matrices not compatible");
545 return;
546 }
547
548 if (this->GetMatrixArray() == a.GetMatrixArray()) {
549 Error("Plus", "this->GetMatrixArray() == a.GetMatrixArray()");
550 return;
551 }
552
553 if (this->GetMatrixArray() == b.GetMatrixArray()) {
554 Error("Plus", "this->GetMatrixArray() == b.GetMatrixArray()");
555 return;
556 }
557 }
558
559 const Element *ap = a.GetMatrixArray();
560 const Element *bp = b.GetMatrixArray();
561 Element *cp = this->GetMatrixArray();
562 const Element *const cp_last = cp + this->fNelems;
563
564 while (cp < cp_last) {
565 *cp = *ap++ + *bp++;
566 cp++;
567 }
568}
569
570////////////////////////////////////////////////////////////////////////////////
571/// General matrix subtraction. Replace this matrix with C such that C = A - B.
572
573template <class Element>
574void TMatrixT<Element>::Minus(const TMatrixT<Element> &a, const TMatrixT<Element> &b)
575{
576 if (gMatrixCheck) {
577 if (!AreCompatible(a, b)) {
578 Error("Minus", "matrices not compatible");
579 return;
580 }
581
582 if (this->GetMatrixArray() == a.GetMatrixArray()) {
583 Error("Minus", "this->GetMatrixArray() == a.GetMatrixArray()");
584 return;
585 }
586
587 if (this->GetMatrixArray() == b.GetMatrixArray()) {
588 Error("Minus", "this->GetMatrixArray() == b.GetMatrixArray()");
589 return;
590 }
591 }
592
593 const Element *ap = a.GetMatrixArray();
594 const Element *bp = b.GetMatrixArray();
595 Element *cp = this->GetMatrixArray();
596 const Element *const cp_last = cp + this->fNelems;
597
598 while (cp < cp_last) {
599 *cp = *ap++ - *bp++;
600 cp++;
601 }
602}
603
604////////////////////////////////////////////////////////////////////////////////
605/// General matrix subtraction. Replace this matrix with C such that C = A - B.
606
607template <class Element>
608void TMatrixT<Element>::Minus(const TMatrixT<Element> &a, const TMatrixTSym<Element> &b)
609{
610 if (gMatrixCheck) {
611 if (!AreCompatible(a, b)) {
612 Error("Minus", "matrices not compatible");
613 return;
614 }
615
616 if (this->GetMatrixArray() == a.GetMatrixArray()) {
617 Error("Minus", "this->GetMatrixArray() == a.GetMatrixArray()");
618 return;
619 }
620
621 if (this->GetMatrixArray() == b.GetMatrixArray()) {
622 Error("Minus", "this->GetMatrixArray() == b.GetMatrixArray()");
623 return;
624 }
625 }
626
627 const Element *ap = a.GetMatrixArray();
628 const Element *bp = b.GetMatrixArray();
629 Element *cp = this->GetMatrixArray();
630 const Element *const cp_last = cp + this->fNelems;
631
632 while (cp < cp_last) {
633 *cp = *ap++ - *bp++;
634 cp++;
635 }
636}
637
638////////////////////////////////////////////////////////////////////////////////
639/// General matrix multiplication. Replace this matrix with C such that C = A * B.
640
641template <class Element>
642void TMatrixT<Element>::Mult(const TMatrixT<Element> &a, const TMatrixT<Element> &b)
643{
644 if (gMatrixCheck) {
645 if (a.GetNcols() != b.GetNrows() || a.GetColLwb() != b.GetRowLwb()) {
646 Error("Mult", "A rows and B columns incompatible");
647 return;
648 }
649
650 if (this->GetMatrixArray() == a.GetMatrixArray()) {
651 Error("Mult", "this->GetMatrixArray() == a.GetMatrixArray()");
652 return;
653 }
654
655 if (this->GetMatrixArray() == b.GetMatrixArray()) {
656 Error("Mult", "this->GetMatrixArray() == b.GetMatrixArray()");
657 return;
658 }
659 }
660
661#ifdef CBLAS
662 const Element *ap = a.GetMatrixArray();
663 const Element *bp = b.GetMatrixArray();
664 Element *cp = this->GetMatrixArray();
665 if (typeid(Element) == typeid(Double_t))
666 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, fNrows, fNcols, a.GetNcols(), 1.0, ap, a.GetNcols(), bp,
667 b.GetNcols(), 1.0, cp, fNcols);
668 else if (typeid(Element) != typeid(Float_t))
669 cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, fNrows, fNcols, a.GetNcols(), 1.0, ap, a.GetNcols(), bp,
670 b.GetNcols(), 1.0, cp, fNcols);
671 else
672 Error("Mult", "type %s not implemented in BLAS library", typeid(Element));
673#else
674 const Int_t na = a.GetNoElements();
675 const Int_t nb = b.GetNoElements();
676 const Int_t ncolsa = a.GetNcols();
677 const Int_t ncolsb = b.GetNcols();
678 const Element *const ap = a.GetMatrixArray();
679 const Element *const bp = b.GetMatrixArray();
680 Element *cp = this->GetMatrixArray();
681
682 AMultB(ap, na, ncolsa, bp, nb, ncolsb, cp);
683#endif
684}
685
686////////////////////////////////////////////////////////////////////////////////
687/// Matrix multiplication, with A symmetric and B general.
688/// Replace this matrix with C such that C = A * B.
689
690template <class Element>
691void TMatrixT<Element>::Mult(const TMatrixTSym<Element> &a, const TMatrixT<Element> &b)
692{
693 if (gMatrixCheck) {
694 R__ASSERT(a.IsValid());
695 R__ASSERT(b.IsValid());
696 if (a.GetNcols() != b.GetNrows() || a.GetColLwb() != b.GetRowLwb()) {
697 Error("Mult", "A rows and B columns incompatible");
698 return;
699 }
700
701 if (this->GetMatrixArray() == a.GetMatrixArray()) {
702 Error("Mult", "this->GetMatrixArray() == a.GetMatrixArray()");
703 return;
704 }
705
706 if (this->GetMatrixArray() == b.GetMatrixArray()) {
707 Error("Mult", "this->GetMatrixArray() == b.GetMatrixArray()");
708 return;
709 }
710 }
711
712#ifdef CBLAS
713 const Element *ap = a.GetMatrixArray();
714 const Element *bp = b.GetMatrixArray();
715 Element *cp = this->GetMatrixArray();
716 if (typeid(Element) == typeid(Double_t))
717 cblas_dsymm(CblasRowMajor, CblasLeft, CblasUpper, fNrows, fNcols, 1.0, ap, a.GetNcols(), bp, b.GetNcols(), 0.0,
718 cp, fNcols);
719 else if (typeid(Element) != typeid(Float_t))
720 cblas_ssymm(CblasRowMajor, CblasLeft, CblasUpper, fNrows, fNcols, 1.0, ap, a.GetNcols(), bp, b.GetNcols(), 0.0,
721 cp, fNcols);
722 else
723 Error("Mult", "type %s not implemented in BLAS library", typeid(Element));
724#else
725 const Int_t na = a.GetNoElements();
726 const Int_t nb = b.GetNoElements();
727 const Int_t ncolsa = a.GetNcols();
728 const Int_t ncolsb = b.GetNcols();
729 const Element *const ap = a.GetMatrixArray();
730 const Element *const bp = b.GetMatrixArray();
731 Element *cp = this->GetMatrixArray();
732
733 AMultB(ap, na, ncolsa, bp, nb, ncolsb, cp);
734
735#endif
736}
737
738////////////////////////////////////////////////////////////////////////////////
739/// Matrix multiplication, with A general and B symmetric.
740/// Replace this matrix with C such that C = A * B.
741
742template <class Element>
743void TMatrixT<Element>::Mult(const TMatrixT<Element> &a, const TMatrixTSym<Element> &b)
744{
745 if (gMatrixCheck) {
746 R__ASSERT(a.IsValid());
747 R__ASSERT(b.IsValid());
748 if (a.GetNcols() != b.GetNrows() || a.GetColLwb() != b.GetRowLwb()) {
749 Error("Mult", "A rows and B columns incompatible");
750 return;
751 }
752
753 if (this->GetMatrixArray() == a.GetMatrixArray()) {
754 Error("Mult", "this->GetMatrixArray() == a.GetMatrixArray()");
755 return;
756 }
757
758 if (this->GetMatrixArray() == b.GetMatrixArray()) {
759 Error("Mult", "this->GetMatrixArray() == b.GetMatrixArray()");
760 return;
761 }
762 }
763
764#ifdef CBLAS
765 const Element *ap = a.GetMatrixArray();
766 const Element *bp = b.GetMatrixArray();
767 Element *cp = this->GetMatrixArray();
768 if (typeid(Element) == typeid(Double_t))
769 cblas_dsymm(CblasRowMajor, CblasRight, CblasUpper, fNrows, fNcols, 1.0, bp, b.GetNcols(), ap, a.GetNcols(), 0.0,
770 cp, fNcols);
771 else if (typeid(Element) != typeid(Float_t))
772 cblas_ssymm(CblasRowMajor, CblasRight, CblasUpper, fNrows, fNcols, 1.0, bp, b.GetNcols(), ap, a.GetNcols(), 0.0,
773 cp, fNcols);
774 else
775 Error("Mult", "type %s not implemented in BLAS library", typeid(Element));
776#else
777 const Int_t na = a.GetNoElements();
778 const Int_t nb = b.GetNoElements();
779 const Int_t ncolsa = a.GetNcols();
780 const Int_t ncolsb = b.GetNcols();
781 const Element *const ap = a.GetMatrixArray();
782 const Element *const bp = b.GetMatrixArray();
783 Element *cp = this->GetMatrixArray();
784
785 AMultB(ap, na, ncolsa, bp, nb, ncolsb, cp);
786#endif
787}
788
789////////////////////////////////////////////////////////////////////////////////
790/// Matrix multiplication, with A symmetric and B symmetric.
791/// (Actually copied for the moment routine for B general)
792/// Replace this matrix with C such that C = A * B.
793
794template <class Element>
795void TMatrixT<Element>::Mult(const TMatrixTSym<Element> &a, const TMatrixTSym<Element> &b)
796{
797 if (gMatrixCheck) {
798 R__ASSERT(a.IsValid());
799 R__ASSERT(b.IsValid());
800 if (a.GetNcols() != b.GetNrows() || a.GetColLwb() != b.GetRowLwb()) {
801 Error("Mult", "A rows and B columns incompatible");
802 return;
803 }
804
805 if (this->GetMatrixArray() == a.GetMatrixArray()) {
806 Error("Mult", "this->GetMatrixArray() == a.GetMatrixArray()");
807 return;
808 }
809
810 if (this->GetMatrixArray() == b.GetMatrixArray()) {
811 Error("Mult", "this->GetMatrixArray() == b.GetMatrixArray()");
812 return;
813 }
814 }
815
816#ifdef CBLAS
817 const Element *ap = a.GetMatrixArray();
818 const Element *bp = b.GetMatrixArray();
819 Element *cp = this->GetMatrixArray();
820 if (typeid(Element) == typeid(Double_t))
821 cblas_dsymm(CblasRowMajor, CblasLeft, CblasUpper, fNrows, fNcols, 1.0, ap, a.GetNcols(), bp, b.GetNcols(), 0.0,
822 cp, fNcols);
823 else if (typeid(Element) != typeid(Float_t))
824 cblas_ssymm(CblasRowMajor, CblasLeft, CblasUpper, fNrows, fNcols, 1.0, ap, a.GetNcols(), bp, b.GetNcols(), 0.0,
825 cp, fNcols);
826 else
827 Error("Mult", "type %s not implemented in BLAS library", typeid(Element));
828#else
829 const Int_t na = a.GetNoElements();
830 const Int_t nb = b.GetNoElements();
831 const Int_t ncolsa = a.GetNcols();
832 const Int_t ncolsb = b.GetNcols();
833 const Element *const ap = a.GetMatrixArray();
834 const Element *const bp = b.GetMatrixArray();
835 Element *cp = this->GetMatrixArray();
836
837 AMultB(ap, na, ncolsa, bp, nb, ncolsb, cp);
838#endif
839}
840
841////////////////////////////////////////////////////////////////////////////////
842/// Replace this matrix with C such that C = A' * B. In other words,
843/// c[i,j] = SUM{ a[k,i] * b[k,j] }.
844
845template <class Element>
846void TMatrixT<Element>::TMult(const TMatrixT<Element> &a, const TMatrixT<Element> &b)
847{
848 if (gMatrixCheck) {
849 R__ASSERT(a.IsValid());
850 R__ASSERT(b.IsValid());
851 if (a.GetNrows() != b.GetNrows() || a.GetRowLwb() != b.GetRowLwb()) {
852 Error("TMult", "A rows and B columns incompatible");
853 return;
854 }
855
856 if (this->GetMatrixArray() == a.GetMatrixArray()) {
857 Error("TMult", "this->GetMatrixArray() == a.GetMatrixArray()");
858 return;
859 }
860
861 if (this->GetMatrixArray() == b.GetMatrixArray()) {
862 Error("TMult", "this->GetMatrixArray() == b.GetMatrixArray()");
863 return;
864 }
865 }
866
867#ifdef CBLAS
868 const Element *ap = a.GetMatrixArray();
869 const Element *bp = b.GetMatrixArray();
870 Element *cp = this->GetMatrixArray();
871 if (typeid(Element) == typeid(Double_t))
872 cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans, this->fNrows, this->fNcols, a.GetNrows(), 1.0, ap,
873 a.GetNcols(), bp, b.GetNcols(), 1.0, cp, this->fNcols);
874 else if (typeid(Element) != typeid(Float_t))
875 cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans, fNrows, fNcols, a.GetNrows(), 1.0, ap, a.GetNcols(), bp,
876 b.GetNcols(), 1.0, cp, fNcols);
877 else
878 Error("TMult", "type %s not implemented in BLAS library", typeid(Element));
879#else
880 const Int_t nb = b.GetNoElements();
881 const Int_t ncolsa = a.GetNcols();
882 const Int_t ncolsb = b.GetNcols();
883 const Element *const ap = a.GetMatrixArray();
884 const Element *const bp = b.GetMatrixArray();
885 Element *cp = this->GetMatrixArray();
886
887 AtMultB(ap, ncolsa, bp, nb, ncolsb, cp);
888#endif
889}
890
891////////////////////////////////////////////////////////////////////////////////
892/// Replace this matrix with C such that C = A' * B. In other words,
893/// c[i,j] = SUM{ a[k,i] * b[k,j] }.
894
895template <class Element>
896void TMatrixT<Element>::TMult(const TMatrixT<Element> &a, const TMatrixTSym<Element> &b)
897{
898 if (gMatrixCheck) {
899 R__ASSERT(a.IsValid());
900 R__ASSERT(b.IsValid());
901 if (a.GetNrows() != b.GetNrows() || a.GetRowLwb() != b.GetRowLwb()) {
902 Error("TMult", "A rows and B columns incompatible");
903 return;
904 }
905
906 if (this->GetMatrixArray() == a.GetMatrixArray()) {
907 Error("TMult", "this->GetMatrixArray() == a.GetMatrixArray()");
908 return;
909 }
910
911 if (this->GetMatrixArray() == b.GetMatrixArray()) {
912 Error("TMult", "this->GetMatrixArray() == b.GetMatrixArray()");
913 return;
914 }
915 }
916
917#ifdef CBLAS
918 const Element *ap = a.GetMatrixArray();
919 const Element *bp = b.GetMatrixArray();
920 Element *cp = this->GetMatrixArray();
921 if (typeid(Element) == typeid(Double_t))
922 cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans, fNrows, fNcols, a.GetNrows(), 1.0, ap, a.GetNcols(), bp,
923 b.GetNcols(), 1.0, cp, fNcols);
924 else if (typeid(Element) != typeid(Float_t))
925 cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans, fNrows, fNcols, a.GetNrows(), 1.0, ap, a.GetNcols(), bp,
926 b.GetNcols(), 1.0, cp, fNcols);
927 else
928 Error("TMult", "type %s not implemented in BLAS library", typeid(Element));
929#else
930 const Int_t nb = b.GetNoElements();
931 const Int_t ncolsa = a.GetNcols();
932 const Int_t ncolsb = b.GetNcols();
933 const Element *const ap = a.GetMatrixArray();
934 const Element *const bp = b.GetMatrixArray();
935 Element *cp = this->GetMatrixArray();
936
937 AtMultB(ap, ncolsa, bp, nb, ncolsb, cp);
938#endif
939}
940
941////////////////////////////////////////////////////////////////////////////////
942/// General matrix multiplication. Replace this matrix with C such that C = A * B^T.
943
944template <class Element>
945void TMatrixT<Element>::MultT(const TMatrixT<Element> &a, const TMatrixT<Element> &b)
946{
947 if (gMatrixCheck) {
948 R__ASSERT(a.IsValid());
949 R__ASSERT(b.IsValid());
950
951 if (a.GetNcols() != b.GetNcols() || a.GetColLwb() != b.GetColLwb()) {
952 Error("MultT", "A rows and B columns incompatible");
953 return;
954 }
955
956 if (this->GetMatrixArray() == a.GetMatrixArray()) {
957 Error("MultT", "this->GetMatrixArray() == a.GetMatrixArray()");
958 return;
959 }
960
961 if (this->GetMatrixArray() == b.GetMatrixArray()) {
962 Error("MultT", "this->GetMatrixArray() == b.GetMatrixArray()");
963 return;
964 }
965 }
966
967#ifdef CBLAS
968 const Element *ap = a.GetMatrixArray();
969 const Element *bp = b.GetMatrixArray();
970 Element *cp = this->GetMatrixArray();
971 if (typeid(Element) == typeid(Double_t))
972 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans, fNrows, fNcols, a.GetNcols(), 1.0, ap, a.GetNcols(), bp,
973 b.GetNcols(), 1.0, cp, fNcols);
974 else if (typeid(Element) != typeid(Float_t))
975 cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasTrans, fNrows, fNcols, a.GetNcols(), 1.0, ap, a.GetNcols(), bp,
976 b.GetNcols(), 1.0, cp, fNcols);
977 else
978 Error("MultT", "type %s not implemented in BLAS library", typeid(Element));
979#else
980 const Int_t na = a.GetNoElements();
981 const Int_t nb = b.GetNoElements();
982 const Int_t ncolsa = a.GetNcols();
983 const Int_t ncolsb = b.GetNcols();
984 const Element *const ap = a.GetMatrixArray();
985 const Element *const bp = b.GetMatrixArray();
986 Element *cp = this->GetMatrixArray();
987
988 AMultBt(ap, na, ncolsa, bp, nb, ncolsb, cp);
989#endif
990}
991
992////////////////////////////////////////////////////////////////////////////////
993/// Matrix multiplication, with A symmetric and B general.
994/// Replace this matrix with C such that C = A * B^T.
995
996template <class Element>
997void TMatrixT<Element>::MultT(const TMatrixTSym<Element> &a, const TMatrixT<Element> &b)
998{
999 if (gMatrixCheck) {
1000 R__ASSERT(a.IsValid());
1001 R__ASSERT(b.IsValid());
1002 if (a.GetNcols() != b.GetNcols() || a.GetColLwb() != b.GetColLwb()) {
1003 Error("MultT", "A rows and B columns incompatible");
1004 return;
1005 }
1006
1007 if (this->GetMatrixArray() == a.GetMatrixArray()) {
1008 Error("MultT", "this->GetMatrixArray() == a.GetMatrixArray()");
1009 return;
1010 }
1011
1012 if (this->GetMatrixArray() == b.GetMatrixArray()) {
1013 Error("MultT", "this->GetMatrixArray() == b.GetMatrixArray()");
1014 return;
1015 }
1016 }
1017
1018#ifdef CBLAS
1019 const Element *ap = a.GetMatrixArray();
1020 const Element *bp = b.GetMatrixArray();
1021 Element *cp = this->GetMatrixArray();
1022 if (typeid(Element) == typeid(Double_t))
1023 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans, this->fNrows, this->fNcols, a.GetNcols(), 1.0, ap,
1024 a.GetNcols(), bp, b.GetNcols(), 1.0, cp, this->fNcols);
1025 else if (typeid(Element) != typeid(Float_t))
1026 cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasTrans, fNrows, fNcols, a.GetNcols(), 1.0, ap, a.GetNcols(), bp,
1027 b.GetNcols(), 1.0, cp, fNcols);
1028 else
1029 Error("MultT", "type %s not implemented in BLAS library", typeid(Element));
1030#else
1031 const Int_t na = a.GetNoElements();
1032 const Int_t nb = b.GetNoElements();
1033 const Int_t ncolsa = a.GetNcols();
1034 const Int_t ncolsb = b.GetNcols();
1035 const Element *const ap = a.GetMatrixArray();
1036 const Element *const bp = b.GetMatrixArray();
1037 Element *cp = this->GetMatrixArray();
1038
1039 AMultBt(ap, na, ncolsa, bp, nb, ncolsb, cp);
1040#endif
1041}
1042
1043////////////////////////////////////////////////////////////////////////////////
1044/// Use the array data to fill the matrix ([row_lwb..row_upb] x [col_lwb..col_upb])
1045
1046template <class Element>
1047TMatrixT<Element> &TMatrixT<Element>::Use(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb, Element *data)
1048{
1049 if (gMatrixCheck) {
1050 if (row_upb < row_lwb) {
1051 Error("Use", "row_upb=%d < row_lwb=%d", row_upb, row_lwb);
1052 return *this;
1053 }
1054 if (col_upb < col_lwb) {
1055 Error("Use", "col_upb=%d < col_lwb=%d", col_upb, col_lwb);
1056 return *this;
1057 }
1058 }
1059
1060 Clear();
1061 this->fNrows = row_upb - row_lwb + 1;
1062 this->fNcols = col_upb - col_lwb + 1;
1063 this->fRowLwb = row_lwb;
1064 this->fColLwb = col_lwb;
1065 this->fNelems = this->fNrows * this->fNcols;
1066 fElements = data;
1067 this->fIsOwner = kFALSE;
1068
1069 return *this;
1070}
1071
1072////////////////////////////////////////////////////////////////////////////////
1073/// Get submatrix [row_lwb..row_upb] x [col_lwb..col_upb]; The indexing range of the
1074/// returned matrix depends on the argument option:
1075///
1076/// option == "S" : return [0..row_upb-row_lwb][0..col_upb-col_lwb] (default)
1077/// else : return [row_lwb..row_upb][col_lwb..col_upb]
1078
1079template <class Element>
1080TMatrixTBase<Element> &TMatrixT<Element>::GetSub(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb,
1081 TMatrixTBase<Element> &target, Option_t *option) const
1082{
1083 if (gMatrixCheck) {
1084 R__ASSERT(this->IsValid());
1085 if (row_lwb < this->fRowLwb || row_lwb > this->fRowLwb + this->fNrows - 1) {
1086 Error("GetSub", "row_lwb out of bounds");
1087 return target;
1088 }
1089 if (col_lwb < this->fColLwb || col_lwb > this->fColLwb + this->fNcols - 1) {
1090 Error("GetSub", "col_lwb out of bounds");
1091 return target;
1092 }
1093 if (row_upb < this->fRowLwb || row_upb > this->fRowLwb + this->fNrows - 1) {
1094 Error("GetSub", "row_upb out of bounds");
1095 return target;
1096 }
1097 if (col_upb < this->fColLwb || col_upb > this->fColLwb + this->fNcols - 1) {
1098 Error("GetSub", "col_upb out of bounds");
1099 return target;
1100 }
1101 if (row_upb < row_lwb || col_upb < col_lwb) {
1102 Error("GetSub", "row_upb < row_lwb || col_upb < col_lwb");
1103 return target;
1104 }
1105 }
1106
1107 TString opt(option);
1108 opt.ToUpper();
1109 const Int_t shift = (opt.Contains("S")) ? 1 : 0;
1110
1111 const Int_t row_lwb_sub = (shift) ? 0 : row_lwb;
1112 const Int_t row_upb_sub = (shift) ? row_upb - row_lwb : row_upb;
1113 const Int_t col_lwb_sub = (shift) ? 0 : col_lwb;
1114 const Int_t col_upb_sub = (shift) ? col_upb - col_lwb : col_upb;
1115
1116 target.ResizeTo(row_lwb_sub, row_upb_sub, col_lwb_sub, col_upb_sub);
1117 const Int_t nrows_sub = row_upb_sub - row_lwb_sub + 1;
1118 const Int_t ncols_sub = col_upb_sub - col_lwb_sub + 1;
1119
1120 if (target.GetRowIndexArray() && target.GetColIndexArray()) {
1121 for (Int_t irow = 0; irow < nrows_sub; irow++) {
1122 for (Int_t icol = 0; icol < ncols_sub; icol++) {
1123 target(irow + row_lwb_sub, icol + col_lwb_sub) = (*this)(row_lwb + irow, col_lwb + icol);
1124 }
1125 }
1126 } else {
1127 const Element *ap = this->GetMatrixArray() + (row_lwb - this->fRowLwb) * this->fNcols + (col_lwb - this->fColLwb);
1128 Element *bp = target.GetMatrixArray();
1129
1130 for (Int_t irow = 0; irow < nrows_sub; irow++) {
1131 const Element *ap_sub = ap;
1132 for (Int_t icol = 0; icol < ncols_sub; icol++) {
1133 *bp++ = *ap_sub++;
1134 }
1135 ap += this->fNcols;
1136 }
1137 }
1138
1139 return target;
1140}
1141
1142////////////////////////////////////////////////////////////////////////////////
1143/// Insert matrix source starting at [row_lwb][col_lwb], thereby overwriting the part
1144/// [row_lwb..row_lwb+nrows_source][col_lwb..col_lwb+ncols_source];
1145
1146template <class Element>
1147TMatrixTBase<Element> &TMatrixT<Element>::SetSub(Int_t row_lwb, Int_t col_lwb, const TMatrixTBase<Element> &source)
1148{
1149 if (gMatrixCheck) {
1150 R__ASSERT(this->IsValid());
1151 R__ASSERT(source.IsValid());
1152
1153 if (row_lwb < this->fRowLwb || row_lwb > this->fRowLwb + this->fNrows - 1) {
1154 Error("SetSub", "row_lwb outof bounds");
1155 return *this;
1156 }
1157 if (col_lwb < this->fColLwb || col_lwb > this->fColLwb + this->fNcols - 1) {
1158 Error("SetSub", "col_lwb outof bounds");
1159 return *this;
1160 }
1161 if (row_lwb + source.GetNrows() > this->fRowLwb + this->fNrows ||
1162 col_lwb + source.GetNcols() > this->fColLwb + this->fNcols) {
1163 Error("SetSub", "source matrix too large");
1164 return *this;
1165 }
1166 }
1167
1168 const Int_t nRows_source = source.GetNrows();
1169 const Int_t nCols_source = source.GetNcols();
1170
1171 if (source.GetRowIndexArray() && source.GetColIndexArray()) {
1172 const Int_t rowlwb_s = source.GetRowLwb();
1173 const Int_t collwb_s = source.GetColLwb();
1174 for (Int_t irow = 0; irow < nRows_source; irow++) {
1175 for (Int_t icol = 0; icol < nCols_source; icol++) {
1176 (*this)(row_lwb + irow, col_lwb + icol) = source(rowlwb_s + irow, collwb_s + icol);
1177 }
1178 }
1179 } else {
1180 const Element *bp = source.GetMatrixArray();
1181 Element *ap = this->GetMatrixArray() + (row_lwb - this->fRowLwb) * this->fNcols + (col_lwb - this->fColLwb);
1182
1183 for (Int_t irow = 0; irow < nRows_source; irow++) {
1184 Element *ap_sub = ap;
1185 for (Int_t icol = 0; icol < nCols_source; icol++) {
1186 *ap_sub++ = *bp++;
1187 }
1188 ap += this->fNcols;
1189 }
1190 }
1191
1192 return *this;
1193}
1194
1195////////////////////////////////////////////////////////////////////////////////
1196/// Set size of the matrix to nrows x ncols
1197/// New dynamic elements are created, the overlapping part of the old ones are
1198/// copied to the new structures, then the old elements are deleted.
1199
1200template <class Element>
1201TMatrixTBase<Element> &TMatrixT<Element>::ResizeTo(Int_t nrows, Int_t ncols, Int_t /*nr_nonzeros*/)
1202{
1203 R__ASSERT(this->IsValid());
1204 if (!this->fIsOwner) {
1205 Error("ResizeTo(Int_t,Int_t)", "Not owner of data array,cannot resize");
1206 return *this;
1207 }
1208
1209 if (this->fNelems > 0) {
1210 if (this->fNrows == nrows && this->fNcols == ncols)
1211 return *this;
1212 else if (nrows == 0 || ncols == 0) {
1213 this->fNrows = nrows;
1214 this->fNcols = ncols;
1215 Clear();
1216 return *this;
1217 }
1218
1219 Element *elements_old = GetMatrixArray();
1220 const Int_t nelems_old = this->fNelems;
1221 const Int_t nrows_old = this->fNrows;
1222 const Int_t ncols_old = this->fNcols;
1223
1224 Allocate(nrows, ncols);
1225 R__ASSERT(this->IsValid());
1226
1227 Element *elements_new = GetMatrixArray();
1228 // new memory should be initialized but be careful not to wipe out the stack
1229 // storage. Initialize all when old or new storage was on the heap
1230 if (this->fNelems > this->kSizeMax || nelems_old > this->kSizeMax)
1231 memset(elements_new, 0, this->fNelems * sizeof(Element));
1232 else if (this->fNelems > nelems_old)
1233 memset(elements_new + nelems_old, 0, (this->fNelems - nelems_old) * sizeof(Element));
1234
1235 // Copy overlap
1236 const Int_t ncols_copy = TMath::Min(this->fNcols, ncols_old);
1237 const Int_t nrows_copy = TMath::Min(this->fNrows, nrows_old);
1238
1239 const Int_t nelems_new = this->fNelems;
1240 if (ncols_old < this->fNcols) {
1241 for (Int_t i = nrows_copy - 1; i >= 0; i--) {
1242 Memcpy_m(elements_new + i * this->fNcols, elements_old + i * ncols_old, ncols_copy, nelems_new, nelems_old);
1243 if (this->fNelems <= this->kSizeMax && nelems_old <= this->kSizeMax)
1244 memset(elements_new + i * this->fNcols + ncols_copy, 0, (this->fNcols - ncols_copy) * sizeof(Element));
1245 }
1246 } else {
1247 for (Int_t i = 0; i < nrows_copy; i++)
1248 Memcpy_m(elements_new + i * this->fNcols, elements_old + i * ncols_old, ncols_copy, nelems_new, nelems_old);
1249 }
1250
1251 Delete_m(nelems_old, elements_old);
1252 } else {
1253 Allocate(nrows, ncols, 0, 0, 1);
1254 }
1255
1256 return *this;
1257}
1258
1259////////////////////////////////////////////////////////////////////////////////
1260/// Set size of the matrix to [row_lwb:row_upb] x [col_lwb:col_upb]
1261/// New dynamic elements are created, the overlapping part of the old ones are
1262/// copied to the new structures, then the old elements are deleted.
1263
1264template <class Element>
1265TMatrixTBase<Element> &
1266TMatrixT<Element>::ResizeTo(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb, Int_t /*nr_nonzeros*/)
1267{
1268 R__ASSERT(this->IsValid());
1269 if (!this->fIsOwner) {
1270 Error("ResizeTo(Int_t,Int_t,Int_t,Int_t)", "Not owner of data array,cannot resize");
1271 return *this;
1272 }
1273
1274 const Int_t new_nrows = row_upb - row_lwb + 1;
1275 const Int_t new_ncols = col_upb - col_lwb + 1;
1276
1277 if (this->fNelems > 0) {
1278
1279 if (this->fNrows == new_nrows && this->fNcols == new_ncols && this->fRowLwb == row_lwb &&
1280 this->fColLwb == col_lwb)
1281 return *this;
1282 else if (new_nrows == 0 || new_ncols == 0) {
1283 this->fNrows = new_nrows;
1284 this->fNcols = new_ncols;
1285 this->fRowLwb = row_lwb;
1286 this->fColLwb = col_lwb;
1287 Clear();
1288 return *this;
1289 }
1290
1291 Element *elements_old = GetMatrixArray();
1292 const Int_t nelems_old = this->fNelems;
1293 const Int_t nrows_old = this->fNrows;
1294 const Int_t ncols_old = this->fNcols;
1295 const Int_t rowLwb_old = this->fRowLwb;
1296 const Int_t colLwb_old = this->fColLwb;
1297
1298 Allocate(new_nrows, new_ncols, row_lwb, col_lwb);
1299 R__ASSERT(this->IsValid());
1300
1301 Element *elements_new = GetMatrixArray();
1302 // new memory should be initialized but be careful not to wipe out the stack
1303 // storage. Initialize all when old or new storage was on the heap
1304 if (this->fNelems > this->kSizeMax || nelems_old > this->kSizeMax)
1305 memset(elements_new, 0, this->fNelems * sizeof(Element));
1306 else if (this->fNelems > nelems_old)
1307 memset(elements_new + nelems_old, 0, (this->fNelems - nelems_old) * sizeof(Element));
1308
1309 // Copy overlap
1310 const Int_t rowLwb_copy = TMath::Max(this->fRowLwb, rowLwb_old);
1311 const Int_t colLwb_copy = TMath::Max(this->fColLwb, colLwb_old);
1312 const Int_t rowUpb_copy = TMath::Min(this->fRowLwb + this->fNrows - 1, rowLwb_old + nrows_old - 1);
1313 const Int_t colUpb_copy = TMath::Min(this->fColLwb + this->fNcols - 1, colLwb_old + ncols_old - 1);
1314
1315 const Int_t nrows_copy = rowUpb_copy - rowLwb_copy + 1;
1316 const Int_t ncols_copy = colUpb_copy - colLwb_copy + 1;
1317
1318 if (nrows_copy > 0 && ncols_copy > 0) {
1319 const Int_t colOldOff = colLwb_copy - colLwb_old;
1320 const Int_t colNewOff = colLwb_copy - this->fColLwb;
1321 if (ncols_old < this->fNcols) {
1322 for (Int_t i = nrows_copy - 1; i >= 0; i--) {
1323 const Int_t iRowOld = rowLwb_copy + i - rowLwb_old;
1324 const Int_t iRowNew = rowLwb_copy + i - this->fRowLwb;
1325 Memcpy_m(elements_new + iRowNew * this->fNcols + colNewOff,
1326 elements_old + iRowOld * ncols_old + colOldOff, ncols_copy, this->fNelems, nelems_old);
1327 if (this->fNelems <= this->kSizeMax && nelems_old <= this->kSizeMax)
1328 memset(elements_new + iRowNew * this->fNcols + colNewOff + ncols_copy, 0,
1329 (this->fNcols - ncols_copy) * sizeof(Element));
1330 }
1331 } else {
1332 for (Int_t i = 0; i < nrows_copy; i++) {
1333 const Int_t iRowOld = rowLwb_copy + i - rowLwb_old;
1334 const Int_t iRowNew = rowLwb_copy + i - this->fRowLwb;
1335 Memcpy_m(elements_new + iRowNew * this->fNcols + colNewOff,
1336 elements_old + iRowOld * ncols_old + colOldOff, ncols_copy, this->fNelems, nelems_old);
1337 }
1338 }
1339 }
1340
1341 Delete_m(nelems_old, elements_old);
1342 } else {
1343 Allocate(new_nrows, new_ncols, row_lwb, col_lwb, 1);
1344 }
1345
1346 return *this;
1347}
1348
1349////////////////////////////////////////////////////////////////////////////////
1350/// Return the matrix determinant
1351
1352template <class Element>
1353Double_t TMatrixT<Element>::Determinant() const
1354{
1355 const TMatrixT<Element> &tmp = *this;
1356 TDecompLU lu(tmp, this->fTol);
1357 Double_t d1, d2;
1358 lu.Det(d1, d2);
1359 return d1 * TMath::Power(2.0, d2);
1360}
1361
1362////////////////////////////////////////////////////////////////////////////////
1363/// Return the matrix determinant as d1,d2 where det = d1*TMath::Power(2.0,d2)
1364
1365template <class Element>
1366void TMatrixT<Element>::Determinant(Double_t &d1, Double_t &d2) const
1367{
1368 const TMatrixT<Element> &tmp = *this;
1369 TDecompLU lu(tmp, Double_t(this->fTol));
1370 lu.Det(d1, d2);
1371}
1372
1373////////////////////////////////////////////////////////////////////////////////
1374/// Invert the matrix and calculate its determinant
1375
1376template <>
1377TMatrixT<Double_t> &TMatrixT<Double_t>::Invert(Double_t *det)
1378{
1379 R__ASSERT(this->IsValid());
1380 TDecompLU::InvertLU(*this, Double_t(fTol), det);
1381 return *this;
1382}
1383
1384////////////////////////////////////////////////////////////////////////////////
1385/// Invert the matrix and calculate its determinant
1386
1387template <class Element>
1388TMatrixT<Element> &TMatrixT<Element>::Invert(Double_t *det)
1389{
1390 TMatrixD tmp(*this);
1391 if (TDecompLU::InvertLU(tmp, Double_t(this->fTol), det))
1392 std::copy(tmp.GetMatrixArray(), tmp.GetMatrixArray() + this->GetNoElements(), this->GetMatrixArray());
1393
1394 return *this;
1395}
1396
1397////////////////////////////////////////////////////////////////////////////////
1398/// Invert the matrix and calculate its determinant, however upto (6x6)
1399/// a fast Cramer inversion is used .
1400
1401template <class Element>
1402TMatrixT<Element> &TMatrixT<Element>::InvertFast(Double_t *det)
1403{
1404 R__ASSERT(this->IsValid());
1405
1406 const Char_t nRows = Char_t(this->GetNrows());
1407 switch (nRows) {
1408 case 1: {
1409 if (this->GetNrows() != this->GetNcols() || this->GetRowLwb() != this->GetColLwb()) {
1410 Error("Invert()", "matrix should be square");
1411 } else {
1412 Element *pM = this->GetMatrixArray();
1413 if (*pM == 0.) {
1414 Error("InvertFast", "matrix is singular");
1415 *det = 0;
1416 } else {
1417 *det = *pM;
1418 *pM = 1.0 / (*pM);
1419 }
1420 }
1421 return *this;
1422 }
1423 case 2: {
1424 TMatrixTCramerInv::Inv2x2<Element>(*this, det);
1425 return *this;
1426 }
1427 case 3: {
1428 TMatrixTCramerInv::Inv3x3<Element>(*this, det);
1429 return *this;
1430 }
1431 case 4: {
1432 TMatrixTCramerInv::Inv4x4<Element>(*this, det);
1433 return *this;
1434 }
1435 case 5: {
1436 TMatrixTCramerInv::Inv5x5<Element>(*this, det);
1437 return *this;
1438 }
1439 case 6: {
1440 TMatrixTCramerInv::Inv6x6<Element>(*this, det);
1441 return *this;
1442 }
1443 default: {
1444 return Invert(det);
1445 }
1446 }
1447}
1448
1449////////////////////////////////////////////////////////////////////////////////
1450/// Transpose matrix source.
1451
1452template <class Element>
1453TMatrixT<Element> &TMatrixT<Element>::Transpose(const TMatrixT<Element> &source)
1454{
1455 R__ASSERT(this->IsValid());
1456 R__ASSERT(source.IsValid());
1457
1458 if (this->GetMatrixArray() == source.GetMatrixArray()) {
1459 Element *ap = this->GetMatrixArray();
1460 if (this->fNrows == this->fNcols && this->fRowLwb == this->fColLwb) {
1461 for (Int_t i = 0; i < this->fNrows; i++) {
1462 const Int_t off_i = i * this->fNrows;
1463 for (Int_t j = i + 1; j < this->fNcols; j++) {
1464 const Int_t off_j = j * this->fNcols;
1465 const Element tmp = ap[off_i + j];
1466 ap[off_i + j] = ap[off_j + i];
1467 ap[off_j + i] = tmp;
1468 }
1469 }
1470 } else {
1471 Element *oldElems = new Element[source.GetNoElements()];
1472 memcpy(oldElems, source.GetMatrixArray(), source.GetNoElements() * sizeof(Element));
1473 const Int_t nrows_old = this->fNrows;
1474 const Int_t ncols_old = this->fNcols;
1475 const Int_t rowlwb_old = this->fRowLwb;
1476 const Int_t collwb_old = this->fColLwb;
1477
1478 this->fNrows = ncols_old;
1479 this->fNcols = nrows_old;
1480 this->fRowLwb = collwb_old;
1481 this->fColLwb = rowlwb_old;
1482 for (Int_t irow = this->fRowLwb; irow < this->fRowLwb + this->fNrows; irow++) {
1483 for (Int_t icol = this->fColLwb; icol < this->fColLwb + this->fNcols; icol++) {
1484 const Int_t off = (icol - collwb_old) * ncols_old;
1485 (*this)(irow, icol) = oldElems[off + irow - rowlwb_old];
1486 }
1487 }
1488 delete[] oldElems;
1489 }
1490 } else {
1491 if (this->fNrows != source.GetNcols() || this->fNcols != source.GetNrows() ||
1492 this->fRowLwb != source.GetColLwb() || this->fColLwb != source.GetRowLwb()) {
1493 Error("Transpose", "matrix has wrong shape");
1494 return *this;
1495 }
1496
1497 const Element *sp1 = source.GetMatrixArray();
1498 const Element *scp = sp1; // Row source pointer
1499 Element *tp = this->GetMatrixArray();
1500 const Element *const tp_last = this->GetMatrixArray() + this->fNelems;
1501
1502 // (This: target) matrix is traversed row-wise way,
1503 // whilst the source matrix is scanned column-wise
1504 while (tp < tp_last) {
1505 const Element *sp2 = scp++;
1506
1507 // Move tp to the next elem in the row and sp to the next elem in the curr col
1508 while (sp2 < sp1 + this->fNelems) {
1509 *tp++ = *sp2;
1510 sp2 += this->fNrows;
1511 }
1512 }
1513 R__ASSERT(tp == tp_last && scp == sp1 + this->fNrows);
1514 }
1515
1516 return *this;
1517}
1518
1519////////////////////////////////////////////////////////////////////////////////
1520/// Perform a rank 1 operation on matrix A:
1521/// A += alpha * v * v^T
1522
1523template <class Element>
1524TMatrixT<Element> &TMatrixT<Element>::Rank1Update(const TVectorT<Element> &v, Element alpha)
1525{
1526 if (gMatrixCheck) {
1527 R__ASSERT(this->IsValid());
1528 R__ASSERT(v.IsValid());
1529 if (v.GetNoElements() < TMath::Max(this->fNrows, this->fNcols)) {
1530 Error("Rank1Update", "vector too short");
1531 return *this;
1532 }
1533 }
1534
1535 const Element *const pv = v.GetMatrixArray();
1536 Element *mp = this->GetMatrixArray();
1537
1538 for (Int_t i = 0; i < this->fNrows; i++) {
1539 const Element tmp = alpha * pv[i];
1540 for (Int_t j = 0; j < this->fNcols; j++)
1541 *mp++ += tmp * pv[j];
1542 }
1543
1544 return *this;
1545}
1546
1547////////////////////////////////////////////////////////////////////////////////
1548/// Perform a rank 1 operation on matrix A:
1549/// A += alpha * v1 * v2^T
1550
1551template <class Element>
1552TMatrixT<Element> &
1553TMatrixT<Element>::Rank1Update(const TVectorT<Element> &v1, const TVectorT<Element> &v2, Element alpha)
1554{
1555 if (gMatrixCheck) {
1556 R__ASSERT(this->IsValid());
1557 R__ASSERT(v1.IsValid());
1558 R__ASSERT(v2.IsValid());
1559 if (v1.GetNoElements() < this->fNrows) {
1560 Error("Rank1Update", "vector v1 too short");
1561 return *this;
1562 }
1563
1564 if (v2.GetNoElements() < this->fNcols) {
1565 Error("Rank1Update", "vector v2 too short");
1566 return *this;
1567 }
1568 }
1569
1570 const Element *const pv1 = v1.GetMatrixArray();
1571 const Element *const pv2 = v2.GetMatrixArray();
1572 Element *mp = this->GetMatrixArray();
1573
1574 for (Int_t i = 0; i < this->fNrows; i++) {
1575 const Element tmp = alpha * pv1[i];
1576 for (Int_t j = 0; j < this->fNcols; j++)
1577 *mp++ += tmp * pv2[j];
1578 }
1579
1580 return *this;
1581}
1582
1583////////////////////////////////////////////////////////////////////////////////
1584/// Calculate scalar v * (*this) * v^T
1585
1586template <class Element>
1587Element TMatrixT<Element>::Similarity(const TVectorT<Element> &v) const
1588{
1589 if (gMatrixCheck) {
1590 R__ASSERT(this->IsValid());
1591 R__ASSERT(v.IsValid());
1592 if (this->fNcols != this->fNrows || this->fColLwb != this->fRowLwb) {
1593 Error("Similarity(const TVectorT &)", "matrix is not square");
1594 return -1.;
1595 }
1596
1597 if (this->fNcols != v.GetNrows() || this->fColLwb != v.GetLwb()) {
1598 Error("Similarity(const TVectorT &)", "vector and matrix incompatible");
1599 return -1.;
1600 }
1601 }
1602
1603 const Element *mp = this->GetMatrixArray(); // Matrix row ptr
1604 const Element *vp = v.GetMatrixArray(); // vector ptr
1605
1606 Element sum1 = 0;
1607 const Element *const vp_first = vp;
1608 const Element *const vp_last = vp + v.GetNrows();
1609 while (vp < vp_last) {
1610 Element sum2 = 0;
1611 for (const Element *sp = vp_first; sp < vp_last;)
1612 sum2 += *mp++ * *sp++;
1613 sum1 += sum2 * *vp++;
1614 }
1615
1616 R__ASSERT(mp == this->GetMatrixArray() + this->GetNoElements());
1617
1618 return sum1;
1619}
1620
1621////////////////////////////////////////////////////////////////////////////////
1622/// Multiply/divide matrix columns by a vector:
1623/// option:
1624/// "D" : b(i,j) = a(i,j)/v(i) i = 0,fNrows-1 (default)
1625/// else : b(i,j) = a(i,j)*v(i)
1626
1627template <class Element>
1628TMatrixT<Element> &TMatrixT<Element>::NormByColumn(const TVectorT<Element> &v, Option_t *option)
1629{
1630 if (gMatrixCheck) {
1631 R__ASSERT(this->IsValid());
1632 R__ASSERT(v.IsValid());
1633 if (v.GetNoElements() < this->fNrows) {
1634 Error("NormByColumn", "vector shorter than matrix column");
1635 return *this;
1636 }
1637 }
1638
1639 TString opt(option);
1640 opt.ToUpper();
1641 const Int_t divide = (opt.Contains("D")) ? 1 : 0;
1642
1643 const Element *pv = v.GetMatrixArray();
1644 Element *mp = this->GetMatrixArray();
1645 const Element *const mp_last = mp + this->fNelems;
1646
1647 if (divide) {
1648 for (; mp < mp_last; pv++) {
1649 for (Int_t j = 0; j < this->fNcols; j++) {
1650 if (*pv != 0.0)
1651 *mp++ /= *pv;
1652 else {
1653 Error("NormbyColumn", "vector element %ld is zero", Long_t(pv - v.GetMatrixArray()));
1654 mp++;
1655 }
1656 }
1657 }
1658 } else {
1659 for (; mp < mp_last; pv++)
1660 for (Int_t j = 0; j < this->fNcols; j++)
1661 *mp++ *= *pv;
1662 }
1663
1664 return *this;
1665}
1666
1667////////////////////////////////////////////////////////////////////////////////
1668/// Multiply/divide matrix rows with a vector:
1669/// option:
1670/// "D" : b(i,j) = a(i,j)/v(j) i = 0,fNcols-1 (default)
1671/// else : b(i,j) = a(i,j)*v(j)
1672
1673template <class Element>
1674TMatrixT<Element> &TMatrixT<Element>::NormByRow(const TVectorT<Element> &v, Option_t *option)
1675{
1676 if (gMatrixCheck) {
1677 R__ASSERT(this->IsValid());
1678 R__ASSERT(v.IsValid());
1679 if (v.GetNoElements() < this->fNcols) {
1680 Error("NormByRow", "vector shorter than matrix column");
1681 return *this;
1682 }
1683 }
1684
1685 TString opt(option);
1686 opt.ToUpper();
1687 const Int_t divide = (opt.Contains("D")) ? 1 : 0;
1688
1689 const Element *pv0 = v.GetMatrixArray();
1690 const Element *pv = pv0;
1691 Element *mp = this->GetMatrixArray();
1692 const Element *const mp_last = mp + this->fNelems;
1693
1694 if (divide) {
1695 for (; mp < mp_last; pv = pv0)
1696 for (Int_t j = 0; j < this->fNcols; j++) {
1697 if (*pv != 0.0)
1698 *mp++ /= *pv++;
1699 else {
1700 Error("NormbyRow", "vector element %ld is zero", Long_t(pv - pv0));
1701 mp++;
1702 }
1703 }
1704 } else {
1705 for (; mp < mp_last; pv = pv0)
1706 for (Int_t j = 0; j < this->fNcols; j++)
1707 *mp++ *= *pv++;
1708 }
1709
1710 return *this;
1711}
1712
1713////////////////////////////////////////////////////////////////////////////////
1714/// Assignment operator
1715
1716template <class Element>
1717TMatrixT<Element> &TMatrixT<Element>::operator=(const TMatrixT<Element> &source)
1718{
1719 if (gMatrixCheck && !AreCompatible(*this, source)) {
1720 Error("operator=(const TMatrixT &)", "matrices not compatible");
1721 return *this;
1722 }
1723
1724 if (this->GetMatrixArray() != source.GetMatrixArray()) {
1725 TObject::operator=(source);
1726 memcpy(fElements, source.GetMatrixArray(), this->fNelems * sizeof(Element));
1727 this->fTol = source.GetTol();
1728 }
1729 return *this;
1730}
1731
1732////////////////////////////////////////////////////////////////////////////////
1733/// Assignment operator
1734
1735template <class Element>
1736TMatrixT<Element> &TMatrixT<Element>::operator=(const TMatrixTSym<Element> &source)
1737{
1738 if (gMatrixCheck && !AreCompatible(*this, source)) {
1739 Error("operator=(const TMatrixTSym &)", "matrices not compatible");
1740 return *this;
1741 }
1742
1743 if (this->GetMatrixArray() != source.GetMatrixArray()) {
1744 TObject::operator=(source);
1745 memcpy(fElements, source.GetMatrixArray(), this->fNelems * sizeof(Element));
1746 this->fTol = source.GetTol();
1747 }
1748 return *this;
1749}
1750
1751////////////////////////////////////////////////////////////////////////////////
1752/// Assignment operator
1753
1754template <class Element>
1755TMatrixT<Element> &TMatrixT<Element>::operator=(const TMatrixTSparse<Element> &source)
1756{
1757 if ((gMatrixCheck && this->GetNrows() != source.GetNrows()) || this->GetNcols() != source.GetNcols() ||
1758 this->GetRowLwb() != source.GetRowLwb() || this->GetColLwb() != source.GetColLwb()) {
1759 Error("operator=(const TMatrixTSparse &", "matrices not compatible");
1760 return *this;
1761 }
1762
1763 if (this->GetMatrixArray() != source.GetMatrixArray()) {
1764 TObject::operator=(source);
1765 memset(fElements, 0, this->fNelems * sizeof(Element));
1766
1767 const Element *const sp = source.GetMatrixArray();
1768 Element *tp = this->GetMatrixArray();
1769
1770 const Int_t *const pRowIndex = source.GetRowIndexArray();
1771 const Int_t *const pColIndex = source.GetColIndexArray();
1772
1773 for (Int_t irow = 0; irow < this->fNrows; irow++) {
1774 const Int_t off = irow * this->fNcols;
1775 const Int_t sIndex = pRowIndex[irow];
1776 const Int_t eIndex = pRowIndex[irow + 1];
1777 for (Int_t index = sIndex; index < eIndex; index++)
1778 tp[off + pColIndex[index]] = sp[index];
1779 }
1780 this->fTol = source.GetTol();
1781 }
1782 return *this;
1783}
1784
1785////////////////////////////////////////////////////////////////////////////////
1786/// Assignment operator
1787
1788template <class Element>
1789TMatrixT<Element> &TMatrixT<Element>::operator=(const TMatrixTLazy<Element> &lazy_constructor)
1790{
1791 R__ASSERT(this->IsValid());
1792
1793 if (lazy_constructor.GetRowUpb() != this->GetRowUpb() || lazy_constructor.GetColUpb() != this->GetColUpb() ||
1794 lazy_constructor.GetRowLwb() != this->GetRowLwb() || lazy_constructor.GetColLwb() != this->GetColLwb()) {
1795 Error("operator=(const TMatrixTLazy&)", "matrix is incompatible with "
1796 "the assigned Lazy matrix");
1797 return *this;
1798 }
1799
1800 lazy_constructor.FillIn(*this);
1801 return *this;
1802}
1803
1804////////////////////////////////////////////////////////////////////////////////
1805/// Assign val to every element of the matrix.
1806
1807template <class Element>
1808TMatrixT<Element> &TMatrixT<Element>::operator=(Element val)
1809{
1810 R__ASSERT(this->IsValid());
1811
1812 Element *ep = this->GetMatrixArray();
1813 const Element *const ep_last = ep + this->fNelems;
1814 while (ep < ep_last)
1815 *ep++ = val;
1816
1817 return *this;
1818}
1819
1820////////////////////////////////////////////////////////////////////////////////
1821/// Add val to every element of the matrix.
1822
1823template <class Element>
1824TMatrixT<Element> &TMatrixT<Element>::operator+=(Element val)
1825{
1826 R__ASSERT(this->IsValid());
1827
1828 Element *ep = this->GetMatrixArray();
1829 const Element *const ep_last = ep + this->fNelems;
1830 while (ep < ep_last)
1831 *ep++ += val;
1832
1833 return *this;
1834}
1835
1836////////////////////////////////////////////////////////////////////////////////
1837/// Subtract val from every element of the matrix.
1838
1839template <class Element>
1840TMatrixT<Element> &TMatrixT<Element>::operator-=(Element val)
1841{
1842 R__ASSERT(this->IsValid());
1843
1844 Element *ep = this->GetMatrixArray();
1845 const Element *const ep_last = ep + this->fNelems;
1846 while (ep < ep_last)
1847 *ep++ -= val;
1848
1849 return *this;
1850}
1851
1852////////////////////////////////////////////////////////////////////////////////
1853/// Multiply every element of the matrix with val.
1854
1855template <class Element>
1856TMatrixT<Element> &TMatrixT<Element>::operator*=(Element val)
1857{
1858 R__ASSERT(this->IsValid());
1859
1860 Element *ep = this->GetMatrixArray();
1861 const Element *const ep_last = ep + this->fNelems;
1862 while (ep < ep_last)
1863 *ep++ *= val;
1864
1865 return *this;
1866}
1867
1868////////////////////////////////////////////////////////////////////////////////
1869/// Add the source matrix.
1870
1871template <class Element>
1872TMatrixT<Element> &TMatrixT<Element>::operator+=(const TMatrixT<Element> &source)
1873{
1874 if (gMatrixCheck && !AreCompatible(*this, source)) {
1875 Error("operator+=(const TMatrixT &)", "matrices not compatible");
1876 return *this;
1877 }
1878
1879 const Element *sp = source.GetMatrixArray();
1880 Element *tp = this->GetMatrixArray();
1881 const Element *const tp_last = tp + this->fNelems;
1882 while (tp < tp_last)
1883 *tp++ += *sp++;
1884
1885 return *this;
1886}
1887
1888////////////////////////////////////////////////////////////////////////////////
1889/// Add the source matrix.
1890
1891template <class Element>
1892TMatrixT<Element> &TMatrixT<Element>::operator+=(const TMatrixTSym<Element> &source)
1893{
1894 if (gMatrixCheck && !AreCompatible(*this, source)) {
1895 Error("operator+=(const TMatrixTSym &)", "matrices not compatible");
1896 return *this;
1897 }
1898
1899 const Element *sp = source.GetMatrixArray();
1900 Element *tp = this->GetMatrixArray();
1901 const Element *const tp_last = tp + this->fNelems;
1902 while (tp < tp_last)
1903 *tp++ += *sp++;
1904
1905 return *this;
1906}
1907
1908////////////////////////////////////////////////////////////////////////////////
1909/// Subtract the source matrix.
1910
1911template <class Element>
1912TMatrixT<Element> &TMatrixT<Element>::operator-=(const TMatrixT<Element> &source)
1913{
1914 if (gMatrixCheck && !AreCompatible(*this, source)) {
1915 Error("operator=-(const TMatrixT &)", "matrices not compatible");
1916 return *this;
1917 }
1918
1919 const Element *sp = source.GetMatrixArray();
1920 Element *tp = this->GetMatrixArray();
1921 const Element *const tp_last = tp + this->fNelems;
1922 while (tp < tp_last)
1923 *tp++ -= *sp++;
1924
1925 return *this;
1926}
1927
1928////////////////////////////////////////////////////////////////////////////////
1929/// Subtract the source matrix.
1930
1931template <class Element>
1932TMatrixT<Element> &TMatrixT<Element>::operator-=(const TMatrixTSym<Element> &source)
1933{
1934 if (gMatrixCheck && !AreCompatible(*this, source)) {
1935 Error("operator=-(const TMatrixTSym &)", "matrices not compatible");
1936 return *this;
1937 }
1938
1939 const Element *sp = source.GetMatrixArray();
1940 Element *tp = this->GetMatrixArray();
1941 const Element *const tp_last = tp + this->fNelems;
1942 while (tp < tp_last)
1943 *tp++ -= *sp++;
1944
1945 return *this;
1946}
1947
1948////////////////////////////////////////////////////////////////////////////////
1949/// Compute target = target * source inplace. Strictly speaking, it can't be
1950/// done inplace, though only the row of the target matrix needs to be saved.
1951/// "Inplace" multiplication is only allowed when the 'source' matrix is square.
1952
1953template <class Element>
1954TMatrixT<Element> &TMatrixT<Element>::operator*=(const TMatrixT<Element> &source)
1955{
1956 if (gMatrixCheck) {
1957 R__ASSERT(this->IsValid());
1958 R__ASSERT(source.IsValid());
1959 if (this->fNcols != source.GetNrows() || this->fColLwb != source.GetRowLwb() ||
1960 this->fNcols != source.GetNcols() || this->fColLwb != source.GetColLwb()) {
1961 Error("operator*=(const TMatrixT &)", "source matrix has wrong shape");
1962 return *this;
1963 }
1964 }
1965
1966 // Check for A *= A;
1967 const Element *sp;
1968 TMatrixT<Element> tmp;
1969 if (this->GetMatrixArray() == source.GetMatrixArray()) {
1970 tmp.ResizeTo(source);
1971 tmp = source;
1972 sp = tmp.GetMatrixArray();
1973 } else
1974 sp = source.GetMatrixArray();
1975
1976 // One row of the old_target matrix
1977 Element work[kWorkMax];
1978 Bool_t isAllocated = kFALSE;
1979 Element *trp = work;
1980 if (this->fNcols > kWorkMax) {
1981 isAllocated = kTRUE;
1982 trp = new Element[this->fNcols];
1983 }
1984
1985 Element *cp = this->GetMatrixArray();
1986 const Element *trp0 = cp; // Pointer to target[i,0];
1987 const Element *const trp0_last = trp0 + this->fNelems;
1988 while (trp0 < trp0_last) {
1989 memcpy(trp, trp0, this->fNcols * sizeof(Element)); // copy the i-th row of target, Start at target[i,0]
1990 for (const Element *scp = sp; scp < sp + this->fNcols;) { // Pointer to the j-th column of source,
1991 // Start scp = source[0,0]
1992 Element cij = 0;
1993 for (Int_t j = 0; j < this->fNcols; j++) {
1994 cij += trp[j] * *scp; // the j-th col of source
1995 scp += this->fNcols;
1996 }
1997 *cp++ = cij;
1998 scp -= source.GetNoElements() - 1; // Set bcp to the (j+1)-th col
1999 }
2000 trp0 += this->fNcols; // Set trp0 to the (i+1)-th row
2001 R__ASSERT(trp0 == cp);
2002 }
2003
2004 R__ASSERT(cp == trp0_last && trp0 == trp0_last);
2005 if (isAllocated)
2006 delete[] trp;
2007
2008 return *this;
2009}
2010
2011////////////////////////////////////////////////////////////////////////////////
2012/// Compute target = target * source inplace. Strictly speaking, it can't be
2013/// done inplace, though only the row of the target matrix needs to be saved.
2014
2015template <class Element>
2016TMatrixT<Element> &TMatrixT<Element>::operator*=(const TMatrixTSym<Element> &source)
2017{
2018 if (gMatrixCheck) {
2019 R__ASSERT(this->IsValid());
2020 R__ASSERT(source.IsValid());
2021 if (this->fNcols != source.GetNrows() || this->fColLwb != source.GetRowLwb()) {
2022 Error("operator*=(const TMatrixTSym &)", "source matrix has wrong shape");
2023 return *this;
2024 }
2025 }
2026
2027 // Check for A *= A;
2028 const Element *sp;
2029 TMatrixT<Element> tmp;
2030 if (this->GetMatrixArray() == source.GetMatrixArray()) {
2031 tmp.ResizeTo(source);
2032 tmp = source;
2033 sp = tmp.GetMatrixArray();
2034 } else
2035 sp = source.GetMatrixArray();
2036
2037 // One row of the old_target matrix
2038 Element work[kWorkMax];
2039 Bool_t isAllocated = kFALSE;
2040 Element *trp = work;
2041 if (this->fNcols > kWorkMax) {
2042 isAllocated = kTRUE;
2043 trp = new Element[this->fNcols];
2044 }
2045
2046 Element *cp = this->GetMatrixArray();
2047 const Element *trp0 = cp; // Pointer to target[i,0];
2048 const Element *const trp0_last = trp0 + this->fNelems;
2049 while (trp0 < trp0_last) {
2050 memcpy(trp, trp0, this->fNcols * sizeof(Element)); // copy the i-th row of target, Start at target[i,0]
2051 for (const Element *scp = sp; scp < sp + this->fNcols;) { // Pointer to the j-th column of source,
2052 // Start scp = source[0,0]
2053 Element cij = 0;
2054 for (Int_t j = 0; j < this->fNcols; j++) {
2055 cij += trp[j] * *scp; // the j-th col of source
2056 scp += this->fNcols;
2057 }
2058 *cp++ = cij;
2059 scp -= source.GetNoElements() - 1; // Set bcp to the (j+1)-th col
2060 }
2061 trp0 += this->fNcols; // Set trp0 to the (i+1)-th row
2062 R__ASSERT(trp0 == cp);
2063 }
2064
2065 R__ASSERT(cp == trp0_last && trp0 == trp0_last);
2066 if (isAllocated)
2067 delete[] trp;
2068
2069 return *this;
2070}
2071
2072////////////////////////////////////////////////////////////////////////////////
2073/// Multiply a matrix row by the diagonal of another matrix
2074/// matrix(i,j) *= diag(j), j=0,fNcols-1
2075
2076template <class Element>
2077TMatrixT<Element> &TMatrixT<Element>::operator*=(const TMatrixTDiag_const<Element> &diag)
2078{
2079 if (gMatrixCheck) {
2080 R__ASSERT(this->IsValid());
2081 R__ASSERT(diag.GetMatrix()->IsValid());
2082 if (this->fNcols != diag.GetNdiags()) {
2083 Error("operator*=(const TMatrixTDiag_const &)", "wrong diagonal length");
2084 return *this;
2085 }
2086 }
2087
2088 Element *mp = this->GetMatrixArray(); // Matrix ptr
2089 const Element *const mp_last = mp + this->fNelems;
2090 const Int_t inc = diag.GetInc();
2091 while (mp < mp_last) {
2092 const Element *dp = diag.GetPtr();
2093 for (Int_t j = 0; j < this->fNcols; j++) {
2094 *mp++ *= *dp;
2095 dp += inc;
2096 }
2097 }
2098
2099 return *this;
2100}
2101
2102////////////////////////////////////////////////////////////////////////////////
2103/// Divide a matrix row by the diagonal of another matrix
2104/// matrix(i,j) /= diag(j)
2105
2106template <class Element>
2107TMatrixT<Element> &TMatrixT<Element>::operator/=(const TMatrixTDiag_const<Element> &diag)
2108{
2109 if (gMatrixCheck) {
2110 R__ASSERT(this->IsValid());
2111 R__ASSERT(diag.GetMatrix()->IsValid());
2112 if (this->fNcols != diag.GetNdiags()) {
2113 Error("operator/=(const TMatrixTDiag_const &)", "wrong diagonal length");
2114 return *this;
2115 }
2116 }
2117
2118 Element *mp = this->GetMatrixArray(); // Matrix ptr
2119 const Element *const mp_last = mp + this->fNelems;
2120 const Int_t inc = diag.GetInc();
2121 while (mp < mp_last) {
2122 const Element *dp = diag.GetPtr();
2123 for (Int_t j = 0; j < this->fNcols; j++) {
2124 if (*dp != 0.0)
2125 *mp++ /= *dp;
2126 else {
2127 Error("operator/=", "%d-diagonal element is zero", j);
2128 mp++;
2129 }
2130 dp += inc;
2131 }
2132 }
2133
2134 return *this;
2135}
2136
2137////////////////////////////////////////////////////////////////////////////////
2138/// Multiply a matrix by the column of another matrix
2139/// matrix(i,j) *= another(i,k) for fixed k
2140
2141template <class Element>
2142TMatrixT<Element> &TMatrixT<Element>::operator*=(const TMatrixTColumn_const<Element> &col)
2143{
2144 const TMatrixTBase<Element> *mt = col.GetMatrix();
2145
2146 if (gMatrixCheck) {
2147 R__ASSERT(this->IsValid());
2148 R__ASSERT(mt->IsValid());
2149 if (this->fNrows != mt->GetNrows()) {
2150 Error("operator*=(const TMatrixTColumn_const &)", "wrong column length");
2151 return *this;
2152 }
2153 }
2154
2155 const Element *const endp = col.GetPtr() + mt->GetNoElements();
2156 Element *mp = this->GetMatrixArray(); // Matrix ptr
2157 const Element *const mp_last = mp + this->fNelems;
2158 const Element *cp = col.GetPtr(); // ptr
2159 const Int_t inc = col.GetInc();
2160 while (mp < mp_last) {
2161 R__ASSERT(cp < endp);
2162 for (Int_t j = 0; j < this->fNcols; j++)
2163 *mp++ *= *cp;
2164 cp += inc;
2165 }
2166
2167 return *this;
2168}
2169
2170////////////////////////////////////////////////////////////////////////////////
2171/// Divide a matrix by the column of another matrix
2172/// matrix(i,j) /= another(i,k) for fixed k
2173
2174template <class Element>
2175TMatrixT<Element> &TMatrixT<Element>::operator/=(const TMatrixTColumn_const<Element> &col)
2176{
2177 const TMatrixTBase<Element> *mt = col.GetMatrix();
2178
2179 if (gMatrixCheck) {
2180 R__ASSERT(this->IsValid());
2181 R__ASSERT(mt->IsValid());
2182 if (this->fNrows != mt->GetNrows()) {
2183 Error("operator/=(const TMatrixTColumn_const &)", "wrong column matrix");
2184 return *this;
2185 }
2186 }
2187
2188 const Element *const endp = col.GetPtr() + mt->GetNoElements();
2189 Element *mp = this->GetMatrixArray(); // Matrix ptr
2190 const Element *const mp_last = mp + this->fNelems;
2191 const Element *cp = col.GetPtr(); // ptr
2192 const Int_t inc = col.GetInc();
2193 while (mp < mp_last) {
2194 R__ASSERT(cp < endp);
2195 if (*cp != 0.0) {
2196 for (Int_t j = 0; j < this->fNcols; j++)
2197 *mp++ /= *cp;
2198 } else {
2199 const Int_t icol = (cp - mt->GetMatrixArray()) / inc;
2200 Error("operator/=", "%d-row of matrix column is zero", icol);
2201 mp += this->fNcols;
2202 }
2203 cp += inc;
2204 }
2205
2206 return *this;
2207}
2208
2209////////////////////////////////////////////////////////////////////////////////
2210/// Multiply a matrix by the row of another matrix
2211/// matrix(i,j) *= another(k,j) for fixed k
2212
2213template <class Element>
2214TMatrixT<Element> &TMatrixT<Element>::operator*=(const TMatrixTRow_const<Element> &row)
2215{
2216 const TMatrixTBase<Element> *mt = row.GetMatrix();
2217
2218 if (gMatrixCheck) {
2219 R__ASSERT(this->IsValid());
2220 R__ASSERT(mt->IsValid());
2221 if (this->fNcols != mt->GetNcols()) {
2222 Error("operator*=(const TMatrixTRow_const &)", "wrong row length");
2223 return *this;
2224 }
2225 }
2226
2227 const Element *const endp = row.GetPtr() + mt->GetNoElements();
2228 Element *mp = this->GetMatrixArray(); // Matrix ptr
2229 const Element *const mp_last = mp + this->fNelems;
2230 const Int_t inc = row.GetInc();
2231 while (mp < mp_last) {
2232 const Element *rp = row.GetPtr(); // Row ptr
2233 for (Int_t j = 0; j < this->fNcols; j++) {
2234 R__ASSERT(rp < endp);
2235 *mp++ *= *rp;
2236 rp += inc;
2237 }
2238 }
2239
2240 return *this;
2241}
2242
2243////////////////////////////////////////////////////////////////////////////////
2244/// Divide a matrix by the row of another matrix
2245/// matrix(i,j) /= another(k,j) for fixed k
2246
2247template <class Element>
2248TMatrixT<Element> &TMatrixT<Element>::operator/=(const TMatrixTRow_const<Element> &row)
2249{
2250 const TMatrixTBase<Element> *mt = row.GetMatrix();
2251 R__ASSERT(this->IsValid());
2252 R__ASSERT(mt->IsValid());
2253
2254 if (this->fNcols != mt->GetNcols()) {
2255 Error("operator/=(const TMatrixTRow_const &)", "wrong row length");
2256 return *this;
2257 }
2258
2259 const Element *const endp = row.GetPtr() + mt->GetNoElements();
2260 Element *mp = this->GetMatrixArray(); // Matrix ptr
2261 const Element *const mp_last = mp + this->fNelems;
2262 const Int_t inc = row.GetInc();
2263 while (mp < mp_last) {
2264 const Element *rp = row.GetPtr(); // Row ptr
2265 for (Int_t j = 0; j < this->fNcols; j++) {
2266 R__ASSERT(rp < endp);
2267 if (*rp != 0.0) {
2268 *mp++ /= *rp;
2269 } else {
2270 Error("operator/=", "%d-col of matrix row is zero", j);
2271 mp++;
2272 }
2273 rp += inc;
2274 }
2275 }
2276
2277 return *this;
2278}
2279
2280////////////////////////////////////////////////////////////////////////////////
2281/// Return a matrix containing the eigen-vectors ordered by descending values
2282/// of Re^2+Im^2 of the complex eigen-values .
2283/// If the matrix is asymmetric, only the real part of the eigen-values is
2284/// returned . For full functionality use TMatrixDEigen .
2285
2286template <class Element>
2287const TMatrixT<Element> TMatrixT<Element>::EigenVectors(TVectorT<Element> &eigenValues) const
2288{
2289 if (!this->IsSymmetric())
2290 Warning("EigenVectors(TVectorT &)", "Only real part of eigen-values will be returned");
2291 TMatrixDEigen eigen(*this);
2292 eigenValues.ResizeTo(this->fNrows);
2293 eigenValues = eigen.GetEigenValuesRe();
2294 return eigen.GetEigenVectors();
2295}
2296
2297////////////////////////////////////////////////////////////////////////////////
2298/// operation this = source1+source2
2299
2300template <class Element>
2307
2308////////////////////////////////////////////////////////////////////////////////
2309/// operation this = source1+source2
2310
2311template <class Element>
2318
2319////////////////////////////////////////////////////////////////////////////////
2320/// operation this = source1+source2
2321
2322template <class Element>
2327
2328////////////////////////////////////////////////////////////////////////////////
2329/// operation this = source+val
2330
2331template <class Element>
2338
2339////////////////////////////////////////////////////////////////////////////////
2340/// operation this = val+source
2341
2342template <class Element>
2347
2348////////////////////////////////////////////////////////////////////////////////
2349/// operation this = source1-source2
2350
2351template <class Element>
2358
2359////////////////////////////////////////////////////////////////////////////////
2360/// operation this = source1-source2
2361
2362template <class Element>
2369
2370////////////////////////////////////////////////////////////////////////////////
2371/// operation this = source1-source2
2372
2373template <class Element>
2378
2379////////////////////////////////////////////////////////////////////////////////
2380/// operation this = source-val
2381
2382template <class Element>
2389
2390////////////////////////////////////////////////////////////////////////////////
2391/// operation this = val-source
2392
2393template <class Element>
2394TMatrixT<Element> TMatrixTAutoloadOps::operator-(Element val, const TMatrixT<Element> &source)
2395{
2396 return Element(-1.0) * operator-(source, val);
2397}
2398
2399////////////////////////////////////////////////////////////////////////////////
2400/// operation this = val*source
2401
2402template <class Element>
2409
2410////////////////////////////////////////////////////////////////////////////////
2411/// operation this = val*source
2412
2413template <class Element>
2418
2419////////////////////////////////////////////////////////////////////////////////
2420/// operation this = source1*source2
2421
2422template <class Element>
2428
2429////////////////////////////////////////////////////////////////////////////////
2430/// operation this = source1*source2
2431
2432template <class Element>
2438
2439////////////////////////////////////////////////////////////////////////////////
2440/// operation this = source1*source2
2441
2442template <class Element>
2448
2449////////////////////////////////////////////////////////////////////////////////
2450/// operation this = source1*source2
2451
2452template <class Element>
2453TMatrixT<Element>
2459
2460////////////////////////////////////////////////////////////////////////////////
2461/// Logical AND
2462
2463template <class Element>
2464TMatrixT<Element> TMatrixTAutoloadOps::operator&&(const TMatrixT<Element> &source1, const TMatrixT<Element> &source2)
2465{
2466 TMatrixT<Element> target;
2467
2468 if (gMatrixCheck && !AreCompatible(source1, source2)) {
2469 Error("operator&&(const TMatrixT&,const TMatrixT&)", "matrices not compatible");
2470 return target;
2471 }
2472
2473 target.ResizeTo(source1);
2474
2475 const Element *sp1 = source1.GetMatrixArray();
2476 const Element *sp2 = source2.GetMatrixArray();
2477 Element *tp = target.GetMatrixArray();
2478 const Element *const tp_last = tp + target.GetNoElements();
2479 while (tp < tp_last)
2480 *tp++ = (*sp1++ != 0.0 && *sp2++ != 0.0);
2481
2482 return target;
2483}
2484
2485////////////////////////////////////////////////////////////////////////////////
2486/// Logical AND
2487
2488template <class Element>
2489TMatrixT<Element> TMatrixTAutoloadOps::operator&&(const TMatrixT<Element> &source1, const TMatrixTSym<Element> &source2)
2490{
2491 TMatrixT<Element> target;
2492
2493 if (gMatrixCheck && !AreCompatible(source1, source2)) {
2494 Error("operator&&(const TMatrixT&,const TMatrixTSym&)", "matrices not compatible");
2495 return target;
2496 }
2497
2498 target.ResizeTo(source1);
2499
2500 const Element *sp1 = source1.GetMatrixArray();
2501 const Element *sp2 = source2.GetMatrixArray();
2502 Element *tp = target.GetMatrixArray();
2503 const Element *const tp_last = tp + target.GetNoElements();
2504 while (tp < tp_last)
2505 *tp++ = (*sp1++ != 0.0 && *sp2++ != 0.0);
2506
2507 return target;
2508}
2509
2510////////////////////////////////////////////////////////////////////////////////
2511/// Logical AND
2512
2513template <class Element>
2518
2519////////////////////////////////////////////////////////////////////////////////
2520/// Logical OR
2521
2522template <class Element>
2523TMatrixT<Element> TMatrixTAutoloadOps::operator||(const TMatrixT<Element> &source1, const TMatrixT<Element> &source2)
2524{
2525 TMatrixT<Element> target;
2526
2527 if (gMatrixCheck && !AreCompatible(source1, source2)) {
2528 Error("operator||(const TMatrixT&,const TMatrixT&)", "matrices not compatible");
2529 return target;
2530 }
2531
2532 target.ResizeTo(source1);
2533
2534 const Element *sp1 = source1.GetMatrixArray();
2535 const Element *sp2 = source2.GetMatrixArray();
2536 Element *tp = target.GetMatrixArray();
2537 const Element *const tp_last = tp + target.GetNoElements();
2538 while (tp < tp_last)
2539 *tp++ = (*sp1++ != 0.0 || *sp2++ != 0.0);
2540
2541 return target;
2542}
2543
2544////////////////////////////////////////////////////////////////////////////////
2545/// Logical OR
2546
2547template <class Element>
2548TMatrixT<Element> TMatrixTAutoloadOps::operator||(const TMatrixT<Element> &source1, const TMatrixTSym<Element> &source2)
2549{
2550 TMatrixT<Element> target;
2551
2552 if (gMatrixCheck && !AreCompatible(source1, source2)) {
2553 Error("operator||(const TMatrixT&,const TMatrixTSym&)", "matrices not compatible");
2554 return target;
2555 }
2556
2557 target.ResizeTo(source1);
2558
2559 const Element *sp1 = source1.GetMatrixArray();
2560 const Element *sp2 = source2.GetMatrixArray();
2561 Element *tp = target.GetMatrixArray();
2562 const Element *const tp_last = tp + target.GetNoElements();
2563 while (tp < tp_last)
2564 *tp++ = (*sp1++ != 0.0 || *sp2++ != 0.0);
2565
2566 return target;
2567}
2568
2569////////////////////////////////////////////////////////////////////////////////
2570/// Logical OR
2571
2572template <class Element>
2577
2578////////////////////////////////////////////////////////////////////////////////
2579/// logical operation source1 > source2
2580
2581template <class Element>
2582TMatrixT<Element> TMatrixTAutoloadOps::operator>(const TMatrixT<Element> &source1, const TMatrixT<Element> &source2)
2583{
2584 TMatrixT<Element> target;
2585
2586 if (gMatrixCheck && !AreCompatible(source1, source2)) {
2587 Error("operator|(const TMatrixT&,const TMatrixT&)", "matrices not compatible");
2588 return target;
2589 }
2590
2591 target.ResizeTo(source1);
2592
2593 const Element *sp1 = source1.GetMatrixArray();
2594 const Element *sp2 = source2.GetMatrixArray();
2595 Element *tp = target.GetMatrixArray();
2596 const Element *const tp_last = tp + target.GetNoElements();
2597 while (tp < tp_last) {
2598 *tp++ = (*sp1) > (*sp2);
2599 sp1++;
2600 sp2++;
2601 }
2602
2603 return target;
2604}
2605
2606////////////////////////////////////////////////////////////////////////////////
2607/// logical operation source1 > source2
2608
2609template <class Element>
2610TMatrixT<Element> TMatrixTAutoloadOps::operator>(const TMatrixT<Element> &source1, const TMatrixTSym<Element> &source2)
2611{
2612 TMatrixT<Element> target;
2613
2614 if (gMatrixCheck && !AreCompatible(source1, source2)) {
2615 Error("operator>(const TMatrixT&,const TMatrixTSym&)", "matrices not compatible");
2616 return target;
2617 }
2618
2619 target.ResizeTo(source1);
2620
2621 const Element *sp1 = source1.GetMatrixArray();
2622 const Element *sp2 = source2.GetMatrixArray();
2623 Element *tp = target.GetMatrixArray();
2624 const Element *const tp_last = tp + target.GetNoElements();
2625 while (tp < tp_last) {
2626 *tp++ = (*sp1) > (*sp2);
2627 sp1++;
2628 sp2++;
2629 }
2630
2631 return target;
2632}
2633
2634////////////////////////////////////////////////////////////////////////////////
2635/// logical operation source1 > source2
2636
2637template <class Element>
2642
2643////////////////////////////////////////////////////////////////////////////////
2644/// logical operation source1 >= source2
2645
2646template <class Element>
2647TMatrixT<Element> TMatrixTAutoloadOps::operator>=(const TMatrixT<Element> &source1, const TMatrixT<Element> &source2)
2648{
2649 TMatrixT<Element> target;
2650
2651 if (gMatrixCheck && !AreCompatible(source1, source2)) {
2652 Error("operator>=(const TMatrixT&,const TMatrixT&)", "matrices not compatible");
2653 return target;
2654 }
2655
2656 target.ResizeTo(source1);
2657
2658 const Element *sp1 = source1.GetMatrixArray();
2659 const Element *sp2 = source2.GetMatrixArray();
2660 Element *tp = target.GetMatrixArray();
2661 const Element *const tp_last = tp + target.GetNoElements();
2662 while (tp < tp_last) {
2663 *tp++ = (*sp1) >= (*sp2);
2664 sp1++;
2665 sp2++;
2666 }
2667
2668 return target;
2669}
2670
2671////////////////////////////////////////////////////////////////////////////////
2672/// logical operation source1 >= source2
2673
2674template <class Element>
2675TMatrixT<Element> TMatrixTAutoloadOps::operator>=(const TMatrixT<Element> &source1, const TMatrixTSym<Element> &source2)
2676{
2677 TMatrixT<Element> target;
2678
2679 if (gMatrixCheck && !AreCompatible(source1, source2)) {
2680 Error("operator>=(const TMatrixT&,const TMatrixTSym&)", "matrices not compatible");
2681 return target;
2682 }
2683
2684 target.ResizeTo(source1);
2685
2686 const Element *sp1 = source1.GetMatrixArray();
2687 const Element *sp2 = source2.GetMatrixArray();
2688 Element *tp = target.GetMatrixArray();
2689 const Element *const tp_last = tp + target.GetNoElements();
2690 while (tp < tp_last) {
2691 *tp++ = (*sp1) >= (*sp2);
2692 sp1++;
2693 sp2++;
2694 }
2695
2696 return target;
2697}
2698
2699////////////////////////////////////////////////////////////////////////////////
2700/// logical operation source1 >= source2
2701
2702template <class Element>
2703TMatrixT<Element> TMatrixTAutoloadOps::operator>=(const TMatrixTSym<Element> &source1, const TMatrixT<Element> &source2)
2704{
2705 return operator<(source2, source1);
2706}
2707
2708////////////////////////////////////////////////////////////////////////////////
2709/// logical operation source1 <= source2
2710
2711template <class Element>
2712TMatrixT<Element> TMatrixTAutoloadOps::operator<=(const TMatrixT<Element> &source1, const TMatrixT<Element> &source2)
2713{
2714 TMatrixT<Element> target;
2715
2716 if (gMatrixCheck && !AreCompatible(source1, source2)) {
2717 Error("operator<=(const TMatrixT&,const TMatrixT&)", "matrices not compatible");
2718 return target;
2719 }
2720
2721 target.ResizeTo(source1);
2722
2723 const Element *sp1 = source1.GetMatrixArray();
2724 const Element *sp2 = source2.GetMatrixArray();
2725 Element *tp = target.GetMatrixArray();
2726 const Element *const tp_last = tp + target.GetNoElements();
2727 while (tp < tp_last) {
2728 *tp++ = (*sp1) <= (*sp2);
2729 sp1++;
2730 sp2++;
2731 }
2732
2733 return target;
2734}
2735
2736////////////////////////////////////////////////////////////////////////////////
2737/// logical operation source1 <= source2
2738
2739template <class Element>
2740TMatrixT<Element> TMatrixTAutoloadOps::operator<=(const TMatrixT<Element> &source1, const TMatrixTSym<Element> &source2)
2741{
2742 TMatrixT<Element> target;
2743
2744 if (gMatrixCheck && !AreCompatible(source1, source2)) {
2745 Error("operator<=(const TMatrixT&,const TMatrixTSym&)", "matrices not compatible");
2746 return target;
2747 }
2748
2749 target.ResizeTo(source1);
2750
2751 const Element *sp1 = source1.GetMatrixArray();
2752 const Element *sp2 = source2.GetMatrixArray();
2753 Element *tp = target.GetMatrixArray();
2754 const Element *const tp_last = tp + target.GetNoElements();
2755 while (tp < tp_last) {
2756 *tp++ = (*sp1) <= (*sp2);
2757 sp1++;
2758 sp2++;
2759 }
2760
2761 return target;
2762}
2763
2764////////////////////////////////////////////////////////////////////////////////
2765/// logical operation source1 <= source2
2766
2767template <class Element>
2768TMatrixT<Element> TMatrixTAutoloadOps::operator<=(const TMatrixTSym<Element> &source1, const TMatrixT<Element> &source2)
2769{
2770 return operator>(source2, source1);
2771}
2772
2773////////////////////////////////////////////////////////////////////////////////
2774/// logical operation source1 < source2
2775
2776template <class Element>
2777TMatrixT<Element> TMatrixTAutoloadOps::operator<(const TMatrixT<Element> &source1, const TMatrixT<Element> &source2)
2778{
2779 TMatrixT<Element> target;
2780
2781 if (gMatrixCheck && !AreCompatible(source1, source2)) {
2782 Error("operator<(const TMatrixT&,const TMatrixT&)", "matrices not compatible");
2783 return target;
2784 }
2785
2786 const Element *sp1 = source1.GetMatrixArray();
2787 const Element *sp2 = source2.GetMatrixArray();
2788 Element *tp = target.GetMatrixArray();
2789 const Element *const tp_last = tp + target.GetNoElements();
2790 while (tp < tp_last) {
2791 *tp++ = (*sp1) < (*sp2);
2792 sp1++;
2793 sp2++;
2794 }
2795
2796 return target;
2797}
2798
2799////////////////////////////////////////////////////////////////////////////////
2800/// logical operation source1 < source2
2801
2802template <class Element>
2803TMatrixT<Element> TMatrixTAutoloadOps::operator<(const TMatrixT<Element> &source1, const TMatrixTSym<Element> &source2)
2804{
2805 TMatrixT<Element> target;
2806
2807 if (gMatrixCheck && !AreCompatible(source1, source2)) {
2808 Error("operator<(const TMatrixT&,const TMatrixTSym&)", "matrices not compatible");
2809 return target;
2810 }
2811
2812 target.ResizeTo(source1);
2813
2814 const Element *sp1 = source1.GetMatrixArray();
2815 const Element *sp2 = source2.GetMatrixArray();
2816 Element *tp = target.GetMatrixArray();
2817 const Element *const tp_last = tp + target.GetNoElements();
2818 while (tp < tp_last) {
2819 *tp++ = (*sp1) < (*sp2);
2820 sp1++;
2821 sp2++;
2822 }
2823
2824 return target;
2825}
2826
2827////////////////////////////////////////////////////////////////////////////////
2828/// logical operation source1 < source2
2829
2830template <class Element>
2831TMatrixT<Element> TMatrixTAutoloadOps::operator<(const TMatrixTSym<Element> &source1, const TMatrixT<Element> &source2)
2832{
2833 return operator>=(source2, source1);
2834}
2835
2836////////////////////////////////////////////////////////////////////////////////
2837/// logical operation source1 != source2
2838
2839template <class Element>
2840TMatrixT<Element> TMatrixTAutoloadOps::operator!=(const TMatrixT<Element> &source1, const TMatrixT<Element> &source2)
2841{
2842 TMatrixT<Element> target;
2843
2844 if (gMatrixCheck && !AreCompatible(source1, source2)) {
2845 Error("operator!=(const TMatrixT&,const TMatrixT&)", "matrices not compatible");
2846 return target;
2847 }
2848
2849 target.ResizeTo(source1);
2850
2851 const Element *sp1 = source1.GetMatrixArray();
2852 const Element *sp2 = source2.GetMatrixArray();
2853 Element *tp = target.GetMatrixArray();
2854 const Element *const tp_last = tp + target.GetNoElements();
2855 while (tp != tp_last) {
2856 *tp++ = (*sp1) != (*sp2);
2857 sp1++;
2858 sp2++;
2859 }
2860
2861 return target;
2862}
2863
2864////////////////////////////////////////////////////////////////////////////////
2865/// logical operation source1 != source2
2866
2867template <class Element>
2868TMatrixT<Element> TMatrixTAutoloadOps::operator!=(const TMatrixT<Element> &source1, const TMatrixTSym<Element> &source2)
2869{
2870 TMatrixT<Element> target;
2871
2872 if (gMatrixCheck && !AreCompatible(source1, source2)) {
2873 Error("operator!=(const TMatrixT&,const TMatrixTSym&)", "matrices not compatible");
2874 return target;
2875 }
2876
2877 target.ResizeTo(source1);
2878
2879 const Element *sp1 = source1.GetMatrixArray();
2880 const Element *sp2 = source2.GetMatrixArray();
2881 Element *tp = target.GetMatrixArray();
2882 const Element *const tp_last = tp + target.GetNoElements();
2883 while (tp != tp_last) {
2884 *tp++ = (*sp1) != (*sp2);
2885 sp1++;
2886 sp2++;
2887 }
2888
2889 return target;
2890}
2891
2892////////////////////////////////////////////////////////////////////////////////
2893/// logical operation source1 != source2
2894
2895template <class Element>
2900
2901/*
2902////////////////////////////////////////////////////////////////////////////////
2903/// logical operation source1 != val
2904
2905template<class Element>
2906TMatrixT<Element> operator!=(const TMatrixT<Element> &source1,Element val)
2907{
2908 TMatrixT<Element> target; target.ResizeTo(source1);
2909
2910 const Element *sp = source1.GetMatrixArray();
2911 Element *tp = target.GetMatrixArray();
2912 const Element * const tp_last = tp+target.GetNoElements();
2913 while (tp != tp_last) {
2914 *tp++ = (*sp != val); sp++;
2915 }
2916
2917 return target;
2918}
2919
2920////////////////////////////////////////////////////////////////////////////////
2921/// logical operation source1 != val
2922
2923template<class Element>
2924TMatrixT<Element> operator!=(Element val,const TMatrixT<Element> &source1)
2925{
2926 return operator!=(source1,val);
2927}
2928*/
2929
2930////////////////////////////////////////////////////////////////////////////////
2931/// Modify addition: target += scalar * source.
2932
2933template <class Element>
2934TMatrixT<Element> &TMatrixTAutoloadOps::Add(TMatrixT<Element> &target, Element scalar, const TMatrixT<Element> &source)
2935{
2936 if (gMatrixCheck && !AreCompatible(target, source)) {
2937 ::Error("Add(TMatrixT &,Element,const TMatrixT &)", "matrices not compatible");
2938 return target;
2939 }
2940
2941 const Element *sp = source.GetMatrixArray();
2942 Element *tp = target.GetMatrixArray();
2943 const Element *ftp = tp + target.GetNoElements();
2944 if (scalar == 0) {
2945 while (tp < ftp)
2946 *tp++ = scalar * (*sp++);
2947 } else if (scalar == 1.) {
2948 while (tp < ftp)
2949 *tp++ = (*sp++);
2950 } else {
2951 while (tp < ftp)
2952 *tp++ += scalar * (*sp++);
2953 }
2954
2955 return target;
2956}
2957
2958////////////////////////////////////////////////////////////////////////////////
2959/// Modify addition: target += scalar * source.
2960
2961template <class Element>
2962TMatrixT<Element> &
2963TMatrixTAutoloadOps::Add(TMatrixT<Element> &target, Element scalar, const TMatrixTSym<Element> &source)
2964{
2965 if (gMatrixCheck && !AreCompatible(target, source)) {
2966 ::Error("Add(TMatrixT &,Element,const TMatrixTSym &)", "matrices not compatible");
2967 return target;
2968 }
2969
2970 const Element *sp = source.GetMatrixArray();
2971 Element *tp = target.GetMatrixArray();
2972 const Element *ftp = tp + target.GetNoElements();
2973 while (tp < ftp)
2974 *tp++ += scalar * (*sp++);
2975
2976 return target;
2977}
2978
2979////////////////////////////////////////////////////////////////////////////////
2980/// Multiply target by the source, element-by-element.
2981
2982template <class Element>
2983TMatrixT<Element> &TMatrixTAutoloadOps::ElementMult(TMatrixT<Element> &target, const TMatrixT<Element> &source)
2984{
2985 if (gMatrixCheck && !AreCompatible(target, source)) {
2986 ::Error("ElementMult(TMatrixT &,const TMatrixT &)", "matrices not compatible");
2987 return target;
2988 }
2989
2990 const Element *sp = source.GetMatrixArray();
2991 Element *tp = target.GetMatrixArray();
2992 const Element *ftp = tp + target.GetNoElements();
2993 while (tp < ftp)
2994 *tp++ *= *sp++;
2995
2996 return target;
2997}
2998
2999////////////////////////////////////////////////////////////////////////////////
3000/// Multiply target by the source, element-by-element.
3001
3002template <class Element>
3003TMatrixT<Element> &TMatrixTAutoloadOps::ElementMult(TMatrixT<Element> &target, const TMatrixTSym<Element> &source)
3004{
3005 if (gMatrixCheck && !AreCompatible(target, source)) {
3006 ::Error("ElementMult(TMatrixT &,const TMatrixTSym &)", "matrices not compatible");
3007 return target;
3008 }
3009
3010 const Element *sp = source.GetMatrixArray();
3011 Element *tp = target.GetMatrixArray();
3012 const Element *ftp = tp + target.GetNoElements();
3013 while (tp < ftp)
3014 *tp++ *= *sp++;
3015
3016 return target;
3017}
3018
3019////////////////////////////////////////////////////////////////////////////////
3020/// Divide target by the source, element-by-element.
3021
3022template <class Element>
3023TMatrixT<Element> &TMatrixTAutoloadOps::ElementDiv(TMatrixT<Element> &target, const TMatrixT<Element> &source)
3024{
3025 if (gMatrixCheck && !AreCompatible(target, source)) {
3026 ::Error("ElementDiv(TMatrixT &,const TMatrixT &)", "matrices not compatible");
3027 return target;
3028 }
3029
3030 const Element *sp = source.GetMatrixArray();
3031 Element *tp = target.GetMatrixArray();
3032 const Element *ftp = tp + target.GetNoElements();
3033 while (tp < ftp) {
3034 if (*sp != 0.0)
3035 *tp++ /= *sp++;
3036 else {
3037 const Int_t irow = (sp - source.GetMatrixArray()) / source.GetNcols();
3038 const Int_t icol = (sp - source.GetMatrixArray()) % source.GetNcols();
3039 Error("ElementDiv", "source (%d,%d) is zero", irow, icol);
3040 tp++;
3041 }
3042 }
3043
3044 return target;
3045}
3046
3047////////////////////////////////////////////////////////////////////////////////
3048/// Multiply target by the source, element-by-element.
3049
3050template <class Element>
3051TMatrixT<Element> &TMatrixTAutoloadOps::ElementDiv(TMatrixT<Element> &target, const TMatrixTSym<Element> &source)
3052{
3053 if (gMatrixCheck && !AreCompatible(target, source)) {
3054 ::Error("ElementDiv(TMatrixT &,const TMatrixTSym &)", "matrices not compatible");
3055 return target;
3056 }
3057
3058 const Element *sp = source.GetMatrixArray();
3059 Element *tp = target.GetMatrixArray();
3060 const Element *ftp = tp + target.GetNoElements();
3061 while (tp < ftp) {
3062 if (*sp != 0.0)
3063 *tp++ /= *sp++;
3064 else {
3065 const Int_t irow = (sp - source.GetMatrixArray()) / source.GetNcols();
3066 const Int_t icol = (sp - source.GetMatrixArray()) % source.GetNcols();
3067 Error("ElementDiv", "source (%d,%d) is zero", irow, icol);
3068 *tp++ = 0.0;
3069 }
3070 }
3071
3072 return target;
3073}
3074
3075////////////////////////////////////////////////////////////////////////////////
3076/// Elementary routine to calculate matrix multiplication A*B
3077
3078template <class Element>
3079void TMatrixTAutoloadOps::AMultB(const Element *const ap, Int_t na, Int_t ncolsa, const Element *const bp, Int_t /*nb*/,
3080 Int_t ncolsb, Element *cp)
3081{
3082 const Int_t M = na / ncolsa;
3083 const Int_t N = ncolsa;
3084 const Int_t P = ncolsb;
3085
3086 if (M <= 12 && N <= 12 && P <= 12) {
3087 for (Int_t i = 0; i < M; ++i) {
3088 for (Int_t j = 0; j < P; ++j) {
3089 Element sum = Element(0);
3090 for (Int_t k = 0; k < N; ++k) {
3091 sum += ap[i * N + k] * bp[k * P + j];
3092 }
3093 cp[i * P + j] = sum;
3094 }
3095 }
3096 return;
3097 }
3098 const Int_t BLOCK = (M >= 192 && N >= 192 && P >= 192) ? 48 : 32;
3099 for (Int_t i0 = 0; i0 < M; i0 += BLOCK) {
3100 for (Int_t j0 = 0; j0 < P; j0 += BLOCK) {
3101 const Int_t i1 = (i0 + BLOCK < M) ? i0 + BLOCK : M;
3102 const Int_t j1 = (j0 + BLOCK < P) ? j0 + BLOCK : P;
3103 for (Int_t i = i0; i < i1; ++i) {
3104 Int_t j = j0;
3105 for (; j <= j1 - 4; j += 4) {
3106 cp[i * P + j + 0] = Element(0);
3107 cp[i * P + j + 1] = Element(0);
3108 cp[i * P + j + 2] = Element(0);
3109 cp[i * P + j + 3] = Element(0);
3110 }
3111 for (; j < j1; ++j)
3112 cp[i * P + j] = Element(0);
3113 }
3114
3115 // ────────────────────── Accumulate over k blocks ──────────────────────
3116 for (Int_t k0 = 0; k0 < N; k0 += BLOCK) {
3117 const Int_t k1 = (k0 + BLOCK < N) ? k0 + BLOCK : N;
3118
3119 for (Int_t i = i0; i < i1; ++i) {
3120 for (Int_t k = k0; k < k1; ++k) {
3121 const Element aik = ap[i * N + k];
3122
3123 Int_t j = j0;
3124 // Main 4-wide accumulation
3125 for (; j <= j1 - 4; j += 4) {
3126 cp[i * P + j + 0] += aik * bp[k * P + j + 0];
3127 cp[i * P + j + 1] += aik * bp[k * P + j + 1];
3128 cp[i * P + j + 2] += aik * bp[k * P + j + 2];
3129 cp[i * P + j + 3] += aik * bp[k * P + j + 3];
3130 }
3131 // Remainder
3132 for (; j < j1; ++j)
3133 cp[i * P + j] += aik * bp[k * P + j];
3134 }
3135 }
3136 }
3137 }
3138 }
3139}
3140
3141////////////////////////////////////////////////////////////////////////////////
3142/// Elementary routine to calculate matrix multiplication A^T*B
3143
3144template <class Element>
3145void TMatrixTAutoloadOps::AtMultB(const Element *const ap, Int_t ncolsa, const Element *const bp, Int_t nb,
3146 Int_t ncolsb, Element *cp)
3147{
3148 const Element *acp0 = ap; // Pointer to A[i,0];
3149 while (acp0 < ap + ncolsa) {
3150 for (const Element *bcp = bp; bcp < bp + ncolsb;) { // Pointer to the j-th column of B, Start bcp = B[0,0]
3151 const Element *acp = acp0; // Pointer to the i-th column of A, reset to A[0,i]
3152 Element cij = 0;
3153 while (bcp < bp + nb) { // Scan the i-th column of A and
3154 cij += *acp * *bcp; // the j-th col of B
3155 acp += ncolsa;
3156 bcp += ncolsb;
3157 }
3158 *cp++ = cij;
3159 bcp -= nb - 1; // Set bcp to the (j+1)-th col
3160 }
3161 acp0++; // Set acp0 to the (i+1)-th col
3162 }
3163}
3164
3165////////////////////////////////////////////////////////////////////////////////
3166/// Elementary routine to calculate matrix multiplication A*B^T
3167
3168template <class Element>
3169void TMatrixTAutoloadOps::AMultBt(const Element *const ap, Int_t na, Int_t ncolsa, const Element *const bp, Int_t nb,
3170 Int_t ncolsb, Element *cp)
3171{
3172 const Element *arp0 = ap; // Pointer to A[i,0];
3173 while (arp0 < ap + na) {
3174 const Element *brp0 = bp; // Pointer to B[j,0];
3175 while (brp0 < bp + nb) {
3176 const Element *arp = arp0; // Pointer to the i-th row of A, reset to A[i,0]
3177 const Element *brp = brp0; // Pointer to the j-th row of B, reset to B[j,0]
3178 Element cij = 0;
3179 while (brp < brp0 + ncolsb) // Scan the i-th row of A and
3180 cij += *arp++ * *brp++; // the j-th row of B
3181 *cp++ = cij;
3182 brp0 += ncolsb; // Set brp0 to the (j+1)-th row
3183 }
3184 arp0 += ncolsa; // Set arp0 to the (i+1)-th row
3185 }
3186}
3187
3188////////////////////////////////////////////////////////////////////////////////
3189/// Stream an object of class TMatrixT.
3190
3191template <class Element>
3192void TMatrixT<Element>::Streamer(TBuffer &R__b)
3193{
3194 if (R__b.IsReading()) {
3195 UInt_t R__s, R__c;
3196 Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
3197 if (R__v > 2) {
3198 Clear();
3199 R__b.ReadClassBuffer(TMatrixT<Element>::Class(), this, R__v, R__s, R__c);
3200 } else if (R__v == 2) { // process old version 2
3201 Clear();
3202 TObject::Streamer(R__b);
3203 this->MakeValid();
3204 R__b >> this->fNrows;
3205 R__b >> this->fNcols;
3206 R__b >> this->fNelems;
3207 R__b >> this->fRowLwb;
3208 R__b >> this->fColLwb;
3209 Char_t isArray;
3210 R__b >> isArray;
3211 if (isArray) {
3212 if (this->fNelems > 0) {
3213 fElements = new Element[this->fNelems];
3214 R__b.ReadFastArray(fElements, this->fNelems);
3215 } else
3216 fElements = nullptr;
3217 }
3218 R__b.CheckByteCount(R__s, R__c, TMatrixT<Element>::IsA());
3219 } else { //====process old versions before automatic schema evolution
3220 TObject::Streamer(R__b);
3221 this->MakeValid();
3222 R__b >> this->fNrows;
3223 R__b >> this->fNcols;
3224 R__b >> this->fRowLwb;
3225 R__b >> this->fColLwb;
3226 this->fNelems = R__b.ReadArray(fElements);
3227 R__b.CheckByteCount(R__s, R__c, TMatrixT<Element>::IsA());
3228 }
3229 // in version <=2 , the matrix was stored column-wise
3230 if (R__v <= 2 && fElements) {
3231 for (Int_t i = 0; i < this->fNrows; i++) {
3232 const Int_t off_i = i * this->fNcols;
3233 for (Int_t j = i; j < this->fNcols; j++) {
3234 const Int_t off_j = j * this->fNrows;
3235 const Element tmp = fElements[off_i + j];
3236 fElements[off_i + j] = fElements[off_j + i];
3237 fElements[off_j + i] = tmp;
3238 }
3239 }
3240 }
3241 if (this->fNelems > 0 && this->fNelems <= this->kSizeMax) {
3242 if (fElements) {
3243 memcpy(fDataStack, fElements, this->fNelems * sizeof(Element));
3244 delete[] fElements;
3245 }
3246 fElements = fDataStack;
3247 } else if (this->fNelems < 0)
3248 this->Invalidate();
3249 } else {
3250 R__b.WriteClassBuffer(TMatrixT<Element>::Class(), this);
3251 }
3252}
3253
3254template class TMatrixT<Float_t>;
3255
3256#include "TMatrixFfwd.h"
3257#include "TMatrixFSymfwd.h"
3258
3259template TMatrixF TMatrixTAutoloadOps::operator+ <Float_t>(const TMatrixF &source1, const TMatrixF &source2);
3260template TMatrixF TMatrixTAutoloadOps::operator+ <Float_t>(const TMatrixF &source1, const TMatrixFSym &source2);
3261template TMatrixF TMatrixTAutoloadOps::operator+ <Float_t>(const TMatrixFSym &source1, const TMatrixF &source2);
3262template TMatrixF TMatrixTAutoloadOps::operator+ <Float_t>(const TMatrixF &source, Float_t val);
3263template TMatrixF TMatrixTAutoloadOps::operator+ <Float_t>(Float_t val, const TMatrixF &source);
3264template TMatrixF TMatrixTAutoloadOps::operator- <Float_t>(const TMatrixF &source1, const TMatrixF &source2);
3265template TMatrixF TMatrixTAutoloadOps::operator- <Float_t>(const TMatrixF &source1, const TMatrixFSym &source2);
3266template TMatrixF TMatrixTAutoloadOps::operator- <Float_t>(const TMatrixFSym &source1, const TMatrixF &source2);
3267template TMatrixF TMatrixTAutoloadOps::operator- <Float_t>(const TMatrixF &source, Float_t val);
3268template TMatrixF TMatrixTAutoloadOps::operator- <Float_t>(Float_t val, const TMatrixF &source);
3269template TMatrixF TMatrixTAutoloadOps::operator* <Float_t>(Float_t val, const TMatrixF &source);
3270template TMatrixF TMatrixTAutoloadOps::operator* <Float_t>(const TMatrixF &source, Float_t val);
3271template TMatrixF TMatrixTAutoloadOps::operator* <Float_t>(const TMatrixF &source1, const TMatrixF &source2);
3272template TMatrixF TMatrixTAutoloadOps::operator* <Float_t>(const TMatrixF &source1, const TMatrixFSym &source2);
3273template TMatrixF TMatrixTAutoloadOps::operator* <Float_t>(const TMatrixFSym &source1, const TMatrixF &source2);
3274template TMatrixF TMatrixTAutoloadOps::operator* <Float_t>(const TMatrixFSym &source1, const TMatrixFSym &source2);
3275template TMatrixF TMatrixTAutoloadOps::operator&& <Float_t>(const TMatrixF &source1, const TMatrixF &source2);
3276template TMatrixF TMatrixTAutoloadOps::operator&& <Float_t>(const TMatrixF &source1, const TMatrixFSym &source2);
3277template TMatrixF TMatrixTAutoloadOps::operator&& <Float_t>(const TMatrixFSym &source1, const TMatrixF &source2);
3278template TMatrixF TMatrixTAutoloadOps::operator|| <Float_t>(const TMatrixF &source1, const TMatrixF &source2);
3279template TMatrixF TMatrixTAutoloadOps::operator|| <Float_t>(const TMatrixF &source1, const TMatrixFSym &source2);
3280template TMatrixF TMatrixTAutoloadOps::operator|| <Float_t>(const TMatrixFSym &source1, const TMatrixF &source2);
3281template TMatrixF TMatrixTAutoloadOps::operator><Float_t>(const TMatrixF & source1, const TMatrixF & source2);
3282template TMatrixF TMatrixTAutoloadOps::operator><Float_t>(const TMatrixF & source1, const TMatrixFSym & source2);
3283template TMatrixF TMatrixTAutoloadOps::operator><Float_t>(const TMatrixFSym & source1, const TMatrixF & source2);
3284template TMatrixF TMatrixTAutoloadOps::operator>= <Float_t>(const TMatrixF &source1, const TMatrixF &source2);
3285template TMatrixF TMatrixTAutoloadOps::operator>= <Float_t>(const TMatrixF &source1, const TMatrixFSym &source2);
3286template TMatrixF TMatrixTAutoloadOps::operator>= <Float_t>(const TMatrixFSym &source1, const TMatrixF &source2);
3287template TMatrixF TMatrixTAutoloadOps::operator<= <Float_t>(const TMatrixF &source1, const TMatrixF &source2);
3288template TMatrixF TMatrixTAutoloadOps::operator<= <Float_t>(const TMatrixF &source1, const TMatrixFSym &source2);
3289template TMatrixF TMatrixTAutoloadOps::operator<= <Float_t>(const TMatrixFSym &source1, const TMatrixF &source2);
3290template TMatrixF TMatrixTAutoloadOps::operator< <Float_t>(const TMatrixF &source1, const TMatrixF &source2);
3291template TMatrixF TMatrixTAutoloadOps::operator< <Float_t>(const TMatrixF &source1, const TMatrixFSym &source2);
3292template TMatrixF TMatrixTAutoloadOps::operator< <Float_t>(const TMatrixFSym &source1, const TMatrixF &source2);
3293template TMatrixF TMatrixTAutoloadOps::operator!= <Float_t>(const TMatrixF &source1, const TMatrixF &source2);
3294template TMatrixF TMatrixTAutoloadOps::operator!= <Float_t>(const TMatrixF &source1, const TMatrixFSym &source2);
3295template TMatrixF TMatrixTAutoloadOps::operator!= <Float_t>(const TMatrixFSym &source1, const TMatrixF &source2);
3296
3297template TMatrixF &TMatrixTAutoloadOps::Add<Float_t>(TMatrixF &target, Float_t scalar, const TMatrixF &source);
3298template TMatrixF &TMatrixTAutoloadOps::Add<Float_t>(TMatrixF &target, Float_t scalar, const TMatrixFSym &source);
3299template TMatrixF &TMatrixTAutoloadOps::ElementMult<Float_t>(TMatrixF &target, const TMatrixF &source);
3300template TMatrixF &TMatrixTAutoloadOps::ElementMult<Float_t>(TMatrixF &target, const TMatrixFSym &source);
3301template TMatrixF &TMatrixTAutoloadOps::ElementDiv<Float_t>(TMatrixF &target, const TMatrixF &source);
3302template TMatrixF &TMatrixTAutoloadOps::ElementDiv<Float_t>(TMatrixF &target, const TMatrixFSym &source);
3303
3304template void TMatrixTAutoloadOps::AMultB<Float_t>(const Float_t *const ap, Int_t na, Int_t ncolsa,
3305 const Float_t *const bp, Int_t nb, Int_t ncolsb, Float_t *cp);
3306template void TMatrixTAutoloadOps::AtMultB<Float_t>(const Float_t *const ap, Int_t ncolsa, const Float_t *const bp,
3307 Int_t nb, Int_t ncolsb, Float_t *cp);
3308template void TMatrixTAutoloadOps::AMultBt<Float_t>(const Float_t *const ap, Int_t na, Int_t ncolsa,
3309 const Float_t *const bp, Int_t nb, Int_t ncolsb, Float_t *cp);
3310
3311#include "TMatrixDfwd.h"
3312#include "TMatrixDSymfwd.h"
3313
3314template class TMatrixT<Double_t>;
3315
3316template TMatrixD TMatrixTAutoloadOps::operator+ <Double_t>(const TMatrixD &source1, const TMatrixD &source2);
3317template TMatrixD TMatrixTAutoloadOps::operator+ <Double_t>(const TMatrixD &source1, const TMatrixDSym &source2);
3318template TMatrixD TMatrixTAutoloadOps::operator+ <Double_t>(const TMatrixDSym &source1, const TMatrixD &source2);
3319template TMatrixD TMatrixTAutoloadOps::operator+ <Double_t>(const TMatrixD &source, Double_t val);
3320template TMatrixD TMatrixTAutoloadOps::operator+ <Double_t>(Double_t val, const TMatrixD &source);
3321template TMatrixD TMatrixTAutoloadOps::operator- <Double_t>(const TMatrixD &source1, const TMatrixD &source2);
3322template TMatrixD TMatrixTAutoloadOps::operator- <Double_t>(const TMatrixD &source1, const TMatrixDSym &source2);
3323template TMatrixD TMatrixTAutoloadOps::operator- <Double_t>(const TMatrixDSym &source1, const TMatrixD &source2);
3324template TMatrixD TMatrixTAutoloadOps::operator- <Double_t>(const TMatrixD &source, Double_t val);
3325template TMatrixD TMatrixTAutoloadOps::operator- <Double_t>(Double_t val, const TMatrixD &source);
3326template TMatrixD TMatrixTAutoloadOps::operator* <Double_t>(Double_t val, const TMatrixD &source);
3327template TMatrixD TMatrixTAutoloadOps::operator* <Double_t>(const TMatrixD &source, Double_t val);
3328template TMatrixD TMatrixTAutoloadOps::operator* <Double_t>(const TMatrixD &source1, const TMatrixD &source2);
3329template TMatrixD TMatrixTAutoloadOps::operator* <Double_t>(const TMatrixD &source1, const TMatrixDSym &source2);
3330template TMatrixD TMatrixTAutoloadOps::operator* <Double_t>(const TMatrixDSym &source1, const TMatrixD &source2);
3331template TMatrixD TMatrixTAutoloadOps::operator* <Double_t>(const TMatrixDSym &source1, const TMatrixDSym &source2);
3332template TMatrixD TMatrixTAutoloadOps::operator&& <Double_t>(const TMatrixD &source1, const TMatrixD &source2);
3333template TMatrixD TMatrixTAutoloadOps::operator&& <Double_t>(const TMatrixD &source1, const TMatrixDSym &source2);
3334template TMatrixD TMatrixTAutoloadOps::operator&& <Double_t>(const TMatrixDSym &source1, const TMatrixD &source2);
3335template TMatrixD TMatrixTAutoloadOps::operator|| <Double_t>(const TMatrixD &source1, const TMatrixD &source2);
3336template TMatrixD TMatrixTAutoloadOps::operator|| <Double_t>(const TMatrixD &source1, const TMatrixDSym &source2);
3337template TMatrixD TMatrixTAutoloadOps::operator|| <Double_t>(const TMatrixDSym &source1, const TMatrixD &source2);
3338template TMatrixD TMatrixTAutoloadOps::operator><Double_t>(const TMatrixD & source1, const TMatrixD & source2);
3339template TMatrixD TMatrixTAutoloadOps::operator><Double_t>(const TMatrixD & source1, const TMatrixDSym & source2);
3340template TMatrixD TMatrixTAutoloadOps::operator><Double_t>(const TMatrixDSym & source1, const TMatrixD & source2);
3341template TMatrixD TMatrixTAutoloadOps::operator>= <Double_t>(const TMatrixD &source1, const TMatrixD &source2);
3342template TMatrixD TMatrixTAutoloadOps::operator>= <Double_t>(const TMatrixD &source1, const TMatrixDSym &source2);
3343template TMatrixD TMatrixTAutoloadOps::operator>= <Double_t>(const TMatrixDSym &source1, const TMatrixD &source2);
3344template TMatrixD TMatrixTAutoloadOps::operator<= <Double_t>(const TMatrixD &source1, const TMatrixD &source2);
3345template TMatrixD TMatrixTAutoloadOps::operator<= <Double_t>(const TMatrixD &source1, const TMatrixDSym &source2);
3346template TMatrixD TMatrixTAutoloadOps::operator<= <Double_t>(const TMatrixDSym &source1, const TMatrixD &source2);
3347template TMatrixD TMatrixTAutoloadOps::operator< <Double_t>(const TMatrixD &source1, const TMatrixD &source2);
3348template TMatrixD TMatrixTAutoloadOps::operator< <Double_t>(const TMatrixD &source1, const TMatrixDSym &source2);
3349template TMatrixD TMatrixTAutoloadOps::operator< <Double_t>(const TMatrixDSym &source1, const TMatrixD &source2);
3350template TMatrixD TMatrixTAutoloadOps::operator!= <Double_t>(const TMatrixD &source1, const TMatrixD &source2);
3351template TMatrixD TMatrixTAutoloadOps::operator!= <Double_t>(const TMatrixD &source1, const TMatrixDSym &source2);
3352template TMatrixD TMatrixTAutoloadOps::operator!= <Double_t>(const TMatrixDSym &source1, const TMatrixD &source2);
3353
3354template TMatrixD &TMatrixTAutoloadOps::Add<Double_t>(TMatrixD &target, Double_t scalar, const TMatrixD &source);
3355template TMatrixD &TMatrixTAutoloadOps::Add<Double_t>(TMatrixD &target, Double_t scalar, const TMatrixDSym &source);
3356template TMatrixD &TMatrixTAutoloadOps::ElementMult<Double_t>(TMatrixD &target, const TMatrixD &source);
3357template TMatrixD &TMatrixTAutoloadOps::ElementMult<Double_t>(TMatrixD &target, const TMatrixDSym &source);
3358template TMatrixD &TMatrixTAutoloadOps::ElementDiv<Double_t>(TMatrixD &target, const TMatrixD &source);
3359template TMatrixD &TMatrixTAutoloadOps::ElementDiv<Double_t>(TMatrixD &target, const TMatrixDSym &source);
3360
3361template void TMatrixTAutoloadOps::AMultB<Double_t>(const Double_t *const ap, Int_t na, Int_t ncolsa,
3362 const Double_t *const bp, Int_t nb, Int_t ncolsb, Double_t *cp);
3363template void TMatrixTAutoloadOps::AtMultB<Double_t>(const Double_t *const ap, Int_t ncolsa, const Double_t *const bp,
3364 Int_t nb, Int_t ncolsb, Double_t *cp);
3365template void TMatrixTAutoloadOps::AMultBt<Double_t>(const Double_t *const ap, Int_t na, Int_t ncolsa,
3366 const Double_t *const bp, Int_t nb, Int_t ncolsb, Double_t *cp);
b
#define b(i)
Definition RSha256.hxx:100
a
#define a(i)
Definition RSha256.hxx:99
size
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
Int_t
int Int_t
Signed integer 4 bytes (int)
Definition RtypesCore.h:59
Version_t
short Version_t
Class version identifier (short)
Definition RtypesCore.h:79
Char_t
char Char_t
Character 1 byte (char)
Definition RtypesCore.h:51
Long_t
long Long_t
Signed long integer 4 bytes (long). Size depends on architecture.
Definition RtypesCore.h:68
Float_t
float Float_t
Float 4 bytes (float)
Definition RtypesCore.h:71
kFALSE
constexpr Bool_t kFALSE
Definition RtypesCore.h:108
Double_t
double Double_t
Double 8 bytes.
Definition RtypesCore.h:73
kTRUE
constexpr Bool_t kTRUE
Definition RtypesCore.h:107
Option_t
const char Option_t
Option string (const char)
Definition RtypesCore.h:80
kPlus
@ kPlus
Definition TAttMarker.h:56
TRangeDynCast
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
Definition TCollection.h:360
operator<=
Bool_t operator<=(const TDatime &d1, const TDatime &d2)
Definition TDatime.h:108
operator>
Bool_t operator>(const TDatime &d1, const TDatime &d2)
Definition TDatime.h:110
operator!=
Bool_t operator!=(const TDatime &d1, const TDatime &d2)
Definition TDatime.h:104
operator>=
Bool_t operator>=(const TDatime &d1, const TDatime &d2)
Definition TDatime.h:112
operator<
Bool_t operator<(const TDatime &d1, const TDatime &d2)
Definition TDatime.h:106
R__ASSERT
#define R__ASSERT(e)
Checks condition e and reports a fatal error if it's false.
Definition TError.h:125
Error
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition TError.cxx:208
Warning
void Warning(const char *location, const char *msgfmt,...)
Use this function in warning situations.
Definition TError.cxx:252
N
#define N
option
Option_t Option_t option
Definition TGWin32VirtualXProxy.cxx:44
data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Definition TGWin32VirtualXProxy.cxx:104
target
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
Definition TGWin32VirtualXProxy.cxx:247
index
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t index
Definition TGWin32VirtualXProxy.cxx:168
gMatrixCheck
R__EXTERN Int_t gMatrixCheck
Definition TMatrixTBase.h:82
operator*
TTime operator*(const TTime &t1, const TTime &t2)
Definition TTime.h:85
operator-
TTime operator-(const TTime &t1, const TTime &t2)
Definition TTime.h:83
TBuffer
Buffer base class used for serializing objects.
Definition TBuffer.h:43
TDecompLU
LU Decomposition class.
Definition TDecompLU.h:24
TDecompLU::InvertLU
static Bool_t InvertLU(TMatrixD &a, Double_t tol, Double_t *det=nullptr)
Calculate matrix inversion through in place forward/backward substitution.
Definition TDecompLU.cxx:774
TMatrixDEigen
TMatrixDEigen.
Definition TMatrixDEigen.h:27
TMatrixTBase
TMatrixTBase.
Definition TMatrixTBase.h:84
TMatrixTBase::SetMatrixArray
virtual TMatrixTBase< Element > & SetMatrixArray(const Element *data, Option_t *option="")
Copy array data to matrix .
Definition TMatrixTBase.cxx:187
TMatrixTSym
TMatrixTSym.
Definition TMatrixTSym.h:36
TMatrixT::GetSub
TMatrixTBase< Element > & GetSub(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb, TMatrixTBase< Element > &target, Option_t *option="S") const override
Get submatrix [row_lwb..row_upb] x [col_lwb..col_upb]; The indexing range of the returned matrix depe...
Definition TMatrixT.cxx:1080
TMatrixT::Rank1Update
TMatrixT< Element > & Rank1Update(const TVectorT< Element > &v, Element alpha=1.0)
Perform a rank 1 operation on matrix A: A += alpha * v * v^T.
Definition TMatrixT.cxx:1524
TMatrixT::Use
TMatrixT< Element > & Use(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb, Element *data)
Use the array data to fill the matrix ([row_lwb..row_upb] x [col_lwb..col_upb])
Definition TMatrixT.cxx:1047
TMatrixT::Delete_m
void Delete_m(Int_t size, Element *&)
Delete data pointer m, if it was assigned on the heap.
Definition TMatrixT.cxx:404
TMatrixT::Streamer
void Streamer(TBuffer &) override
Stream an object of class TMatrixT.
Definition TMatrixT.cxx:3192
TMatrixT::operator*=
TMatrixT< Element > & operator*=(Element val)
Multiply every element of the matrix with val.
Definition TMatrixT.cxx:1856
TMatrixT::operator=
TMatrixT< Element > & operator=(const TMatrixT< Element > &source)
Assignment operator.
Definition TMatrixT.cxx:1717
TMatrixT::NormByRow
TMatrixT< Element > & NormByRow(const TVectorT< Element > &v, Option_t *option="D")
Multiply/divide matrix rows with a vector: option: "D" : b(i,j) = a(i,j)/v(j) i = 0,...
Definition TMatrixT.cxx:1674
TMatrixT::EMatrixCreatorsOp2
EMatrixCreatorsOp2
Definition TMatrixT.h:60
TMatrixT::Minus
void Minus(const TMatrixT< Element > &a, const TMatrixT< Element > &b)
General matrix subtraction. Replace this matrix with C such that C = A - B.
Definition TMatrixT.cxx:574
TMatrixT::Similarity
Element Similarity(const TVectorT< Element > &v) const
Calculate scalar v * (*this) * v^T.
Definition TMatrixT.cxx:1587
TMatrixT::Plus
void Plus(const TMatrixT< Element > &a, const TMatrixT< Element > &b)
General matrix summation. Replace this matrix with C such that C = A + B.
Definition TMatrixT.cxx:506
TMatrixT::operator-=
TMatrixT< Element > & operator-=(Element val)
Subtract val from every element of the matrix.
Definition TMatrixT.cxx:1840
TMatrixT::Determinant
Double_t Determinant() const override
Return the matrix determinant.
Definition TMatrixT.cxx:1353
TMatrixT::Transpose
TMatrixT< Element > & Transpose(const TMatrixT< Element > &source)
Transpose matrix source.
Definition TMatrixT.cxx:1453
TMatrixT::MultT
void MultT(const TMatrixT< Element > &a, const TMatrixT< Element > &b)
General matrix multiplication. Replace this matrix with C such that C = A * B^T.
Definition TMatrixT.cxx:945
TMatrixT::InvertFast
TMatrixT< Element > & InvertFast(Double_t *det=nullptr)
Invert the matrix and calculate its determinant, however upto (6x6) a fast Cramer inversion is used .
Definition TMatrixT.cxx:1402
TMatrixT::SetSub
TMatrixTBase< Element > & SetSub(Int_t row_lwb, Int_t col_lwb, const TMatrixTBase< Element > &source) override
Insert matrix source starting at [row_lwb][col_lwb], thereby overwriting the part [row_lwb....
Definition TMatrixT.cxx:1147
TMatrixT::New_m
Element * New_m(Int_t size)
Return data pointer .
Definition TMatrixT.cxx:418
TMatrixT::EMatrixCreatorsOp1
EMatrixCreatorsOp1
Definition TMatrixT.h:59
TMatrixT::ResizeTo
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...
Definition TMatrixT.cxx:1201
TMatrixT::Invert
TMatrixT< Element > & Invert(Double_t *det=nullptr)
Invert the matrix and calculate its determinant.
Definition TMatrixT.cxx:1388
TMatrixT::Allocate
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.
Definition TMatrixT.cxx:462
TMatrixT::TMatrixT
TMatrixT()
Definition TMatrixT.h:62
TMatrixT::operator+=
TMatrixT< Element > & operator+=(Element val)
Add val to every element of the matrix.
Definition TMatrixT.cxx:1824
TMatrixT::NormByColumn
TMatrixT< Element > & NormByColumn(const TVectorT< Element > &v, Option_t *option="D")
Multiply/divide matrix columns by a vector: option: "D" : b(i,j) = a(i,j)/v(i) i = 0,...
Definition TMatrixT.cxx:1628
TMatrixT::TMult
void TMult(const TMatrixT< Element > &a, const TMatrixT< Element > &b)
Replace this matrix with C such that C = A' * B.
Definition TMatrixT.cxx:846
TMatrixT::EigenVectors
const TMatrixT< Element > EigenVectors(TVectorT< Element > &eigenValues) const
Return a matrix containing the eigen-vectors ordered by descending values of Re^2+Im^2 of the complex...
Definition TMatrixT.cxx:2287
TMatrixT::Mult
void Mult(const TMatrixT< Element > &a, const TMatrixT< Element > &b)
General matrix multiplication. Replace this matrix with C such that C = A * B.
Definition TMatrixT.cxx:642
TMatrixT::Memcpy_m
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 .
Definition TMatrixT.cxx:437
TMatrixT::operator/=
TMatrixT< Element > & operator/=(const TMatrixTDiag_const< Element > &diag)
Divide a matrix row by the diagonal of another matrix matrix(i,j) /= diag(j)
Definition TMatrixT.cxx:2107
TObject::operator=
TObject & operator=(const TObject &rhs) noexcept
TObject assignment operator.
Definition TObject.h:299
TObject::Streamer
virtual void Streamer(TBuffer &)
Stream an object of class TObject.
Definition TObject.cxx:989
TString
Basic string class.
Definition TString.h:138
TString::ToUpper
void ToUpper()
Change string to upper case.
Definition TString.cxx:1202
TString::Contains
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition TString.h:641
ROOT::operator+
ELogLevel operator+(ELogLevel severity, int offset)
Definition RLogger.hxx:42
TMath::Max
Short_t Max(Short_t a, Short_t b)
Returns the largest of a and b.
Definition TMathBase.h:249
TMath::Power
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition TMath.h:732
TMath::Min
Short_t Min(Short_t a, Short_t b)
Returns the smallest of a and b.
Definition TMathBase.h:197
TMatrixTAutoloadOps::operator!=
TMatrixT< Element > operator!=(const TMatrixT< Element > &source1, const TMatrixT< Element > &source2)
logical operation source1 != source2
Definition TMatrixT.cxx:2840
TMatrixTAutoloadOps::AMultBt
void AMultBt(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^T.
Definition TMatrixT.cxx:3169
TMatrixTAutoloadOps::AMultB
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.
Definition TMatrixT.cxx:3079
TMatrixTAutoloadOps::operator+
TMatrixT< Element > operator+(const TMatrixT< Element > &source1, const TMatrixT< Element > &source2)
operation this = source1+source2
Definition TMatrixT.cxx:2301
TMatrixTAutoloadOps::AtMultB
void AtMultB(const Element *const ap, Int_t ncolsa, const Element *const bp, Int_t nb, Int_t ncolsb, Element *cp)
Elementary routine to calculate matrix multiplication A^T*B.
Definition TMatrixT.cxx:3145
TMatrixTAutoloadOps::operator>=
TMatrixT< Element > operator>=(const TMatrixT< Element > &source1, const TMatrixT< Element > &source2)
logical operation source1 >= source2
Definition TMatrixT.cxx:2647
TMatrixTAutoloadOps::ElementMult
TMatrixT< Element > & ElementMult(TMatrixT< Element > &target, const TMatrixT< Element > &source)
Multiply target by the source, element-by-element.
Definition TMatrixT.cxx:2983
TMatrixTAutoloadOps::operator||
TMatrixT< Element > operator||(const TMatrixT< Element > &source1, const TMatrixT< Element > &source2)
Logical OR.
Definition TMatrixT.cxx:2523
TMatrixTAutoloadOps::operator<
TMatrixT< Element > operator<(const TMatrixT< Element > &source1, const TMatrixT< Element > &source2)
logical operation source1 < source2
Definition TMatrixT.cxx:2777
TMatrixTAutoloadOps::operator>
TMatrixT< Element > operator>(const TMatrixT< Element > &source1, const TMatrixT< Element > &source2)
logical operation source1 > source2
Definition TMatrixT.cxx:2582
TMatrixTAutoloadOps::Add
TMatrixT< Element > & Add(TMatrixT< Element > &target, Element scalar, const TMatrixT< Element > &source)
Modify addition: target += scalar * source.
Definition TMatrixT.cxx:2934
TMatrixTAutoloadOps::ElementDiv
TMatrixT< Element > & ElementDiv(TMatrixT< Element > &target, const TMatrixT< Element > &source)
Divide target by the source, element-by-element.
Definition TMatrixT.cxx:3023
TMatrixTAutoloadOps::Mult
Element1 Mult(const TVectorT< Element1 > &v1, const TMatrixT< Element2 > &m, const TVectorT< Element3 > &v2)
Perform v1 * M * v2, a scalar result.
Definition TVectorT.cxx:1541
TMatrixTAutoloadOps::operator-
TMatrixT< Element > operator-(const TMatrixT< Element > &source1, const TMatrixT< Element > &source2)
operation this = source1-source2
Definition TMatrixT.cxx:2352
TMatrixTAutoloadOps::operator<=
TMatrixT< Element > operator<=(const TMatrixT< Element > &source1, const TMatrixT< Element > &source2)
logical operation source1 <= source2
Definition TMatrixT.cxx:2712
TMatrixTAutoloadOps::operator&&
TMatrixT< Element > operator&&(const TMatrixT< Element > &source1, const TMatrixT< Element > &source2)
Logical AND.
Definition TMatrixT.cxx:2464
TMatrixTAutoloadOps::operator*
TMatrixT< Element > operator*(Element val, const TMatrixT< Element > &source)
operation this = val*source
Definition TMatrixT.cxx:2403
TMatrixTAutoloadOps::AreCompatible
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 .
Definition TMatrixTBase.cxx:887
v2
@ v2
Definition rootcling_impl.cxx:3686
v
@ v
Definition rootcling_impl.cxx:3683
v1
@ v1
Definition rootcling_impl.cxx:3685
m
TMarker m
Definition textangle.C:8
sum
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2339