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