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