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