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