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