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