Logo ROOT   6.16/01
Reference Guide
TMatrixTSym.cxx
Go to the documentation of this file.
1// @(#)root/matrix:$Id$
2// Authors: Fons Rademakers, Eddy Offermann Nov 2003
3
4/*************************************************************************
5 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6 * All rights reserved. *
7 * *
8 * For the licensing terms see $ROOTSYS/LICENSE. *
9 * For the list of contributors see $ROOTSYS/README/CREDITS. *
10 *************************************************************************/
11
12/** \class TMatrixTSym
13 \ingroup Matrix
14
15 TMatrixTSym
16
17 Template class of a symmetric matrix in the linear algebra package
18
19 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"
29#include "TDecompLU.h"
30#include "TMatrixDSymEigen.h"
31#include "TClass.h"
32#include "TMath.h"
33
35
36
37////////////////////////////////////////////////////////////////////////////////
38
39template<class Element>
42 Allocate(no_rows,no_rows,0,0,1);
44
45////////////////////////////////////////////////////////////////////////////////
46
47template<class Element>
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
63template<class Element>
64TMatrixTSym<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
76template<class Element>
77TMatrixTSym<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
89template<class Element>
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
102template<class Element>
104{
105 R__ASSERT(prototype.IsValid());
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;
125 case kInverted:
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;
143 default:
144 Error("TMatrixTSym(EMatrixCreatorOp1,const TMatrixTSym)",
145 "operation %d not yet implemented", op);
146 }
147}
148
149////////////////////////////////////////////////////////////////////////////////
150
151template<class Element>
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;
162 default:
163 Error("TMatrixTSym(EMatrixCreatorOp1,const TMatrixT)",
164 "operation %d not yet implemented", op);
167
168////////////////////////////////////////////////////////////////////////////////
170template<class Element>
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
198template<class Element>
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
213template<class Element>
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
227template<class Element>
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
245template<class Element>
246Int_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
272template<class Element>
273void 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;
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
310template<class Element>
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
344template<class Element>
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
379template<class Element>
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
427template<class Element>
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
473template<class Element>
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
501template<class Element>
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
567template<class Element>
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
634template<class Element>
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
684template<class Element>
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
741template<class Element>
743{
745 if (!this->IsSymmetric()) {
746 Error("SetMatrixArray","Matrix is not symmetric after Set");
747 }
748
749 return *this;
750}
751
752////////////////////////////////////////////////////////////////////////////////
753
754template<class Element>
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
769template<class Element>
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
839template<class Element>
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
933template<class Element>
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
945template<class Element>
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
959template<class Element>
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
977template<class Element>
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
1039template<class Element>
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
1061template<class Element>
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
1097template<class Element>
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
1188template<class Element>
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
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
1290template<class Element>
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
1325template<class Element>
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
1412template<class Element>
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
1429template<class Element>
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
1448template<class Element>
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
1464template<class Element>
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
1480template<class Element>
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
1496template<class Element>
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
1512template<class Element>
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
1532template<class Element>
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
1551template<class Element>
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
1578template<class Element>
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
1606template<class Element>
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
1637template<class Element>
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
1679template<class Element>
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
1692template<class Element>
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
1702template<class Element>
1704{
1705 TMatrixTSym<Element> target(source1);
1706 target += source2;
1707 return target;
1708}
1709
1710////////////////////////////////////////////////////////////////////////////////
1711
1712template<class Element>
1714{
1715 TMatrixTSym<Element> target(source1);
1716 target += val;
1717 return target;
1718}
1719
1720////////////////////////////////////////////////////////////////////////////////
1721
1722template<class Element>
1724{
1725 return operator+(source1,val);
1726}
1727
1728////////////////////////////////////////////////////////////////////////////////
1729
1730template<class Element>
1732{
1733 TMatrixTSym<Element> target(source1);
1734 target -= source2;
1735 return target;
1736}
1737
1738////////////////////////////////////////////////////////////////////////////////
1739
1740template<class Element>
1742{
1743 TMatrixTSym<Element> target(source1);
1744 target -= val;
1745 return target;
1746}
1747
1748////////////////////////////////////////////////////////////////////////////////
1749
1750template<class Element>
1752{
1753 return Element(-1.0)*operator-(source1,val);
1754}
1755
1756////////////////////////////////////////////////////////////////////////////////
1757
1758template<class Element>
1760{
1761 TMatrixTSym<Element> target(source1);
1762 target *= val;
1763 return target;
1764}
1765
1766////////////////////////////////////////////////////////////////////////////////
1767
1768template<class Element>
1770{
1771 return operator*(source1,val);
1772}
1773
1774////////////////////////////////////////////////////////////////////////////////
1775/// Logical AND
1776
1777template<class Element>
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
1802template<class Element>
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
1827template<class Element>
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
1853template<class Element>
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
1879template<class Element>
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
1905template<class Element>
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
1931template<class Element>
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
1964template<class Element>
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
1996template<class Element>
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
2035template<class Element>
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 {
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
2070template class TMatrixTSym<Float_t>;
2071
2072template Bool_t operator== <Float_t>(const TMatrixFSym &source1,const TMatrixFSym &source2);
2073template TMatrixFSym operator+ <Float_t>(const TMatrixFSym &source1,const TMatrixFSym &source2);
2074template TMatrixFSym operator+ <Float_t>(const TMatrixFSym &source1, Float_t val);
2075template TMatrixFSym operator+ <Float_t>( Float_t val ,const TMatrixFSym &source2);
2076template TMatrixFSym operator- <Float_t>(const TMatrixFSym &source1,const TMatrixFSym &source2);
2077template TMatrixFSym operator- <Float_t>(const TMatrixFSym &source1, Float_t val);
2078template TMatrixFSym operator- <Float_t>( Float_t val ,const TMatrixFSym &source2);
2079template TMatrixFSym operator* <Float_t>(const TMatrixFSym &source, Float_t val );
2080template TMatrixFSym operator* <Float_t>( Float_t val, const TMatrixFSym &source );
2081template TMatrixFSym operator&& <Float_t>(const TMatrixFSym &source1,const TMatrixFSym &source2);
2082template TMatrixFSym operator|| <Float_t>(const TMatrixFSym &source1,const TMatrixFSym &source2);
2083template TMatrixFSym operator> <Float_t>(const TMatrixFSym &source1,const TMatrixFSym &source2);
2084template TMatrixFSym operator>= <Float_t>(const TMatrixFSym &source1,const TMatrixFSym &source2);
2085template TMatrixFSym operator<= <Float_t>(const TMatrixFSym &source1,const TMatrixFSym &source2);
2086template TMatrixFSym operator< <Float_t>(const TMatrixFSym &source1,const TMatrixFSym &source2);
2087
2088template TMatrixFSym &Add <Float_t>(TMatrixFSym &target, Float_t scalar,const TMatrixFSym &source);
2091
2092#include "TMatrixDSymfwd.h"
2093
2094template class TMatrixTSym<Double_t>;
2095
2096template Bool_t operator== <Double_t>(const TMatrixDSym &source1,const TMatrixDSym &source2);
2097template TMatrixDSym operator+ <Double_t>(const TMatrixDSym &source1,const TMatrixDSym &source2);
2098template TMatrixDSym operator+ <Double_t>(const TMatrixDSym &source1, Double_t val);
2099template TMatrixDSym operator+ <Double_t>( Double_t val ,const TMatrixDSym &source2);
2100template TMatrixDSym operator- <Double_t>(const TMatrixDSym &source1,const TMatrixDSym &source2);
2101template TMatrixDSym operator- <Double_t>(const TMatrixDSym &source1, Double_t val);
2102template TMatrixDSym operator- <Double_t>( Double_t val ,const TMatrixDSym &source2);
2103template TMatrixDSym operator* <Double_t>(const TMatrixDSym &source, Double_t val );
2104template TMatrixDSym operator* <Double_t>( Double_t val, const TMatrixDSym &source );
2105template TMatrixDSym operator&& <Double_t>(const TMatrixDSym &source1,const TMatrixDSym &source2);
2106template TMatrixDSym operator|| <Double_t>(const TMatrixDSym &source1,const TMatrixDSym &source2);
2107template TMatrixDSym operator> <Double_t>(const TMatrixDSym &source1,const TMatrixDSym &source2);
2108template TMatrixDSym operator>= <Double_t>(const TMatrixDSym &source1,const TMatrixDSym &source2);
2109template TMatrixDSym operator<= <Double_t>(const TMatrixDSym &source1,const TMatrixDSym &source2);
2110template TMatrixDSym operator< <Double_t>(const TMatrixDSym &source1,const TMatrixDSym &source2);
2111
2112template TMatrixDSym &Add <Double_t>(TMatrixDSym &target, Double_t scalar,const TMatrixDSym &source);
SVector< double, 2 > v
Definition: Dict.h:5
#define b(i)
Definition: RSha256.hxx:100
static double p1(double t, double a, double b)
static double p2(double t, double a, double b, double c)
static Int_t init()
int Int_t
Definition: RtypesCore.h:41
short Version_t
Definition: RtypesCore.h:61
char Char_t
Definition: RtypesCore.h:29
unsigned int UInt_t
Definition: RtypesCore.h:42
const Bool_t kFALSE
Definition: RtypesCore.h:88
bool Bool_t
Definition: RtypesCore.h:59
double Double_t
Definition: RtypesCore.h:55
float Float_t
Definition: RtypesCore.h:53
const Bool_t kTRUE
Definition: RtypesCore.h:87
const char Option_t
Definition: RtypesCore.h:62
#define templateClassImp(name)
Definition: Rtypes.h:407
@ kPlus
Definition: TAttMarker.h:46
#define R__ASSERT(e)
Definition: TError.h:96
void Error(const char *location, const char *msgfmt,...)
R__EXTERN Int_t gMatrixCheck
Definition: TMatrixTBase.h:81
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 .
template TMatrixDSym & ElementMult< Double_t >(TMatrixDSym &target, const TMatrixDSym &source)
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, Element val)
TMatrixTSym< Element > operator<(const TMatrixTSym< Element > &source1, const TMatrixTSym< Element > &source2)
source1 < source2
template TMatrixFSym & ElementMult< Float_t >(TMatrixFSym &target, const TMatrixFSym &source)
TMatrixTSym< Element > operator+(const TMatrixTSym< Element > &source1, const TMatrixTSym< Element > &source2)
template TMatrixFSym & Add< Float_t >(TMatrixFSym &target, Float_t scalar, const TMatrixFSym &source)
template TMatrixFSym & ElementDiv< Float_t >(TMatrixFSym &target, const TMatrixFSym &source)
template TMatrixDSym & Add< Double_t >(TMatrixDSym &target, Double_t scalar, const TMatrixDSym &source)
TMatrixTSym< Element > operator-(const TMatrixTSym< Element > &source1, const TMatrixTSym< Element > &source2)
template TMatrixFSym operator<=< Float_t >(const TMatrixFSym &source1, const TMatrixFSym &source2)
TMatrixTSym< Element > & Add(TMatrixTSym< Element > &target, Element scalar, const TMatrixTSym< Element > &source)
Modify addition: target += scalar * source.
TMatrixTSym< Element > & ElementMult(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
TMatrixTSym< Element > operator>=(const TMatrixTSym< Element > &source1, const TMatrixTSym< Element > &source2)
source1 >= source2
TMatrixTSym< Element > operator<=(const TMatrixTSym< Element > &source1, const TMatrixTSym< Element > &source2)
source1 <= source2
template TMatrixFSym operator<<Float_t >(const TMatrixFSym &source1, const TMatrixFSym &source2)
template TMatrixDSym & ElementDiv< Double_t >(TMatrixDSym &target, const TMatrixDSym &source)
TMatrixTSym< Element > operator||(const TMatrixTSym< Element > &source1, const TMatrixTSym< Element > &source2)
Logical Or.
TMatrixTSym< Element > operator&&(const TMatrixTSym< Element > &source1, const TMatrixTSym< Element > &source2)
Logical AND.
template TMatrixDSym operator<<Double_t >(const TMatrixDSym &source1, const TMatrixDSym &source2)
Bool_t operator==(const TMatrixTSym< Element > &m1, const TMatrixTSym< Element > &m2)
Check to see if two matrices are identical.
template TMatrixDSym operator<=< Double_t >(const TMatrixDSym &source1, const TMatrixDSym &source2)
Double_t Drand(Double_t &ix)
Random number generator [0....1] with seed ix.
Buffer base class used for serializing objects.
Definition: TBuffer.h:40
virtual Int_t ReadClassBuffer(const TClass *cl, void *pointer, const TClass *onfile_class=0)=0
virtual Version_t ReadVersion(UInt_t *start=0, UInt_t *bcnt=0, const TClass *cl=0)=0
virtual void ReadFastArray(Bool_t *b, Int_t n)=0
Bool_t IsReading() const
Definition: TBuffer.h:83
virtual Int_t WriteClassBuffer(const TClass *cl, void *pointer)=0
virtual void WriteFastArray(const Bool_t *b, Int_t n)=0
LU Decomposition class.
Definition: TDecompLU.h:24
virtual void Det(Double_t &d1, Double_t &d2)
Calculate determinant det = d1*TMath::Power(2.,d2)
Definition: TDecompLU.cxx:509
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
virtual void Operation(Element &element) const =0
virtual void Operation(Element &element) const =0
TMatrixDSymEigen.
const TVectorD & GetEigenValues() const
const TMatrixD & GetEigenVectors() const
Linear Algebra Package.
Definition: TMatrixTBase.h:85
virtual const Element * GetMatrixArray() const =0
Int_t GetNrows() const
Definition: TMatrixTBase.h:124
virtual const Int_t * GetRowIndexArray() const =0
virtual const Int_t * GetColIndexArray() const =0
virtual TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t nr_nonzeros=-1)=0
Int_t GetRowLwb() const
Definition: TMatrixTBase.h:122
Int_t GetColLwb() const
Definition: TMatrixTBase.h:125
Int_t GetNoElements() const
Definition: TMatrixTBase.h:128
Bool_t IsValid() const
Definition: TMatrixTBase.h:147
Int_t GetNcols() const
Definition: TMatrixTBase.h:127
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 Bool_t IsSymmetric() const
Check whether matrix is symmetric.
virtual TMatrixTBase< Element > & SetMatrixArray(const Element *data, Option_t *option="")
Copy array data to matrix .
Int_t GetRowUpb() const
Definition: TMatrixTLazy.h:111
Int_t GetRowLwb() const
Definition: TMatrixTLazy.h:110
virtual void FillIn(TMatrixTSym< Element > &m) const =0
TMatrixTSym.
Definition: TMatrixTSym.h:34
TMatrixTSym< Element > & operator*=(Element val)
Multiply every element of the matrix with val.
virtual const Element * GetMatrixArray() const
Definition: TMatrixTSym.h:189
TMatrixTSym< Element > & operator+=(Element val)
Add val to every element of the matrix.
virtual Double_t Determinant() const
virtual TMatrixTSym< Element > & RandomizePD(Element alpha, Element beta, Double_t &seed)
randomize matrix element values but keep matrix symmetric positive definite
virtual const Int_t * GetColIndexArray() const
Definition: TMatrixTSym.h:86
TMatrixTSym< Element > & Transpose(const TMatrixTSym< Element > &source)
Transpose a matrix.
Element * New_m(Int_t size)
return data pointer .
TMatrixTSym< Element > & Use(Int_t row_lwb, Int_t row_upb, Element *data)
TMatrixTSym< Element > & Rank1Update(const TVectorT< Element > &v, Element alpha=1.0)
Perform a rank 1 operation on the matrix: A += alpha * v * v^T.
void Minus(const TMatrixTSym< Element > &a, const TMatrixTSym< Element > &b)
Symmetric matrix summation. Create a matrix C such that C = A + B.
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.
TMatrixTSym< Element > & Invert(Double_t *det=0)
Invert the matrix and calculate its determinant Notice that the LU decomposition is used instead of B...
void TMult(const TMatrixT< Element > &a)
Create a matrix C such that C = A' * A.
void Delete_m(Int_t size, Element *&)
delete data pointer m, if it was assigned on the heap
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.
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 TMatrixTBase< Element > & SetMatrixArray(const Element *data, Option_t *option="")
Copy array data to matrix .
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 .
const TMatrixT< Element > EigenVectors(TVectorT< Element > &eigenValues) const
Return a matrix containing the eigen-vectors ordered by descending eigen-values.
TMatrixTSym< Element > & operator-=(Element val)
Subtract val from every element of the matrix.
void Plus(const TMatrixTSym< Element > &a, const TMatrixTSym< Element > &b)
Symmetric matrix summation. Create a matrix C such that C = A + B.
virtual const Int_t * GetRowIndexArray() const
Definition: TMatrixTSym.h:84
TMatrixTBase< Element > & Apply(const TElementActionT< Element > &action)
Apply action to each matrix element.
TMatrixTSym< Element > & SimilarityT(const TMatrixT< Element > &n)
Calculate B^T * (*this) * B , final matrix will be (ncolsb x ncolsb) It is more efficient than applyi...
virtual TMatrixTBase< Element > & Randomize(Element alpha, Element beta, Double_t &seed)
randomize matrix element values but keep matrix symmetric
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...
TMatrixTSym< Element > & operator=(const TMatrixTSym< Element > &source)
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....
Element * fElements
data container
Definition: TMatrixTSym.h:39
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...
TMatrixT.
Definition: TMatrixT.h:39
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:1053
virtual const Element * GetMatrixArray() const
Definition: TMatrixT.h:222
void TMult(const TMatrixT< Element > &a, const TMatrixT< Element > &b)
Create a matrix C such that C = A' * B.
Definition: TMatrixT.cxx:852
TObject & operator=(const TObject &rhs)
TObject assignment operator.
Definition: TObject.h:271
Basic string class.
Definition: TString.h:131
void ToUpper()
Change string to upper case.
Definition: TString.cxx:1113
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:619
TVectorT.
Definition: TVectorT.h:27
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Definition: TVectorT.cxx:292
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.
double beta(double x, double y)
Calculates the beta function.
TCppObject_t Allocate(TCppType_t type)
Definition: Cppyy.cxx:271
void AMultB(int n, int m, int k, const double *A, const double *B, double *C)
Definition: cblas.cxx:11
int Invert(LASymMatrix &)
Definition: LaInverse.cxx:22
static constexpr double m2
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:212
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:723
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:180
auto * m
Definition: textangle.C:8
auto * a
Definition: textangle.C:12
REAL epsilon
Definition: triangle.c:617