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