Logo ROOT   6.18/05
Reference Guide
TMatrixTUtils.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 TMatrixTUtils
13 \ingroup Matrix
14
15 Matrix utility classes.
16
17 Templates of utility classes in the Linear Algebra Package.
18 The following classes are defined here:
19
20 Different matrix views without copying data elements :
21~~~
22 TMatrixTRow_const TMatrixTRow
23 TMatrixTColumn_const TMatrixTColumn
24 TMatrixTDiag_const TMatrixTDiag
25 TMatrixTFlat_const TMatrixTFlat
26 TMatrixTSub_const TMatrixTSub
27 TMatrixTSparseRow_const TMatrixTSparseRow
28 TMatrixTSparseDiag_const TMatrixTSparseDiag
29
30 TElementActionT
31 TElementPosActionT
32~~~
33*/
34
35#include "TMatrixTUtils.h"
36#include "TMatrixT.h"
37#include "TMatrixTSym.h"
38#include "TMatrixTSparse.h"
39#include "TMath.h"
40#include "TVectorT.h"
41
42////////////////////////////////////////////////////////////////////////////////
43/// Constructor with row "row" of matrix
44
45template<class Element>
47{
48 R__ASSERT(matrix.IsValid());
49
50 fRowInd = row-matrix.GetRowLwb();
51 if (fRowInd >= matrix.GetNrows() || fRowInd < 0) {
52 Error("TMatrixTRow_const(const TMatrixT<Element> &,Int_t)","row index out of bounds");
53 fMatrix = 0;
54 fPtr = 0;
55 fInc = 0;
56 return;
57 }
58
59 fMatrix = &matrix;
60 fPtr = matrix.GetMatrixArray()+fRowInd*matrix.GetNcols();
61 fInc = 1;
62}
63
64////////////////////////////////////////////////////////////////////////////////
65/// Constructor with row "row" of symmetric matrix
66
67template<class Element>
69{
70 R__ASSERT(matrix.IsValid());
71
72 fRowInd = row-matrix.GetRowLwb();
73 if (fRowInd >= matrix.GetNrows() || fRowInd < 0) {
74 Error("TMatrixTRow_const(const TMatrixTSym &,Int_t)","row index out of bounds");
75 fMatrix = 0;
76 fPtr = 0;
77 fInc = 0;
78 return;
79 }
80
81 fMatrix = &matrix;
82 fPtr = matrix.GetMatrixArray()+fRowInd*matrix.GetNcols();
83 fInc = 1;
84}
85
86////////////////////////////////////////////////////////////////////////////////
87/// Constructor with row "row" of symmetric matrix
88
89template<class Element>
91 :TMatrixTRow_const<Element>(matrix,row)
92{
93}
94
95////////////////////////////////////////////////////////////////////////////////
96/// Constructor with row "row" of symmetric matrix
97
98template<class Element>
100 :TMatrixTRow_const<Element>(matrix,row)
101{
102}
103
104////////////////////////////////////////////////////////////////////////////////
105/// Copy constructor
106
107template<class Element>
109{
110 *this = mr;
111}
112
113////////////////////////////////////////////////////////////////////////////////
114/// Assign val to every element of the matrix row.
115
116template<class Element>
118{
119 R__ASSERT(this->fMatrix->IsValid());
120 Element *rp = const_cast<Element *>(this->fPtr);
121 for ( ; rp < this->fPtr+this->fMatrix->GetNcols(); rp += this->fInc)
122 *rp = val;
123}
124
125template<class Element>
126void TMatrixTRow<Element>::operator=(std::initializer_list<Element> l)
127{
128 R__ASSERT(this->fMatrix->IsValid());
129 Element *rp = const_cast<Element *>(this->fPtr);
130 auto litr = l.begin();
131 for ( ; rp < this->fPtr+this->fMatrix->GetNcols() && litr != l.end(); rp += this->fInc)
132 *rp = *litr++;
133}
134
135////////////////////////////////////////////////////////////////////////////////
136/// Add val to every element of the matrix row.
137
138template<class Element>
140{
141 R__ASSERT(this->fMatrix->IsValid());
142 Element *rp = const_cast<Element *>(this->fPtr);
143 for ( ; rp < this->fPtr+this->fMatrix->GetNcols(); rp += this->fInc)
144 *rp += val;
145}
146
147////////////////////////////////////////////////////////////////////////////////
148/// Multiply every element of the matrix row with val.
149
150template<class Element>
152{
153 R__ASSERT(this->fMatrix->IsValid());
154 Element *rp = const_cast<Element *>(this->fPtr);
155 for ( ; rp < this->fPtr + this->fMatrix->GetNcols(); rp += this->fInc)
156 *rp *= val;
157}
158
159////////////////////////////////////////////////////////////////////////////////
160/// Assignment operator
161
162template<class Element>
164{
165 const TMatrixTBase<Element> *mt = mr.GetMatrix();
166 if (this->fMatrix->GetMatrixArray() == mt->GetMatrixArray() && this->fRowInd == mr.GetRowIndex()) return;
167
168 R__ASSERT(this->fMatrix->IsValid());
169 R__ASSERT(mt->IsValid());
170
171 if (this->fMatrix->GetNcols() != mt->GetNcols() || this->fMatrix->GetColLwb() != mt->GetColLwb()) {
172 Error("operator=(const TMatrixTRow_const &)", "matrix rows not compatible");
173 return;
174 }
175
176 Element *rp1 = const_cast<Element *>(this->fPtr);
177 const Element *rp2 = mr.GetPtr();
178 for ( ; rp1 < this->fPtr+this->fMatrix->GetNcols(); rp1 += this->fInc,rp2 += mr.GetInc())
179 *rp1 = *rp2;
180}
181
182////////////////////////////////////////////////////////////////////////////////
183/// Assign a vector to a matrix row. The vector is considered row-vector
184/// to allow the assignment in the strict sense.
185
186template<class Element>
188{
189 R__ASSERT(this->fMatrix->IsValid());
190 R__ASSERT(vec.IsValid());
191
192 if (this->fMatrix->GetColLwb() != vec.GetLwb() || this->fMatrix->GetNcols() != vec.GetNrows()) {
193 Error("operator=(const TVectorT &)","vector length != matrix-row length");
194 return;
195 }
196
197 Element *rp = const_cast<Element *>(this->fPtr);
198 const Element *vp = vec.GetMatrixArray();
199 for ( ; rp < this->fPtr+this->fMatrix->GetNcols(); rp += this->fInc)
200 *rp = *vp++;
201}
202
203////////////////////////////////////////////////////////////////////////////////
204/// Add to every element of the matrix row the corresponding element of row r.
205
206template<class Element>
208{
209 const TMatrixTBase<Element> *mt = r.GetMatrix();
210
211 R__ASSERT(this->fMatrix->IsValid());
212 R__ASSERT(mt->IsValid());
213
214 if (this->fMatrix->GetColLwb() != mt->GetColLwb() || this->fMatrix->GetNcols() != mt->GetNcols()) {
215 Error("operator+=(const TMatrixTRow_const &)","different row lengths");
216 return;
217 }
218
219 Element *rp1 = const_cast<Element *>(this->fPtr);
220 const Element *rp2 = r.GetPtr();
221 for ( ; rp1 < this->fPtr+this->fMatrix->GetNcols(); rp1 += this->fInc,rp2 += r.GetInc())
222 *rp1 += *rp2;
223}
224
225////////////////////////////////////////////////////////////////////////////////
226/// Multiply every element of the matrix row with the
227/// corresponding element of row r.
228
229template<class Element>
231{
232 const TMatrixTBase<Element> *mt = r.GetMatrix();
233
234 R__ASSERT(this->fMatrix->IsValid());
235 R__ASSERT(mt->IsValid());
236
237 if (this->fMatrix->GetColLwb() != mt->GetColLwb() || this->fMatrix->GetNcols() != mt->GetNcols()) {
238 Error("operator*=(const TMatrixTRow_const &)","different row lengths");
239 return;
240 }
241
242 Element *rp1 = const_cast<Element *>(this->fPtr);
243 const Element *rp2 = r.GetPtr();
244 for ( ; rp1 < this->fPtr+this->fMatrix->GetNcols(); rp1 += this->fInc,rp2 += r.GetInc())
245 *rp1 *= *rp2;
246}
247
248////////////////////////////////////////////////////////////////////////////////
249/// Constructor with column "col" of matrix
250
251template<class Element>
253{
254 R__ASSERT(matrix.IsValid());
255
256 this->fColInd = col-matrix.GetColLwb();
257 if (this->fColInd >= matrix.GetNcols() || this->fColInd < 0) {
258 Error("TMatrixTColumn_const(const TMatrixT &,Int_t)","column index out of bounds");
259 fMatrix = 0;
260 fPtr = 0;
261 fInc = 0;
262 return;
263 }
264
265 fMatrix = &matrix;
266 fPtr = matrix.GetMatrixArray()+fColInd;
267 fInc = matrix.GetNcols();
268}
269
270////////////////////////////////////////////////////////////////////////////////
271/// Constructor with column "col" of matrix
272
273template<class Element>
275{
276 R__ASSERT(matrix.IsValid());
277
278 fColInd = col-matrix.GetColLwb();
279 if (fColInd >= matrix.GetNcols() || fColInd < 0) {
280 Error("TMatrixTColumn_const(const TMatrixTSym &,Int_t)","column index out of bounds");
281 fMatrix = 0;
282 fPtr = 0;
283 fInc = 0;
284 return;
285 }
286
287 fMatrix = &matrix;
288 fPtr = matrix.GetMatrixArray()+fColInd;
289 fInc = matrix.GetNcols();
290}
291
292////////////////////////////////////////////////////////////////////////////////
293/// Constructor with column "col" of matrix
294
295template<class Element>
297 :TMatrixTColumn_const<Element>(matrix,col)
298{
299}
300
301////////////////////////////////////////////////////////////////////////////////
302/// Constructor with column "col" of matrix
303
304template<class Element>
306 :TMatrixTColumn_const<Element>(matrix,col)
307{
308}
309
310////////////////////////////////////////////////////////////////////////////////
311/// Copy constructor
312
313template<class Element>
315{
316 *this = mc;
317}
318
319////////////////////////////////////////////////////////////////////////////////
320/// Assign val to every element of the matrix column.
321
322template<class Element>
324{
325 R__ASSERT(this->fMatrix->IsValid());
326 Element *cp = const_cast<Element *>(this->fPtr);
327 for ( ; cp < this->fPtr+this->fMatrix->GetNoElements(); cp += this->fInc)
328 *cp = val;
329}
330
331////////////////////////////////////////////////////////////////////////////////
332/// Assign element of the matrix column using given initializer list
333
334template<class Element>
335void TMatrixTColumn<Element>::operator=(std::initializer_list<Element> l)
336{
337 R__ASSERT(this->fMatrix->IsValid());
338 Element *rp = const_cast<Element *>(this->fPtr);
339 auto litr = l.begin();
340 for ( ; rp < this->fPtr+this->fMatrix->GetNoElements() && litr != l.end(); rp += this->fInc)
341 *rp = *litr++;
342}
343
344////////////////////////////////////////////////////////////////////////////////
345/// Add val to every element of the matrix column.
346
347template<class Element>
349{
350 R__ASSERT(this->fMatrix->IsValid());
351 Element *cp = const_cast<Element *>(this->fPtr);
352 for ( ; cp < this->fPtr+this->fMatrix->GetNoElements(); cp += this->fInc)
353 *cp += val;
354}
355
356////////////////////////////////////////////////////////////////////////////////
357/// Multiply every element of the matrix column with val.
358
359template<class Element>
361{
362 R__ASSERT(this->fMatrix->IsValid());
363 Element *cp = const_cast<Element *>(this->fPtr);
364 for ( ; cp < this->fPtr+this->fMatrix->GetNoElements(); cp += this->fInc)
365 *cp *= val;
366}
367
368////////////////////////////////////////////////////////////////////////////////
369/// Assignment operator
370
371template<class Element>
373{
374 const TMatrixTBase<Element> *mt = mc.GetMatrix();
375 if (this->fMatrix->GetMatrixArray() == mt->GetMatrixArray() && this->fColInd == mc.GetColIndex()) return;
376
377 R__ASSERT(this->fMatrix->IsValid());
378 R__ASSERT(mt->IsValid());
379
380 if (this->fMatrix->GetNrows() != mt->GetNrows() || this->fMatrix->GetRowLwb() != mt->GetRowLwb()) {
381 Error("operator=(const TMatrixTColumn_const &)", "matrix columns not compatible");
382 return;
383 }
384
385 Element *cp1 = const_cast<Element *>(this->fPtr);
386 const Element *cp2 = mc.GetPtr();
387 for ( ; cp1 < this->fPtr+this->fMatrix->GetNoElements(); cp1 += this->fInc,cp2 += mc.GetInc())
388 *cp1 = *cp2;
389}
390
391////////////////////////////////////////////////////////////////////////////////
392/// Assign a vector to a matrix column.
393
394template<class Element>
396{
397 R__ASSERT(this->fMatrix->IsValid());
398 R__ASSERT(vec.IsValid());
399
400 if (this->fMatrix->GetRowLwb() != vec.GetLwb() || this->fMatrix->GetNrows() != vec.GetNrows()) {
401 Error("operator=(const TVectorT &)","vector length != matrix-column length");
402 return;
403 }
404
405 Element *cp = const_cast<Element *>(this->fPtr);
406 const Element *vp = vec.GetMatrixArray();
407 for ( ; cp < this->fPtr+this->fMatrix->GetNoElements(); cp += this->fInc)
408 *cp = *vp++;
409
410 R__ASSERT(vp == vec.GetMatrixArray()+vec.GetNrows());
411}
412
413////////////////////////////////////////////////////////////////////////////////
414/// Add to every element of the matrix row the corresponding element of row mc.
415
416template<class Element>
418{
419 const TMatrixTBase<Element> *mt = mc.GetMatrix();
420
421 R__ASSERT(this->fMatrix->IsValid());
422 R__ASSERT(mt->IsValid());
423
424 if (this->fMatrix->GetRowLwb() != mt->GetRowLwb() || this->fMatrix->GetNrows() != mt->GetNrows()) {
425 Error("operator+=(const TMatrixTColumn_const &)","different row lengths");
426 return;
427 }
428
429 Element *cp1 = const_cast<Element *>(this->fPtr);
430 const Element *cp2 = mc.GetPtr();
431 for ( ; cp1 < this->fPtr+this->fMatrix->GetNoElements(); cp1 += this->fInc,cp2 += mc.GetInc())
432 *cp1 += *cp2;
433}
434
435////////////////////////////////////////////////////////////////////////////////
436/// Multiply every element of the matrix column with the
437/// corresponding element of column mc.
438
439template<class Element>
441{
442 const TMatrixTBase<Element> *mt = mc.GetMatrix();
443
444 R__ASSERT(this->fMatrix->IsValid());
445 R__ASSERT(mt->IsValid());
446
447 if (this->fMatrix->GetRowLwb() != mt->GetRowLwb() || this->fMatrix->GetNrows() != mt->GetNrows()) {
448 Error("operator*=(const TMatrixTColumn_const &)","different row lengths");
449 return;
450 }
451
452 Element *cp1 = const_cast<Element *>(this->fPtr);
453 const Element *cp2 = mc.GetPtr();
454 for ( ; cp1 < this->fPtr+this->fMatrix->GetNoElements(); cp1 += this->fInc,cp2 += mc.GetInc())
455 *cp1 *= *cp2;
456}
457
458////////////////////////////////////////////////////////////////////////////////
459/// Constructor
460
461template<class Element>
463{
464 R__ASSERT(matrix.IsValid());
465
466 fMatrix = &matrix;
467 fNdiag = TMath::Min(matrix.GetNrows(),matrix.GetNcols());
468 fPtr = matrix.GetMatrixArray();
469 fInc = matrix.GetNcols()+1;
470}
471
472////////////////////////////////////////////////////////////////////////////////
473/// Constructor
474
475template<class Element>
477{
478 R__ASSERT(matrix.IsValid());
479
480 fMatrix = &matrix;
481 fNdiag = TMath::Min(matrix.GetNrows(),matrix.GetNcols());
482 fPtr = matrix.GetMatrixArray();
483 fInc = matrix.GetNcols()+1;
484}
485
486////////////////////////////////////////////////////////////////////////////////
487/// Constructor
488
489template<class Element>
491 :TMatrixTDiag_const<Element>(matrix)
492{
493}
494
495////////////////////////////////////////////////////////////////////////////////
496/// Constructor
497
498template<class Element>
500 :TMatrixTDiag_const<Element>(matrix)
501{
502}
503
504////////////////////////////////////////////////////////////////////////////////
505/// Copy constructor
506
507template<class Element>
509{
510 *this = md;
511}
512
513////////////////////////////////////////////////////////////////////////////////
514/// Assign val to every element of the matrix diagonal.
515
516template<class Element>
518{
519 R__ASSERT(this->fMatrix->IsValid());
520 Element *dp = const_cast<Element *>(this->fPtr);
521 for (Int_t i = 0; i < this->fNdiag; i++, dp += this->fInc)
522 *dp = val;
523}
524
525////////////////////////////////////////////////////////////////////////////////
526/// Assign val to every element of the matrix diagonal.
527
528template<class Element>
530{
531 R__ASSERT(this->fMatrix->IsValid());
532 Element *dp = const_cast<Element *>(this->fPtr);
533 for (Int_t i = 0; i < this->fNdiag; i++, dp += this->fInc)
534 *dp += val;
535}
536
537////////////////////////////////////////////////////////////////////////////////
538/// Assign val to every element of the matrix diagonal.
539
540template<class Element>
542{
543 R__ASSERT(this->fMatrix->IsValid());
544 Element *dp = const_cast<Element *>(this->fPtr);
545 for (Int_t i = 0; i < this->fNdiag; i++, dp += this->fInc)
546 *dp *= val;
547}
548
549////////////////////////////////////////////////////////////////////////////////
550/// Assignment operator
551
552template<class Element>
554{
555 const TMatrixTBase<Element> *mt = md.GetMatrix();
556 if (this->fMatrix == mt) return;
557
558 R__ASSERT(this->fMatrix->IsValid());
559 R__ASSERT(mt->IsValid());
560
561 if (this->GetNdiags() != md.GetNdiags()) {
562 Error("operator=(const TMatrixTDiag_const &)","diagonals not compatible");
563 return;
564 }
565
566 Element *dp1 = const_cast<Element *>(this->fPtr);
567 const Element *dp2 = md.GetPtr();
568 for (Int_t i = 0; i < this->fNdiag; i++, dp1 += this->fInc, dp2 += md.GetInc())
569 *dp1 = *dp2;
570}
571
572////////////////////////////////////////////////////////////////////////////////
573/// Assign a vector to the matrix diagonal.
574
575template<class Element>
577{
578 R__ASSERT(this->fMatrix->IsValid());
579 R__ASSERT(vec.IsValid());
580
581 if (this->fNdiag != vec.GetNrows()) {
582 Error("operator=(const TVectorT &)","vector length != matrix-diagonal length");
583 return;
584 }
585
586 Element *dp = const_cast<Element *>(this->fPtr);
587 const Element *vp = vec.GetMatrixArray();
588 for ( ; vp < vec.GetMatrixArray()+vec.GetNrows(); dp += this->fInc)
589 *dp = *vp++;
590}
591
592////////////////////////////////////////////////////////////////////////////////
593/// Add to every element of the matrix diagonal the
594/// corresponding element of diagonal md.
595
596template<class Element>
598{
599 const TMatrixTBase<Element> *mt = md.GetMatrix();
600
601 R__ASSERT(this->fMatrix->IsValid());
602 R__ASSERT(mt->IsValid());
603 if (this->fNdiag != md.GetNdiags()) {
604 Error("operator=(const TMatrixTDiag_const &)","matrix-diagonal's different length");
605 return;
606 }
607
608 Element *dp1 = const_cast<Element *>(this->fPtr);
609 const Element *dp2 = md.GetPtr();
610 for (Int_t i = 0; i < this->fNdiag; i++, dp1 += this->fInc, dp2 += md.GetInc())
611 *dp1 += *dp2;
612}
613
614////////////////////////////////////////////////////////////////////////////////
615/// Multiply every element of the matrix diagonal with the
616/// corresponding element of diagonal md.
617
618template<class Element>
620{
621 const TMatrixTBase<Element> *mt = md.GetMatrix();
622
623 R__ASSERT(this->fMatrix->IsValid());
624 R__ASSERT(mt->IsValid());
625 if (this->fNdiag != md.GetNdiags()) {
626 Error("operator*=(const TMatrixTDiag_const &)","matrix-diagonal's different length");
627 return;
628 }
629
630 Element *dp1 = const_cast<Element *>(this->fPtr);
631 const Element *dp2 = md.GetPtr();
632 for (Int_t i = 0; i < this->fNdiag; i++, dp1 += this->fInc, dp2 += md.GetInc())
633 *dp1 *= *dp2;
634}
635
636////////////////////////////////////////////////////////////////////////////////
637/// Constructor
638
639template<class Element>
641{
642 R__ASSERT(matrix.IsValid());
643
644 fMatrix = &matrix;
645 fPtr = matrix.GetMatrixArray();
646 fNelems = matrix.GetNoElements();
647}
648
649////////////////////////////////////////////////////////////////////////////////
650/// Constructor
651
652template<class Element>
654{
655 R__ASSERT(matrix.IsValid());
656
657 fMatrix = &matrix;
658 fPtr = matrix.GetMatrixArray();
659 fNelems = matrix.GetNoElements();
660}
661
662////////////////////////////////////////////////////////////////////////////////
663/// Constructor
664
665template<class Element>
667 :TMatrixTFlat_const<Element>(matrix)
668{
669}
670
671////////////////////////////////////////////////////////////////////////////////
672/// Constructor
673
674template<class Element>
676 :TMatrixTFlat_const<Element>(matrix)
677{
678}
679
680////////////////////////////////////////////////////////////////////////////////
681/// Copy constructor
682
683template<class Element>
685{
686 *this = mf;
687}
688
689////////////////////////////////////////////////////////////////////////////////
690/// Assign val to every element of the matrix.
691
692template<class Element>
694{
695 R__ASSERT(this->fMatrix->IsValid());
696 Element *fp = const_cast<Element *>(this->fPtr);
697 while (fp < this->fPtr+this->fMatrix->GetNoElements())
698 *fp++ = val;
699}
700
701////////////////////////////////////////////////////////////////////////////////
702/// Add val to every element of the matrix.
703
704template<class Element>
706{
707 R__ASSERT(this->fMatrix->IsValid());
708 Element *fp = const_cast<Element *>(this->fPtr);
709 while (fp < this->fPtr+this->fMatrix->GetNoElements())
710 *fp++ += val;
711}
712
713////////////////////////////////////////////////////////////////////////////////
714/// Multiply every element of the matrix with val.
715
716template<class Element>
718{
719 R__ASSERT(this->fMatrix->IsValid());
720 Element *fp = const_cast<Element *>(this->fPtr);
721 while (fp < this->fPtr+this->fMatrix->GetNoElements())
722 *fp++ *= val;
723}
724
725////////////////////////////////////////////////////////////////////////////////
726/// Assignment operator
727
728template<class Element>
730{
731 const TMatrixTBase<Element> *mt = mf.GetMatrix();
732 if (this->fMatrix->GetMatrixArray() == mt->GetMatrixArray()) return;
733
734 R__ASSERT(this->fMatrix->IsValid());
735 R__ASSERT(mt->IsValid());
736 if (this->fMatrix->GetNoElements() != mt->GetNoElements()) {
737 Error("operator=(const TMatrixTFlat_const &)","matrix lengths different");
738 return;
739 }
740
741 Element *fp1 = const_cast<Element *>(this->fPtr);
742 const Element *fp2 = mf.GetPtr();
743 while (fp1 < this->fPtr+this->fMatrix->GetNoElements())
744 *fp1++ = *fp2++;
745}
746
747////////////////////////////////////////////////////////////////////////////////
748/// Assign a vector to the matrix. The matrix is traversed row-wise
749
750template<class Element>
752{
753 R__ASSERT(vec.IsValid());
754
755 if (this->fMatrix->GetNoElements() != vec.GetNrows()) {
756 Error("operator=(const TVectorT &)","vector length != # matrix-elements");
757 return;
758 }
759
760 Element *fp = const_cast<Element *>(this->fPtr);
761 const Element *vp = vec.GetMatrixArray();
762 while (fp < this->fPtr+this->fMatrix->GetNoElements())
763 *fp++ = *vp++;
764}
765
766////////////////////////////////////////////////////////////////////////////////
767/// Add to every element of the matrix the corresponding element of matrix mf.
768
769template<class Element>
771{
772 const TMatrixTBase<Element> *mt = mf.GetMatrix();
773
774 R__ASSERT(this->fMatrix->IsValid());
775 R__ASSERT(mt->IsValid());
776 if (this->fMatrix->GetNoElements() != mt->GetNoElements()) {
777 Error("operator+=(const TMatrixTFlat_const &)","matrices lengths different");
778 return;
779 }
780
781 Element *fp1 = const_cast<Element *>(this->fPtr);
782 const Element *fp2 = mf.GetPtr();
783 while (fp1 < this->fPtr + this->fMatrix->GetNoElements())
784 *fp1++ += *fp2++;
785}
786
787////////////////////////////////////////////////////////////////////////////////
788/// Multiply every element of the matrix with the corresponding element of diagonal mf.
789
790template<class Element>
792{
793 const TMatrixTBase<Element> *mt = mf.GetMatrix();
794
795 R__ASSERT(this->fMatrix->IsValid());
796 R__ASSERT(mt->IsValid());
797 if (this->fMatrix->GetNoElements() != mt->GetNoElements()) {
798 Error("operator*=(const TMatrixTFlat_const &)","matrices lengths different");
799 return;
800 }
801
802 Element *fp1 = const_cast<Element *>(this->fPtr);
803 const Element *fp2 = mf.GetPtr();
804 while (fp1 < this->fPtr + this->fMatrix->GetNoElements())
805 *fp1++ *= *fp2++;
806}
807
808////////////////////////////////////////////////////////////////////////////////
809/// make a reference to submatrix [row_lwbs..row_upbs][col_lwbs..col_upbs];
810/// The indexing range of the reference is
811/// [0..row_upbs-row_lwbs+1][0..col_upb-col_lwbs+1] (default)
812
813template<class Element>
815 Int_t col_lwbs,Int_t col_upbs)
816{
817 R__ASSERT(matrix.IsValid());
818
819 fRowOff = 0;
820 fColOff = 0;
821 fNrowsSub = 0;
822 fNcolsSub = 0;
823 fMatrix = &matrix;
824
825 if (row_upbs < row_lwbs) {
826 Error("TMatrixTSub_const","Request sub-matrix with row_upbs(%d) < row_lwbs(%d)",row_upbs,row_lwbs);
827 return;
828 }
829 if (col_upbs < col_lwbs) {
830 Error("TMatrixTSub_const","Request sub-matrix with col_upbs(%d) < col_lwbs(%d)",col_upbs,col_lwbs);
831 return;
832 }
833
834 const Int_t rowLwb = matrix.GetRowLwb();
835 const Int_t rowUpb = matrix.GetRowUpb();
836 const Int_t colLwb = matrix.GetColLwb();
837 const Int_t colUpb = matrix.GetColUpb();
838
839 if (row_lwbs < rowLwb || row_lwbs > rowUpb) {
840 Error("TMatrixTSub_const","Request row_lwbs sub-matrix(%d) outside matrix range of %d - %d",row_lwbs,rowLwb,rowUpb);
841 return;
842 }
843 if (col_lwbs < colLwb || col_lwbs > colUpb) {
844 Error("TMatrixTSub_const","Request col_lwbs sub-matrix(%d) outside matrix range of %d - %d",col_lwbs,colLwb,colUpb);
845 return;
846 }
847 if (row_upbs < rowLwb || row_upbs > rowUpb) {
848 Error("TMatrixTSub_const","Request row_upbs sub-matrix(%d) outside matrix range of %d - %d",row_upbs,rowLwb,rowUpb);
849 return;
850 }
851 if (col_upbs < colLwb || col_upbs > colUpb) {
852 Error("TMatrixTSub_const","Request col_upbs sub-matrix(%d) outside matrix range of %d - %d",col_upbs,colLwb,colUpb);
853 return;
854 }
855
856 fRowOff = row_lwbs-rowLwb;
857 fColOff = col_lwbs-colLwb;
858 fNrowsSub = row_upbs-row_lwbs+1;
859 fNcolsSub = col_upbs-col_lwbs+1;
860}
861
862////////////////////////////////////////////////////////////////////////////////
863/// make a reference to submatrix [row_lwbs..row_upbs][col_lwbs..col_upbs];
864/// The indexing range of the reference is
865/// [0..row_upbs-row_lwbs+1][0..col_upb-col_lwbs+1] (default)
866
867template<class Element>
869 Int_t col_lwbs,Int_t col_upbs)
870{
871 R__ASSERT(matrix.IsValid());
872
873 fRowOff = 0;
874 fColOff = 0;
875 fNrowsSub = 0;
876 fNcolsSub = 0;
877 fMatrix = &matrix;
878
879 if (row_upbs < row_lwbs) {
880 Error("TMatrixTSub_const","Request sub-matrix with row_upbs(%d) < row_lwbs(%d)",row_upbs,row_lwbs);
881 return;
882 }
883 if (col_upbs < col_lwbs) {
884 Error("TMatrixTSub_const","Request sub-matrix with col_upbs(%d) < col_lwbs(%d)",col_upbs,col_lwbs);
885 return;
886 }
887
888 const Int_t rowLwb = matrix.GetRowLwb();
889 const Int_t rowUpb = matrix.GetRowUpb();
890 const Int_t colLwb = matrix.GetColLwb();
891 const Int_t colUpb = matrix.GetColUpb();
892
893 if (row_lwbs < rowLwb || row_lwbs > rowUpb) {
894 Error("TMatrixTSub_const","Request row_lwbs sub-matrix(%d) outside matrix range of %d - %d",row_lwbs,rowLwb,rowUpb);
895 return;
896 }
897 if (col_lwbs < colLwb || col_lwbs > colUpb) {
898 Error("TMatrixTSub_const","Request col_lwbs sub-matrix(%d) outside matrix range of %d - %d",col_lwbs,colLwb,colUpb);
899 return;
900 }
901 if (row_upbs < rowLwb || row_upbs > rowUpb) {
902 Error("TMatrixTSub_const","Request row_upbs sub-matrix(%d) outside matrix range of %d - %d",row_upbs,rowLwb,rowUpb);
903 return;
904 }
905 if (col_upbs < colLwb || col_upbs > colUpb) {
906 Error("TMatrixTSub_const","Request col_upbs sub-matrix(%d) outside matrix range of %d - %d",col_upbs,colLwb,colUpb);
907 return;
908 }
909
910 fRowOff = row_lwbs-rowLwb;
911 fColOff = col_lwbs-colLwb;
912 fNrowsSub = row_upbs-row_lwbs+1;
913 fNcolsSub = col_upbs-col_lwbs+1;
914}
915
916////////////////////////////////////////////////////////////////////////////////
917/// Constructor
918
919template<class Element>
921 Int_t col_lwbs,Int_t col_upbs)
922 :TMatrixTSub_const<Element>(matrix,row_lwbs,row_upbs,col_lwbs,col_upbs)
923{
924}
925
926////////////////////////////////////////////////////////////////////////////////
927/// Constructor
928
929template<class Element>
931 Int_t col_lwbs,Int_t col_upbs)
932 :TMatrixTSub_const<Element>(matrix,row_lwbs,row_upbs,col_lwbs,col_upbs)
933{
934}
935
936////////////////////////////////////////////////////////////////////////////////
937/// Copy constructor
938
939template<class Element>
941{
942 *this = ms;
943}
944
945////////////////////////////////////////////////////////////////////////////////
946/// Perform a rank 1 operation on the matrix:
947/// A += alpha * v * v^T
948
949template<class Element>
951{
952 R__ASSERT(this->fMatrix->IsValid());
953 R__ASSERT(v.IsValid());
954
955 if (v.GetNoElements() < TMath::Max(this->fNrowsSub,this->fNcolsSub)) {
956 Error("Rank1Update","vector too short");
957 return;
958 }
959
960 const Element * const pv = v.GetMatrixArray();
961 Element *mp = (const_cast<TMatrixTBase<Element> *>(this->fMatrix))->GetMatrixArray();
962
963 const Int_t ncols = this->fMatrix->GetNcols();
964 for (Int_t irow = 0; irow < this->fNrowsSub; irow++) {
965 const Int_t off = (irow+this->fRowOff)*ncols+this->fColOff;
966 const Element tmp = alpha*pv[irow];
967 for (Int_t icol = 0; icol < this->fNcolsSub; icol++)
968 mp[off+icol] += tmp*pv[icol];
969 }
970}
971
972////////////////////////////////////////////////////////////////////////////////
973/// Assign val to every element of the sub matrix.
974
975template<class Element>
977{
978 R__ASSERT(this->fMatrix->IsValid());
979
980 Element *p = (const_cast<TMatrixTBase<Element> *>(this->fMatrix))->GetMatrixArray();
981 const Int_t ncols = this->fMatrix->GetNcols();
982 for (Int_t irow = 0; irow < this->fNrowsSub; irow++) {
983 const Int_t off = (irow+this->fRowOff)*ncols+this->fColOff;
984 for (Int_t icol = 0; icol < this->fNcolsSub; icol++)
985 p[off+icol] = val;
986 }
987}
988
989////////////////////////////////////////////////////////////////////////////////
990/// Add val to every element of the sub matrix.
991
992template<class Element>
994{
995 R__ASSERT(this->fMatrix->IsValid());
996
997 Element *p = (const_cast<TMatrixTBase<Element> *>(this->fMatrix))->GetMatrixArray();
998 const Int_t ncols = this->fMatrix->GetNcols();
999 for (Int_t irow = 0; irow < this->fNrowsSub; irow++) {
1000 const Int_t off = (irow+this->fRowOff)*ncols+this->fColOff;
1001 for (Int_t icol = 0; icol < this->fNcolsSub; icol++)
1002 p[off+icol] += val;
1003 }
1004}
1005
1006////////////////////////////////////////////////////////////////////////////////
1007/// Multiply every element of the sub matrix by val .
1008
1009template<class Element>
1011{
1012 R__ASSERT(this->fMatrix->IsValid());
1013
1014 Element *p = (const_cast<TMatrixTBase<Element> *>(this->fMatrix))->GetMatrixArray();
1015 const Int_t ncols = this->fMatrix->GetNcols();
1016 for (Int_t irow = 0; irow < this->fNrowsSub; irow++) {
1017 const Int_t off = (irow+this->fRowOff)*ncols+this->fColOff;
1018 for (Int_t icol = 0; icol < this->fNcolsSub; icol++)
1019 p[off+icol] *= val;
1020 }
1021}
1022
1023////////////////////////////////////////////////////////////////////////////////
1024/// Assignment operator
1025
1026template<class Element>
1028{
1029 const TMatrixTBase<Element> *mt = ms.GetMatrix();
1030
1031 R__ASSERT(this->fMatrix->IsValid());
1032 R__ASSERT(mt->IsValid());
1033
1034 if (this->fMatrix == mt &&
1035 (this->GetNrows() == ms.GetNrows () && this->GetNcols() == ms.GetNcols () &&
1036 this->GetRowOff() == ms.GetRowOff() && this->GetColOff() == ms.GetColOff()) )
1037 return;
1038
1039 if (this->GetNrows() != ms.GetNrows() || this->GetNcols() != ms.GetNcols()) {
1040 Error("operator=(const TMatrixTSub_const &)","sub matrices have different size");
1041 return;
1042 }
1043
1044 const Int_t rowOff2 = ms.GetRowOff();
1045 const Int_t colOff2 = ms.GetColOff();
1046
1047 Bool_t overlap = (this->fMatrix == mt) &&
1048 ( (rowOff2 >= this->fRowOff && rowOff2 < this->fRowOff+this->fNrowsSub) ||
1049 (colOff2 >= this->fColOff && colOff2 < this->fColOff+this->fNcolsSub) );
1050
1051 Element *p1 = const_cast<Element *>(this->fMatrix->GetMatrixArray());
1052 if (!overlap) {
1053 const Element *p2 = mt->GetMatrixArray();
1054
1055 const Int_t ncols1 = this->fMatrix->GetNcols();
1056 const Int_t ncols2 = mt->GetNcols();
1057 for (Int_t irow = 0; irow < this->fNrowsSub; irow++) {
1058 const Int_t off1 = (irow+this->fRowOff)*ncols1+this->fColOff;
1059 const Int_t off2 = (irow+rowOff2)*ncols2+colOff2;
1060 for (Int_t icol = 0; icol < this->fNcolsSub; icol++)
1061 p1[off1+icol] = p2[off2+icol];
1062 }
1063 } else {
1064 const Int_t row_lwbs = rowOff2+mt->GetRowLwb();
1065 const Int_t row_upbs = row_lwbs+this->fNrowsSub-1;
1066 const Int_t col_lwbs = colOff2+mt->GetColLwb();
1067 const Int_t col_upbs = col_lwbs+this->fNcolsSub-1;
1068 TMatrixT<Element> tmp; mt->GetSub(row_lwbs,row_upbs,col_lwbs,col_upbs,tmp);
1069 const Element *p2 = tmp.GetMatrixArray();
1070
1071 const Int_t ncols1 = this->fMatrix->GetNcols();
1072 const Int_t ncols2 = tmp.GetNcols();
1073 for (Int_t irow = 0; irow < this->fNrowsSub; irow++) {
1074 const Int_t off1 = (irow+this->fRowOff)*ncols1+this->fColOff;
1075 const Int_t off2 = irow*ncols2;
1076 for (Int_t icol = 0; icol < this->fNcolsSub; icol++)
1077 p1[off1+icol] = p2[off2+icol];
1078 }
1079 }
1080}
1081
1082////////////////////////////////////////////////////////////////////////////////
1083/// Assignment operator
1084
1085template<class Element>
1087{
1088 R__ASSERT(this->fMatrix->IsValid());
1089 R__ASSERT(m.IsValid());
1090
1091 if (this->fMatrix->GetMatrixArray() == m.GetMatrixArray()) return;
1092
1093 if (this->fNrowsSub != m.GetNrows() || this->fNcolsSub != m.GetNcols()) {
1094 Error("operator=(const TMatrixTBase<Element> &)","sub matrices and matrix have different size");
1095 return;
1096 }
1097 const Int_t row_lwbs = this->fRowOff+this->fMatrix->GetRowLwb();
1098 const Int_t col_lwbs = this->fColOff+this->fMatrix->GetColLwb();
1099 (const_cast<TMatrixTBase<Element> *>(this->fMatrix))->SetSub(row_lwbs,col_lwbs,m);
1100}
1101
1102////////////////////////////////////////////////////////////////////////////////
1103/// Add to every element of the submatrix the corresponding element of submatrix ms.
1104
1105template<class Element>
1107{
1108 const TMatrixTBase<Element> *mt = ms.GetMatrix();
1109
1110 R__ASSERT(this->fMatrix->IsValid());
1111 R__ASSERT(mt->IsValid());
1112
1113 if (this->GetNrows() != ms.GetNrows() || this->GetNcols() != ms.GetNcols()) {
1114 Error("operator+=(const TMatrixTSub_const &)","sub matrices have different size");
1115 return;
1116 }
1117
1118 const Int_t rowOff2 = ms.GetRowOff();
1119 const Int_t colOff2 = ms.GetColOff();
1120
1121 Bool_t overlap = (this->fMatrix == mt) &&
1122 ( (rowOff2 >= this->fRowOff && rowOff2 < this->fRowOff+this->fNrowsSub) ||
1123 (colOff2 >= this->fColOff && colOff2 < this->fColOff+this->fNcolsSub) );
1124
1125 Element *p1 = const_cast<Element *>(this->fMatrix->GetMatrixArray());
1126 if (!overlap) {
1127 const Element *p2 = mt->GetMatrixArray();
1128
1129 const Int_t ncols1 = this->fMatrix->GetNcols();
1130 const Int_t ncols2 = mt->GetNcols();
1131 for (Int_t irow = 0; irow < this->fNrowsSub; irow++) {
1132 const Int_t off1 = (irow+this->fRowOff)*ncols1+this->fColOff;
1133 const Int_t off2 = (irow+rowOff2)*ncols2+colOff2;
1134 for (Int_t icol = 0; icol < this->fNcolsSub; icol++)
1135 p1[off1+icol] += p2[off2+icol];
1136 }
1137 } else {
1138 const Int_t row_lwbs = rowOff2+mt->GetRowLwb();
1139 const Int_t row_upbs = row_lwbs+this->fNrowsSub-1;
1140 const Int_t col_lwbs = colOff2+mt->GetColLwb();
1141 const Int_t col_upbs = col_lwbs+this->fNcolsSub-1;
1142 TMatrixT<Element> tmp; mt->GetSub(row_lwbs,row_upbs,col_lwbs,col_upbs,tmp);
1143 const Element *p2 = tmp.GetMatrixArray();
1144
1145 const Int_t ncols1 = this->fMatrix->GetNcols();
1146 const Int_t ncols2 = tmp.GetNcols();
1147 for (Int_t irow = 0; irow < this->fNrowsSub; irow++) {
1148 const Int_t off1 = (irow+this->fRowOff)*ncols1+this->fColOff;
1149 const Int_t off2 = irow*ncols2;
1150 for (Int_t icol = 0; icol < this->fNcolsSub; icol++)
1151 p1[off1+icol] += p2[off2+icol];
1152 }
1153 }
1154}
1155
1156////////////////////////////////////////////////////////////////////////////////
1157/// Multiply submatrix with submatrix ms.
1158
1159template<class Element>
1161{
1162 if (this->fNcolsSub != ms.GetNrows() || this->fNcolsSub != ms.GetNcols()) {
1163 Error("operator*=(const TMatrixTSub_const &)","source sub matrix has wrong shape");
1164 return;
1165 }
1166
1167 const TMatrixTBase<Element> *source = ms.GetMatrix();
1168
1169 TMatrixT<Element> source_sub;
1170 {
1171 const Int_t row_lwbs = ms.GetRowOff()+source->GetRowLwb();
1172 const Int_t row_upbs = row_lwbs+this->fNrowsSub-1;
1173 const Int_t col_lwbs = ms.GetColOff()+source->GetColLwb();
1174 const Int_t col_upbs = col_lwbs+this->fNcolsSub-1;
1175 source->GetSub(row_lwbs,row_upbs,col_lwbs,col_upbs,source_sub);
1176 }
1177
1178 const Element *sp = source_sub.GetMatrixArray();
1179 const Int_t ncols = this->fMatrix->GetNcols();
1180
1181 // One row of the old_target matrix
1182 Element work[kWorkMax];
1183 Bool_t isAllocated = kFALSE;
1184 Element *trp = work;
1185 if (this->fNcolsSub > kWorkMax) {
1186 isAllocated = kTRUE;
1187 trp = new Element[this->fNcolsSub];
1188 }
1189
1190 Element *cp = const_cast<Element *>(this->fMatrix->GetMatrixArray())+this->fRowOff*ncols+this->fColOff;
1191 const Element *trp0 = cp; // Pointer to target[i,0];
1192 const Element * const trp0_last = trp0+this->fNrowsSub*ncols;
1193 while (trp0 < trp0_last) {
1194 memcpy(trp,trp0,this->fNcolsSub*sizeof(Element)); // copy the i-th row of target, Start at target[i,0]
1195 for (const Element *scp = sp; scp < sp+this->fNcolsSub; ) { // Pointer to the j-th column of source,
1196 // Start scp = source[0,0]
1197 Element cij = 0;
1198 for (Int_t j = 0; j < this->fNcolsSub; j++) {
1199 cij += trp[j] * *scp; // the j-th col of source
1200 scp += this->fNcolsSub;
1201 }
1202 *cp++ = cij;
1203 scp -= source_sub.GetNoElements()-1; // Set bcp to the (j+1)-th col
1204 }
1205 cp += ncols-this->fNcolsSub;
1206 trp0 += ncols; // Set trp0 to the (i+1)-th row
1207 R__ASSERT(trp0 == cp);
1208 }
1209
1210 R__ASSERT(cp == trp0_last && trp0 == trp0_last);
1211 if (isAllocated)
1212 delete [] trp;
1213}
1214
1215////////////////////////////////////////////////////////////////////////////////
1216/// Add to every element of the submatrix the corresponding element of matrix mt.
1217
1218template<class Element>
1220{
1221 R__ASSERT(this->fMatrix->IsValid());
1222 R__ASSERT(mt.IsValid());
1223
1224 if (this->GetNrows() != mt.GetNrows() || this->GetNcols() != mt.GetNcols()) {
1225 Error("operator+=(const TMatrixTBase<Element> &)","sub matrix and matrix have different size");
1226 return;
1227 }
1228
1229 Element *p1 = const_cast<Element *>(this->fMatrix->GetMatrixArray());
1230 const Element *p2 = mt.GetMatrixArray();
1231
1232 const Int_t ncols1 = this->fMatrix->GetNcols();
1233 const Int_t ncols2 = mt.GetNcols();
1234 for (Int_t irow = 0; irow < this->fNrowsSub; irow++) {
1235 const Int_t off1 = (irow+this->fRowOff)*ncols1+this->fColOff;
1236 const Int_t off2 = irow*ncols2;
1237 for (Int_t icol = 0; icol < this->fNcolsSub; icol++)
1238 p1[off1+icol] += p2[off2+icol];
1239 }
1240}
1241
1242////////////////////////////////////////////////////////////////////////////////
1243/// Multiply submatrix with matrix source.
1244
1245template<class Element>
1247{
1248 if (this->fNcolsSub != source.GetNrows() || this->fNcolsSub != source.GetNcols()) {
1249 Error("operator*=(const TMatrixT<Element> &)","source matrix has wrong shape");
1250 return;
1251 }
1252
1253 // Check for A *= A;
1254 const Element *sp;
1256 if (this->fMatrix->GetMatrixArray() == source.GetMatrixArray()) {
1257 tmp.ResizeTo(source);
1258 tmp = source;
1259 sp = tmp.GetMatrixArray();
1260 }
1261 else
1262 sp = source.GetMatrixArray();
1263
1264 const Int_t ncols = this->fMatrix->GetNcols();
1265
1266 // One row of the old_target matrix
1267 Element work[kWorkMax];
1268 Bool_t isAllocated = kFALSE;
1269 Element *trp = work;
1270 if (this->fNcolsSub > kWorkMax) {
1271 isAllocated = kTRUE;
1272 trp = new Element[this->fNcolsSub];
1273 }
1274
1275 Element *cp = const_cast<Element *>(this->fMatrix->GetMatrixArray())+this->fRowOff*ncols+this->fColOff;
1276 const Element *trp0 = cp; // Pointer to target[i,0];
1277 const Element * const trp0_last = trp0+this->fNrowsSub*ncols;
1278 while (trp0 < trp0_last) {
1279 memcpy(trp,trp0,this->fNcolsSub*sizeof(Element)); // copy the i-th row of target, Start at target[i,0]
1280 for (const Element *scp = sp; scp < sp+this->fNcolsSub; ) { // Pointer to the j-th column of source,
1281 // Start scp = source[0,0]
1282 Element cij = 0;
1283 for (Int_t j = 0; j < this->fNcolsSub; j++) {
1284 cij += trp[j] * *scp; // the j-th col of source
1285 scp += this->fNcolsSub;
1286 }
1287 *cp++ = cij;
1288 scp -= source.GetNoElements()-1; // Set bcp to the (j+1)-th col
1289 }
1290 cp += ncols-this->fNcolsSub;
1291 trp0 += ncols; // Set trp0 to the (i+1)-th row
1292 R__ASSERT(trp0 == cp);
1293 }
1294
1295 R__ASSERT(cp == trp0_last && trp0 == trp0_last);
1296 if (isAllocated)
1297 delete [] trp;
1298}
1299
1300////////////////////////////////////////////////////////////////////////////////
1301/// Multiply submatrix with matrix source.
1302
1303template<class Element>
1305{
1306 if (this->fNcolsSub != source.GetNrows() || this->fNcolsSub != source.GetNcols()) {
1307 Error("operator*=(const TMatrixTSym<Element> &)","source matrix has wrong shape");
1308 return;
1309 }
1310
1311 // Check for A *= A;
1312 const Element *sp;
1314 if (this->fMatrix->GetMatrixArray() == source.GetMatrixArray()) {
1315 tmp.ResizeTo(source);
1316 tmp = source;
1317 sp = tmp.GetMatrixArray();
1318 }
1319 else
1320 sp = source.GetMatrixArray();
1321
1322 const Int_t ncols = this->fMatrix->GetNcols();
1323
1324 // One row of the old_target matrix
1325 Element work[kWorkMax];
1326 Bool_t isAllocated = kFALSE;
1327 Element *trp = work;
1328 if (this->fNcolsSub > kWorkMax) {
1329 isAllocated = kTRUE;
1330 trp = new Element[this->fNcolsSub];
1331 }
1332
1333 Element *cp = const_cast<Element *>(this->fMatrix->GetMatrixArray())+this->fRowOff*ncols+this->fColOff;
1334 const Element *trp0 = cp; // Pointer to target[i,0];
1335 const Element * const trp0_last = trp0+this->fNrowsSub*ncols;
1336 while (trp0 < trp0_last) {
1337 memcpy(trp,trp0,this->fNcolsSub*sizeof(Element)); // copy the i-th row of target, Start at target[i,0]
1338 for (const Element *scp = sp; scp < sp+this->fNcolsSub; ) { // Pointer to the j-th column of source,
1339 // Start scp = source[0,0]
1340 Element cij = 0;
1341 for (Int_t j = 0; j < this->fNcolsSub; j++) {
1342 cij += trp[j] * *scp; // the j-th col of source
1343 scp += this->fNcolsSub;
1344 }
1345 *cp++ = cij;
1346 scp -= source.GetNoElements()-1; // Set bcp to the (j+1)-th col
1347 }
1348 cp += ncols-this->fNcolsSub;
1349 trp0 += ncols; // Set trp0 to the (i+1)-th row
1350 R__ASSERT(trp0 == cp);
1351 }
1352
1353 R__ASSERT(cp == trp0_last && trp0 == trp0_last);
1354 if (isAllocated)
1355 delete [] trp;
1356}
1357
1358////////////////////////////////////////////////////////////////////////////////
1359/// Constructor with row "row" of matrix
1360
1361template<class Element>
1363{
1364 R__ASSERT(matrix.IsValid());
1365
1366 fRowInd = row-matrix.GetRowLwb();
1367 if (fRowInd >= matrix.GetNrows() || fRowInd < 0) {
1368 Error("TMatrixTSparseRow_const(const TMatrixTSparse &,Int_t)","row index out of bounds");
1369 fMatrix = 0;
1370 fNindex = 0;
1371 fColPtr = 0;
1372 fDataPtr = 0;
1373 return;
1374 }
1375
1376 const Int_t sIndex = matrix.GetRowIndexArray()[fRowInd];
1377 const Int_t eIndex = matrix.GetRowIndexArray()[fRowInd+1];
1378 fMatrix = &matrix;
1379 fNindex = eIndex-sIndex;
1380 fColPtr = matrix.GetColIndexArray()+sIndex;
1381 fDataPtr = matrix.GetMatrixArray()+sIndex;
1382}
1383
1384////////////////////////////////////////////////////////////////////////////////
1385
1386template<class Element>
1388{
1389 if (!fMatrix) return TMatrixTBase<Element>::NaNValue();
1390 R__ASSERT(fMatrix->IsValid());
1391 const Int_t acoln = i-fMatrix->GetColLwb();
1392 if (acoln < fMatrix->GetNcols() && acoln >= 0) {
1393 const Int_t index = TMath::BinarySearch(fNindex,fColPtr,acoln);
1394 if (index >= 0 && fColPtr[index] == acoln) return fDataPtr[index];
1395 else return 0.0;
1396 } else {
1397 Error("operator()","Request col(%d) outside matrix range of %d - %d",
1398 i,fMatrix->GetColLwb(),fMatrix->GetColLwb()+fMatrix->GetNcols());
1400 }
1401 }
1402
1403////////////////////////////////////////////////////////////////////////////////
1404/// Constructor with row "row" of matrix
1405
1406template<class Element>
1408 : TMatrixTSparseRow_const<Element>(matrix,row)
1409{
1410}
1411
1412////////////////////////////////////////////////////////////////////////////////
1413/// Copy constructor
1414
1415template<class Element>
1417 : TMatrixTSparseRow_const<Element>(mr)
1418{
1419 *this = mr;
1420}
1421
1422////////////////////////////////////////////////////////////////////////////////
1423
1424template<class Element>
1426{
1427 if (!this->fMatrix) return TMatrixTBase<Element>::NaNValue();
1428 R__ASSERT(this->fMatrix->IsValid());
1429 const Int_t acoln = i-this->fMatrix->GetColLwb();
1430 if (acoln < this->fMatrix->GetNcols() && acoln >= 0) {
1431 const Int_t index = TMath::BinarySearch(this->fNindex,this->fColPtr,acoln);
1432 if (index >= 0 && this->fColPtr[index] == acoln) return this->fDataPtr[index];
1433 else return 0.0;
1434 } else {
1435 Error("operator()","Request col(%d) outside matrix range of %d - %d",
1436 i,this->fMatrix->GetColLwb(),this->fMatrix->GetColLwb()+this->fMatrix->GetNcols());
1438 }
1439}
1440
1441////////////////////////////////////////////////////////////////////////////////
1442/// operator() : pick element row(i)
1443
1444template<class Element>
1446{
1447 if (!this->fMatrix) return TMatrixTBase<Element>::NaNValue();
1448 R__ASSERT(this->fMatrix->IsValid());
1449
1450 const Int_t acoln = i-this->fMatrix->GetColLwb();
1451 if (acoln >= this->fMatrix->GetNcols() || acoln < 0) {
1452 Error("operator()(Int_t","Requested element %d outside range : %d - %d",i,
1453 this->fMatrix->GetColLwb(),this->fMatrix->GetColLwb()+this->fMatrix->GetNcols());
1455 }
1456
1457 Int_t index = TMath::BinarySearch(this->fNindex,this->fColPtr,acoln);
1458 if (index >= 0 && this->fColPtr[index] == acoln)
1459 return (const_cast<Element*>(this->fDataPtr))[index];
1460 else {
1461 TMatrixTSparse<Element> *mt = const_cast<TMatrixTSparse<Element> *>(this->fMatrix);
1462 const Int_t row = this->fRowInd+mt->GetRowLwb();
1463 Element val = 0.;
1464 mt->InsertRow(row,i,&val,1);
1465 const Int_t sIndex = mt->GetRowIndexArray()[this->fRowInd];
1466 const Int_t eIndex = mt->GetRowIndexArray()[this->fRowInd+1];
1467 this->fNindex = eIndex-sIndex;
1468 this->fColPtr = mt->GetColIndexArray()+sIndex;
1469 this->fDataPtr = mt->GetMatrixArray()+sIndex;
1470 index = TMath::BinarySearch(this->fNindex,this->fColPtr,acoln);
1471 if (index >= 0 && this->fColPtr[index] == acoln)
1472 return (const_cast<Element*>(this->fDataPtr))[index];
1473 else {
1474 Error("operator()(Int_t","Insert row failed");
1476 }
1477 }
1478}
1479
1480////////////////////////////////////////////////////////////////////////////////
1481/// Assign val to every non-zero (!) element of the matrix row.
1482
1483template<class Element>
1485{
1486 R__ASSERT(this->fMatrix->IsValid());
1487 Element *rp = const_cast<Element *>(this->fDataPtr);
1488 for ( ; rp < this->fDataPtr+this->fNindex; rp++)
1489 *rp = val;
1490}
1491
1492////////////////////////////////////////////////////////////////////////////////
1493/// Add val to every non-zero (!) element of the matrix row.
1494
1495template<class Element>
1497{
1498 R__ASSERT(this->fMatrix->IsValid());
1499 Element *rp = const_cast<Element *>(this->fDataPtr);
1500 for ( ; rp < this->fDataPtr+this->fNindex; rp++)
1501 *rp += val;
1502}
1503
1504////////////////////////////////////////////////////////////////////////////////
1505/// Multiply every element of the matrix row by val.
1506
1507template<class Element>
1509{
1510 R__ASSERT(this->fMatrix->IsValid());
1511 Element *rp = const_cast<Element *>(this->fDataPtr);
1512 for ( ; rp < this->fDataPtr+this->fNindex; rp++)
1513 *rp *= val;
1514}
1515
1516////////////////////////////////////////////////////////////////////////////////
1517/// Assignment operator
1518
1519template<class Element>
1521{
1522 const TMatrixTBase<Element> *mt = mr.GetMatrix();
1523 if (this->fMatrix == mt) return;
1524
1525 R__ASSERT(this->fMatrix->IsValid());
1526 R__ASSERT(mt->IsValid());
1527 if (this->fMatrix->GetColLwb() != mt->GetColLwb() || this->fMatrix->GetNcols() != mt->GetNcols()) {
1528 Error("operator=(const TMatrixTSparseRow_const &)","matrix rows not compatible");
1529 return;
1530 }
1531
1532 const Int_t ncols = this->fMatrix->GetNcols();
1533 const Int_t row1 = this->fRowInd+this->fMatrix->GetRowLwb();
1534 const Int_t row2 = mr.GetRowIndex()+mt->GetRowLwb();
1535 const Int_t col = this->fMatrix->GetColLwb();
1536
1537 TVectorT<Element> v(ncols);
1538 mt->ExtractRow(row2,col,v.GetMatrixArray());
1539 const_cast<TMatrixTSparse<Element> *>(this->fMatrix)->InsertRow(row1,col,v.GetMatrixArray());
1540
1541 const Int_t sIndex = this->fMatrix->GetRowIndexArray()[this->fRowInd];
1542 const Int_t eIndex = this->fMatrix->GetRowIndexArray()[this->fRowInd+1];
1543 this->fNindex = eIndex-sIndex;
1544 this->fColPtr = this->fMatrix->GetColIndexArray()+sIndex;
1545 this->fDataPtr = this->fMatrix->GetMatrixArray()+sIndex;
1546}
1547
1548////////////////////////////////////////////////////////////////////////////////
1549/// Assign a vector to a matrix row. The vector is considered row-vector
1550/// to allow the assignment in the strict sense.
1551
1552template<class Element>
1554{
1555 R__ASSERT(this->fMatrix->IsValid());
1556 R__ASSERT(vec.IsValid());
1557
1558 if (this->fMatrix->GetColLwb() != vec.GetLwb() || this->fMatrix->GetNcols() != vec.GetNrows()) {
1559 Error("operator=(const TVectorT &)","vector length != matrix-row length");
1560 return;
1561 }
1562
1563 const Element *vp = vec.GetMatrixArray();
1564 const Int_t row = this->fRowInd+this->fMatrix->GetRowLwb();
1565 const Int_t col = this->fMatrix->GetColLwb();
1566 const_cast<TMatrixTSparse<Element> *>(this->fMatrix)->InsertRow(row,col,vp,vec.GetNrows());
1567
1568 const Int_t sIndex = this->fMatrix->GetRowIndexArray()[this->fRowInd];
1569 const Int_t eIndex = this->fMatrix->GetRowIndexArray()[this->fRowInd+1];
1570 this->fNindex = eIndex-sIndex;
1571 this->fColPtr = this->fMatrix->GetColIndexArray()+sIndex;
1572 this->fDataPtr = this->fMatrix->GetMatrixArray()+sIndex;
1573}
1574
1575////////////////////////////////////////////////////////////////////////////////
1576/// Add to every element of the matrix row the corresponding element of row r.
1577
1578template<class Element>
1580{
1581 const TMatrixTBase<Element> *mt = r.GetMatrix();
1582
1583 R__ASSERT(this->fMatrix->IsValid());
1584 R__ASSERT(mt->IsValid());
1585 if (this->fMatrix->GetColLwb() != mt->GetColLwb() || this->fMatrix->GetNcols() != mt->GetNcols()) {
1586 Error("operator+=(const TMatrixTRow_const &)","different row lengths");
1587 return;
1588 }
1589
1590 const Int_t ncols = this->fMatrix->GetNcols();
1591 const Int_t row1 = this->fRowInd+this->fMatrix->GetRowLwb();
1592 const Int_t row2 = r.GetRowIndex()+mt->GetRowLwb();
1593 const Int_t col = this->fMatrix->GetColLwb();
1594
1595 TVectorT<Element> v1(ncols);
1596 TVectorT<Element> v2(ncols);
1597 this->fMatrix->ExtractRow(row1,col,v1.GetMatrixArray());
1598 mt ->ExtractRow(row2,col,v2.GetMatrixArray());
1599 v1 += v2;
1600 const_cast<TMatrixTSparse<Element> *>(this->fMatrix)->InsertRow(row1,col,v1.GetMatrixArray());
1601
1602 const Int_t sIndex = this->fMatrix->GetRowIndexArray()[this->fRowInd];
1603 const Int_t eIndex = this->fMatrix->GetRowIndexArray()[this->fRowInd+1];
1604 this->fNindex = eIndex-sIndex;
1605 this->fColPtr = this->fMatrix->GetColIndexArray()+sIndex;
1606 this->fDataPtr = this->fMatrix->GetMatrixArray()+sIndex;
1607}
1608
1609////////////////////////////////////////////////////////////////////////////////
1610/// Multiply every element of the matrix row with the
1611/// corresponding element of row r.
1612
1613template<class Element>
1615{
1616 const TMatrixTBase<Element> *mt = r.GetMatrix();
1617
1618 R__ASSERT(this->fMatrix->IsValid());
1619 R__ASSERT(mt->IsValid());
1620 if (this->fMatrix->GetColLwb() != mt->GetColLwb() || this->fMatrix->GetNcols() != mt->GetNcols()) {
1621 Error("operator+=(const TMatrixTRow_const &)","different row lengths");
1622 return;
1623 }
1624
1625 const Int_t ncols = this->fMatrix->GetNcols();
1626 const Int_t row1 = r.GetRowIndex()+mt->GetRowLwb();
1627 const Int_t row2 = r.GetRowIndex()+mt->GetRowLwb();
1628 const Int_t col = this->fMatrix->GetColLwb();
1629
1630 TVectorT<Element> v1(ncols);
1631 TVectorT<Element> v2(ncols);
1632 this->fMatrix->ExtractRow(row1,col,v1.GetMatrixArray());
1633 mt ->ExtractRow(row2,col,v2.GetMatrixArray());
1634
1635 ElementMult(v1,v2);
1636 const_cast<TMatrixTSparse<Element> *>(this->fMatrix)->InsertRow(row1,col,v1.GetMatrixArray());
1637
1638 const Int_t sIndex = this->fMatrix->GetRowIndexArray()[this->fRowInd];
1639 const Int_t eIndex = this->fMatrix->GetRowIndexArray()[this->fRowInd+1];
1640 this->fNindex = eIndex-sIndex;
1641 this->fColPtr = this->fMatrix->GetColIndexArray()+sIndex;
1642 this->fDataPtr = this->fMatrix->GetMatrixArray()+sIndex;
1643}
1644
1645////////////////////////////////////////////////////////////////////////////////
1646/// Constructor
1647
1648template<class Element>
1650{
1651 R__ASSERT(matrix.IsValid());
1652
1653 fMatrix = &matrix;
1654 fNdiag = TMath::Min(matrix.GetNrows(),matrix.GetNcols());
1655 fDataPtr = matrix.GetMatrixArray();
1656}
1657
1658////////////////////////////////////////////////////////////////////////////////
1659
1660template<class Element>
1662{
1663 R__ASSERT(fMatrix->IsValid());
1664 if (i < fNdiag && i >= 0) {
1665 const Int_t * const pR = fMatrix->GetRowIndexArray();
1666 const Int_t * const pC = fMatrix->GetColIndexArray();
1667 const Element * const pD = fMatrix->GetMatrixArray();
1668 const Int_t sIndex = pR[i];
1669 const Int_t eIndex = pR[i+1];
1670 const Int_t index = TMath::BinarySearch(eIndex-sIndex,pC+sIndex,i)+sIndex;
1671 if (index >= sIndex && pC[index] == i) return pD[index];
1672 else return 0.0;
1673 } else {
1674 Error("operator()","Request diagonal(%d) outside matrix range of 0 - %d",i,fNdiag);
1675 return 0.0;
1676 }
1677 return 0.0;
1678}
1679
1680////////////////////////////////////////////////////////////////////////////////
1681/// Constructor
1682
1683template<class Element>
1685 :TMatrixTSparseDiag_const<Element>(matrix)
1686{
1687}
1688
1689////////////////////////////////////////////////////////////////////////////////
1690/// Constructor
1691
1692template<class Element>
1694 : TMatrixTSparseDiag_const<Element>(md)
1695{
1696 *this = md;
1697}
1698
1699////////////////////////////////////////////////////////////////////////////////
1700
1701template<class Element>
1703{
1704 R__ASSERT(this->fMatrix->IsValid());
1705 if (i < this->fNdiag && i >= 0) {
1706 const Int_t * const pR = this->fMatrix->GetRowIndexArray();
1707 const Int_t * const pC = this->fMatrix->GetColIndexArray();
1708 const Element * const pD = this->fMatrix->GetMatrixArray();
1709 const Int_t sIndex = pR[i];
1710 const Int_t eIndex = pR[i+1];
1711 const Int_t index = TMath::BinarySearch(eIndex-sIndex,pC+sIndex,i)+sIndex;
1712 if (index >= sIndex && pC[index] == i) return pD[index];
1713 else return 0.0;
1714 } else {
1715 Error("operator()","Request diagonal(%d) outside matrix range of 0 - %d",i,this->fNdiag);
1716 return 0.0;
1717 }
1718 return 0.0;
1719 }
1720
1721////////////////////////////////////////////////////////////////////////////////
1722/// operator() : pick element diag(i)
1723
1724template<class Element>
1726{
1727 R__ASSERT(this->fMatrix->IsValid());
1728
1729 if (i < 0 || i >= this->fNdiag) {
1730 Error("operator()(Int_t","Requested element %d outside range : 0 - %d",i,this->fNdiag);
1731 return (const_cast<Element*>(this->fDataPtr))[0];
1732 }
1733
1734 TMatrixTSparse<Element> *mt = const_cast<TMatrixTSparse<Element> *>(this->fMatrix);
1735 const Int_t *pR = mt->GetRowIndexArray();
1736 const Int_t *pC = mt->GetColIndexArray();
1737 Int_t sIndex = pR[i];
1738 Int_t eIndex = pR[i+1];
1739 Int_t index = TMath::BinarySearch(eIndex-sIndex,pC+sIndex,i)+sIndex;
1740 if (index >= sIndex && pC[index] == i)
1741 return (const_cast<Element*>(this->fDataPtr))[index];
1742 else {
1743 const Int_t row = i+mt->GetRowLwb();
1744 const Int_t col = i+mt->GetColLwb();
1745 Element val = 0.;
1746 mt->InsertRow(row,col,&val,1);
1747 this->fDataPtr = mt->GetMatrixArray();
1748 pR = mt->GetRowIndexArray();
1749 pC = mt->GetColIndexArray();
1750 sIndex = pR[i];
1751 eIndex = pR[i+1];
1752 index = TMath::BinarySearch(eIndex-sIndex,pC+sIndex,i)+sIndex;
1753 if (index >= sIndex && pC[index] == i)
1754 return (const_cast<Element*>(this->fDataPtr))[index];
1755 else {
1756 Error("operator()(Int_t","Insert row failed");
1757 return (const_cast<Element*>(this->fDataPtr))[0];
1758 }
1759 }
1760}
1761
1762////////////////////////////////////////////////////////////////////////////////
1763/// Assign val to every element of the matrix diagonal.
1764
1765template<class Element>
1767{
1768 R__ASSERT(this->fMatrix->IsValid());
1769 for (Int_t i = 0; i < this->fNdiag; i++)
1770 (*this)(i) = val;
1771}
1772
1773////////////////////////////////////////////////////////////////////////////////
1774/// Add val to every element of the matrix diagonal.
1775
1776template<class Element>
1778{
1779 R__ASSERT(this->fMatrix->IsValid());
1780 for (Int_t i = 0; i < this->fNdiag; i++)
1781 (*this)(i) += val;
1782}
1783
1784////////////////////////////////////////////////////////////////////////////////
1785/// Multiply every element of the matrix diagonal by val.
1786
1787template<class Element>
1789{
1790 R__ASSERT(this->fMatrix->IsValid());
1791 for (Int_t i = 0; i < this->fNdiag; i++)
1792 (*this)(i) *= val;
1793}
1794
1795////////////////////////////////////////////////////////////////////////////////
1796/// Assignment operator
1797
1798template<class Element>
1800{
1801 const TMatrixTBase<Element> *mt = md.GetMatrix();
1802 if (this->fMatrix == mt) return;
1803
1804 R__ASSERT(this->fMatrix->IsValid());
1805 R__ASSERT(mt->IsValid());
1806 if (this->fNdiag != md.GetNdiags()) {
1807 Error("operator=(const TMatrixTSparseDiag_const &)","matrix-diagonal's different length");
1808 return;
1809 }
1810
1811 for (Int_t i = 0; i < this->fNdiag; i++)
1812 (*this)(i) = md(i);
1813}
1814
1815////////////////////////////////////////////////////////////////////////////////
1816/// Assign a vector to the matrix diagonal.
1817
1818template<class Element>
1820{
1821 R__ASSERT(this->fMatrix->IsValid());
1822 R__ASSERT(vec.IsValid());
1823
1824 if (this->fNdiag != vec.GetNrows()) {
1825 Error("operator=(const TVectorT &)","vector length != matrix-diagonal length");
1826 return;
1827 }
1828
1829 const Element *vp = vec.GetMatrixArray();
1830 for (Int_t i = 0; i < this->fNdiag; i++)
1831 (*this)(i) = vp[i];
1832}
1833
1834////////////////////////////////////////////////////////////////////////////////
1835/// Add to every element of the matrix diagonal the
1836/// corresponding element of diagonal md.
1837
1838template<class Element>
1840{
1841 const TMatrixTBase<Element> *mt = md.GetMatrix();
1842
1843 R__ASSERT(this->fMatrix->IsValid());
1844 R__ASSERT(mt->IsValid());
1845 if (this->fNdiag != md.GetNdiags()) {
1846 Error("operator+=(const TMatrixTSparseDiag_const &)","matrix-diagonal's different length");
1847 return;
1848 }
1849
1850 for (Int_t i = 0; i < this->fNdiag; i++)
1851 (*this)(i) += md(i);
1852}
1853
1854////////////////////////////////////////////////////////////////////////////////
1855/// Multiply every element of the matrix diagonal with the
1856/// corresponding element of diagonal md.
1857
1858template<class Element>
1860{
1861 const TMatrixTBase<Element> *mt = md.GetMatrix();
1862
1863 R__ASSERT(this->fMatrix->IsValid());
1864 R__ASSERT(mt->IsValid());
1865 if (this->fNdiag != md.GetNdiags()) {
1866 Error("operator*=(const TMatrixTSparseDiag_const &)","matrix-diagonal's different length");
1867 return;
1868 }
1869
1870 for (Int_t i = 0; i < this->fNdiag; i++)
1871 (*this)(i) *= md(i);
1872}
1873
1874////////////////////////////////////////////////////////////////////////////////
1875/// Random number generator [0....1] with seed ix
1876
1878{
1879 const Double_t a = 16807.0;
1880 const Double_t b15 = 32768.0;
1881 const Double_t b16 = 65536.0;
1882 const Double_t p = 2147483647.0;
1883 Double_t xhi = ix/b16;
1884 Int_t xhiint = (Int_t) xhi;
1885 xhi = xhiint;
1886 Double_t xalo = (ix-xhi*b16)*a;
1887
1888 Double_t leftlo = xalo/b16;
1889 Int_t leftloint = (int) leftlo;
1890 leftlo = leftloint;
1891 Double_t fhi = xhi*a+leftlo;
1892 Double_t k = fhi/b15;
1893 Int_t kint = (Int_t) k;
1894 k = kint;
1895 ix = (((xalo-leftlo*b16)-p)+(fhi-k*b15)*b16)+k;
1896 if (ix < 0.0) ix = ix+p;
1897
1898 return (ix*4.656612875e-10);
1899}
1900
1901template class TMatrixTRow_const <Float_t>;
1902template class TMatrixTColumn_const <Float_t>;
1903template class TMatrixTDiag_const <Float_t>;
1904template class TMatrixTFlat_const <Float_t>;
1905template class TMatrixTSub_const <Float_t>;
1908template class TMatrixTRow <Float_t>;
1909template class TMatrixTColumn <Float_t>;
1910template class TMatrixTDiag <Float_t>;
1911template class TMatrixTFlat <Float_t>;
1912template class TMatrixTSub <Float_t>;
1913template class TMatrixTSparseRow <Float_t>;
1914template class TMatrixTSparseDiag <Float_t>;
1915template class TElementActionT <Float_t>;
1916template class TElementPosActionT <Float_t>;
1917
1918template class TMatrixTRow_const <Double_t>;
1919template class TMatrixTColumn_const <Double_t>;
1920template class TMatrixTDiag_const <Double_t>;
1921template class TMatrixTFlat_const <Double_t>;
1922template class TMatrixTSub_const <Double_t>;
1925template class TMatrixTRow <Double_t>;
1926template class TMatrixTColumn <Double_t>;
1927template class TMatrixTDiag <Double_t>;
1928template class TMatrixTFlat <Double_t>;
1929template class TMatrixTSub <Double_t>;
1930template class TMatrixTSparseRow <Double_t>;
1931template class TMatrixTSparseDiag <Double_t>;
1932template class TElementActionT <Double_t>;
1933template class TElementPosActionT <Double_t>;
SVector< double, 2 > v
Definition: Dict.h:5
ROOT::R::TRInterface & r
Definition: Object.C:4
int Int_t
Definition: RtypesCore.h:41
const Bool_t kFALSE
Definition: RtypesCore.h:88
bool Bool_t
Definition: RtypesCore.h:59
double Double_t
Definition: RtypesCore.h:55
const Bool_t kTRUE
Definition: RtypesCore.h:87
#define R__ASSERT(e)
Definition: TError.h:96
void Error(const char *location, const char *msgfmt,...)
Double_t Drand(Double_t &ix)
Random number generator [0....1] with seed ix.
TMatrixT< Element > & ElementMult(TMatrixT< Element > &target, const TMatrixT< Element > &source)
Multiply target by the source, element-by-element.
Definition: TMatrixT.cxx:2981
Linear Algebra Package.
Definition: TMatrixTBase.h:85
virtual const Element * GetMatrixArray() const =0
Int_t GetNrows() const
Definition: TMatrixTBase.h:124
virtual TMatrixTBase< Element > & GetSub(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb, TMatrixTBase< Element > &target, Option_t *option="S") const =0
Int_t GetRowUpb() const
Definition: TMatrixTBase.h:123
Int_t GetRowLwb() const
Definition: TMatrixTBase.h:122
Int_t GetColLwb() const
Definition: TMatrixTBase.h:125
Int_t GetNoElements() const
Definition: TMatrixTBase.h:128
static Element & NaNValue()
Int_t GetColUpb() const
Definition: TMatrixTBase.h:126
Bool_t IsValid() const
Definition: TMatrixTBase.h:147
Int_t GetNcols() const
Definition: TMatrixTBase.h:127
virtual void ExtractRow(Int_t row, Int_t col, Element *v, Int_t n=-1) const
Store in array v, n matrix elements of row rown starting at column coln.
Int_t GetColIndex() const
const TMatrixTBase< Element > * GetMatrix() const
Int_t GetInc() const
const Element * GetPtr() const
void operator+=(Element val)
Add val to every element of the matrix column.
void operator=(Element val)
void Assign(Element val)
Assign val to every element of the matrix column.
void operator*=(Element val)
Multiply every element of the matrix column with val.
const Element * GetPtr() const
Int_t GetNdiags() const
const TMatrixTBase< Element > * GetMatrix() const
Int_t GetInc() const
void operator+=(Element val)
Assign val to every element of the matrix diagonal.
void operator=(Element val)
Assign val to every element of the matrix diagonal.
void operator*=(Element val)
Assign val to every element of the matrix diagonal.
const TMatrixTBase< Element > * GetMatrix() const
const Element * GetPtr() const
void operator+=(Element val)
Add val to every element of the matrix.
void operator*=(Element val)
Multiply every element of the matrix with val.
void operator=(Element val)
Assign val to every element of the matrix.
const Element * GetPtr() const
Int_t GetInc() const
Int_t GetRowIndex() const
const TMatrixTBase< Element > * GetMatrix() const
void Assign(Element val)
Assign val to every element of the matrix row.
void operator*=(Element val)
Multiply every element of the matrix row with val.
void operator+=(Element val)
Add val to every element of the matrix row.
void operator=(std::initializer_list< Element > l)
Element operator()(Int_t i) const
const TMatrixTBase< Element > * GetMatrix() const
Element operator()(Int_t i) const
void operator=(Element val)
Assign val to every element of the matrix diagonal.
void operator*=(Element val)
Multiply every element of the matrix diagonal by val.
void operator+=(Element val)
Add val to every element of the matrix diagonal.
Int_t GetRowIndex() const
Element operator()(Int_t i) const
const TMatrixTBase< Element > * GetMatrix() const
void operator+=(Element val)
Add val to every non-zero (!) element of the matrix row.
void operator=(Element val)
Assign val to every non-zero (!) element of the matrix row.
void operator*=(Element val)
Multiply every element of the matrix row by val.
Element operator()(Int_t i) const
TMatrixTSparse.
virtual const Int_t * GetRowIndexArray() const
virtual TMatrixTBase< Element > & InsertRow(Int_t row, Int_t col, const Element *v, Int_t n=-1)
Insert in row rown, n elements of array v at column coln.
virtual const Element * GetMatrixArray() const
virtual const Int_t * GetColIndexArray() const
void Rank1Update(const TVectorT< Element > &vec, Element alpha=1.0)
Perform a rank 1 operation on the matrix: A += alpha * v * v^T.
void operator*=(Element val)
Multiply every element of the sub matrix by val .
void operator+=(Element val)
Add val to every element of the sub matrix.
void operator=(Element val)
Assign val to every element of the sub matrix.
TMatrixTSym.
Definition: TMatrixTSym.h:34
virtual const Element * GetMatrixArray() const
Definition: TMatrixTSym.h:189
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
virtual TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t=-1)
Set size of the matrix to nrows x ncols New dynamic elements are created, the overlapping part of the...
Definition: TMatrixT.cxx:1210
virtual const Element * GetMatrixArray() const
Definition: TMatrixT.h:222
TVectorT.
Definition: TVectorT.h:27
Int_t GetNrows() const
Definition: TVectorT.h:75
Bool_t IsValid() const
Definition: TVectorT.h:83
Int_t GetLwb() const
Definition: TVectorT.h:73
Element * GetMatrixArray()
Definition: TVectorT.h:78
static constexpr double ms
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:212
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:180
Long64_t BinarySearch(Long64_t n, const T *array, T value)
Definition: TMathBase.h:278
auto * m
Definition: textangle.C:8
auto * l
Definition: textangle.C:4
auto * a
Definition: textangle.C:12