Logo ROOT   master
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 
62 namespace ROOT {
63 
64 namespace Math {
65 
66 
67 
68 //==============================================================================
69 // Constructors
70 //==============================================================================
71 template <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
78 template <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 
92 template <class T, unsigned int D1, unsigned int D2, class R>
94  fRep = rhs.fRep;
95 }
96 
97 
98 template <class T, unsigned int D1, unsigned int D2, class R>
99 template <class R2>
101  operator=(rhs);
102 }
103 
104 
105 template <class T, unsigned int D1, unsigned int D2, class R>
106 template <class A, class R2>
108  operator=(rhs);
109 }
110 
111 
112 //=============================================================================
113 // New Constructors from STL interfaces
114 //=============================================================================
115 template <class T, unsigned int D1, unsigned int D2, class R>
116 template <class InputIterator>
117 SMatrix<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 
123 template <class T, unsigned int D1, unsigned int D2, class R>
124 template <class InputIterator>
125 SMatrix<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 //==============================================================================
137 template <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 
143 template <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 
153 template <class T, unsigned int D1, unsigned int D2, class R>
154 template <class M>
156  fRep = rhs.fRep;
157  return *this;
158 }
159 
160 template <class T, unsigned int D1, unsigned int D2, class R>
162  fRep = rhs.fRep;
163  return *this;
164 }
165 
166 template <class T, unsigned int D1, unsigned int D2, class R>
167 template <class A, class R2>
169 
171  return *this;
172 }
173 
174 //=============================================================================
175 // assign from an identity
176 template <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 //=============================================================================
196 template <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 
205 template <class T, unsigned int D1, unsigned int D2, class R>
206 template <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 
214 
215 template <class T, unsigned int D1, unsigned int D2, class R>
216 template <class A, class R2>
218  // self-addition with an expression
220  return *this;
221 }
222 
223 
224 //==============================================================================
225 // operator-=
226 //==============================================================================
227 template <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 
236 template <class T, unsigned int D1, unsigned int D2, class R>
237 template <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 }
244 
245 
246 template <class T, unsigned int D1, unsigned int D2, class R>
247 template <class A, class R2>
249  // self-subtraction with an expression
251  return *this;
252 }
253 
254 //==============================================================================
255 // operator*=
256 //==============================================================================
257 template <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 
266 template <class T, unsigned int D1, unsigned int D2, class R>
267 template <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 }
273 
274 template <class T, unsigned int D1, unsigned int D2, class R>
275 template <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 //==============================================================================
286 template <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 //==============================================================================
298 template <class T, unsigned int D1, unsigned int D2, class R>
299 bool 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 
307 template <class T, unsigned int D1, unsigned int D2, class R>
308 template <class R2>
310  return fRep == rhs.fRep;
311 }
312 
313 template <class T, unsigned int D1, unsigned int D2, class R>
314 template <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 //==============================================================================
326 template <class T, unsigned int D1, unsigned int D2, class R>
327 inline bool SMatrix<T,D1,D2,R>::operator!=(const T& rhs) const {
328  return !operator==(rhs);
329 }
330 
331 template <class T, unsigned int D1, unsigned int D2, class R>
332 inline bool SMatrix<T,D1,D2,R>::operator!=(const SMatrix<T,D1,D2,R>& rhs) const {
333  return !operator==(rhs);
334 }
335 
336 template <class T, unsigned int D1, unsigned int D2, class R>
337 template <class A, class R2>
338 inline 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 //==============================================================================
346 template <class T, unsigned int D1, unsigned int D2, class R>
347 bool 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 
355 template <class T, unsigned int D1, unsigned int D2, class R>
356 template <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;
363 }
364 
365 template <class T, unsigned int D1, unsigned int D2, class R>
366 template <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 }
374 
375 //==============================================================================
376 // operator< (element wise comparison)
377 //==============================================================================
378 template <class T, unsigned int D1, unsigned int D2, class R>
379 bool 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 
387 template <class T, unsigned int D1, unsigned int D2, class R>
388 template <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 
397 template <class T, unsigned int D1, unsigned int D2, class R>
398 template <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 //==============================================================================
411 template <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
418 template <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
428 template <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
435 template <class T, unsigned int D1, unsigned int D2, class R>
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
445 template <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 
451 template <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 
460 
461 
462 //==============================================================================
463 // det
464 //==============================================================================
465 template <class T, unsigned int D1, unsigned int D2, class R>
466 inline 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 }
472 template <class T, unsigned int D1, unsigned int D2, class R>
473 inline 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 //==============================================================================
482 template <class T, unsigned int D1, unsigned int D2, class R>
483 template <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 //==============================================================================
499 template <class T, unsigned int D1, unsigned int D2, class R>
500 template <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 //==============================================================================
516 template <class T, unsigned int D1, unsigned int D2, class R>
517 template <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;
528 }
529 
530 //==============================================================================
531 // place_in_col
532 //==============================================================================
533 template <class T, unsigned int D1, unsigned int D2, class R>
534 template <class A, unsigned int D>
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 //==============================================================================
550 template <class T, unsigned int D1, unsigned int D2, class R>
551 template <unsigned int D3, unsigned int D4, class R2>
553  unsigned int row,
554  unsigned int col) {
556  return *this;
557 }
558 
559 //==============================================================================
560 // place_at
561 //==============================================================================
562 template <class T, unsigned int D1, unsigned int D2, class R>
563 template <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 //==============================================================================
574 template <class T, unsigned int D1, unsigned int D2, class R>
575 SVector<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 //==============================================================================
589 template <class T, unsigned int D1, unsigned int D2, class R>
590 SVector<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 //==============================================================================
602 template <class T, unsigned int D1, unsigned int D2, class R>
603 std::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 //==============================================================================
626 template <class T, unsigned int D1, unsigned int D2, class R>
627 inline T SMatrix<T,D1,D2,R>::apply(unsigned int i) const { return fRep[i]; }
628 
629 template <class T, unsigned int D1, unsigned int D2, class R>
630 inline const T* SMatrix<T,D1,D2,R>::Array() const { return fRep.Array(); }
631 
632 template <class T, unsigned int D1, unsigned int D2, class R>
633 inline T* SMatrix<T,D1,D2,R>::Array() { return fRep.Array(); }
634 
635 //==============================================================================
636 // Operators
637 //==============================================================================
638 template <class T, unsigned int D1, unsigned int D2, class R>
639 inline const T& SMatrix<T,D1,D2,R>::operator()(unsigned int i, unsigned int j) const {
640  return fRep(i,j);
641 }
642 
643 template <class T, unsigned int D1, unsigned int D2, class R>
644 inline 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 //==============================================================================
652 template <class T, unsigned int D1, unsigned int D2, class R>
653 inline 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 
659 template <class T, unsigned int D1, unsigned int D2, class R>
660 inline 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 //==============================================================================
669 template <class T, unsigned int D1, unsigned int D2, class R>
671  return fRep.Array();
672 }
673 
674 template <class T, unsigned int D1, unsigned int D2, class R>
676  return fRep.Array() + R::kSize;
677 }
678 
679 template <class T, unsigned int D1, unsigned int D2, class R>
680 inline const T * SMatrix<T,D1,D2,R>::begin() const {
681  return fRep.Array();
682 }
683 
684 template <class T, unsigned int D1, unsigned int D2, class R>
685 inline const T * SMatrix<T,D1,D2,R>::end() const {
686  return fRep.Array() + R::kSize;
687 }
688 
689 
690 template <class T, unsigned int D1, unsigned int D2, class R>
691 template <class InputIterator>
692 void 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 
697 template <class T, unsigned int D1, unsigned int D2, class R>
698 template <class InputIterator>
699 void 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 //==============================================================================
710 template <class T, unsigned int D1, unsigned int D2, class R>
711 template <class SubVector>
712 SubVector 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 
726 template <class T, unsigned int D1, unsigned int D2, class R>
727 template <class SubVector>
728 SubVector 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
743 template <class T, unsigned int D1, unsigned int D2, class R>
744 template <class SubMatrix>
745 SubMatrix 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
754 template <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
768 template <class T, unsigned int D1, unsigned int D2, class R>
769 template <class Vector>
770 void 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
783 template <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
795 template <class T, unsigned int D1, unsigned int D2, class R>
796 #ifndef UNSUPPORTED_TEMPLATE_EXPRESSION
797 SVector<T, D1 * (D2 +1)/2 > SMatrix<T,D1,D2,R>::UpperBlock( ) const {
798 #else
799 template <class SubVector>
800 SubVector 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
824 template <class T, unsigned int D1, unsigned int D2, class R>
825 #ifndef UNSUPPORTED_TEMPLATE_EXPRESSION
826 SVector<T, D1 * (D2 +1)/2 > SMatrix<T,D1,D2,R>::LowerBlock( ) const {
827 #else
828 template <class SubVector>
829 SubVector 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
856 template <class T, unsigned int D1, unsigned int D2, class R>
857 #ifndef UNSUPPORTED_TEMPLATE_EXPRESSION
858 SMatrix<T,D1,D2,R>::SMatrix(const SVector<T, D1*(D2+1)/2 > & v, bool lower ) {
859 #else
860 template <unsigned int N>
861 SMatrix<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 
894 template <class T, unsigned int D1, unsigned int D2, class R>
895 bool 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
SMatrix< T, D1, D2, R > InverseChol(int &ifail) const
Invert of a symmetric positive defined Matrix using Choleski decomposition.
Definition: SMatrix.icc:452
SMatrix< T, D1, D2, R > & operator+=(const T &rhs)
addition with a scalar
Definition: SMatrix.icc:197
SMatrix< T, D1, D2, R > & operator*=(const T &rhs)
multiplication with a scalar
Definition: SMatrix.icc:258
void SetDiagonal(const Vector &v)
Set the diagonal elements from a Vector Require that vector implements kSize since a check (staticall...
Definition: SMatrix.icc:770
SMatrix< T, D1, D2, R > & operator/=(const T &rhs)
division with a scalar
Definition: SMatrix.icc:287
bool operator!=(const T &rhs) const
element wise comparison
Definition: SMatrix.icc:327
Returns the available number of logical cores.
Definition: RNumpyDS.hxx:30
bool operator<(const T &rhs) const
element wise comparison
Definition: SMatrix.icc:379
static bool Dinv(MatrixRep &)
Definition: Dinv.h:317
double T(double x)
Definition: ChebyshevPol.h:34
SMatrix< T, D1, D2, R > & operator-=(const T &rhs)
subtraction with a scalar
Definition: SMatrix.icc:228
static bool Dinv(MatrixRep &rhs)
Definition: Dinv.h:148
static void Evaluate(SMatrix< T, D1, D2, R1 > &lhs, const Expr< A, T, D1, D2, R2 > &rhs)
Definition: HelperOps.h:197
SMatrix< T, D1, D2, R > & operator=(const M &rhs)
Assign from another compatible matrix.
Definition: SMatrix.icc:155
const T * Array() const
return read-only pointer to internal array
Definition: SMatrix.icc:630
#define N
bool Invert()
Invert a square Matrix ( this method changes the current matrix).
Definition: SMatrix.icc:412
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
bool InvertFast()
Fast Invertion of a square Matrix ( this method changes the current matrix).
Definition: SMatrix.icc:429
bool InvertChol()
Invertion of a symmetric positive defined Matrix using Choleski decomposition.
Definition: SMatrix.icc:446
#define STATIC_CHECK(expr, msg)
Definition: StaticCheck.h:56
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
SMatrix< T, D1, D2, R > InverseFast(int &ifail) const
Invert a square Matrix and returns a new matrix.
Definition: SMatrix.icc:436
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
SubVector SubRow(unsigned int therow, unsigned int col0=0) const
return a slice of therow as a vector starting at the colum value col0 until col0+N, where N is the size of the vector (SubVector::kSize ) Condition col0+N <= D2
Definition: SMatrix.icc:712
bool Det(T &det)
determinant of square Matrix via Dfact.
Definition: SMatrix.icc:466
SMatrix: a generic fixed size D1 x D2 Matrix class.
SMatrix()
Default constructor:
Definition: SMatrix.icc:72
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
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
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
R fRep
Matrix Storage Object containing matrix data.
Definition: SMatrix.h:709
std::ostream & Print(std::ostream &os) const
Print: used by operator<<()
Definition: SMatrix.icc:603
iterator begin()
STL iterator interface.
Definition: SMatrix.icc:670
SVector< T, D1 > Diagonal() const
return diagonal elements of a matrix as a Vector.
Definition: SMatrix.icc:755
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
T apply(unsigned int i) const
access the parse tree. Index starts from zero
Definition: SVector.icc:537
bool Det2(T &det) const
determinant of square Matrix via Dfact.
Definition: SMatrix.icc:473
SMatrix< T, D1, D2, R > Inverse(int &ifail) const
Invert a square Matrix and returns a new matrix.
Definition: SMatrix.icc:419
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
bool operator>(const T &rhs) const
element wise comparison
Definition: SMatrix.icc:347
static void Evaluate(SMatrix< T, D1, D2, R1 > &lhs, const Expr< A, T, D1, D2, R2 > &rhs)
Definition: HelperOps.h:281
static void Evaluate(SMatrix< T, D1, D2, R > &lhs, Iterator begin, Iterator end, bool triang, bool lower, bool check=true)
Definition: HelperOps.h:519
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
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
T Trace() const
return the trace of a matrix Sum of the diagonal elements
Definition: SMatrix.icc:784
static bool Dinv(MatrixRep &rhs)
matrix inversion for a generic square matrix using LU factorization (code originally from CERNLIB and...
Definition: Dinv.h:75
Namespace for new Math classes and functions.
Binding & operator=(OUT(*fun)(void))
bool operator==(const T &rhs) const
element wise comparison
Definition: SMatrix.icc:299
Expression wrapper class for Vector objects.
Definition: Expression.h:64
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
Bool_t operator==(const TDatime &d1, const TDatime &d2)
Definition: TDatime.h:102
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
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
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
static bool Dfact(MatRepStd< T, n, idim > &rhs, T &det)
Definition: Dfact.h:53
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
void SetElements(InputIterator begin, InputIterator end, bool triang=false, bool lower=true)
Set matrix elements with STL iterator interface.
Definition: SMatrix.icc:692
T apply(unsigned int i) const
Definition: Expression.h:150
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
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
SVector: a generic fixed size Vector class.
T apply(unsigned int i) const
Definition: Expression.h:77
iterator end()
STL iterator interface.
Definition: SMatrix.icc:675