Logo ROOT  

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