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