Logo ROOT   master
Reference Guide
TMatrixTSym.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 TMatrixTSym
13  \ingroup Matrix
14 
15  TMatrixTSym
16 
17  Template class of a symmetric matrix in the linear algebra package.
18 
19  See \ref MatrixPage for the documentation of the linear algebra package.
20 
21  Note that in this implementation both matrix element m[i][j] and
22  m[j][i] are updated and stored in memory . However, when making the
23  object persistent only the upper right triangle is stored .
24 
25 */
26 
27 #include "TMatrixTSym.h"
28 #include "TBuffer.h"
29 #include "TMatrixTLazy.h"
30 #include "TMatrixTSymCramerInv.h"
31 #include "TDecompLU.h"
32 #include "TMatrixDSymEigen.h"
33 #include "TClass.h"
34 #include "TMath.h"
35 
37 
38 
39 ////////////////////////////////////////////////////////////////////////////////
40 
41 template<class Element>
43 {
44  Allocate(no_rows,no_rows,0,0,1);
45 }
46 
47 ////////////////////////////////////////////////////////////////////////////////
48 
49 template<class Element>
51 {
52  const Int_t no_rows = row_upb-row_lwb+1;
53  Allocate(no_rows,no_rows,row_lwb,row_lwb,1);
54 }
55 
56 ////////////////////////////////////////////////////////////////////////////////
57 /// option=
58 /// - "F": array elements contains the matrix stored column-wise
59 /// like in Fortran, so a[i,j] = elements[i+no_rows*j],
60 /// - else it is supposed that array elements are stored row-wise
61 /// a[i,j] = elements[i*no_cols+j]
62 ///
63 /// array elements are copied
64 
65 template<class Element>
66 TMatrixTSym<Element>::TMatrixTSym(Int_t no_rows,const Element *elements,Option_t *option)
67 {
68  Allocate(no_rows,no_rows);
69  SetMatrixArray(elements,option);
70  if (!this->IsSymmetric()) {
71  Error("TMatrixTSym(Int_t,Element*,Option_t*)","matrix not symmetric");
72  }
73 }
74 
75 ////////////////////////////////////////////////////////////////////////////////
76 /// array elements are copied
77 
78 template<class Element>
79 TMatrixTSym<Element>::TMatrixTSym(Int_t row_lwb,Int_t row_upb,const Element *elements,Option_t *option)
80 {
81  const Int_t no_rows = row_upb-row_lwb+1;
82  Allocate(no_rows,no_rows,row_lwb,row_lwb);
83  SetMatrixArray(elements,option);
84  if (!this->IsSymmetric()) {
85  Error("TMatrixTSym(Int_t,Int_t,Element*,Option_t*)","matrix not symmetric");
86  }
87 }
88 
89 ////////////////////////////////////////////////////////////////////////////////
90 
91 template<class Element>
93 {
94  R__ASSERT(another.IsValid());
95  Allocate(another.GetNrows(),another.GetNcols(),another.GetRowLwb(),another.GetColLwb());
96  *this = another;
97 }
98 
99 ////////////////////////////////////////////////////////////////////////////////
100 /// Create a matrix applying a specific operation to the prototype.
101 /// Example: TMatrixTSym<Element> a(10,12); ...; TMatrixTSym<Element> b(TMatrixT::kTransposed, a);
102 /// Supported operations are: kZero, kUnit, kTransposed, kInverted and kAtA.
103 
104 template<class Element>
106 {
107  R__ASSERT(prototype.IsValid());
108 
109  switch(op) {
110  case kZero:
111  Allocate(prototype.GetNrows(),prototype.GetNcols(),
112  prototype.GetRowLwb(),prototype.GetColLwb(),1);
113  break;
114 
115  case kUnit:
116  Allocate(prototype.GetNrows(),prototype.GetNcols(),
117  prototype.GetRowLwb(),prototype.GetColLwb(),1);
118  this->UnitMatrix();
119  break;
120 
121  case kTransposed:
122  Allocate(prototype.GetNcols(), prototype.GetNrows(),
123  prototype.GetColLwb(),prototype.GetRowLwb());
124  Transpose(prototype);
125  break;
126 
127  case kInverted:
128  {
129  Allocate(prototype.GetNrows(),prototype.GetNcols(),
130  prototype.GetRowLwb(),prototype.GetColLwb(),1);
131  *this = prototype;
132  // Since the user can not control the tolerance of this newly created matrix
133  // we put it to the smallest possible number
134  const Element oldTol = this->SetTol(std::numeric_limits<Element>::min());
135  this->Invert();
136  this->SetTol(oldTol);
137  break;
138  }
139 
140  case kAtA:
141  Allocate(prototype.GetNcols(),prototype.GetNcols(),prototype.GetColLwb(),prototype.GetColLwb(),1);
142  TMult(prototype);
143  break;
144 
145  default:
146  Error("TMatrixTSym(EMatrixCreatorOp1,const TMatrixTSym)",
147  "operation %d not yet implemented", op);
148  }
149 }
150 
151 ////////////////////////////////////////////////////////////////////////////////
152 
153 template<class Element>
155 {
156  R__ASSERT(prototype.IsValid());
157 
158  switch(op) {
159  case kAtA:
160  Allocate(prototype.GetNcols(),prototype.GetNcols(),prototype.GetColLwb(),prototype.GetColLwb(),1);
161  TMult(prototype);
162  break;
163 
164  default:
165  Error("TMatrixTSym(EMatrixCreatorOp1,const TMatrixT)",
166  "operation %d not yet implemented", op);
167  }
168 }
169 
170 ////////////////////////////////////////////////////////////////////////////////
171 
172 template<class Element>
174 {
175  R__ASSERT(a.IsValid());
176  R__ASSERT(b.IsValid());
177 
178  switch(op) {
179  case kPlus:
180  {
181  Allocate(a.GetNcols(),a.GetNcols(),a.GetColLwb(),a.GetColLwb(),1);
182  Plus(a,b);
183  break;
184  }
185 
186  case kMinus:
187  {
188  Allocate(a.GetNcols(),a.GetNcols(),a.GetColLwb(),a.GetColLwb(),1);
189  Minus(a,b);
190  break;
191  }
192 
193  default:
194  Error("TMatrixTSym(EMatrixCreatorOp2)", "operation %d not yet implemented", op);
195  }
196 }
197 
198 ////////////////////////////////////////////////////////////////////////////////
199 
200 template<class Element>
202 {
203  Allocate(lazy_constructor.GetRowUpb()-lazy_constructor.GetRowLwb()+1,
204  lazy_constructor.GetRowUpb()-lazy_constructor.GetRowLwb()+1,
205  lazy_constructor.GetRowLwb(),lazy_constructor.GetRowLwb(),1);
206  lazy_constructor.FillIn(*this);
207  if (!this->IsSymmetric()) {
208  Error("TMatrixTSym(TMatrixTSymLazy)","matrix not symmetric");
209  }
210 }
211 
212 ////////////////////////////////////////////////////////////////////////////////
213 /// delete data pointer m, if it was assigned on the heap
214 
215 template<class Element>
217 {
218  if (m) {
219  if (size > this->kSizeMax)
220  delete [] m;
221  m = 0;
222  }
223 }
224 
225 ////////////////////////////////////////////////////////////////////////////////
226 /// return data pointer . if requested size <= kSizeMax, assign pointer
227 /// to the stack space
228 
229 template<class Element>
231 {
232  if (size == 0) return 0;
233  else {
234  if ( size <= this->kSizeMax )
235  return fDataStack;
236  else {
237  Element *heap = new Element[size];
238  return heap;
239  }
240  }
241 }
242 
243 ////////////////////////////////////////////////////////////////////////////////
244 /// copy copySize doubles from *oldp to *newp . However take care of the
245 /// situation where both pointers are assigned to the same stack space
246 
247 template<class Element>
248 Int_t TMatrixTSym<Element>::Memcpy_m(Element *newp,const Element *oldp,Int_t copySize,
249  Int_t newSize,Int_t oldSize)
250 {
251  if (copySize == 0 || oldp == newp)
252  return 0;
253  else {
254  if ( newSize <= this->kSizeMax && oldSize <= this->kSizeMax ) {
255  // both pointers are inside fDataStack, be careful with copy direction !
256  if (newp > oldp) {
257  for (Int_t i = copySize-1; i >= 0; i--)
258  newp[i] = oldp[i];
259  } else {
260  for (Int_t i = 0; i < copySize; i++)
261  newp[i] = oldp[i];
262  }
263  }
264  else
265  memcpy(newp,oldp,copySize*sizeof(Element));
266  }
267  return 0;
268 }
269 
270 ////////////////////////////////////////////////////////////////////////////////
271 /// Allocate new matrix. Arguments are number of rows, columns, row
272 /// lowerbound (0 default) and column lowerbound (0 default).
273 
274 template<class Element>
275 void TMatrixTSym<Element>::Allocate(Int_t no_rows,Int_t no_cols,Int_t row_lwb,Int_t col_lwb,
276  Int_t init,Int_t /*nr_nonzeros*/)
277 {
278  this->fIsOwner = kTRUE;
280  this->fNrows = 0;
281  this->fNcols = 0;
282  this->fRowLwb = 0;
283  this->fColLwb = 0;
284  this->fNelems = 0;
285  fElements = 0;
286 
287  if (no_rows < 0 || no_cols < 0)
288  {
289  Error("Allocate","no_rows=%d no_cols=%d",no_rows,no_cols);
290  this->Invalidate();
291  return;
292  }
293 
294  this->MakeValid();
295  this->fNrows = no_rows;
296  this->fNcols = no_cols;
297  this->fRowLwb = row_lwb;
298  this->fColLwb = col_lwb;
299  this->fNelems = this->fNrows*this->fNcols;
300 
301  if (this->fNelems > 0) {
302  fElements = New_m(this->fNelems);
303  if (init)
304  memset(fElements,0,this->fNelems*sizeof(Element));
305  } else
306  fElements = 0;
307 }
308 
309 ////////////////////////////////////////////////////////////////////////////////
310 /// Symmetric matrix summation. Create a matrix C such that C = A + B.
311 
312 template<class Element>
314 {
315  if (gMatrixCheck) {
316  if (!AreCompatible(a,b)) {
317  Error("Plus","matrices not compatible");
318  return;
319  }
320 
321  if (this->GetMatrixArray() == a.GetMatrixArray()) {
322  Error("Plus","this->GetMatrixArray() == a.GetMatrixArray()");
323  return;
324  }
325 
326  if (this->GetMatrixArray() == b.GetMatrixArray()) {
327  Error("Plus","this->GetMatrixArray() == b.GetMatrixArray()");
328  return;
329  }
330  }
331 
332  const Element * ap = a.GetMatrixArray();
333  const Element * bp = b.GetMatrixArray();
334  Element * cp = this->GetMatrixArray();
335  const Element * const cp_last = cp+this->fNelems;
336 
337  while (cp < cp_last) {
338  *cp = *ap++ + *bp++;
339  cp++;
340  }
341 }
342 
343 ////////////////////////////////////////////////////////////////////////////////
344 /// Symmetric matrix summation. Create a matrix C such that C = A + B.
345 
346 template<class Element>
348 {
349  if (gMatrixCheck) {
350  if (!AreCompatible(a,b)) {
351  Error("Minus","matrices not compatible");
352  return;
353  }
354 
355  if (this->GetMatrixArray() == a.GetMatrixArray()) {
356  Error("Minus","this->GetMatrixArray() == a.GetMatrixArray()");
357  return;
358  }
359 
360  if (this->GetMatrixArray() == b.GetMatrixArray()) {
361  Error("Minus","this->GetMatrixArray() == b.GetMatrixArray()");
362  return;
363  }
364  }
365 
366  const Element * ap = a.GetMatrixArray();
367  const Element * bp = b.GetMatrixArray();
368  Element * cp = this->GetMatrixArray();
369  const Element * const cp_last = cp+this->fNelems;
370 
371  while (cp < cp_last) {
372  *cp = *ap++ - *bp++;
373  cp++;
374  }
375 }
376 
377 ////////////////////////////////////////////////////////////////////////////////
378 /// Create a matrix C such that C = A' * A. In other words,
379 /// c[i,j] = SUM{ a[k,i] * a[k,j] }.
380 
381 template<class Element>
383 {
384  R__ASSERT(a.IsValid());
385 
386 #ifdef CBLAS
387  const Element *ap = a.GetMatrixArray();
388  Element *cp = this->GetMatrixArray();
389  if (typeid(Element) == typeid(Double_t))
390  cblas_dgemm (CblasRowMajor,CblasTrans,CblasNoTrans,this->fNrows,this->fNcols,a.GetNrows(),
391  1.0,ap,a.GetNcols(),ap,a.GetNcols(),1.0,cp,this->fNcols);
392  else if (typeid(Element) != typeid(Float_t))
393  cblas_sgemm (CblasRowMajor,CblasTrans,CblasNoTrans,fNrows,fNcols,a.GetNrows(),
394  1.0,ap,a.GetNcols(),ap,a.GetNcols(),1.0,cp,fNcols);
395  else
396  Error("TMult","type %s not implemented in BLAS library",typeid(Element));
397 #else
398  const Int_t nb = a.GetNoElements();
399  const Int_t ncolsa = a.GetNcols();
400  const Int_t ncolsb = ncolsa;
401  const Element * const ap = a.GetMatrixArray();
402  const Element * const bp = ap;
403  Element * cp = this->GetMatrixArray();
404 
405  const Element *acp0 = ap; // Pointer to A[i,0];
406  while (acp0 < ap+a.GetNcols()) {
407  for (const Element *bcp = bp; bcp < bp+ncolsb; ) { // Pointer to the j-th column of A, Start bcp = A[0,0]
408  const Element *acp = acp0; // Pointer to the i-th column of A, reset to A[0,i]
409  Element cij = 0;
410  while (bcp < bp+nb) { // Scan the i-th column of A and
411  cij += *acp * *bcp; // the j-th col of A
412  acp += ncolsa;
413  bcp += ncolsb;
414  }
415  *cp++ = cij;
416  bcp -= nb-1; // Set bcp to the (j+1)-th col
417  }
418  acp0++; // Set acp0 to the (i+1)-th col
419  }
420 
421  R__ASSERT(cp == this->GetMatrixArray()+this->fNelems && acp0 == ap+ncolsa);
422 #endif
423 }
424 
425 ////////////////////////////////////////////////////////////////////////////////
426 /// Matrix multiplication, with A symmetric
427 /// Create a matrix C such that C = A' * A = A * A = A * A'
428 
429 template<class Element>
431 {
432  R__ASSERT(a.IsValid());
433 
434 #ifdef CBLAS
435  const Element *ap = a.GetMatrixArray();
436  Element *cp = this->GetMatrixArray();
437  if (typeid(Element) == typeid(Double_t))
438  cblas_dsymm (CblasRowMajor,CblasLeft,CblasUpper,this->fNrows,this->fNcols,1.0,
439  ap,a.GetNcols(),ap,a.GetNcols(),0.0,cp,this->fNcols);
440  else if (typeid(Element) != typeid(Float_t))
441  cblas_ssymm (CblasRowMajor,CblasLeft,CblasUpper,fNrows,fNcols,1.0,
442  ap1,a.GetNcols(),ap,a.GetNcols(),0.0,cp,fNcols);
443  else
444  Error("TMult","type %s not implemented in BLAS library",typeid(Element));
445 #else
446  const Int_t nb = a.GetNoElements();
447  const Int_t ncolsa = a.GetNcols();
448  const Int_t ncolsb = ncolsa;
449  const Element * const ap = a.GetMatrixArray();
450  const Element * const bp = ap;
451  Element * cp = this->GetMatrixArray();
452 
453  const Element *acp0 = ap; // Pointer to A[i,0];
454  while (acp0 < ap+a.GetNcols()) {
455  for (const Element *bcp = bp; bcp < bp+ncolsb; ) { // Pointer to the j-th column of A, Start bcp = A[0,0]
456  const Element *acp = acp0; // Pointer to the i-th column of A, reset to A[0,i]
457  Element cij = 0;
458  while (bcp < bp+nb) { // Scan the i-th column of A and
459  cij += *acp * *bcp; // the j-th col of A
460  acp += ncolsa;
461  bcp += ncolsb;
462  }
463  *cp++ = cij;
464  bcp -= nb-1; // Set bcp to the (j+1)-th col
465  }
466  acp0++; // Set acp0 to the (i+1)-th col
467  }
468 
469  R__ASSERT(cp == this->GetMatrixArray()+this->fNelems && acp0 == ap+ncolsa);
470 #endif
471 }
472 
473 ////////////////////////////////////////////////////////////////////////////////
474 
475 template<class Element>
477 {
478  if (gMatrixCheck && row_upb < row_lwb)
479  {
480  Error("Use","row_upb=%d < row_lwb=%d",row_upb,row_lwb);
481  return *this;
482  }
483 
484  this->Clear();
485  this->fNrows = row_upb-row_lwb+1;
486  this->fNcols = this->fNrows;
487  this->fRowLwb = row_lwb;
488  this->fColLwb = row_lwb;
489  this->fNelems = this->fNrows*this->fNcols;
490  fElements = data;
491  this->fIsOwner = kFALSE;
492 
493  return *this;
494 }
495 
496 ////////////////////////////////////////////////////////////////////////////////
497 /// Get submatrix [row_lwb..row_upb][row_lwb..row_upb]; The indexing range of the
498 /// returned matrix depends on the argument option:
499 ///
500 /// option == "S" : return [0..row_upb-row_lwb+1][0..row_upb-row_lwb+1] (default)
501 /// else : return [row_lwb..row_upb][row_lwb..row_upb]
502 
503 template<class Element>
505  TMatrixTSym<Element> &target,Option_t *option) const
506 {
507  if (gMatrixCheck) {
508  R__ASSERT(this->IsValid());
509  if (row_lwb < this->fRowLwb || row_lwb > this->fRowLwb+this->fNrows-1) {
510  Error("GetSub","row_lwb out of bounds");
511  return target;
512  }
513  if (row_upb < this->fRowLwb || row_upb > this->fRowLwb+this->fNrows-1) {
514  Error("GetSub","row_upb out of bounds");
515  return target;
516  }
517  if (row_upb < row_lwb) {
518  Error("GetSub","row_upb < row_lwb");
519  return target;
520  }
521  }
522 
523  TString opt(option);
524  opt.ToUpper();
525  const Int_t shift = (opt.Contains("S")) ? 1 : 0;
526 
527  Int_t row_lwb_sub;
528  Int_t row_upb_sub;
529  if (shift) {
530  row_lwb_sub = 0;
531  row_upb_sub = row_upb-row_lwb;
532  } else {
533  row_lwb_sub = row_lwb;
534  row_upb_sub = row_upb;
535  }
536 
537  target.ResizeTo(row_lwb_sub,row_upb_sub,row_lwb_sub,row_upb_sub);
538  const Int_t nrows_sub = row_upb_sub-row_lwb_sub+1;
539 
540  if (target.GetRowIndexArray() && target.GetColIndexArray()) {
541  for (Int_t irow = 0; irow < nrows_sub; irow++) {
542  for (Int_t icol = 0; icol < nrows_sub; icol++) {
543  target(irow+row_lwb_sub,icol+row_lwb_sub) = (*this)(row_lwb+irow,row_lwb+icol);
544  }
545  }
546  } else {
547  const Element *ap = this->GetMatrixArray()+(row_lwb-this->fRowLwb)*this->fNrows+(row_lwb-this->fRowLwb);
548  Element *bp = target.GetMatrixArray();
549 
550  for (Int_t irow = 0; irow < nrows_sub; irow++) {
551  const Element *ap_sub = ap;
552  for (Int_t icol = 0; icol < nrows_sub; icol++) {
553  *bp++ = *ap_sub++;
554  }
555  ap += this->fNrows;
556  }
557  }
558 
559  return target;
560 }
561 
562 ////////////////////////////////////////////////////////////////////////////////
563 /// Get submatrix [row_lwb..row_upb][col_lwb..col_upb]; The indexing range of the
564 /// returned matrix depends on the argument option:
565 ///
566 /// option == "S" : return [0..row_upb-row_lwb+1][0..col_upb-col_lwb+1] (default)
567 /// else : return [row_lwb..row_upb][col_lwb..col_upb]
568 
569 template<class Element>
571  TMatrixTBase<Element> &target,Option_t *option) const
572 {
573  if (gMatrixCheck) {
574  R__ASSERT(this->IsValid());
575  if (row_lwb < this->fRowLwb || row_lwb > this->fRowLwb+this->fNrows-1) {
576  Error("GetSub","row_lwb out of bounds");
577  return target;
578  }
579  if (col_lwb < this->fColLwb || col_lwb > this->fColLwb+this->fNcols-1) {
580  Error("GetSub","col_lwb out of bounds");
581  return target;
582  }
583  if (row_upb < this->fRowLwb || row_upb > this->fRowLwb+this->fNrows-1) {
584  Error("GetSub","row_upb out of bounds");
585  return target;
586  }
587  if (col_upb < this->fColLwb || col_upb > this->fColLwb+this->fNcols-1) {
588  Error("GetSub","col_upb out of bounds");
589  return target;
590  }
591  if (row_upb < row_lwb || col_upb < col_lwb) {
592  Error("GetSub","row_upb < row_lwb || col_upb < col_lwb");
593  return target;
594  }
595  }
596 
597  TString opt(option);
598  opt.ToUpper();
599  const Int_t shift = (opt.Contains("S")) ? 1 : 0;
600 
601  const Int_t row_lwb_sub = (shift) ? 0 : row_lwb;
602  const Int_t row_upb_sub = (shift) ? row_upb-row_lwb : row_upb;
603  const Int_t col_lwb_sub = (shift) ? 0 : col_lwb;
604  const Int_t col_upb_sub = (shift) ? col_upb-col_lwb : col_upb;
605 
606  target.ResizeTo(row_lwb_sub,row_upb_sub,col_lwb_sub,col_upb_sub);
607  const Int_t nrows_sub = row_upb_sub-row_lwb_sub+1;
608  const Int_t ncols_sub = col_upb_sub-col_lwb_sub+1;
609 
610  if (target.GetRowIndexArray() && target.GetColIndexArray()) {
611  for (Int_t irow = 0; irow < nrows_sub; irow++) {
612  for (Int_t icol = 0; icol < ncols_sub; icol++) {
613  target(irow+row_lwb_sub,icol+col_lwb_sub) = (*this)(row_lwb+irow,col_lwb+icol);
614  }
615  }
616  } else {
617  const Element *ap = this->GetMatrixArray()+(row_lwb-this->fRowLwb)*this->fNcols+(col_lwb-this->fColLwb);
618  Element *bp = target.GetMatrixArray();
619 
620  for (Int_t irow = 0; irow < nrows_sub; irow++) {
621  const Element *ap_sub = ap;
622  for (Int_t icol = 0; icol < ncols_sub; icol++) {
623  *bp++ = *ap_sub++;
624  }
625  ap += this->fNcols;
626  }
627  }
628 
629  return target;
630 }
631 
632 ////////////////////////////////////////////////////////////////////////////////
633 /// Insert matrix source starting at [row_lwb][row_lwb], thereby overwriting the part
634 /// [row_lwb..row_lwb+nrows_source][row_lwb..row_lwb+nrows_source];
635 
636 template<class Element>
638 {
639  if (gMatrixCheck) {
640  R__ASSERT(this->IsValid());
641  R__ASSERT(source.IsValid());
642 
643  if (!source.IsSymmetric()) {
644  Error("SetSub","source matrix is not symmetric");
645  return *this;
646  }
647  if (row_lwb < this->fRowLwb || row_lwb > this->fRowLwb+this->fNrows-1) {
648  Error("SetSub","row_lwb outof bounds");
649  return *this;
650  }
651  if (row_lwb+source.GetNrows() > this->fRowLwb+this->fNrows) {
652  Error("SetSub","source matrix too large");
653  return *this;
654  }
655  }
656 
657  const Int_t nRows_source = source.GetNrows();
658 
659  if (source.GetRowIndexArray() && source.GetColIndexArray()) {
660  const Int_t rowlwb_s = source.GetRowLwb();
661  for (Int_t irow = 0; irow < nRows_source; irow++) {
662  for (Int_t icol = 0; icol < nRows_source; icol++) {
663  (*this)(row_lwb+irow,row_lwb+icol) = source(rowlwb_s+irow,rowlwb_s+icol);
664  }
665  }
666  } else {
667  const Element *bp = source.GetMatrixArray();
668  Element *ap = this->GetMatrixArray()+(row_lwb-this->fRowLwb)*this->fNrows+(row_lwb-this->fRowLwb);
669 
670  for (Int_t irow = 0; irow < nRows_source; irow++) {
671  Element *ap_sub = ap;
672  for (Int_t icol = 0; icol < nRows_source; icol++) {
673  *ap_sub++ = *bp++;
674  }
675  ap += this->fNrows;
676  }
677  }
678 
679  return *this;
680 }
681 
682 ////////////////////////////////////////////////////////////////////////////////
683 /// Insert matrix source starting at [row_lwb][col_lwb] in a symmetric fashion, thereby overwriting the part
684 /// [row_lwb..row_lwb+nrows_source][row_lwb..row_lwb+nrows_source];
685 
686 template<class Element>
688 {
689  if (gMatrixCheck) {
690  R__ASSERT(this->IsValid());
691  R__ASSERT(source.IsValid());
692 
693  if (row_lwb < this->fRowLwb || row_lwb > this->fRowLwb+this->fNrows-1) {
694  Error("SetSub","row_lwb out of bounds");
695  return *this;
696  }
697  if (col_lwb < this->fColLwb || col_lwb > this->fColLwb+this->fNcols-1) {
698  Error("SetSub","col_lwb out of bounds");
699  return *this;
700  }
701 
702  if (row_lwb+source.GetNrows() > this->fRowLwb+this->fNrows || col_lwb+source.GetNcols() > this->fRowLwb+this->fNrows) {
703  Error("SetSub","source matrix too large");
704  return *this;
705  }
706  if (col_lwb+source.GetNcols() > this->fRowLwb+this->fNrows || row_lwb+source.GetNrows() > this->fRowLwb+this->fNrows) {
707  Error("SetSub","source matrix too large");
708  return *this;
709  }
710  }
711 
712  const Int_t nRows_source = source.GetNrows();
713  const Int_t nCols_source = source.GetNcols();
714 
715  const Int_t rowlwb_s = source.GetRowLwb();
716  const Int_t collwb_s = source.GetColLwb();
717  if (row_lwb >= col_lwb) {
718  // lower triangle
719  Int_t irow;
720  for (irow = 0; irow < nRows_source; irow++) {
721  for (Int_t icol = 0; col_lwb+icol <= row_lwb+irow &&
722  icol < nCols_source; icol++) {
723  (*this)(row_lwb+irow,col_lwb+icol) = source(irow+rowlwb_s,icol+collwb_s);
724  }
725  }
726 
727  // upper triangle
728  for (irow = 0; irow < nCols_source; irow++) {
729  for (Int_t icol = nRows_source-1; row_lwb+icol > irow+col_lwb &&
730  icol >= 0; icol--) {
731  (*this)(col_lwb+irow,row_lwb+icol) = source(icol+rowlwb_s,irow+collwb_s);
732  }
733  }
734  } else {
735 
736  }
737 
738  return *this;
739 }
740 
741 ////////////////////////////////////////////////////////////////////////////////
742 
743 template<class Element>
745 {
747  if (!this->IsSymmetric()) {
748  Error("SetMatrixArray","Matrix is not symmetric after Set");
749  }
750 
751  return *this;
752 }
753 
754 ////////////////////////////////////////////////////////////////////////////////
755 
756 template<class Element>
758 {
759  if (row_shift != col_shift) {
760  Error("Shift","row_shift != col_shift");
761  return *this;
762  }
763  return TMatrixTBase<Element>::Shift(row_shift,col_shift);
764 }
765 
766 ////////////////////////////////////////////////////////////////////////////////
767 /// Set size of the matrix to nrows x ncols
768 /// New dynamic elements are created, the overlapping part of the old ones are
769 /// copied to the new structures, then the old elements are deleted.
770 
771 template<class Element>
773 {
774  R__ASSERT(this->IsValid());
775  if (!this->fIsOwner) {
776  Error("ResizeTo(Int_t,Int_t)","Not owner of data array,cannot resize");
777  return *this;
778  }
779 
780  if (nrows != ncols) {
781  Error("ResizeTo(Int_t,Int_t)","nrows != ncols");
782  return *this;
783  }
784 
785  if (this->fNelems > 0) {
786  if (this->fNrows == nrows && this->fNcols == ncols)
787  return *this;
788  else if (nrows == 0 || ncols == 0) {
789  this->fNrows = nrows; this->fNcols = ncols;
790  this->Clear();
791  return *this;
792  }
793 
794  Element *elements_old = GetMatrixArray();
795  const Int_t nelems_old = this->fNelems;
796  const Int_t nrows_old = this->fNrows;
797  const Int_t ncols_old = this->fNcols;
798 
799  Allocate(nrows,ncols);
800  R__ASSERT(this->IsValid());
801 
802  Element *elements_new = GetMatrixArray();
803  // new memory should be initialized but be careful not to wipe out the stack
804  // storage. Initialize all when old or new storage was on the heap
805  if (this->fNelems > this->kSizeMax || nelems_old > this->kSizeMax)
806  memset(elements_new,0,this->fNelems*sizeof(Element));
807  else if (this->fNelems > nelems_old)
808  memset(elements_new+nelems_old,0,(this->fNelems-nelems_old)*sizeof(Element));
809 
810  // Copy overlap
811  const Int_t ncols_copy = TMath::Min(this->fNcols,ncols_old);
812  const Int_t nrows_copy = TMath::Min(this->fNrows,nrows_old);
813 
814  const Int_t nelems_new = this->fNelems;
815  if (ncols_old < this->fNcols) {
816  for (Int_t i = nrows_copy-1; i >= 0; i--) {
817  Memcpy_m(elements_new+i*this->fNcols,elements_old+i*ncols_old,ncols_copy,
818  nelems_new,nelems_old);
819  if (this->fNelems <= this->kSizeMax && nelems_old <= this->kSizeMax)
820  memset(elements_new+i*this->fNcols+ncols_copy,0,(this->fNcols-ncols_copy)*sizeof(Element));
821  }
822  } else {
823  for (Int_t i = 0; i < nrows_copy; i++)
824  Memcpy_m(elements_new+i*this->fNcols,elements_old+i*ncols_old,ncols_copy,
825  nelems_new,nelems_old);
826  }
827 
828  Delete_m(nelems_old,elements_old);
829  } else {
830  Allocate(nrows,ncols,0,0,1);
831  }
832 
833  return *this;
834 }
835 
836 ////////////////////////////////////////////////////////////////////////////////
837 /// Set size of the matrix to [row_lwb:row_upb] x [col_lwb:col_upb]
838 /// New dynamic elements are created, the overlapping part of the old ones are
839 /// copied to the new structures, then the old elements are deleted.
840 
841 template<class Element>
843  Int_t /*nr_nonzeros*/)
844 {
845  R__ASSERT(this->IsValid());
846  if (!this->fIsOwner) {
847  Error("ResizeTo(Int_t,Int_t,Int_t,Int_t)","Not owner of data array,cannot resize");
848  return *this;
849  }
850 
851  if (row_lwb != col_lwb) {
852  Error("ResizeTo(Int_t,Int_t,Int_t,Int_t)","row_lwb != col_lwb");
853  return *this;
854  }
855  if (row_upb != col_upb) {
856  Error("ResizeTo(Int_t,Int_t,Int_t,Int_t)","row_upb != col_upb");
857  return *this;
858  }
859 
860  const Int_t new_nrows = row_upb-row_lwb+1;
861  const Int_t new_ncols = col_upb-col_lwb+1;
862 
863  if (this->fNelems > 0) {
864 
865  if (this->fNrows == new_nrows && this->fNcols == new_ncols &&
866  this->fRowLwb == row_lwb && this->fColLwb == col_lwb)
867  return *this;
868  else if (new_nrows == 0 || new_ncols == 0) {
869  this->fNrows = new_nrows; this->fNcols = new_ncols;
870  this->fRowLwb = row_lwb; this->fColLwb = col_lwb;
871  this->Clear();
872  return *this;
873  }
874 
875  Element *elements_old = GetMatrixArray();
876  const Int_t nelems_old = this->fNelems;
877  const Int_t nrows_old = this->fNrows;
878  const Int_t ncols_old = this->fNcols;
879  const Int_t rowLwb_old = this->fRowLwb;
880  const Int_t colLwb_old = this->fColLwb;
881 
882  Allocate(new_nrows,new_ncols,row_lwb,col_lwb);
883  R__ASSERT(this->IsValid());
884 
885  Element *elements_new = GetMatrixArray();
886  // new memory should be initialized but be careful ot to wipe out the stack
887  // storage. Initialize all when old or new storage was on the heap
888  if (this->fNelems > this->kSizeMax || nelems_old > this->kSizeMax)
889  memset(elements_new,0,this->fNelems*sizeof(Element));
890  else if (this->fNelems > nelems_old)
891  memset(elements_new+nelems_old,0,(this->fNelems-nelems_old)*sizeof(Element));
892 
893  // Copy overlap
894  const Int_t rowLwb_copy = TMath::Max(this->fRowLwb,rowLwb_old);
895  const Int_t colLwb_copy = TMath::Max(this->fColLwb,colLwb_old);
896  const Int_t rowUpb_copy = TMath::Min(this->fRowLwb+this->fNrows-1,rowLwb_old+nrows_old-1);
897  const Int_t colUpb_copy = TMath::Min(this->fColLwb+this->fNcols-1,colLwb_old+ncols_old-1);
898 
899  const Int_t nrows_copy = rowUpb_copy-rowLwb_copy+1;
900  const Int_t ncols_copy = colUpb_copy-colLwb_copy+1;
901 
902  if (nrows_copy > 0 && ncols_copy > 0) {
903  const Int_t colOldOff = colLwb_copy-colLwb_old;
904  const Int_t colNewOff = colLwb_copy-this->fColLwb;
905  if (ncols_old < this->fNcols) {
906  for (Int_t i = nrows_copy-1; i >= 0; i--) {
907  const Int_t iRowOld = rowLwb_copy+i-rowLwb_old;
908  const Int_t iRowNew = rowLwb_copy+i-this->fRowLwb;
909  Memcpy_m(elements_new+iRowNew*this->fNcols+colNewOff,
910  elements_old+iRowOld*ncols_old+colOldOff,ncols_copy,this->fNelems,nelems_old);
911  if (this->fNelems <= this->kSizeMax && nelems_old <= this->kSizeMax)
912  memset(elements_new+iRowNew*this->fNcols+colNewOff+ncols_copy,0,
913  (this->fNcols-ncols_copy)*sizeof(Element));
914  }
915  } else {
916  for (Int_t i = 0; i < nrows_copy; i++) {
917  const Int_t iRowOld = rowLwb_copy+i-rowLwb_old;
918  const Int_t iRowNew = rowLwb_copy+i-this->fRowLwb;
919  Memcpy_m(elements_new+iRowNew*this->fNcols+colNewOff,
920  elements_old+iRowOld*ncols_old+colOldOff,ncols_copy,this->fNelems,nelems_old);
921  }
922  }
923  }
924 
925  Delete_m(nelems_old,elements_old);
926  } else {
927  Allocate(new_nrows,new_ncols,row_lwb,col_lwb,1);
928  }
929 
930  return *this;
931 }
932 
933 ////////////////////////////////////////////////////////////////////////////////
934 
935 template<class Element>
937 {
938  const TMatrixT<Element> &tmp = *this;
939  TDecompLU lu(tmp,this->fTol);
940  Double_t d1,d2;
941  lu.Det(d1,d2);
942  return d1*TMath::Power(2.0,d2);
943 }
944 
945 ////////////////////////////////////////////////////////////////////////////////
946 
947 template<class Element>
949 {
950  const TMatrixT<Element> &tmp = *this;
951  TDecompLU lu(tmp,this->fTol);
952  lu.Det(d1,d2);
953 }
954 
955 ////////////////////////////////////////////////////////////////////////////////
956 /// Invert the matrix and calculate its determinant
957 /// Notice that the LU decomposition is used instead of Bunch-Kaufman
958 /// Bunch-Kaufman guarantees a symmetric inverted matrix but is slower than LU .
959 /// The user can access Bunch-Kaufman through the TDecompBK class .
960 
961 template<class Element>
963 {
964  R__ASSERT(this->IsValid());
965  TMatrixD tmp(*this);
966  if (TDecompLU::InvertLU(tmp,Double_t(this->fTol),det)) {
967  const Double_t *p1 = tmp.GetMatrixArray();
968  Element *p2 = this->GetMatrixArray();
969  for (Int_t i = 0; i < this->GetNoElements(); i++)
970  p2[i] = p1[i];
971  }
972 
973  return *this;
974 }
975 
976 ////////////////////////////////////////////////////////////////////////////////
977 /// Invert the matrix and calculate its determinant
978 
979 template<class Element>
981 {
982  R__ASSERT(this->IsValid());
983 
984  const Char_t nRows = Char_t(this->GetNrows());
985  switch (nRows) {
986  case 1:
987  {
988  Element *pM = this->GetMatrixArray();
989  if (*pM == 0.) {
990  Error("InvertFast","matrix is singular");
991  *det = 0;
992  } else {
993  *det = *pM;
994  *pM = 1.0/(*pM);
995  }
996  return *this;
997  }
998  case 2:
999  {
1000  TMatrixTSymCramerInv::Inv2x2<Element>(*this,det);
1001  return *this;
1002  }
1003  case 3:
1004  {
1005  TMatrixTSymCramerInv::Inv3x3<Element>(*this,det);
1006  return *this;
1007  }
1008  case 4:
1009  {
1010  TMatrixTSymCramerInv::Inv4x4<Element>(*this,det);
1011  return *this;
1012  }
1013  case 5:
1014  {
1015  TMatrixTSymCramerInv::Inv5x5<Element>(*this,det);
1016  return *this;
1017  }
1018  case 6:
1019  {
1020  TMatrixTSymCramerInv::Inv6x6<Element>(*this,det);
1021  return *this;
1022  }
1023 
1024  default:
1025  {
1026  TMatrixD tmp(*this);
1027  if (TDecompLU::InvertLU(tmp,Double_t(this->fTol),det)) {
1028  const Double_t *p1 = tmp.GetMatrixArray();
1029  Element *p2 = this->GetMatrixArray();
1030  for (Int_t i = 0; i < this->GetNoElements(); i++)
1031  p2[i] = p1[i];
1032  }
1033  return *this;
1034  }
1035  }
1036 }
1037 
1038 ////////////////////////////////////////////////////////////////////////////////
1039 /// Transpose a matrix.
1040 
1041 template<class Element>
1043 {
1044  if (gMatrixCheck) {
1045  R__ASSERT(this->IsValid());
1046  R__ASSERT(source.IsValid());
1047 
1048  if (this->fNrows != source.GetNcols() || this->fRowLwb != source.GetColLwb())
1049  {
1050  Error("Transpose","matrix has wrong shape");
1051  return *this;
1052  }
1053  }
1054 
1055  *this = source;
1056  return *this;
1057 }
1058 
1059 ////////////////////////////////////////////////////////////////////////////////
1060 /// Perform a rank 1 operation on the matrix:
1061 /// A += alpha * v * v^T
1062 
1063 template<class Element>
1065 {
1066  if (gMatrixCheck) {
1067  R__ASSERT(this->IsValid());
1068  R__ASSERT(v.IsValid());
1069  if (v.GetNoElements() < this->fNrows) {
1070  Error("Rank1Update","vector too short");
1071  return *this;
1072  }
1073  }
1074 
1075  const Element * const pv = v.GetMatrixArray();
1076  Element *trp = this->GetMatrixArray(); // pointer to UR part and diagonal, traverse row-wise
1077  Element *tcp = trp; // pointer to LL part, traverse col-wise
1078  for (Int_t i = 0; i < this->fNrows; i++) {
1079  trp += i; // point to [i,i]
1080  tcp += i*this->fNcols; // point to [i,i]
1081  const Element tmp = alpha*pv[i];
1082  for (Int_t j = i; j < this->fNcols; j++) {
1083  if (j > i) *tcp += tmp*pv[j];
1084  *trp++ += tmp*pv[j];
1085  tcp += this->fNcols;
1086  }
1087  tcp -= this->fNelems-1; // point to [0,i]
1088  }
1089 
1090  return *this;
1091 }
1092 
1093 ////////////////////////////////////////////////////////////////////////////////
1094 /// Calculate B * (*this) * B^T , final matrix will be (nrowsb x nrowsb)
1095 /// This is a similarity transform when B is orthogonal . It is more
1096 /// efficient than applying the actual multiplication because this
1097 /// routine realizes that the final matrix is symmetric .
1098 
1099 template<class Element>
1101 {
1102  if (gMatrixCheck) {
1103  R__ASSERT(this->IsValid());
1104  R__ASSERT(b.IsValid());
1105  if (this->fNcols != b.GetNcols() || this->fColLwb != b.GetColLwb()) {
1106  Error("Similarity(const TMatrixT &)","matrices incompatible");
1107  return *this;
1108  }
1109  }
1110 
1111  const Int_t ncolsa = this->fNcols;
1112  const Int_t nb = b.GetNoElements();
1113  const Int_t nrowsb = b.GetNrows();
1114  const Int_t ncolsb = b.GetNcols();
1115 
1116  const Element * const bp = b.GetMatrixArray();
1117 
1118  Element work[kWorkMax];
1119  Bool_t isAllocated = kFALSE;
1120  Element *bap = work;
1121  if (nrowsb*ncolsa > kWorkMax) {
1122  isAllocated = kTRUE;
1123  bap = new Element[nrowsb*ncolsa];
1124  }
1125 
1126  AMultB(bp,nb,ncolsb,this->fElements,this->fNelems,this->fNcols,bap);
1127 
1128  if (nrowsb != this->fNrows)
1129  this->ResizeTo(nrowsb,nrowsb);
1130 
1131 #ifdef CBLAS
1132  Element *cp = this->GetMatrixArray();
1133  if (typeid(Element) == typeid(Double_t))
1134  cblas_dgemm (CblasRowMajor,CblasNoTrans,CblasTrans,this->fNrows,this->fNcols,ba.GetNcols(),
1135  1.0,bap,ba.GetNcols(),bp,b.GetNcols(),1.0,cp,this->fNcols);
1136  else if (typeid(Element) != typeid(Float_t))
1137  cblas_sgemm (CblasRowMajor,CblasNoTrans,CblasTrans,this->fNrows,this->fNcols,ba.GetNcols(),
1138  1.0,bap,ba.GetNcols(),bp,b.GetNcols(),1.0,cp,this->fNcols);
1139  else
1140  Error("Similarity","type %s not implemented in BLAS library",typeid(Element));
1141 #else
1142  const Int_t nba = nrowsb*ncolsa;
1143  const Int_t ncolsba = ncolsa;
1144  const Element * bi1p = bp;
1145  Element * cp = this->GetMatrixArray();
1146  Element * const cp0 = cp;
1147 
1148  Int_t ishift = 0;
1149  const Element *barp0 = bap;
1150  while (barp0 < bap+nba) {
1151  const Element *brp0 = bi1p;
1152  while (brp0 < bp+nb) {
1153  const Element *barp = barp0;
1154  const Element *brp = brp0;
1155  Element cij = 0;
1156  while (brp < brp0+ncolsb)
1157  cij += *barp++ * *brp++;
1158  *cp++ = cij;
1159  brp0 += ncolsb;
1160  }
1161  barp0 += ncolsba;
1162  bi1p += ncolsb;
1163  cp += ++ishift;
1164  }
1165 
1166  R__ASSERT(cp == cp0+this->fNelems+ishift && barp0 == bap+nba);
1167 
1168  cp = cp0;
1169  for (Int_t irow = 0; irow < this->fNrows; irow++) {
1170  const Int_t rowOff1 = irow*this->fNrows;
1171  for (Int_t icol = 0; icol < irow; icol++) {
1172  const Int_t rowOff2 = icol*this->fNrows;
1173  cp[rowOff1+icol] = cp[rowOff2+irow];
1174  }
1175  }
1176 #endif
1177 
1178  if (isAllocated)
1179  delete [] bap;
1180 
1181  return *this;
1182 }
1183 
1184 ////////////////////////////////////////////////////////////////////////////////
1185 /// Calculate B * (*this) * B^T , final matrix will be (nrowsb x nrowsb)
1186 /// This is a similarity transform when B is orthogonal . It is more
1187 /// efficient than applying the actual multiplication because this
1188 /// routine realizes that the final matrix is symmetric .
1189 
1190 template<class Element>
1192 {
1193  if (gMatrixCheck) {
1194  R__ASSERT(this->IsValid());
1195  R__ASSERT(b.IsValid());
1196  if (this->fNcols != b.GetNcols() || this->fColLwb != b.GetColLwb()) {
1197  Error("Similarity(const TMatrixTSym &)","matrices incompatible");
1198  return *this;
1199  }
1200  }
1201 
1202 #ifdef CBLAS
1203  const Int_t nrowsb = b.GetNrows();
1204  const Int_t ncolsa = this->GetNcols();
1205 
1206  Element work[kWorkMax];
1207  Bool_t isAllocated = kFALSE;
1208  Element *abtp = work;
1209  if (this->fNcols > kWorkMax) {
1210  isAllocated = kTRUE;
1211  abtp = new Element[this->fNcols];
1212  }
1213 
1215 
1216  const Element *bp = b.GetMatrixArray();
1217  Element *cp = this->GetMatrixArray();
1218  if (typeid(Element) == typeid(Double_t))
1219  cblas_dsymm (CblasRowMajor,CblasLeft,CblasUpper,this->fNrows,this->fNcols,1.0,
1220  bp,b.GetNcols(),abtp,abt.GetNcols(),0.0,cp,this->fNcols);
1221  else if (typeid(Element) != typeid(Float_t))
1222  cblas_ssymm (CblasRowMajor,CblasLeft,CblasUpper,this->fNrows,this->fNcols,1.0,
1223  bp,b.GetNcols(),abtp,abt.GetNcols(),0.0,cp,this->fNcols);
1224  else
1225  Error("Similarity","type %s not implemented in BLAS library",typeid(Element));
1226 
1227  if (isAllocated)
1228  delete [] abtp;
1229 #else
1230  const Int_t ncolsa = this->GetNcols();
1231  const Int_t nb = b.GetNoElements();
1232  const Int_t nrowsb = b.GetNrows();
1233  const Int_t ncolsb = b.GetNcols();
1234 
1235  const Element * const bp = b.GetMatrixArray();
1236 
1237  Element work[kWorkMax];
1238  Bool_t isAllocated = kFALSE;
1239  Element *bap = work;
1240  if (nrowsb*ncolsa > kWorkMax) {
1241  isAllocated = kTRUE;
1242  bap = new Element[nrowsb*ncolsa];
1243  }
1244 
1245  AMultB(bp,nb,ncolsb,this->fElements,this->fNelems,this->fNcols,bap);
1246 
1247  const Int_t nba = nrowsb*ncolsa;
1248  const Int_t ncolsba = ncolsa;
1249  const Element * bi1p = bp;
1250  Element * cp = this->GetMatrixArray();
1251  Element * const cp0 = cp;
1252 
1253  Int_t ishift = 0;
1254  const Element *barp0 = bap;
1255  while (barp0 < bap+nba) {
1256  const Element *brp0 = bi1p;
1257  while (brp0 < bp+nb) {
1258  const Element *barp = barp0;
1259  const Element *brp = brp0;
1260  Element cij = 0;
1261  while (brp < brp0+ncolsb)
1262  cij += *barp++ * *brp++;
1263  *cp++ = cij;
1264  brp0 += ncolsb;
1265  }
1266  barp0 += ncolsba;
1267  bi1p += ncolsb;
1268  cp += ++ishift;
1269  }
1270 
1271  R__ASSERT(cp == cp0+this->fNelems+ishift && barp0 == bap+nba);
1272 
1273  cp = cp0;
1274  for (Int_t irow = 0; irow < this->fNrows; irow++) {
1275  const Int_t rowOff1 = irow*this->fNrows;
1276  for (Int_t icol = 0; icol < irow; icol++) {
1277  const Int_t rowOff2 = icol*this->fNrows;
1278  cp[rowOff1+icol] = cp[rowOff2+irow];
1279  }
1280  }
1281 
1282  if (isAllocated)
1283  delete [] bap;
1284 #endif
1285 
1286  return *this;
1287 }
1288 
1289 ////////////////////////////////////////////////////////////////////////////////
1290 /// Calculate scalar v * (*this) * v^T
1291 
1292 template<class Element>
1294 {
1295  if (gMatrixCheck) {
1296  R__ASSERT(this->IsValid());
1297  R__ASSERT(v.IsValid());
1298  if (this->fNcols != v.GetNrows() || this->fColLwb != v.GetLwb()) {
1299  Error("Similarity(const TVectorT &)","vector and matrix incompatible");
1300  return -1.;
1301  }
1302  }
1303 
1304  const Element *mp = this->GetMatrixArray(); // Matrix row ptr
1305  const Element *vp = v.GetMatrixArray(); // vector ptr
1306 
1307  Element sum1 = 0;
1308  const Element * const vp_first = vp;
1309  const Element * const vp_last = vp+v.GetNrows();
1310  while (vp < vp_last) {
1311  Element sum2 = 0;
1312  for (const Element *sp = vp_first; sp < vp_last; )
1313  sum2 += *mp++ * *sp++;
1314  sum1 += sum2 * *vp++;
1315  }
1316 
1317  R__ASSERT(mp == this->GetMatrixArray()+this->GetNoElements());
1318 
1319  return sum1;
1320 }
1321 
1322 ////////////////////////////////////////////////////////////////////////////////
1323 /// Calculate B^T * (*this) * B , final matrix will be (ncolsb x ncolsb)
1324 /// It is more efficient than applying the actual multiplication because this
1325 /// routine realizes that the final matrix is symmetric .
1326 
1327 template<class Element>
1329 {
1330  if (gMatrixCheck) {
1331  R__ASSERT(this->IsValid());
1332  R__ASSERT(b.IsValid());
1333  if (this->fNrows != b.GetNrows() || this->fRowLwb != b.GetRowLwb()) {
1334  Error("SimilarityT(const TMatrixT &)","matrices incompatible");
1335  return *this;
1336  }
1337  }
1338 
1339  const Int_t ncolsb = b.GetNcols();
1340  const Int_t ncolsa = this->GetNcols();
1341 
1342  Element work[kWorkMax];
1343  Bool_t isAllocated = kFALSE;
1344  Element *btap = work;
1345  if (ncolsb*ncolsa > kWorkMax) {
1346  isAllocated = kTRUE;
1347  btap = new Element[ncolsb*ncolsa];
1348  }
1349 
1350  TMatrixT<Element> bta; bta.Use(ncolsb,ncolsa,btap);
1351  bta.TMult(b,*this);
1352 
1353  if (ncolsb != this->fNcols)
1354  this->ResizeTo(ncolsb,ncolsb);
1355 
1356 #ifdef CBLAS
1357  const Element *bp = b.GetMatrixArray();
1358  Element *cp = this->GetMatrixArray();
1359  if (typeid(Element) == typeid(Double_t))
1360  cblas_dgemm (CblasRowMajor,CblasNoTrans,CblasNoTrans,this->fNrows,this->fNcols,bta.GetNcols(),
1361  1.0,btap,bta.GetNcols(),bp,b.GetNcols(),1.0,cp,this->fNcols);
1362  else if (typeid(Element) != typeid(Float_t))
1363  cblas_sgemm (CblasRowMajor,CblasNoTrans,CblasNoTrans,this->fNrows,this->fNcols,bta.GetNcols(),
1364  1.0,btap,bta.GetNcols(),bp,b.GetNcols(),1.0,cp,this->fNcols);
1365  else
1366  Error("similarityT","type %s not implemented in BLAS library",typeid(Element));
1367 #else
1368  const Int_t nbta = bta.GetNoElements();
1369  const Int_t nb = b.GetNoElements();
1370  const Int_t ncolsbta = bta.GetNcols();
1371  const Element * const bp = b.GetMatrixArray();
1372  Element * cp = this->GetMatrixArray();
1373  Element * const cp0 = cp;
1374 
1375  Int_t ishift = 0;
1376  const Element *btarp0 = btap; // Pointer to A[i,0];
1377  const Element *bcp0 = bp;
1378  while (btarp0 < btap+nbta) {
1379  for (const Element *bcp = bcp0; bcp < bp+ncolsb; ) { // Pointer to the j-th column of B, Start bcp = B[0,0]
1380  const Element *btarp = btarp0; // Pointer to the i-th row of A, reset to A[i,0]
1381  Element cij = 0;
1382  while (bcp < bp+nb) { // Scan the i-th row of A and
1383  cij += *btarp++ * *bcp; // the j-th col of B
1384  bcp += ncolsb;
1385  }
1386  *cp++ = cij;
1387  bcp -= nb-1; // Set bcp to the (j+1)-th col
1388  }
1389  btarp0 += ncolsbta; // Set ap to the (i+1)-th row
1390  bcp0++;
1391  cp += ++ishift;
1392  }
1393 
1394  R__ASSERT(cp == cp0+this->fNelems+ishift && btarp0 == btap+nbta);
1395 
1396  cp = cp0;
1397  for (Int_t irow = 0; irow < this->fNrows; irow++) {
1398  const Int_t rowOff1 = irow*this->fNrows;
1399  for (Int_t icol = 0; icol < irow; icol++) {
1400  const Int_t rowOff2 = icol*this->fNrows;
1401  cp[rowOff1+icol] = cp[rowOff2+irow];
1402  }
1403  }
1404 #endif
1405 
1406  if (isAllocated)
1407  delete [] btap;
1408 
1409  return *this;
1410 }
1411 
1412 ////////////////////////////////////////////////////////////////////////////////
1413 
1414 template<class Element>
1416 {
1417  if (gMatrixCheck && !AreCompatible(*this,source)) {
1418  Error("operator=","matrices not compatible");
1419  return *this;
1420  }
1421 
1422  if (this->GetMatrixArray() != source.GetMatrixArray()) {
1423  TObject::operator=(source);
1424  memcpy(this->GetMatrixArray(),source.fElements,this->fNelems*sizeof(Element));
1425  }
1426  return *this;
1427 }
1428 
1429 ////////////////////////////////////////////////////////////////////////////////
1430 
1431 template<class Element>
1433 {
1434  R__ASSERT(this->IsValid());
1435 
1436  if (lazy_constructor.fRowUpb != this->GetRowUpb() ||
1437  lazy_constructor.fRowLwb != this->GetRowLwb()) {
1438  Error("operator=(const TMatrixTSymLazy&)", "matrix is incompatible with "
1439  "the assigned Lazy matrix");
1440  return *this;
1441  }
1442 
1443  lazy_constructor.FillIn(*this);
1444  return *this;
1445 }
1446 
1447 ////////////////////////////////////////////////////////////////////////////////
1448 /// Assign val to every element of the matrix.
1449 
1450 template<class Element>
1452 {
1453  R__ASSERT(this->IsValid());
1454 
1455  Element *ep = fElements;
1456  const Element * const ep_last = ep+this->fNelems;
1457  while (ep < ep_last)
1458  *ep++ = val;
1459 
1460  return *this;
1461 }
1462 
1463 ////////////////////////////////////////////////////////////////////////////////
1464 /// Add val to every element of the matrix.
1465 
1466 template<class Element>
1468 {
1469  R__ASSERT(this->IsValid());
1470 
1471  Element *ep = fElements;
1472  const Element * const ep_last = ep+this->fNelems;
1473  while (ep < ep_last)
1474  *ep++ += val;
1475 
1476  return *this;
1477 }
1478 
1479 ////////////////////////////////////////////////////////////////////////////////
1480 /// Subtract val from every element of the matrix.
1481 
1482 template<class Element>
1484 {
1485  R__ASSERT(this->IsValid());
1486 
1487  Element *ep = fElements;
1488  const Element * const ep_last = ep+this->fNelems;
1489  while (ep < ep_last)
1490  *ep++ -= val;
1491 
1492  return *this;
1493 }
1494 
1495 ////////////////////////////////////////////////////////////////////////////////
1496 /// Multiply every element of the matrix with val.
1497 
1498 template<class Element>
1500 {
1501  R__ASSERT(this->IsValid());
1502 
1503  Element *ep = fElements;
1504  const Element * const ep_last = ep+this->fNelems;
1505  while (ep < ep_last)
1506  *ep++ *= val;
1507 
1508  return *this;
1509 }
1510 
1511 ////////////////////////////////////////////////////////////////////////////////
1512 /// Add the source matrix.
1513 
1514 template<class Element>
1516 {
1517  if (gMatrixCheck && !AreCompatible(*this,source)) {
1518  Error("operator+=","matrices not compatible");
1519  return *this;
1520  }
1521 
1522  const Element *sp = source.GetMatrixArray();
1523  Element *tp = this->GetMatrixArray();
1524  const Element * const tp_last = tp+this->fNelems;
1525  while (tp < tp_last)
1526  *tp++ += *sp++;
1527 
1528  return *this;
1529 }
1530 
1531 ////////////////////////////////////////////////////////////////////////////////
1532 /// Subtract the source matrix.
1533 
1534 template<class Element>
1536 {
1537  if (gMatrixCheck && !AreCompatible(*this,source)) {
1538  Error("operator-=","matrices not compatible");
1539  return *this;
1540  }
1541 
1542  const Element *sp = source.GetMatrixArray();
1543  Element *tp = this->GetMatrixArray();
1544  const Element * const tp_last = tp+this->fNelems;
1545  while (tp < tp_last)
1546  *tp++ -= *sp++;
1547 
1548  return *this;
1549 }
1550 
1551 ////////////////////////////////////////////////////////////////////////////////
1552 
1553 template<class Element>
1555 {
1556  R__ASSERT(this->IsValid());
1557 
1558  Element val = 0;
1559  Element *trp = this->GetMatrixArray(); // pointer to UR part and diagonal, traverse row-wise
1560  Element *tcp = trp; // pointer to LL part, traverse col-wise
1561  for (Int_t i = 0; i < this->fNrows; i++) {
1562  trp += i; // point to [i,i]
1563  tcp += i*this->fNcols; // point to [i,i]
1564  for (Int_t j = i; j < this->fNcols; j++) {
1565  action.Operation(val);
1566  if (j > i) *tcp = val;
1567  *trp++ = val;
1568  tcp += this->fNcols;
1569  }
1570  tcp -= this->fNelems-1; // point to [0,i]
1571  }
1572 
1573  return *this;
1574 }
1575 
1576 ////////////////////////////////////////////////////////////////////////////////
1577 /// Apply action to each element of the matrix. To action the location
1578 /// of the current element is passed.
1579 
1580 template<class Element>
1582 {
1583  R__ASSERT(this->IsValid());
1584 
1585  Element val = 0;
1586  Element *trp = this->GetMatrixArray(); // pointer to UR part and diagonal, traverse row-wise
1587  Element *tcp = trp; // pointer to LL part, traverse col-wise
1588  for (Int_t i = 0; i < this->fNrows; i++) {
1589  action.fI = i+this->fRowLwb;
1590  trp += i; // point to [i,i]
1591  tcp += i*this->fNcols; // point to [i,i]
1592  for (Int_t j = i; j < this->fNcols; j++) {
1593  action.fJ = j+this->fColLwb;
1594  action.Operation(val);
1595  if (j > i) *tcp = val;
1596  *trp++ = val;
1597  tcp += this->fNcols;
1598  }
1599  tcp -= this->fNelems-1; // point to [0,i]
1600  }
1601 
1602  return *this;
1603 }
1604 
1605 ////////////////////////////////////////////////////////////////////////////////
1606 /// randomize matrix element values but keep matrix symmetric
1607 
1608 template<class Element>
1610 {
1611  if (gMatrixCheck) {
1612  R__ASSERT(this->IsValid());
1613  if (this->fNrows != this->fNcols || this->fRowLwb != this->fColLwb) {
1614  Error("Randomize(Element,Element,Element &","matrix should be square");
1615  return *this;
1616  }
1617  }
1618 
1619  const Element scale = beta-alpha;
1620  const Element shift = alpha/scale;
1621 
1622  Element *ep = GetMatrixArray();
1623  for (Int_t i = 0; i < this->fNrows; i++) {
1624  const Int_t off = i*this->fNcols;
1625  for (Int_t j = 0; j <= i; j++) {
1626  ep[off+j] = scale*(Drand(seed)+shift);
1627  if (i != j) {
1628  ep[j*this->fNcols+i] = ep[off+j];
1629  }
1630  }
1631  }
1632 
1633  return *this;
1634 }
1635 
1636 ////////////////////////////////////////////////////////////////////////////////
1637 /// randomize matrix element values but keep matrix symmetric positive definite
1638 
1639 template<class Element>
1641 {
1642  if (gMatrixCheck) {
1643  R__ASSERT(this->IsValid());
1644  if (this->fNrows != this->fNcols || this->fRowLwb != this->fColLwb) {
1645  Error("RandomizeSym(Element,Element,Element &","matrix should be square");
1646  return *this;
1647  }
1648  }
1649 
1650  const Element scale = beta-alpha;
1651  const Element shift = alpha/scale;
1652 
1653  Element *ep = GetMatrixArray();
1654  Int_t i;
1655  for (i = 0; i < this->fNrows; i++) {
1656  const Int_t off = i*this->fNcols;
1657  for (Int_t j = 0; j <= i; j++)
1658  ep[off+j] = scale*(Drand(seed)+shift);
1659  }
1660 
1661  for (i = this->fNrows-1; i >= 0; i--) {
1662  const Int_t off1 = i*this->fNcols;
1663  for (Int_t j = i; j >= 0; j--) {
1664  const Int_t off2 = j*this->fNcols;
1665  ep[off1+j] *= ep[off2+j];
1666  for (Int_t k = j-1; k >= 0; k--) {
1667  ep[off1+j] += ep[off1+k]*ep[off2+k];
1668  }
1669  if (i != j)
1670  ep[off2+i] = ep[off1+j];
1671  }
1672  }
1673 
1674  return *this;
1675 }
1676 
1677 ////////////////////////////////////////////////////////////////////////////////
1678 /// Return a matrix containing the eigen-vectors ordered by descending eigen-values.
1679 /// For full functionality use TMatrixDSymEigen .
1680 
1681 template<class Element>
1683 {
1684  TMatrixDSym tmp = *this;
1685  TMatrixDSymEigen eigen(tmp);
1686  eigenValues.ResizeTo(this->fNrows);
1687  eigenValues = eigen.GetEigenValues();
1688  return eigen.GetEigenVectors();
1689 }
1690 
1691 ////////////////////////////////////////////////////////////////////////////////
1692 /// Check to see if two matrices are identical.
1693 
1694 template<class Element>
1696 {
1697  if (!AreCompatible(m1,m2)) return kFALSE;
1698  return (memcmp(m1.GetMatrixArray(),m2.GetMatrixArray(),
1699  m1.GetNoElements()*sizeof(Element)) == 0);
1700 }
1701 
1702 ////////////////////////////////////////////////////////////////////////////////
1703 
1704 template<class Element>
1706 {
1707  TMatrixTSym<Element> target(source1);
1708  target += source2;
1709  return target;
1710 }
1711 
1712 ////////////////////////////////////////////////////////////////////////////////
1713 
1714 template<class Element>
1716 {
1717  TMatrixTSym<Element> target(source1);
1718  target += val;
1719  return target;
1720 }
1721 
1722 ////////////////////////////////////////////////////////////////////////////////
1723 
1724 template<class Element>
1726 {
1727  return operator+(source1,val);
1728 }
1729 
1730 ////////////////////////////////////////////////////////////////////////////////
1731 
1732 template<class Element>
1734 {
1735  TMatrixTSym<Element> target(source1);
1736  target -= source2;
1737  return target;
1738 }
1739 
1740 ////////////////////////////////////////////////////////////////////////////////
1741 
1742 template<class Element>
1744 {
1745  TMatrixTSym<Element> target(source1);
1746  target -= val;
1747  return target;
1748 }
1749 
1750 ////////////////////////////////////////////////////////////////////////////////
1751 
1752 template<class Element>
1754 {
1755  return Element(-1.0)*operator-(source1,val);
1756 }
1757 
1758 ////////////////////////////////////////////////////////////////////////////////
1759 
1760 template<class Element>
1762 {
1763  TMatrixTSym<Element> target(source1);
1764  target *= val;
1765  return target;
1766 }
1767 
1768 ////////////////////////////////////////////////////////////////////////////////
1769 
1770 template<class Element>
1772 {
1773  return operator*(source1,val);
1774 }
1775 
1776 ////////////////////////////////////////////////////////////////////////////////
1777 /// Logical AND
1778 
1779 template<class Element>
1781 {
1782  TMatrixTSym<Element> target;
1783 
1784  if (gMatrixCheck && !AreCompatible(source1,source2)) {
1785  Error("operator&&(const TMatrixTSym&,const TMatrixTSym&)","matrices not compatible");
1786  return target;
1787  }
1788 
1789  target.ResizeTo(source1);
1790 
1791  const Element *sp1 = source1.GetMatrixArray();
1792  const Element *sp2 = source2.GetMatrixArray();
1793  Element *tp = target.GetMatrixArray();
1794  const Element * const tp_last = tp+target.GetNoElements();
1795  while (tp < tp_last)
1796  *tp++ = (*sp1++ != 0.0 && *sp2++ != 0.0);
1797 
1798  return target;
1799 }
1800 
1801 ////////////////////////////////////////////////////////////////////////////////
1802 /// Logical Or
1803 
1804 template<class Element>
1806 {
1807  TMatrixTSym<Element> target;
1808 
1809  if (gMatrixCheck && !AreCompatible(source1,source2)) {
1810  Error("operator||(const TMatrixTSym&,const TMatrixTSym&)","matrices not compatible");
1811  return target;
1812  }
1813 
1814  target.ResizeTo(source1);
1815 
1816  const Element *sp1 = source1.GetMatrixArray();
1817  const Element *sp2 = source2.GetMatrixArray();
1818  Element *tp = target.GetMatrixArray();
1819  const Element * const tp_last = tp+target.GetNoElements();
1820  while (tp < tp_last)
1821  *tp++ = (*sp1++ != 0.0 || *sp2++ != 0.0);
1822 
1823  return target;
1824 }
1825 
1826 ////////////////////////////////////////////////////////////////////////////////
1827 /// source1 > source2
1828 
1829 template<class Element>
1831 {
1832  TMatrixTSym<Element> target;
1833 
1834  if (gMatrixCheck && !AreCompatible(source1,source2)) {
1835  Error("operator>(const TMatrixTSym&,const TMatrixTSym&)","matrices not compatible");
1836  return target;
1837  }
1838 
1839  target.ResizeTo(source1);
1840 
1841  const Element *sp1 = source1.GetMatrixArray();
1842  const Element *sp2 = source2.GetMatrixArray();
1843  Element *tp = target.GetMatrixArray();
1844  const Element * const tp_last = tp+target.GetNoElements();
1845  while (tp < tp_last) {
1846  *tp++ = (*sp1) > (*sp2); sp1++; sp2++;
1847  }
1848 
1849  return target;
1850 }
1851 
1852 ////////////////////////////////////////////////////////////////////////////////
1853 /// source1 >= source2
1854 
1855 template<class Element>
1857 {
1858  TMatrixTSym<Element> target;
1859 
1860  if (gMatrixCheck && !AreCompatible(source1,source2)) {
1861  Error("operator>=(const TMatrixTSym&,const TMatrixTSym&)","matrices not compatible");
1862  return target;
1863  }
1864 
1865  target.ResizeTo(source1);
1866 
1867  const Element *sp1 = source1.GetMatrixArray();
1868  const Element *sp2 = source2.GetMatrixArray();
1869  Element *tp = target.GetMatrixArray();
1870  const Element * const tp_last = tp+target.GetNoElements();
1871  while (tp < tp_last) {
1872  *tp++ = (*sp1) >= (*sp2); sp1++; sp2++;
1873  }
1874 
1875  return target;
1876 }
1877 
1878 ////////////////////////////////////////////////////////////////////////////////
1879 /// source1 <= source2
1880 
1881 template<class Element>
1882 TMatrixTSym<Element> operator<=(const TMatrixTSym<Element> &source1,const TMatrixTSym<Element> &source2)
1883 {
1884  TMatrixTSym<Element> target;
1885 
1886  if (gMatrixCheck && !AreCompatible(source1,source2)) {
1887  Error("operator<=(const TMatrixTSym&,const TMatrixTSym&)","matrices not compatible");
1888  return target;
1889  }
1890 
1891  target.ResizeTo(source1);
1892 
1893  const Element *sp1 = source1.GetMatrixArray();
1894  const Element *sp2 = source2.GetMatrixArray();
1895  Element *tp = target.GetMatrixArray();
1896  const Element * const tp_last = tp+target.GetNoElements();
1897  while (tp < tp_last) {
1898  *tp++ = (*sp1) <= (*sp2); sp1++; sp2++;
1899  }
1900 
1901  return target;
1902 }
1903 
1904 ////////////////////////////////////////////////////////////////////////////////
1905 /// source1 < source2
1906 
1907 template<class Element>
1908 TMatrixTSym<Element> operator<(const TMatrixTSym<Element> &source1,const TMatrixTSym<Element> &source2)
1909 {
1910  TMatrixTSym<Element> target;
1911 
1912  if (gMatrixCheck && !AreCompatible(source1,source2)) {
1913  Error("operator<(const TMatrixTSym&,const TMatrixTSym&)","matrices not compatible");
1914  return target;
1915  }
1916 
1917  target.ResizeTo(source1);
1918 
1919  const Element *sp1 = source1.GetMatrixArray();
1920  const Element *sp2 = source2.GetMatrixArray();
1921  Element *tp = target.GetMatrixArray();
1922  const Element * const tp_last = tp+target.GetNoElements();
1923  while (tp < tp_last) {
1924  *tp++ = (*sp1) < (*sp2); sp1++; sp2++;
1925  }
1926 
1927  return target;
1928 }
1929 
1930 ////////////////////////////////////////////////////////////////////////////////
1931 /// Modify addition: target += scalar * source.
1932 
1933 template<class Element>
1935 {
1936  if (gMatrixCheck && !AreCompatible(target,source)) {
1937  ::Error("Add","matrices not compatible");
1938  return target;
1939  }
1940 
1941  const Int_t nrows = target.GetNrows();
1942  const Int_t ncols = target.GetNcols();
1943  const Int_t nelems = target.GetNoElements();
1944  const Element *sp = source.GetMatrixArray();
1945  Element *trp = target.GetMatrixArray(); // pointer to UR part and diagonal, traverse row-wise
1946  Element *tcp = target.GetMatrixArray(); // pointer to LL part, traverse col-wise
1947  for (Int_t i = 0; i < nrows; i++) {
1948  sp += i;
1949  trp += i; // point to [i,i]
1950  tcp += i*ncols; // point to [i,i]
1951  for (Int_t j = i; j < ncols; j++) {
1952  const Element tmp = scalar * *sp++;
1953  if (j > i) *tcp += tmp;
1954  *trp++ += tmp;
1955  tcp += ncols;
1956  }
1957  tcp -= nelems-1; // point to [0,i]
1958  }
1959 
1960  return target;
1961 }
1962 
1963 ////////////////////////////////////////////////////////////////////////////////
1964 /// Multiply target by the source, element-by-element.
1965 
1966 template<class Element>
1968 {
1969  if (gMatrixCheck && !AreCompatible(target,source)) {
1970  ::Error("ElementMult","matrices not compatible");
1971  return target;
1972  }
1973 
1974  const Int_t nrows = target.GetNrows();
1975  const Int_t ncols = target.GetNcols();
1976  const Int_t nelems = target.GetNoElements();
1977  const Element *sp = source.GetMatrixArray();
1978  Element *trp = target.GetMatrixArray(); // pointer to UR part and diagonal, traverse row-wise
1979  Element *tcp = target.GetMatrixArray(); // pointer to LL part, traverse col-wise
1980  for (Int_t i = 0; i < nrows; i++) {
1981  sp += i;
1982  trp += i; // point to [i,i]
1983  tcp += i*ncols; // point to [i,i]
1984  for (Int_t j = i; j < ncols; j++) {
1985  if (j > i) *tcp *= *sp;
1986  *trp++ *= *sp++;
1987  tcp += ncols;
1988  }
1989  tcp -= nelems-1; // point to [0,i]
1990  }
1991 
1992  return target;
1993 }
1994 
1995 ////////////////////////////////////////////////////////////////////////////////
1996 /// Multiply target by the source, element-by-element.
1997 
1998 template<class Element>
2000 {
2001  if (gMatrixCheck && !AreCompatible(target,source)) {
2002  ::Error("ElementDiv","matrices not compatible");
2003  return target;
2004  }
2005 
2006  const Int_t nrows = target.GetNrows();
2007  const Int_t ncols = target.GetNcols();
2008  const Int_t nelems = target.GetNoElements();
2009  const Element *sp = source.GetMatrixArray();
2010  Element *trp = target.GetMatrixArray(); // pointer to UR part and diagonal, traverse row-wise
2011  Element *tcp = target.GetMatrixArray(); // pointer to LL part, traverse col-wise
2012  for (Int_t i = 0; i < nrows; i++) {
2013  sp += i;
2014  trp += i; // point to [i,i]
2015  tcp += i*ncols; // point to [i,i]
2016  for (Int_t j = i; j < ncols; j++) {
2017  if (*sp != 0.0) {
2018  if (j > i) *tcp /= *sp;
2019  *trp++ /= *sp++;
2020  } else {
2021  const Int_t irow = (sp-source.GetMatrixArray())/source.GetNcols();
2022  const Int_t icol = (sp-source.GetMatrixArray())%source.GetNcols();
2023  Error("ElementDiv","source (%d,%d) is zero",irow,icol);
2024  trp++;
2025  }
2026  tcp += ncols;
2027  }
2028  tcp -= nelems-1; // point to [0,i]
2029  }
2030 
2031  return target;
2032 }
2033 
2034 ////////////////////////////////////////////////////////////////////////////////
2035 /// Stream an object of class TMatrixTSym.
2036 
2037 template<class Element>
2039 {
2040  if (R__b.IsReading()) {
2041  UInt_t R__s, R__c;
2042  Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
2043  Clear();
2044  R__b.ReadClassBuffer(TMatrixTBase<Element>::Class(),this,R__v,R__s,R__c);
2045  fElements = new Element[this->fNelems];
2046  Int_t i;
2047  for (i = 0; i < this->fNrows; i++) {
2048  R__b.ReadFastArray(fElements+i*this->fNcols+i,this->fNcols-i);
2049  }
2050  // copy to Lower left triangle
2051  for (i = 0; i < this->fNrows; i++) {
2052  for (Int_t j = 0; j < i; j++) {
2053  fElements[i*this->fNcols+j] = fElements[j*this->fNrows+i];
2054  }
2055  }
2056  if (this->fNelems <= this->kSizeMax) {
2057  memcpy(fDataStack,fElements,this->fNelems*sizeof(Element));
2058  delete [] fElements;
2059  fElements = fDataStack;
2060  }
2061  } else {
2063  // Only write the Upper right triangle
2064  for (Int_t i = 0; i < this->fNrows; i++) {
2065  R__b.WriteFastArray(fElements+i*this->fNcols+i,this->fNcols-i);
2066  }
2067  }
2068 }
2069 
2070 #include "TMatrixFSymfwd.h"
2071 
2072 template class TMatrixTSym<Float_t>;
2073 
2074 template Bool_t operator== <Float_t>(const TMatrixFSym &source1,const TMatrixFSym &source2);
2075 template TMatrixFSym operator+ <Float_t>(const TMatrixFSym &source1,const TMatrixFSym &source2);
2076 template TMatrixFSym operator+ <Float_t>(const TMatrixFSym &source1, Float_t val);
2077 template TMatrixFSym operator+ <Float_t>( Float_t val ,const TMatrixFSym &source2);
2078 template TMatrixFSym operator- <Float_t>(const TMatrixFSym &source1,const TMatrixFSym &source2);
2079 template TMatrixFSym operator- <Float_t>(const TMatrixFSym &source1, Float_t val);
2080 template TMatrixFSym operator- <Float_t>( Float_t val ,const TMatrixFSym &source2);
2081 template TMatrixFSym operator* <Float_t>(const TMatrixFSym &source, Float_t val );
2082 template TMatrixFSym operator* <Float_t>( Float_t val, const TMatrixFSym &source );
2083 template TMatrixFSym operator&& <Float_t>(const TMatrixFSym &source1,const TMatrixFSym &source2);
2084 template TMatrixFSym operator|| <Float_t>(const TMatrixFSym &source1,const TMatrixFSym &source2);
2085 template TMatrixFSym operator> <Float_t>(const TMatrixFSym &source1,const TMatrixFSym &source2);
2086 template TMatrixFSym operator>= <Float_t>(const TMatrixFSym &source1,const TMatrixFSym &source2);
2087 template TMatrixFSym operator<= <Float_t>(const TMatrixFSym &source1,const TMatrixFSym &source2);
2088 template TMatrixFSym operator< <Float_t>(const TMatrixFSym &source1,const TMatrixFSym &source2);
2089 
2090 template TMatrixFSym &Add <Float_t>(TMatrixFSym &target, Float_t scalar,const TMatrixFSym &source);
2091 template TMatrixFSym &ElementMult<Float_t>(TMatrixFSym &target,const TMatrixFSym &source);
2092 template TMatrixFSym &ElementDiv <Float_t>(TMatrixFSym &target,const TMatrixFSym &source);
2093 
2094 #include "TMatrixDSymfwd.h"
2095 
2096 template class TMatrixTSym<Double_t>;
2097 
2098 template Bool_t operator== <Double_t>(const TMatrixDSym &source1,const TMatrixDSym &source2);
2099 template TMatrixDSym operator+ <Double_t>(const TMatrixDSym &source1,const TMatrixDSym &source2);
2100 template TMatrixDSym operator+ <Double_t>(const TMatrixDSym &source1, Double_t val);
2101 template TMatrixDSym operator+ <Double_t>( Double_t val ,const TMatrixDSym &source2);
2102 template TMatrixDSym operator- <Double_t>(const TMatrixDSym &source1,const TMatrixDSym &source2);
2103 template TMatrixDSym operator- <Double_t>(const TMatrixDSym &source1, Double_t val);
2104 template TMatrixDSym operator- <Double_t>( Double_t val ,const TMatrixDSym &source2);
2105 template TMatrixDSym operator* <Double_t>(const TMatrixDSym &source, Double_t val );
2106 template TMatrixDSym operator* <Double_t>( Double_t val, const TMatrixDSym &source );
2107 template TMatrixDSym operator&& <Double_t>(const TMatrixDSym &source1,const TMatrixDSym &source2);
2108 template TMatrixDSym operator|| <Double_t>(const TMatrixDSym &source1,const TMatrixDSym &source2);
2109 template TMatrixDSym operator> <Double_t>(const TMatrixDSym &source1,const TMatrixDSym &source2);
2110 template TMatrixDSym operator>= <Double_t>(const TMatrixDSym &source1,const TMatrixDSym &source2);
2111 template TMatrixDSym operator<= <Double_t>(const TMatrixDSym &source1,const TMatrixDSym &source2);
2112 template TMatrixDSym operator< <Double_t>(const TMatrixDSym &source1,const TMatrixDSym &source2);
2113 
2114 template TMatrixDSym &Add <Double_t>(TMatrixDSym &target, Double_t scalar,const TMatrixDSym &source);
2115 template TMatrixDSym &ElementMult<Double_t>(TMatrixDSym &target,const TMatrixDSym &source);
2116 template TMatrixDSym &ElementDiv <Double_t>(TMatrixDSym &target,const TMatrixDSym &source);
TMatrixTSym< Element > & SimilarityT(const TMatrixT< Element > &n)
Calculate B^T * (*this) * B , final matrix will be (ncolsb x ncolsb) It is more efficient than applyi...
virtual const Element * GetMatrixArray() const =0
virtual Bool_t IsSymmetric() const
Check whether matrix is symmetric.
virtual TMatrixTBase< Element > & Shift(Int_t row_shift, Int_t col_shift)
Shift the row index by adding row_shift and the column index by adding col_shift, respectively...
Bool_t IsReading() const
Definition: TBuffer.h:85
virtual TMatrixTBase< Element > & Shift(Int_t row_shift, Int_t col_shift)
Shift the row index by adding row_shift and the column index by adding col_shift, respectively...
template TMatrixFSym operator<=< Float_t >(const TMatrixFSym &source1, const TMatrixFSym &source2)
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
virtual Int_t WriteClassBuffer(const TClass *cl, void *pointer)=0
TMatrixTSym< Element > & ElementMult(TMatrixTSym< Element > &target, const TMatrixTSym< Element > &source)
Multiply target by the source, element-by-element.
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...
TMatrixTSym< Element > & Use(Int_t row_lwb, Int_t row_upb, Element *data)
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 .
short Version_t
Definition: RtypesCore.h:63
TMatrixTSym< Element > & Add(TMatrixTSym< Element > &target, Element scalar, const TMatrixTSym< Element > &source)
Modify addition: target += scalar * 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.
float Float_t
Definition: RtypesCore.h:55
TMatrixTSym< Element > & operator=(const TMatrixTSym< Element > &source)
const char Option_t
Definition: RtypesCore.h:64
virtual const Element * GetMatrixArray() const
Definition: TMatrixT.h:222
Int_t GetNcols() const
Definition: TMatrixTBase.h:127
TMatrixTSym< Element > & ElementDiv(TMatrixTSym< Element > &target, const TMatrixTSym< Element > &source)
Multiply target by the source, element-by-element.
TMatrixTSym< Element > operator>=(const TMatrixTSym< Element > &source1, const TMatrixTSym< Element > &source2)
source1 >= source2
TVectorT.
Definition: TMatrixTBase.h:79
void ToUpper()
Change string to upper case.
Definition: TString.cxx:1138
Buffer base class used for serializing objects.
Definition: TBuffer.h:42
void TMult(const TMatrixT< Element > &a)
Create a matrix C such that C = A&#39; * A.
#define R__ASSERT(e)
Definition: TError.h:96
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.
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
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
bool Bool_t
Definition: RtypesCore.h:61
virtual void Det(Double_t &d1, Double_t &d2)
Calculate determinant det = d1*TMath::Power(2.,d2)
Definition: TDecompLU.cxx:509
virtual const Element * GetMatrixArray() const
Definition: TMatrixTSym.h:189
TMatrixTSym< Element > & operator*=(Element val)
Multiply every element of the matrix with val.
virtual TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t nr_nonzeros=-1)=0
Int_t GetRowUpb() const
Definition: TMatrixTLazy.h:111
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:725
TMatrixTSym< Element > & Rank1Update(const TVectorT< Element > &v, Element alpha=1.0)
Perform a rank 1 operation on the matrix: A += alpha * v * v^T.
double beta(double x, double y)
Calculates the beta function.
TMatrixDSymEigen.
TMatrixT.
Definition: TMatrixDfwd.h:22
Int_t GetNoElements() const
Definition: TMatrixTBase.h:128
LU Decomposition class.
Definition: TDecompLU.h:23
template TMatrixDSym & ElementDiv< Double_t >(TMatrixDSym &target, const TMatrixDSym &source)
#define templateClassImp(name)
Definition: Rtypes.h:405
TMatrixTSym< Element > & operator+=(Element val)
Add val to every element of the matrix.
Int_t GetColLwb() const
Definition: TMatrixTBase.h:125
virtual void Operation(Element &element) const =0
void AMultB(int n, int m, int k, const double *A, const double *B, double *C)
Definition: cblas.cxx:11
template TMatrixDSym & ElementMult< Double_t >(TMatrixDSym &target, const TMatrixDSym &source)
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
virtual TMatrixTBase< Element > & SetMatrixArray(const Element *data, Option_t *option="")
Copy array data to matrix .
virtual void FillIn(TMatrixTSym< Element > &m) const =0
void Error(const char *location, const char *msgfmt,...)
void Minus(const TMatrixTSym< Element > &a, const TMatrixTSym< Element > &b)
Symmetric matrix summation. Create a matrix C such that C = A + B.
TMatrixTSym< Element > operator>(const TMatrixTSym< Element > &source1, const TMatrixTSym< Element > &source2)
source1 > source2
virtual const Int_t * GetColIndexArray() const
Definition: TMatrixTSym.h:86
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 .
TMatrixTSym< Element > operator &&(const TMatrixTSym< Element > &source1, const TMatrixTSym< Element > &source2)
Logical AND.
virtual Double_t Determinant() const
static constexpr double m2
TMatrixTBase< Element > & Apply(const TElementActionT< Element > &action)
Apply action to each matrix element.
auto * a
Definition: textangle.C:12
template TMatrixFSym & Add< Float_t >(TMatrixFSym &target, Float_t scalar, const TMatrixFSym &source)
template TMatrixDSym operator<=< Double_t >(const TMatrixDSym &source1, const TMatrixDSym &source2)
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)
template TMatrixFSym & ElementDiv< Float_t >(TMatrixFSym &target, const TMatrixFSym &source)
Bool_t operator==(const TMatrixTSym< Element > &m1, const TMatrixTSym< Element > &m2)
Check to see if two matrices are identical.
R__EXTERN Int_t gMatrixCheck
Definition: TMatrixTBase.h:81
virtual TMatrixTBase< Element > & Randomize(Element alpha, Element beta, Double_t &seed)
randomize matrix element values but keep matrix symmetric
virtual void ReadFastArray(Bool_t *b, Int_t n)=0
virtual void WriteFastArray(const Bool_t *b, Int_t n)=0
const TMatrixD & GetEigenVectors() const
TMatrixTBase.
REAL epsilon
Definition: triangle.c:617
TMatrixTSym< Element > & GetSub(Int_t row_lwb, Int_t row_upb, TMatrixTSym< Element > &target, Option_t *option="S") const
Get submatrix [row_lwb..row_upb][row_lwb..row_upb]; The indexing range of the returned matrix depends...
Int_t GetNrows() const
Definition: TMatrixTBase.h:124
const Bool_t kFALSE
Definition: RtypesCore.h:90
TMatrixTSym< Element > & Transpose(const TMatrixTSym< Element > &source)
Transpose a matrix.
virtual Int_t ReadClassBuffer(const TClass *cl, void *pointer, const TClass *onfile_class=0)=0
const TVectorD & GetEigenValues() const
void Plus(const TMatrixTSym< Element > &a, const TMatrixTSym< Element > &b)
Symmetric matrix summation. Create a matrix C such that C = A + B.
double Double_t
Definition: RtypesCore.h:57
TMatrixTSym< Element > & operator-=(Element val)
Subtract val from every element of the matrix.
virtual TMatrixTSym< Element > & RandomizePD(Element alpha, Element beta, Double_t &seed)
randomize matrix element values but keep matrix symmetric positive definite
TMatrixTSym< Element > & InvertFast(Double_t *det=0)
Invert the matrix and calculate its determinant.
TMatrixTSym< Element > & SetSub(Int_t row_lwb, const TMatrixTBase< Element > &source)
Insert matrix source starting at [row_lwb][row_lwb], thereby overwriting the part [row_lwb...
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:619
TMatrixTSym< Element > & Similarity(const TMatrixT< Element > &n)
Calculate B * (*this) * B^T , final matrix will be (nrowsb x nrowsb) This is a similarity transform w...
void Delete_m(Int_t size, Element *&)
delete data pointer m, if it was assigned on the heap
virtual void Operation(Element &element) const =0
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
TMatrixTSym< Element > & Invert(Double_t *det=0)
Invert the matrix and calculate its determinant Notice that the LU decomposition is used instead of B...
TMatrixTSym< Element > operator||(const TMatrixTSym< Element > &source1, const TMatrixTSym< Element > &source2)
Logical Or.
char Char_t
Definition: RtypesCore.h:31
Int_t GetRowLwb() const
Definition: TMatrixTLazy.h:110
TMatrixTSym< Element > operator+(const TMatrixTSym< Element > &source1, const TMatrixTSym< Element > &source2)
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:212
template TMatrixDSym & Add< Double_t >(TMatrixDSym &target, Double_t scalar, const TMatrixDSym &source)
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
const TMatrixT< Element > EigenVectors(TVectorT< Element > &eigenValues) const
Return a matrix containing the eigen-vectors ordered by descending eigen-values.
Element * New_m(Int_t size)
return data pointer .
TMatrixTSym< Element > operator*(const TMatrixTSym< Element > &source1, Element val)
template TMatrixFSym & ElementMult< Float_t >(TMatrixFSym &target, const TMatrixFSym &source)
const Bool_t kTRUE
Definition: RtypesCore.h:89
Element * fElements
data container
Definition: TMatrixTSym.h:39
Double_t Drand(Double_t &ix)
Random number generator [0....1] with seed ix.
TMatrixTSym< Element > operator-(const TMatrixTSym< Element > &source1, const TMatrixTSym< Element > &source2)
virtual const Int_t * GetRowIndexArray() const
Definition: TMatrixTSym.h:84
Bool_t IsValid() const
Definition: TMatrixTBase.h:147
virtual const Int_t * GetColIndexArray() const =0
virtual Version_t ReadVersion(UInt_t *start=0, UInt_t *bcnt=0, const TClass *cl=0)=0
virtual const Int_t * GetRowIndexArray() const =0