Logo ROOT  
Reference Guide
SMatrix.icc
Go to the documentation of this file.
1// @(#)root/smatrix:$Id$
2// Authors: T. Glebe, L. Moneta 2005
3
4#ifndef ROOT_Math_SMatrix_icc
5#define ROOT_Math_SMatrix_icc
6// ********************************************************************
7//
8// source:
9//
10// type: source code
11//
12// created: 21. Mar 2001
13//
14// author: Thorsten Glebe
15// HERA-B Collaboration
16// Max-Planck-Institut fuer Kernphysik
17// Saupfercheckweg 1
18// 69117 Heidelberg
19// Germany
20// E-mail: T.Glebe@mpi-hd.mpg.de
21//
22// Description: SMatrix implementation file
23//
24// changes:
25// 21 Mar 2001 (TG) creation
26// 26 Mar 2001 (TG) place_in_row(), place_in_col() added
27// 03 Apr 2001 (TG) invert() added
28// 07 Apr 2001 (TG) CTOR from SVertex (dyadic product) added
29// 09 Apr 2001 (TG) CTOR from array added
30// 25 Mai 2001 (TG) row(), col() added
31// 11 Jul 2001 (TG) added #include Functions.hh
32// 11 Jan 2002 (TG) added operator==(), operator!=()
33// 14 Jan 2002 (TG) added more operator==(), operator!=(), operator>(), operator<()
34//
35// ********************************************************************
36
37#ifndef ROOT_Math_SMatrix
38#error "Do not use SMatrix.icc directly. #include \"Math/SMatrix.h\" instead."
39#endif // ROOT_Math_SMatrix
40
41#include <iostream>
42#include <iomanip>
43#include <assert.h>
44//#ifndef ROOT_Math_Dsinv
45//#include "Math/Dsinv.h"
46//#endif
47//#include "Math/Dsinv_array.h"
48//#include "Math/Dsfact.h"
49
50#include "Math/Dfact.h"
51#include "Math/Dinv.h"
52#include "Math/Functions.h"
53#include "Math/HelperOps.h"
54#include "Math/StaticCheck.h"
55
56
57
58
59
60
61
62namespace ROOT {
63
64namespace Math {
65
66
67
68//==============================================================================
69// Constructors
70//==============================================================================
71template <class T, unsigned int D1, unsigned int D2, class R>
73 // operator=(0); // if operator=(T ) is defined
74 for(unsigned int i=0; i<R::kSize; ++i) fRep.Array()[i] = 0;
75}
76
77//identity
78template <class T, unsigned int D1, unsigned int D2, class R>
80 for(unsigned int i=0; i<R::kSize; ++i)
81 fRep.Array()[i] = 0;
82 if (D1 <= D2) {
83 for(unsigned int i=0; i<D1; ++i)
84 fRep[i*D2+i] = 1;
85 }
86 else {
87 for(unsigned int i=0; i<D2; ++i)
88 fRep[i*D2+i] = 1;
89 }
90}
91
92template <class T, unsigned int D1, unsigned int D2, class R>
94 fRep = rhs.fRep;
95}
96
97
98template <class T, unsigned int D1, unsigned int D2, class R>
99template <class R2>
101 operator=(rhs);
102}
103
104
105template <class T, unsigned int D1, unsigned int D2, class R>
106template <class A, class R2>
108 operator=(rhs);
109}
110
111
112//=============================================================================
113// New Constructors from STL interfaces
114//=============================================================================
115template <class T, unsigned int D1, unsigned int D2, class R>
116template <class InputIterator>
117SMatrix<T,D1,D2,R>::SMatrix(InputIterator ibegin, InputIterator iend, bool triang, bool lower) {
118 // assume iterator size == matrix size
119 for(unsigned int i=0; i<R::kSize; ++i) fRep.Array()[i] = 0;
120 AssignItr<T,D1,D2,R>::Evaluate(*this,ibegin,iend,triang,lower);
121}
122
123template <class T, unsigned int D1, unsigned int D2, class R>
124template <class InputIterator>
125SMatrix<T,D1,D2,R>::SMatrix(InputIterator ibegin, unsigned int size, bool triang, bool lower) {
126 // assume iterator size <= matrix size (no check needed in AssignItr)
127 assert( size <= R::kSize);
128 for(unsigned int i=0; i<R::kSize; ++i) fRep.Array()[i] = 0;
129 AssignItr<T,D1,D2,R>::Evaluate(*this,ibegin,ibegin+size,triang,lower,false);
130}
131
132
133//==============================================================================
134// Assignment and operator= for scalar types for matrices of size 1
135// compiles only for matrices of size 1
136//==============================================================================
137template <class T, unsigned int D1, unsigned int D2, class R>
139 STATIC_CHECK( kSize == 1,SMatrix_NOT_of_size_1 );
140 fRep[0] = rhs;
141}
142
143template <class T, unsigned int D1, unsigned int D2, class R>
145 STATIC_CHECK( kSize == 1,SMatrix_NOT_of_size_1 );
146 fRep[0] = rhs;
147 return *this;
148}
149
150//=============================================================================
151//=============================================================================
152
153template <class T, unsigned int D1, unsigned int D2, class R>
154template <class M>
156 fRep = rhs.fRep;
157 return *this;
158}
159
160template <class T, unsigned int D1, unsigned int D2, class R>
162 fRep = rhs.fRep;
163 return *this;
164}
165
166template <class T, unsigned int D1, unsigned int D2, class R>
167template <class A, class R2>
169
171 return *this;
172}
173
174//=============================================================================
175// assign from an identity
176template <class T, unsigned int D1, unsigned int D2, class R>
178 for(unsigned int i=0; i<R::kSize; ++i)
179 fRep.Array()[i] = 0;
180 if (D1 <= D2) {
181 for(unsigned int i=0; i<D1; ++i)
182 fRep[i*D2+i] = 1;
183 }
184 else {
185 for(unsigned int i=0; i<D2; ++i)
186 fRep[i*D2+i] = 1;
187 }
188 return *this;
189}
190
191
192
193//=============================================================================
194// operator+=
195//=============================================================================
196template <class T, unsigned int D1, unsigned int D2, class R>
198 // self-addition with a scalar value
199 for(unsigned int i=0; i<R::kSize; ++i) {
200 fRep.Array()[i] += rhs;
201 }
202 return *this;
203}
204
205template <class T, unsigned int D1, unsigned int D2, class R>
206template <class R2>
208 // self-addition with another matrix (any representation)
209 // use operator+= of the representation object
210 fRep += rhs.fRep;
211 return *this;
212}
213
215template <class T, unsigned int D1, unsigned int D2, class R>
216template <class A, class R2>
218 // self-addition with an expression
220 return *this;
221}
222
223
224//==============================================================================
225// operator-=
226//==============================================================================
227template <class T, unsigned int D1, unsigned int D2, class R>
229 // self-subtraction with a scalar value
230 for(unsigned int i=0; i<R::kSize; ++i) {
231 fRep.Array()[i] -= rhs;
232 }
233 return *this;
234}
235
236template <class T, unsigned int D1, unsigned int D2, class R>
237template <class R2>
239 // self-subtraction with another matrix (any representation)
240 // use operator-= of the representation object
241 fRep -= rhs.fRep;
242 return *this;
243}
245
246template <class T, unsigned int D1, unsigned int D2, class R>
247template <class A, class R2>
249 // self-subtraction with an expression
251 return *this;
252}
253
254//==============================================================================
255// operator*=
256//==============================================================================
257template <class T, unsigned int D1, unsigned int D2, class R>
259 // case of multiplication with a scalar
260 for(unsigned int i=0; i<R::kSize; ++i) {
261 fRep.Array()[i] *= rhs;
262 }
263 return *this;
264}
265
266template <class T, unsigned int D1, unsigned int D2, class R>
267template <class R2>
269 // self-multiplication with another matrix (will work only for square matrices)
270 // a temporary is needed and will be created automatically to store intermediate result
271 return operator=(*this * rhs);
272}
274template <class T, unsigned int D1, unsigned int D2, class R>
275template <class A, class R2>
277 // self-multiplication with an expression (will work only for square matrices)
278 // a temporary is needed and will be created automatically to store intermediate result
279 return operator=(*this * rhs);
280}
281
282
283//==============================================================================
284// operator/= (only for scalar values)
285//==============================================================================
286template <class T, unsigned int D1, unsigned int D2, class R>
288 // division with a scalar
289 for(unsigned int i=0; i<R::kSize; ++i) {
290 fRep.Array()[i] /= rhs;
291 }
292 return *this;
293}
294
295//==============================================================================
296// operator== (element wise comparison)
297//==============================================================================
298template <class T, unsigned int D1, unsigned int D2, class R>
299bool SMatrix<T,D1,D2,R>::operator==(const T& rhs) const {
300 bool rc = true;
301 for(unsigned int i=0; i<R::kSize; ++i) {
302 rc = rc && (fRep.Array()[i] == rhs);
303 }
304 return rc;
305}
306
307template <class T, unsigned int D1, unsigned int D2, class R>
308template <class R2>
310 return fRep == rhs.fRep;
311}
312
313template <class T, unsigned int D1, unsigned int D2, class R>
314template <class A, class R2>
316 bool rc = true;
317 for(unsigned int i=0; i<D1*D2; ++i) {
318 rc = rc && (fRep[i] == rhs.apply(i));
319 }
320 return rc;
321}
322
323//==============================================================================
324// operator!= (element wise comparison)
325//==============================================================================
326template <class T, unsigned int D1, unsigned int D2, class R>
327inline bool SMatrix<T,D1,D2,R>::operator!=(const T& rhs) const {
328 return !operator==(rhs);
329}
330
331template <class T, unsigned int D1, unsigned int D2, class R>
333 return !operator==(rhs);
335
336template <class T, unsigned int D1, unsigned int D2, class R>
337template <class A, class R2>
339 return !operator==(rhs);
340}
342
343//==============================================================================
344// operator> (element wise comparison)
345//==============================================================================
346template <class T, unsigned int D1, unsigned int D2, class R>
347bool SMatrix<T,D1,D2,R>::operator>(const T& rhs) const {
348 bool rc = true;
349 for(unsigned int i=0; i<D1*D2; ++i) {
350 rc = rc && (fRep[i] > rhs);
351 }
352 return rc;
353}
354
355template <class T, unsigned int D1, unsigned int D2, class R>
356template <class R2>
358 bool rc = true;
359 for(unsigned int i=0; i<D1*D2; ++i) {
360 rc = rc && (fRep[i] > rhs.fRep[i]);
361 }
362 return rc;
364
365template <class T, unsigned int D1, unsigned int D2, class R>
366template <class A, class R2>
368 bool rc = true;
369 for(unsigned int i=0; i<D1*D2; ++i) {
370 rc = rc && (fRep[i] > rhs.apply(i));
371 }
372 return rc;
373}
375//==============================================================================
376// operator< (element wise comparison)
377//==============================================================================
378template <class T, unsigned int D1, unsigned int D2, class R>
379bool SMatrix<T,D1,D2,R>::operator<(const T& rhs) const {
380 bool rc = true;
381 for(unsigned int i=0; i<D1*D2; ++i) {
382 rc = rc && (fRep[i] < rhs);
383 }
384 return rc;
385}
386
387template <class T, unsigned int D1, unsigned int D2, class R>
388template <class R2>
390 bool rc = true;
391 for(unsigned int i=0; i<D1*D2; ++i) {
392 rc = rc && (fRep[i] < rhs.fRep[i]);
393 }
394 return rc;
395}
396
397template <class T, unsigned int D1, unsigned int D2, class R>
398template <class A, class R2>
400 bool rc = true;
401 for(unsigned int i=0; i<D1*D2; ++i) {
402 rc = rc && (fRep[i] < rhs.apply(i));
403 }
404 return rc;
405}
406
407
408//==============================================================================
409// invert
410//==============================================================================
411template <class T, unsigned int D1, unsigned int D2, class R>
413 STATIC_CHECK( D1 == D2,SMatrix_not_square);
414 return Inverter<D1,D2>::Dinv((*this).fRep);
415}
416
417// invert returning a new matrix
418template <class T, unsigned int D1, unsigned int D2, class R>
420 SMatrix<T,D1,D2,R> tmp(*this);
421 bool ok = tmp.Invert();
422 ifail = 0;
423 if (!ok) ifail = 1;
424 return tmp;
425}
426
427// fast inversion
428template <class T, unsigned int D1, unsigned int D2, class R>
430 STATIC_CHECK( D1 == D2,SMatrix_not_square);
431 return FastInverter<D1,D2>::Dinv((*this).fRep);
432}
433
434// fast inversion returning a new matrix
435template <class T, unsigned int D1, unsigned int D2, class R>
438 bool ok = tmp.InvertFast();
439 ifail = 0;
440 if (!ok) ifail = 1;
441 return tmp;
442}
444// Choleski inversion for symmetric and positive defined matrices
445template <class T, unsigned int D1, unsigned int D2, class R>
447 STATIC_CHECK( D1 == D2,SMatrix_not_square);
448 return CholInverter<D1>::Dinv((*this).fRep);
449}
450
451template <class T, unsigned int D1, unsigned int D2, class R>
453 SMatrix<T,D1,D2,R > tmp(*this);
454 bool ok = tmp.InvertChol();
455 ifail = 0;
456 if (!ok) ifail = 1;
457 return tmp;
458}
459
461
462//==============================================================================
463// det
464//==============================================================================
465template <class T, unsigned int D1, unsigned int D2, class R>
466inline bool SMatrix<T,D1,D2,R>::Det(T& det) {
467 STATIC_CHECK( D1 == D2,SMatrix_not_square);
468 // return Dfact<SMatrix<T,D1,D1>, D1, D1>(*this,det);
469 //return Dfact<R, D1, D1>(fRep, det);
470 return Determinant<D1,D2>::Dfact(fRep, det);
471}
472template <class T, unsigned int D1, unsigned int D2, class R>
473inline bool SMatrix<T,D1,D2,R>::Det2(T& det) const {
474 SMatrix<T,D1,D2,R> tmp(*this);
475 return tmp.Det(det);
476}
477
478
479//==============================================================================
480// place_in_row
481//==============================================================================
482template <class T, unsigned int D1, unsigned int D2, class R>
483template <unsigned int D>
485 unsigned int row,
486 unsigned int col) {
487
488 assert(col+D <= D2);
489
490 for(unsigned int i=row*D2+col, j=0; j<D; ++i, ++j) {
491 fRep[i] = rhs.apply(j);
492 }
493 return *this;
494}
495
496//==============================================================================
497// place_in_row
498//==============================================================================
499template <class T, unsigned int D1, unsigned int D2, class R>
500template <class A, unsigned int D>
502 unsigned int row,
503 unsigned int col) {
504
505 assert(col+D <= D2);
506
507 for(unsigned int i=row*D2+col, j=0; j<D; ++i, ++j) {
508 fRep[i] = rhs.apply(j);
509 }
510 return *this;
511}
512
513//==============================================================================
514// place_in_col
515//==============================================================================
516template <class T, unsigned int D1, unsigned int D2, class R>
517template <unsigned int D>
519 unsigned int row,
520 unsigned int col) {
521
522 assert(row+D <= D1);
523
524 for(unsigned int i=row*D2+col, j=0; j<D; i+=D2, ++j) {
525 fRep[i] = rhs.apply(j);
526 }
527 return *this;
529
530//==============================================================================
531// place_in_col
532//==============================================================================
533template <class T, unsigned int D1, unsigned int D2, class R>
534template <class A, unsigned int D>
536 unsigned int row,
537 unsigned int col) {
538
539 assert(row+D <= D1);
541 for(unsigned int i=row*D2+col, j=0; j<D; i+=D2, ++j) {
542 fRep[i] = rhs.apply(j);
543 }
544 return *this;
545}
546
547//==============================================================================
548// place_at
549//==============================================================================
550template <class T, unsigned int D1, unsigned int D2, class R>
551template <unsigned int D3, unsigned int D4, class R2>
553 unsigned int row,
554 unsigned int col) {
556 return *this;
557}
559//==============================================================================
560// place_at
561//==============================================================================
562template <class T, unsigned int D1, unsigned int D2, class R>
563template <class A, unsigned int D3, unsigned int D4, class R2>
565 unsigned int row,
566 unsigned int col) {
568 return *this;
569}
570
571//==============================================================================
572// row
573//==============================================================================
574template <class T, unsigned int D1, unsigned int D2, class R>
575SVector<T,D2> SMatrix<T,D1,D2,R>::Row(unsigned int therow) const {
576
577 const unsigned int offset = therow*D2;
579 /*static*/ SVector<T,D2> tmp;
580 for(unsigned int i=0; i<D2; ++i) {
581 tmp[i] = fRep[offset+i];
582 }
583 return tmp;
584}
585
586//==============================================================================
587// col
588//==============================================================================
589template <class T, unsigned int D1, unsigned int D2, class R>
590SVector<T,D1> SMatrix<T,D1,D2,R>::Col(unsigned int thecol) const {
591
592 /*static */ SVector<T,D1> tmp;
593 for(unsigned int i=0; i<D1; ++i) {
594 tmp[i] = fRep[thecol+i*D2];
595 }
596 return tmp;
597}
599//==============================================================================
600// print
601//==============================================================================
602template <class T, unsigned int D1, unsigned int D2, class R>
603std::ostream& SMatrix<T,D1,D2,R>::Print(std::ostream& os) const {
604 const std::ios_base::fmtflags prevFmt = os.setf(std::ios::right,std::ios::adjustfield);
605 // os.setf(ios::fixed);
606
607 os << "[ ";
608 for (unsigned int i=0; i < D1; ++i) {
609 for (unsigned int j=0; j < D2; ++j) {
610 os << std::setw(12) << fRep[i*D2+j];
611 if ((!((j+1)%12)) && (j < D2-1))
612 os << std::endl << " ...";
613 }
614 if (i != D1 - 1)
615 os << std::endl << " ";
616 }
617 os << " ]";
619 if (prevFmt != os.flags() ) os.setf(prevFmt, std::ios::adjustfield);
620 return os;
621}
622
623//==============================================================================
624// Access functions
625//==============================================================================
626template <class T, unsigned int D1, unsigned int D2, class R>
627inline T SMatrix<T,D1,D2,R>::apply(unsigned int i) const { return fRep[i]; }
628
629template <class T, unsigned int D1, unsigned int D2, class R>
630inline const T* SMatrix<T,D1,D2,R>::Array() const { return fRep.Array(); }
631
632template <class T, unsigned int D1, unsigned int D2, class R>
633inline T* SMatrix<T,D1,D2,R>::Array() { return fRep.Array(); }
635//==============================================================================
636// Operators
637//==============================================================================
638template <class T, unsigned int D1, unsigned int D2, class R>
639inline const T& SMatrix<T,D1,D2,R>::operator()(unsigned int i, unsigned int j) const {
640 return fRep(i,j);
641}
642
643template <class T, unsigned int D1, unsigned int D2, class R>
644inline T& SMatrix<T,D1,D2,R>::operator()(unsigned int i, unsigned int j) {
645 return fRep(i,j);
646}
647
649//==============================================================================
650// Element access with At()
651//==============================================================================
652template <class T, unsigned int D1, unsigned int D2, class R>
653inline const T& SMatrix<T,D1,D2,R>::At(unsigned int i, unsigned int j) const {
654 assert(i < D1);
655 assert(j < D2);
656 return fRep(i,j);
657}
658
659template <class T, unsigned int D1, unsigned int D2, class R>
660inline T& SMatrix<T,D1,D2,R>::At(unsigned int i, unsigned int j) {
661 assert(i < D1);
662 assert(j < D2);
663 return fRep(i,j);
664}
665
666//==============================================================================
667// STL interface
668//==============================================================================
669template <class T, unsigned int D1, unsigned int D2, class R>
671 return fRep.Array();
672}
673
674template <class T, unsigned int D1, unsigned int D2, class R>
676 return fRep.Array() + R::kSize;
677}
678
679template <class T, unsigned int D1, unsigned int D2, class R>
680inline const T * SMatrix<T,D1,D2,R>::begin() const {
681 return fRep.Array();
682}
683
684template <class T, unsigned int D1, unsigned int D2, class R>
685inline const T * SMatrix<T,D1,D2,R>::end() const {
686 return fRep.Array() + R::kSize;
687}
688
689
690template <class T, unsigned int D1, unsigned int D2, class R>
691template <class InputIterator>
692void SMatrix<T,D1,D2,R>::SetElements(InputIterator ibegin, InputIterator iend, bool triang, bool lower) {
693 // assume iterator size == matrix size when filling full matrix
694 AssignItr<T,D1,D2,R>::Evaluate(*this,ibegin,iend,triang,lower);
695}
696
697template <class T, unsigned int D1, unsigned int D2, class R>
698template <class InputIterator>
699void SMatrix<T,D1,D2,R>::SetElements(InputIterator ibegin, unsigned int size, bool triang, bool lower) {
700 // assume iterator size <= matrix size (no check to be done in AssignItr)
701 assert( size <= R::kSize);
702 AssignItr<T,D1,D2,R>::Evaluate(*this,ibegin,ibegin+size,triang,lower,false);
703}
704
705
706
707//==============================================================================
708// SubMatrices and slices of columns and rows
709//==============================================================================
710template <class T, unsigned int D1, unsigned int D2, class R>
711template <class SubVector>
712SubVector SMatrix<T,D1,D2,R>::SubRow(unsigned int therow, unsigned int col0 ) const {
713
714 const unsigned int offset = therow*D2 + col0;
715
716 STATIC_CHECK( SubVector::kSize <= D2,SVector_dimension_too_small);
717 assert(col0 + SubVector::kSize <= D2);
718
719 SubVector tmp;
720 for(unsigned int i=0; i<SubVector::kSize; ++i) {
721 tmp[i] = fRep[offset+i];
722 }
723 return tmp;
724}
725
726template <class T, unsigned int D1, unsigned int D2, class R>
727template <class SubVector>
728SubVector SMatrix<T,D1,D2,R>::SubCol(unsigned int thecol, unsigned int row0 ) const {
729
730 const unsigned int offset = thecol + row0*D1;
731
732 STATIC_CHECK( SubVector::kSize <= D1,SVector_dimension_too_small);
733 assert(row0 + SubVector::kSize <= D1);
734
735 SubVector tmp;
736 for(unsigned int i=0; i<SubVector::kSize; ++i) {
737 tmp[i] = fRep[offset+i*D1];
738 }
739 return tmp;
740}
741
742// submatrix
743template <class T, unsigned int D1, unsigned int D2, class R>
744template <class SubMatrix>
745SubMatrix SMatrix<T,D1,D2,R>::Sub(unsigned int row0, unsigned int col0) const {
746
747 SubMatrix tmp;
749 (tmp,*this,row0,col0);
750 return tmp;
751}
752
753//diagonal
754template <class T, unsigned int D1, unsigned int D2, class R>
756
757 // only for squared matrices
758 STATIC_CHECK( D1 == D2,SMatrix_NOT_square );
759
760 SVector<T,D1> tmp;
761 for(unsigned int i=0; i<D1; ++i) {
762 tmp[i] = fRep[ i*D2 + i];
763 }
764 return tmp;
765}
766
767//set diagonal
768template <class T, unsigned int D1, unsigned int D2, class R>
769template <class Vector>
770void SMatrix<T,D1,D2,R>::SetDiagonal( const Vector & v) {
771
772 // check size that size of vector is correct
773 STATIC_CHECK( ( ( D1 <= D2) && Vector::kSize == D1 ) ||
774 ( ( D2 < D1 ) && Vector::kSize == D2 ), SVector_size_NOT_correct );
775
776
777 for(unsigned int i=0; i<Vector::kSize; ++i) {
778 fRep[ i*D2 + i] = v[i];
779 }
780}
781
782// matrix trace
783template <class T, unsigned int D1, unsigned int D2, class R>
785 unsigned int diagSize = D1;
786 if (D2 < D1) diagSize = D2;
787 T trace = 0;
788 for(unsigned int i=0; i< diagSize; ++i) {
789 trace += fRep[ i*D2 + i] ;
790 }
791 return trace;
792}
793
794//upper block
795template <class T, unsigned int D1, unsigned int D2, class R>
796#ifndef UNSUPPORTED_TEMPLATE_EXPRESSION
797SVector<T, D1 * (D2 +1)/2 > SMatrix<T,D1,D2,R>::UpperBlock( ) const {
798#else
799template <class SubVector>
800SubVector SMatrix<T,D1,D2,R>::UpperBlock( ) const {
801#endif
802 // only for squared matrices
803 STATIC_CHECK( D1 == D2,SMatrix_NOT_square );
804
805#ifndef UNSUPPORTED_TEMPLATE_EXPRESSION
806 SVector<T, D1 * (D2 +1)/2 > tmp;
807#else
808 // N must be equal D1 *(D1 +1)/2
809 STATIC_CHECK( SubVector::kSize == D1*(D1+1)/2,SVector_Wrong_Size );
810 SubVector tmp;
811#endif
812
813 int k = 0;
814 for(unsigned int i=0; i<D1; ++i) {
815 for(unsigned int j=i; j<D2; ++j) {
816 tmp[k] = fRep[ i*D2 + j];
817 k++;
818 }
819 }
820 return tmp;
821}
822
823//lower block
824template <class T, unsigned int D1, unsigned int D2, class R>
825#ifndef UNSUPPORTED_TEMPLATE_EXPRESSION
826SVector<T, D1 * (D2 +1)/2 > SMatrix<T,D1,D2,R>::LowerBlock( ) const {
827#else
828template <class SubVector>
829SubVector SMatrix<T,D1,D2,R>::LowerBlock( ) const {
830#endif
831
832 // only for squared matrices
833 STATIC_CHECK( D1 == D2,SMatrix_NOT_square );
834
835#ifndef UNSUPPORTED_TEMPLATE_EXPRESSION
836 SVector<T, D1 * (D2 +1)/2 > tmp;
837#else
838 // N must be equal D1 *(D1 +1)/2
839 STATIC_CHECK( SubVector::kSize == D1*(D1+1)/2,SVector_Wrong_Size );
840 SubVector tmp;
841#endif
842
843 int k = 0;
844 for(unsigned int i=0; i<D1; ++i) {
845 for(unsigned int j=0; j<=i; ++j) {
846 tmp[k] = fRep[ i*D2 + j];
847 k++;
848 }
849 }
850 return tmp;
851}
852
853/// construct from upper/lower block
854
855//lower block
856template <class T, unsigned int D1, unsigned int D2, class R>
857#ifndef UNSUPPORTED_TEMPLATE_EXPRESSION
858SMatrix<T,D1,D2,R>::SMatrix(const SVector<T, D1*(D2+1)/2 > & v, bool lower ) {
859#else
860template <unsigned int N>
861SMatrix<T,D1,D2,R>::SMatrix(const SVector<T, N > & v, bool lower ) {
862#endif
863
864 // only for squared matrices
865 STATIC_CHECK( D1 == D2,SMatrix_NOT_square );
866
867#ifdef UNSUPPORTED_TEMPLATE_EXPRESSION
868 STATIC_CHECK( N == D1*(D1+1)/2,SVector_Wrong_Size );
869#endif
870
871 int k = 0;
872 if (lower) {
873 // case of lower block
874 for(unsigned int i=0; i<D1; ++i) {
875 for(unsigned int j=0; j<=i; ++j) {
876 fRep[ i*D2 + j] = v[k];
877 if ( i != j) fRep[ j*D2 + i] = v[k];
878 k++;
879 }
880 }
881 } else {
882 // case of upper block
883 for(unsigned int i=0; i<D1; ++i) {
884 for(unsigned int j=i; j<D2; ++j) {
885 fRep[ i*D2 + j] = v[k];
886 if ( i != j) fRep[ j*D2 + i] = v[k];
887 k++;
888 }
889 }
890 }
891}
892
893
894template <class T, unsigned int D1, unsigned int D2, class R>
895bool SMatrix<T,D1,D2,R>::IsInUse( const T * p) const {
896 return p == fRep.Array();
897}
898
899
900
901
902 } // namespace Math
903
904} // namespace ROOT
905
906
907
908#endif
#define STATIC_CHECK(expr, msg)
Definition: StaticCheck.h:56
Bool_t operator>(const TDatime &d1, const TDatime &d2)
Definition: TDatime.h:110
Bool_t operator!=(const TDatime &d1, const TDatime &d2)
Definition: TDatime.h:104
Bool_t operator<(const TDatime &d1, const TDatime &d2)
Definition: TDatime.h:106
Bool_t operator==(const TDatime &d1, const TDatime &d2)
Definition: TDatime.h:102
#define N
TRObject operator()(const T1 &t1) const
Binding & operator=(OUT(*fun)(void))
std::string & operator+=(std::string &left, const TString &right)
Definition: TString.h:473
@ kSize
Definition: TStructNode.h:26
Detrminant for a general squared matrix Function to compute the determinant from a square matrix ( ) ...
Definition: Dfact.h:49
T apply(unsigned int i) const
Definition: Expression.h:150
Fast Matrix Inverter class Class to specialize calls to Dinv.
Definition: Dinv.h:144
Matrix Inverter class Class to specialize calls to Dinv.
Definition: Dinv.h:69
SMatrix: a generic fixed size D1 x D2 Matrix class.
Definition: SMatrix.h:124
SMatrix()
Default constructor:
Definition: SMatrix.icc:72
bool Det(T &det)
determinant of square Matrix via Dfact.
Definition: SMatrix.icc:466
R fRep
Matrix Storage Object containing matrix data.
Definition: SMatrix.h:709
bool operator==(const T &rhs) const
element wise comparison
Definition: SMatrix.icc:299
bool InvertChol()
Invertion of a symmetric positive defined Matrix using Choleski decomposition.
Definition: SMatrix.icc:446
bool Invert()
Invert a square Matrix ( this method changes the current matrix).
Definition: SMatrix.icc:412
const T * Array() const
return read-only pointer to internal array
Definition: SMatrix.icc:630
SVector: a generic fixed size Vector class.
Definition: SVector.h:75
T apply(unsigned int i) const
access the parse tree. Index starts from zero
Definition: SVector.icc:537
Expression wrapper class for Vector objects.
Definition: Expression.h:64
T apply(unsigned int i) const
Definition: Expression.h:77
Namespace for new Math classes and functions.
double T(double x)
Definition: ChebyshevPol.h:34
void Print(std::ostream &os, const OptionType &opt)
int Invert(LASymMatrix &)
Definition: LaInverse.cxx:22
ABObj< sym, MatrixInverse< sym, ABObj< sym, LASymMatrix, double >, double >, double > Inverse(const ABObj< sym, LASymMatrix, double > &obj)
LAPACK Algebra functions specialize the Invert function for LASymMatrix.
Definition: LaInverse.h:27
VSD Structures.
Definition: StringConv.hxx:21
const char * Array
Structure for assignment to a general matrix from iterator.
Definition: HelperOps.h:517
Structure to assign from an expression based to general matrix to general matrix.
Definition: HelperOps.h:49
Evaluate the expression performing a -= operation Need to check whether creating a temporary object w...
Definition: HelperOps.h:280
Structure to deal when a submatrix is placed in a matrix.
Definition: HelperOps.h:362
Evaluate the expression performing a += operation Need to check whether creating a temporary object w...
Definition: HelperOps.h:196
Structure for getting sub matrices We have different cases according to the matrix representations.
Definition: HelperOps.h:462