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