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