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