Logo ROOT  

Reference Guide
 
Loading...
Searching...
No Matches
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>
72SMatrix<T,D1,D2,R>::SMatrix() {
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>
79SMatrix<T,D1,D2,R>::SMatrix( SMatrixIdentity ) {
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>
96
97
98template <class T, unsigned int D1, unsigned int D2, class R>
99template <class R2>
103
104
105template <class T, unsigned int D1, unsigned int D2, class R>
106template <class A, class R2>
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>
138SMatrix<T,D1,D2,R>::SMatrix(const T& rhs) {
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>
144SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::operator=(const T& rhs) {
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>
155SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::operator=(const M& rhs) {
156 fRep = rhs.fRep;
157 return *this;
158}
159
160template <class T, unsigned int D1, unsigned int D2, class R>
161SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::operator=(const SMatrix<T,D1,D2,R>& rhs) {
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>
168SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::operator=(const Expr<A,T,D1,D2,R2>& rhs) {
169
170 Assign<T,D1,D2,A,R,R2>::Evaluate(*this, rhs);
171 return *this;
172}
173
174//=============================================================================
175// assign from an identity
176template <class T, unsigned int D1, unsigned int D2, class R>
177SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::operator= ( SMatrixIdentity ) {
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>
197SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::operator+=(const T& rhs) {
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>
207SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::operator+=(const SMatrix<T,D1,D2,R2>& rhs) {
208 // self-addition with another matrix (any representation)
209 // use operator+= of the representation object
210 fRep += rhs.fRep;
211 return *this;
212}
213
214
215template <class T, unsigned int D1, unsigned int D2, class R>
216template <class A, class R2>
217SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::operator+=(const Expr<A,T,D1,D2,R2>& rhs) {
218 // self-addition with an expression
219 PlusEquals<T,D1,D2,A,R,R2>::Evaluate(*this, rhs);
220 return *this;
221}
222
223
224//==============================================================================
225// operator-=
226//==============================================================================
227template <class T, unsigned int D1, unsigned int D2, class R>
228SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::operator-=(const T& rhs) {
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>
238SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::operator-=(const SMatrix<T,D1,D2,R2>& rhs) {
239 // self-subtraction with another matrix (any representation)
240 // use operator-= of the representation object
241 fRep -= rhs.fRep;
242 return *this;
243}
244
245
246template <class T, unsigned int D1, unsigned int D2, class R>
247template <class A, class R2>
248SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::operator-=(const Expr<A,T,D1,D2,R2>& rhs) {
249 // self-subtraction with an expression
250 MinusEquals<T,D1,D2,A,R,R2>::Evaluate(*this, rhs);
251 return *this;
252}
253
254//==============================================================================
255// operator*=
256//==============================================================================
257template <class T, unsigned int D1, unsigned int D2, class R>
258SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::operator*=(const T& rhs) {
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>
268SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::operator*=(const SMatrix<T,D1,D2,R2>& rhs) {
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}
273
274template <class T, unsigned int D1, unsigned int D2, class R>
275template <class A, class R2>
276SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::operator*=(const Expr<A,T,D1,D2,R2>& rhs) {
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>
287SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::operator/=(const T& rhs) {
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>
309bool SMatrix<T,D1,D2,R>::operator==(const SMatrix<T,D1,D2,R2>& rhs) const {
310 return fRep == rhs.fRep;
311}
312
313template <class T, unsigned int D1, unsigned int D2, class R>
314template <class A, class R2>
315bool SMatrix<T,D1,D2,R>::operator==(const Expr<A,T,D1,D2,R2>& rhs) const {
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>
332inline bool SMatrix<T,D1,D2,R>::operator!=(const SMatrix<T,D1,D2,R>& rhs) const {
333 return !operator==(rhs);
334}
335
336template <class T, unsigned int D1, unsigned int D2, class R>
337template <class A, class R2>
338inline bool SMatrix<T,D1,D2,R>::operator!=(const Expr<A,T,D1,D2,R2>& rhs) const {
339 return !operator==(rhs);
340}
341
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>
357bool SMatrix<T,D1,D2,R>::operator>(const SMatrix<T,D1,D2,R2>& rhs) const {
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;
363}
364
365template <class T, unsigned int D1, unsigned int D2, class R>
366template <class A, class R2>
367bool SMatrix<T,D1,D2,R>::operator>(const Expr<A,T,D1,D2,R2>& rhs) const {
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}
374
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>
389bool SMatrix<T,D1,D2,R>::operator<(const SMatrix<T,D1,D2,R2>& rhs) const {
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>
399bool SMatrix<T,D1,D2,R>::operator<(const Expr<A,T,D1,D2,R2>& rhs) const {
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>
412inline bool SMatrix<T,D1,D2,R>::Invert() {
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>
419inline SMatrix<T,D1,D2,R> SMatrix<T,D1,D2,R>::Inverse(int & ifail) const {
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>
429inline bool SMatrix<T,D1,D2,R>::InvertFast() {
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>
436inline SMatrix<T,D1,D2,R> SMatrix<T,D1,D2,R>::InverseFast(int & ifail) const {
437 SMatrix<T,D1,D2,R> tmp(*this);
438 bool ok = tmp.InvertFast();
439 ifail = 0;
440 if (!ok) ifail = 1;
441 return tmp;
442}
443
444// Choleski inversion for symmetric and positive defined matrices
445template <class T, unsigned int D1, unsigned int D2, class R>
446inline bool SMatrix<T,D1,D2, R>::InvertChol() {
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>
452inline SMatrix<T,D1,D2, R> SMatrix<T,D1,D2, R>::InverseChol(int &ifail) const {
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
460
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>
484SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::Place_in_row(const SVector<T,D>& rhs,
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>
501SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::Place_in_row(const VecExpr<A,T,D>& rhs,
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>
518SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::Place_in_col(const SVector<T,D>& rhs,
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;
528}
529
530//==============================================================================
531// place_in_col
532//==============================================================================
533template <class T, unsigned int D1, unsigned int D2, class R>
534template <class A, unsigned int D>
535SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::Place_in_col(const VecExpr<A,T,D>& rhs,
536 unsigned int row,
537 unsigned int col) {
538
539 assert(row+D <= D1);
540
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>
552SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::Place_at(const SMatrix<T,D3,D4,R2>& rhs,
553 unsigned int row,
554 unsigned int col) {
555 PlaceMatrix<T,D1,D2,D3,D4,R,R2>::Evaluate(*this,rhs,row,col);
556 return *this;
557}
558
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>
564SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::Place_at(const Expr<A,T,D3,D4,R2>& rhs,
565 unsigned int row,
566 unsigned int col) {
567 PlaceExpr<T,D1,D2,D3,D4,A,R,R2>::Evaluate(*this,rhs,row,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;
578
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}
598
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 << " ]";
618
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(); }
634
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
648
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>
670inline T * SMatrix<T,D1,D2,R>::begin() {
671 return fRep.Array();
672}
673
674template <class T, unsigned int D1, unsigned int D2, class R>
675inline T * SMatrix<T,D1,D2,R>::end() {
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>
752
753//diagonal
754template <class T, unsigned int D1, unsigned int D2, class R>
755SVector<T,D1> SMatrix<T,D1,D2,R>::Diagonal( ) const {
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>
784T SMatrix<T,D1,D2,R>::Trace( ) const {
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
size
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
STATIC_CHECK
#define STATIC_CHECK(expr, msg)
Definition StaticCheck.h:56
TRangeDynCast
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
Definition TCollection.h:358
N
#define N
p
winID h TVirtualViewer3D TVirtualGLPainter p
Definition TGWin32VirtualGLProxy.cxx:51
offset
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h offset
Definition TGWin32VirtualXProxy.cxx:245
operator=
void operator=(const TProof &)
kSize
@ kSize
Definition TStructNode.h:26
ROOT::Math::CholInverter::Dinv
static bool Dinv(MatrixRep &)
Definition Dinv.h:317
ROOT::Math::Determinant::Dfact
static bool Dfact(MatRepStd< T, n, idim > &rhs, T &det)
Definition Dfact.h:53
ROOT::Math::FastInverter::Dinv
static bool Dinv(MatrixRep &rhs)
Definition Dinv.h:148
ROOT::Math::Inverter::Dinv
static bool Dinv(MatrixRep &rhs)
matrix inversion for a generic square matrix using LU factorization (code originally from CERNLIB and...
Definition Dinv.h:75
ROOT::Math::SMatrix
SMatrix: a generic fixed size D1 x D2 Matrix class.
Definition SMatrix.h:101
ROOT::Math::SMatrix::Col
SVector< T, D1 > Col(unsigned int thecol) const
return a full Matrix column as a vector (copy the content in a new vector)
Definition SMatrix.icc:590
ROOT::Math::SMatrix::SMatrix
SMatrix()
Default constructor:
Definition SMatrix.icc:72
ROOT::Math::SMatrix::operator-=
SMatrix< T, D1, D2, R > & operator-=(const T &rhs)
subtraction with a scalar
Definition SMatrix.icc:228
ROOT::Math::SMatrix::apply
T apply(unsigned int i) const
access the parse tree with the index starting from zero and following the C convention for the order ...
Definition SMatrix.icc:627
ROOT::Math::SMatrix::Det2
bool Det2(T &det) const
determinant of square Matrix via Dfact.
Definition SMatrix.icc:473
ROOT::Math::SMatrix::UpperBlock
SVector< T, D1 *(D2+1)/2 > UpperBlock() const
return the upper Triangular block of the matrices (including the diagonal) as a vector of sizes N = D...
Definition SMatrix.icc:797
ROOT::Math::SMatrix::end
iterator end()
STL iterator interface.
Definition SMatrix.icc:675
ROOT::Math::SMatrix::At
const T & At(unsigned int i, unsigned int j) const
read only access to matrix element, with indices starting from 0.
Definition SMatrix.icc:653
ROOT::Math::SMatrix::Det
bool Det(T &det)
determinant of square Matrix via Dfact.
Definition SMatrix.icc:466
ROOT::Math::SMatrix::SubCol
SubVector SubCol(unsigned int thecol, unsigned int row0=0) const
return a slice of the column as a vector starting at the row value row0 until row0+Dsub.
Definition SMatrix.icc:728
ROOT::Math::SMatrix::operator>
bool operator>(const T &rhs) const
element wise comparison
Definition SMatrix.icc:347
ROOT::Math::SMatrix::Print
std::ostream & Print(std::ostream &os) const
Print: used by operator<<()
Definition SMatrix.icc:603
ROOT::Math::SMatrix::operator=
SMatrix< T, D1, D2, R > & operator=(const M &rhs)
Assign from another compatible matrix.
Definition SMatrix.icc:155
ROOT::Math::SMatrix::Place_at
SMatrix< T, D1, D2, R > & Place_at(const SMatrix< T, D3, D4, R2 > &rhs, unsigned int row, unsigned int col)
place a matrix in this matrix
Definition SMatrix.icc:552
ROOT::Math::SMatrix::Inverse
SMatrix< T, D1, D2, R > Inverse(int &ifail) const
Invert a square Matrix and returns a new matrix.
Definition SMatrix.icc:419
ROOT::Math::SMatrix::Sub
SubMatrix Sub(unsigned int row0, unsigned int col0) const
return a submatrix with the upper left corner at the values (row0, col0) and with sizes N1,...
Definition SMatrix.icc:745
ROOT::Math::SMatrix::operator*=
SMatrix< T, D1, D2, R > & operator*=(const T &rhs)
multiplication with a scalar
Definition SMatrix.icc:258
ROOT::Math::SMatrix::operator<
bool operator<(const T &rhs) const
element wise comparison
Definition SMatrix.icc:379
ROOT::Math::SMatrix::begin
iterator begin()
STL iterator interface.
Definition SMatrix.icc:670
ROOT::Math::SMatrix::SetElements
void SetElements(InputIterator begin, InputIterator end, bool triang=false, bool lower=true)
Set matrix elements with STL iterator interface.
Definition SMatrix.icc:692
ROOT::Math::SMatrix::InverseFast
SMatrix< T, D1, D2, R > InverseFast(int &ifail) const
Invert a square Matrix and returns a new matrix.
Definition SMatrix.icc:436
ROOT::Math::SMatrix::operator==
bool operator==(const T &rhs) const
element wise comparison
Definition SMatrix.icc:299
ROOT::Math::SMatrix::Place_in_row
SMatrix< T, D1, D2, R > & Place_in_row(const SVector< T, D > &rhs, unsigned int row, unsigned int col)
place a vector in a Matrix row
Definition SMatrix.icc:484
ROOT::Math::SMatrix::Diagonal
SVector< T, D1 > Diagonal() const
return diagonal elements of a matrix as a Vector.
Definition SMatrix.icc:755
ROOT::Math::SMatrix::InverseChol
SMatrix< T, D1, D2, R > InverseChol(int &ifail) const
Invert of a symmetric positive defined Matrix using Choleski decomposition.
Definition SMatrix.icc:452
ROOT::Math::SMatrix::SetDiagonal
void SetDiagonal(const Vector &v)
Set the diagonal elements from a Vector Require that vector implements SVector::kSize since a check (...
Definition SMatrix.icc:770
ROOT::Math::SMatrix::operator!=
bool operator!=(const T &rhs) const
element wise comparison
Definition SMatrix.icc:327
ROOT::Math::SMatrix::operator+=
SMatrix< T, D1, D2, R > & operator+=(const T &rhs)
addition with a scalar
Definition SMatrix.icc:197
ROOT::Math::SMatrix::InvertChol
bool InvertChol()
Inversion of a symmetric positive defined Matrix using Choleski decomposition.
Definition SMatrix.icc:446
ROOT::Math::SMatrix::Invert
bool Invert()
Invert a square Matrix ( this method changes the current matrix).
Definition SMatrix.icc:412
ROOT::Math::SMatrix::IsInUse
bool IsInUse(const T *p) const
Function to check if a matrix is sharing same memory location of the passed pointer This function is ...
Definition SMatrix.icc:895
ROOT::Math::SMatrix::operator/=
SMatrix< T, D1, D2, R > & operator/=(const T &rhs)
division with a scalar
Definition SMatrix.icc:287
ROOT::Math::SMatrix::InvertFast
bool InvertFast()
Fast inversion of a square Matrix ( this method changes the current matrix).
Definition SMatrix.icc:429
ROOT::Math::SMatrix::LowerBlock
SVector< T, D1 *(D2+1)/2 > LowerBlock() const
return the lower Triangular block of the matrices (including the diagonal) as a vector of sizes N = D...
Definition SMatrix.icc:826
ROOT::Math::SMatrix::Trace
T Trace() const
return the trace of a matrix Sum of the diagonal elements
Definition SMatrix.icc:784
ROOT::Math::SMatrix::SubRow
SubVector SubRow(unsigned int therow, unsigned int col0=0) const
return a slice of therow as a vector starting at the column value col0 until col0+N,...
Definition SMatrix.icc:712
ROOT::Math::SMatrix::operator()
const T & operator()(unsigned int i, unsigned int j) const
read only access to matrix element, with indices starting from 0
Definition SMatrix.icc:639
ROOT::Math::SMatrix::Row
SVector< T, D2 > Row(unsigned int therow) const
return a full Matrix row as a vector (copy the content in a new vector)
Definition SMatrix.icc:575
ROOT::Math::SMatrix::Array
const T * Array() const
return read-only pointer to internal array
Definition SMatrix.icc:630
ROOT::Math::SMatrix::Place_in_col
SMatrix< T, D1, D2, R > & Place_in_col(const SVector< T, D > &rhs, unsigned int row, unsigned int col)
place a vector in a Matrix column
Definition SMatrix.icc:518
ROOT::Math::SVector
SVector: a generic fixed size Vector class.
Definition SVector.h:75
Math
Namespace for new Math classes and functions.
ROOT::Internal::operator==
bool operator==(const RConcurrentHashColl::HashValue &lhs, const RConcurrentHashColl::HashValue &rhs)
Definition RConcurrentHashColl.hxx:70
ROOT
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
Definition EExecutionPolicy.hxx:4
v
@ v
Definition rootcling_impl.cxx:3699
ROOT::Math::AssignItr::Evaluate
static void Evaluate(SMatrix< T, D1, D2, R > &lhs, Iterator begin, Iterator end, bool triang, bool lower, bool check=true)
Definition HelperOps.h:519
ROOT::Math::Assign::Evaluate
static void Evaluate(SMatrix< T, D1, D2, R1 > &lhs, const Expr< A, T, D1, D2, R2 > &rhs)
Evaluate the expression from general to general matrices.
Definition HelperOps.h:55
ROOT::Math::MinusEquals::Evaluate
static void Evaluate(SMatrix< T, D1, D2, R1 > &lhs, const Expr< A, T, D1, D2, R2 > &rhs)
Definition HelperOps.h:281
ROOT::Math::PlaceExpr::Evaluate
static void Evaluate(SMatrix< T, D1, D2, R1 > &lhs, const Expr< A, T, D3, D4, R2 > &rhs, unsigned int row, unsigned int col)
Definition HelperOps.h:380
ROOT::Math::PlaceMatrix::Evaluate
static void Evaluate(SMatrix< T, D1, D2, R1 > &lhs, const SMatrix< T, D3, D4, R2 > &rhs, unsigned int row, unsigned int col)
Definition HelperOps.h:363
ROOT::Math::PlusEquals::Evaluate
static void Evaluate(SMatrix< T, D1, D2, R1 > &lhs, const Expr< A, T, D1, D2, R2 > &rhs)
Definition HelperOps.h:197
ROOT::Math::RetrieveMatrix::Evaluate
static void Evaluate(SMatrix< T, D1, D2, R1 > &lhs, const SMatrix< T, D3, D4, R2 > &rhs, unsigned int row, unsigned int col)
Definition HelperOps.h:463