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