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