Logo ROOT   6.13/01
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>
161 template <class A, class R2>
163 
165  return *this;
166 }
167 
168 //=============================================================================
169 // assign from an identity
170 template <class T, unsigned int D1, unsigned int D2, class R>
172  for(unsigned int i=0; i<R::kSize; ++i)
173  fRep.Array()[i] = 0;
174  if (D1 <= D2) {
175  for(unsigned int i=0; i<D1; ++i)
176  fRep[i*D2+i] = 1;
177  }
178  else {
179  for(unsigned int i=0; i<D2; ++i)
180  fRep[i*D2+i] = 1;
181  }
182  return *this;
183 }
184 
185 
186 
187 //=============================================================================
188 // operator+=
189 //=============================================================================
190 template <class T, unsigned int D1, unsigned int D2, class R>
192  // self-addition with a scalar value
193  for(unsigned int i=0; i<R::kSize; ++i) {
194  fRep.Array()[i] += rhs;
195  }
196  return *this;
197 }
198 
199 template <class T, unsigned int D1, unsigned int D2, class R>
200 template <class R2>
202  // self-addition with another matrix (any representation)
203  // use operator+= of the representation object
204  fRep += rhs.fRep;
205  return *this;
206 }
207 
208 
209 template <class T, unsigned int D1, unsigned int D2, class R>
210 template <class A, class R2>
212  // self-addition with an expression
214  return *this;
215 }
216 
217 
218 //==============================================================================
219 // operator-=
220 //==============================================================================
221 template <class T, unsigned int D1, unsigned int D2, class R>
223  // self-subtraction with a scalar value
224  for(unsigned int i=0; i<R::kSize; ++i) {
225  fRep.Array()[i] -= rhs;
226  }
227  return *this;
228 }
229 
230 template <class T, unsigned int D1, unsigned int D2, class R>
231 template <class R2>
233  // self-subtraction with another matrix (any representation)
234  // use operator-= of the representation object
235  fRep -= rhs.fRep;
236  return *this;
237 }
238 
239 
240 template <class T, unsigned int D1, unsigned int D2, class R>
241 template <class A, class R2>
243  // self-subtraction with an expression
245  return *this;
246 }
247 
248 //==============================================================================
249 // operator*=
250 //==============================================================================
251 template <class T, unsigned int D1, unsigned int D2, class R>
253  // case of multiplication with a scalar
254  for(unsigned int i=0; i<R::kSize; ++i) {
255  fRep.Array()[i] *= rhs;
256  }
257  return *this;
258 }
259 
260 template <class T, unsigned int D1, unsigned int D2, class R>
261 template <class R2>
263  // self-multiplication with another matrix (will work only for square matrices)
264  // a temporary is needed and will be created automatically to store intermediate result
265  return operator=(*this * rhs);
266 }
267 
268 template <class T, unsigned int D1, unsigned int D2, class R>
269 template <class A, class R2>
271  // self-multiplication with an expression (will work only for square matrices)
272  // a temporary is needed and will be created automatically to store intermediate result
273  return operator=(*this * rhs);
274 }
275 
276 
277 //==============================================================================
278 // operator/= (only for scalar values)
279 //==============================================================================
280 template <class T, unsigned int D1, unsigned int D2, class R>
282  // division with a scalar
283  for(unsigned int i=0; i<R::kSize; ++i) {
284  fRep.Array()[i] /= rhs;
285  }
286  return *this;
287 }
288 
289 //==============================================================================
290 // operator== (element wise comparison)
291 //==============================================================================
292 template <class T, unsigned int D1, unsigned int D2, class R>
293 bool SMatrix<T,D1,D2,R>::operator==(const T& rhs) const {
294  bool rc = true;
295  for(unsigned int i=0; i<R::kSize; ++i) {
296  rc = rc && (fRep.Array()[i] == rhs);
297  }
298  return rc;
299 }
300 
301 template <class T, unsigned int D1, unsigned int D2, class R>
302 template <class R2>
304  return fRep == rhs.fRep;
305 }
306 
307 template <class T, unsigned int D1, unsigned int D2, class R>
308 template <class A, class R2>
310  bool rc = true;
311  for(unsigned int i=0; i<D1*D2; ++i) {
312  rc = rc && (fRep[i] == rhs.apply(i));
313  }
314  return rc;
315 }
316 
317 //==============================================================================
318 // operator!= (element wise comparison)
319 //==============================================================================
320 template <class T, unsigned int D1, unsigned int D2, class R>
321 inline bool SMatrix<T,D1,D2,R>::operator!=(const T& rhs) const {
322  return !operator==(rhs);
323 }
324 
325 template <class T, unsigned int D1, unsigned int D2, class R>
326 inline bool SMatrix<T,D1,D2,R>::operator!=(const SMatrix<T,D1,D2,R>& rhs) const {
327  return !operator==(rhs);
328 }
329 
330 template <class T, unsigned int D1, unsigned int D2, class R>
331 template <class A, class R2>
332 inline bool SMatrix<T,D1,D2,R>::operator!=(const Expr<A,T,D1,D2,R2>& rhs) const {
333  return !operator==(rhs);
334 }
335 
336 
337 //==============================================================================
338 // operator> (element wise comparison)
339 //==============================================================================
340 template <class T, unsigned int D1, unsigned int D2, class R>
341 bool SMatrix<T,D1,D2,R>::operator>(const T& rhs) const {
342  bool rc = true;
343  for(unsigned int i=0; i<D1*D2; ++i) {
344  rc = rc && (fRep[i] > rhs);
345  }
346  return rc;
347 }
348 
349 template <class T, unsigned int D1, unsigned int D2, class R>
350 template <class R2>
352  bool rc = true;
353  for(unsigned int i=0; i<D1*D2; ++i) {
354  rc = rc && (fRep[i] > rhs.fRep[i]);
355  }
356  return rc;
357 }
358 
359 template <class T, unsigned int D1, unsigned int D2, class R>
360 template <class A, class R2>
362  bool rc = true;
363  for(unsigned int i=0; i<D1*D2; ++i) {
364  rc = rc && (fRep[i] > rhs.apply(i));
365  }
366  return rc;
367 }
368 
369 //==============================================================================
370 // operator< (element wise comparison)
371 //==============================================================================
372 template <class T, unsigned int D1, unsigned int D2, class R>
373 bool SMatrix<T,D1,D2,R>::operator<(const T& rhs) const {
374  bool rc = true;
375  for(unsigned int i=0; i<D1*D2; ++i) {
376  rc = rc && (fRep[i] < rhs);
377  }
378  return rc;
379 }
380 
381 template <class T, unsigned int D1, unsigned int D2, class R>
382 template <class R2>
384  bool rc = true;
385  for(unsigned int i=0; i<D1*D2; ++i) {
386  rc = rc && (fRep[i] < rhs.fRep[i]);
387  }
388  return rc;
389 }
390 
391 template <class T, unsigned int D1, unsigned int D2, class R>
392 template <class A, class R2>
394  bool rc = true;
395  for(unsigned int i=0; i<D1*D2; ++i) {
396  rc = rc && (fRep[i] < rhs.apply(i));
397  }
398  return rc;
399 }
400 
401 
402 //==============================================================================
403 // invert
404 //==============================================================================
405 template <class T, unsigned int D1, unsigned int D2, class R>
407  STATIC_CHECK( D1 == D2,SMatrix_not_square);
408  return Inverter<D1,D2>::Dinv((*this).fRep);
409 }
410 
411 // invert returning a new matrix
412 template <class T, unsigned int D1, unsigned int D2, class R>
414  SMatrix<T,D1,D2,R> tmp(*this);
415  bool ok = tmp.Invert();
416  ifail = 0;
417  if (!ok) ifail = 1;
418  return tmp;
419 }
420 
421 // fast inversion
422 template <class T, unsigned int D1, unsigned int D2, class R>
424  STATIC_CHECK( D1 == D2,SMatrix_not_square);
425  return FastInverter<D1,D2>::Dinv((*this).fRep);
426 }
427 
428 // fast inversion returning a new matrix
429 template <class T, unsigned int D1, unsigned int D2, class R>
431  SMatrix<T,D1,D2,R> tmp(*this);
432  bool ok = tmp.InvertFast();
433  ifail = 0;
434  if (!ok) ifail = 1;
435  return tmp;
436 }
437 
438 // Choleski inversion for symmetric and positive defined matrices
439 template <class T, unsigned int D1, unsigned int D2, class R>
441  STATIC_CHECK( D1 == D2,SMatrix_not_square);
442  return CholInverter<D1>::Dinv((*this).fRep);
443 }
444 
445 template <class T, unsigned int D1, unsigned int D2, class R>
447  SMatrix<T,D1,D2,R > tmp(*this);
448  bool ok = tmp.InvertChol();
449  ifail = 0;
450  if (!ok) ifail = 1;
451  return tmp;
452 }
453 
454 
455 
456 //==============================================================================
457 // det
458 //==============================================================================
459 template <class T, unsigned int D1, unsigned int D2, class R>
460 inline bool SMatrix<T,D1,D2,R>::Det(T& det) {
461  STATIC_CHECK( D1 == D2,SMatrix_not_square);
462  // return Dfact<SMatrix<T,D1,D1>, D1, D1>(*this,det);
463  //return Dfact<R, D1, D1>(fRep, det);
464  return Determinant<D1,D2>::Dfact(fRep, det);
465 }
466 template <class T, unsigned int D1, unsigned int D2, class R>
467 inline bool SMatrix<T,D1,D2,R>::Det2(T& det) const {
468  SMatrix<T,D1,D2,R> tmp(*this);
469  return tmp.Det(det);
470 }
471 
472 
473 //==============================================================================
474 // place_in_row
475 //==============================================================================
476 template <class T, unsigned int D1, unsigned int D2, class R>
477 template <unsigned int D>
479  unsigned int row,
480  unsigned int col) {
481 
482  assert(col+D <= D2);
483 
484  for(unsigned int i=row*D2+col, j=0; j<D; ++i, ++j) {
485  fRep[i] = rhs.apply(j);
486  }
487  return *this;
488 }
489 
490 //==============================================================================
491 // place_in_row
492 //==============================================================================
493 template <class T, unsigned int D1, unsigned int D2, class R>
494 template <class A, unsigned int D>
496  unsigned int row,
497  unsigned int col) {
498 
499  assert(col+D <= D2);
500 
501  for(unsigned int i=row*D2+col, j=0; j<D; ++i, ++j) {
502  fRep[i] = rhs.apply(j);
503  }
504  return *this;
505 }
506 
507 //==============================================================================
508 // place_in_col
509 //==============================================================================
510 template <class T, unsigned int D1, unsigned int D2, class R>
511 template <unsigned int D>
513  unsigned int row,
514  unsigned int col) {
515 
516  assert(row+D <= D1);
517 
518  for(unsigned int i=row*D2+col, j=0; j<D; i+=D2, ++j) {
519  fRep[i] = rhs.apply(j);
520  }
521  return *this;
522 }
523 
524 //==============================================================================
525 // place_in_col
526 //==============================================================================
527 template <class T, unsigned int D1, unsigned int D2, class R>
528 template <class A, unsigned int D>
530  unsigned int row,
531  unsigned int col) {
532 
533  assert(row+D <= D1);
534 
535  for(unsigned int i=row*D2+col, j=0; j<D; i+=D2, ++j) {
536  fRep[i] = rhs.apply(j);
537  }
538  return *this;
539 }
540 
541 //==============================================================================
542 // place_at
543 //==============================================================================
544 template <class T, unsigned int D1, unsigned int D2, class R>
545 template <unsigned int D3, unsigned int D4, class R2>
547  unsigned int row,
548  unsigned int col) {
550  return *this;
551 }
552 
553 //==============================================================================
554 // place_at
555 //==============================================================================
556 template <class T, unsigned int D1, unsigned int D2, class R>
557 template <class A, unsigned int D3, unsigned int D4, class R2>
559  unsigned int row,
560  unsigned int col) {
562  return *this;
563 }
564 
565 //==============================================================================
566 // row
567 //==============================================================================
568 template <class T, unsigned int D1, unsigned int D2, class R>
569 SVector<T,D2> SMatrix<T,D1,D2,R>::Row(unsigned int therow) const {
570 
571  const unsigned int offset = therow*D2;
572 
573  /*static*/ SVector<T,D2> tmp;
574  for(unsigned int i=0; i<D2; ++i) {
575  tmp[i] = fRep[offset+i];
576  }
577  return tmp;
578 }
579 
580 //==============================================================================
581 // col
582 //==============================================================================
583 template <class T, unsigned int D1, unsigned int D2, class R>
584 SVector<T,D1> SMatrix<T,D1,D2,R>::Col(unsigned int thecol) const {
585 
586  /*static */ SVector<T,D1> tmp;
587  for(unsigned int i=0; i<D1; ++i) {
588  tmp[i] = fRep[thecol+i*D2];
589  }
590  return tmp;
591 }
592 
593 //==============================================================================
594 // print
595 //==============================================================================
596 template <class T, unsigned int D1, unsigned int D2, class R>
597 std::ostream& SMatrix<T,D1,D2,R>::Print(std::ostream& os) const {
598  const std::ios_base::fmtflags prevFmt = os.setf(std::ios::right,std::ios::adjustfield);
599  // os.setf(ios::fixed);
600 
601  os << "[ ";
602  for (unsigned int i=0; i < D1; ++i) {
603  for (unsigned int j=0; j < D2; ++j) {
604  os << std::setw(12) << fRep[i*D2+j];
605  if ((!((j+1)%12)) && (j < D2-1))
606  os << std::endl << " ...";
607  }
608  if (i != D1 - 1)
609  os << std::endl << " ";
610  }
611  os << " ]";
612 
613  if (prevFmt != os.flags() ) os.setf(prevFmt, std::ios::adjustfield);
614  return os;
615 }
616 
617 //==============================================================================
618 // Access functions
619 //==============================================================================
620 template <class T, unsigned int D1, unsigned int D2, class R>
621 inline T SMatrix<T,D1,D2,R>::apply(unsigned int i) const { return fRep[i]; }
622 
623 template <class T, unsigned int D1, unsigned int D2, class R>
624 inline const T* SMatrix<T,D1,D2,R>::Array() const { return fRep.Array(); }
625 
626 template <class T, unsigned int D1, unsigned int D2, class R>
627 inline T* SMatrix<T,D1,D2,R>::Array() { return fRep.Array(); }
628 
629 //==============================================================================
630 // Operators
631 //==============================================================================
632 template <class T, unsigned int D1, unsigned int D2, class R>
633 inline const T& SMatrix<T,D1,D2,R>::operator()(unsigned int i, unsigned int j) const {
634  return fRep(i,j);
635 }
636 
637 template <class T, unsigned int D1, unsigned int D2, class R>
638 inline T& SMatrix<T,D1,D2,R>::operator()(unsigned int i, unsigned int j) {
639  return fRep(i,j);
640 }
641 
642 
643 //==============================================================================
644 // Element access with At()
645 //==============================================================================
646 template <class T, unsigned int D1, unsigned int D2, class R>
647 inline const T& SMatrix<T,D1,D2,R>::At(unsigned int i, unsigned int j) const {
648  assert(i < D1);
649  assert(j < D2);
650  return fRep(i,j);
651 }
652 
653 template <class T, unsigned int D1, unsigned int D2, class R>
654 inline T& SMatrix<T,D1,D2,R>::At(unsigned int i, unsigned int j) {
655  assert(i < D1);
656  assert(j < D2);
657  return fRep(i,j);
658 }
659 
660 //==============================================================================
661 // STL interface
662 //==============================================================================
663 template <class T, unsigned int D1, unsigned int D2, class R>
665  return fRep.Array();
666 }
667 
668 template <class T, unsigned int D1, unsigned int D2, class R>
670  return fRep.Array() + R::kSize;
671 }
672 
673 template <class T, unsigned int D1, unsigned int D2, class R>
674 inline const T * SMatrix<T,D1,D2,R>::begin() const {
675  return fRep.Array();
676 }
677 
678 template <class T, unsigned int D1, unsigned int D2, class R>
679 inline const T * SMatrix<T,D1,D2,R>::end() const {
680  return fRep.Array() + R::kSize;
681 }
682 
683 
684 template <class T, unsigned int D1, unsigned int D2, class R>
685 template <class InputIterator>
686 void SMatrix<T,D1,D2,R>::SetElements(InputIterator ibegin, InputIterator iend, bool triang, bool lower) {
687  // assume iterator size == matrix size when filling full matrix
688  AssignItr<T,D1,D2,R>::Evaluate(*this,ibegin,iend,triang,lower);
689 }
690 
691 template <class T, unsigned int D1, unsigned int D2, class R>
692 template <class InputIterator>
693 void SMatrix<T,D1,D2,R>::SetElements(InputIterator ibegin, unsigned int size, bool triang, bool lower) {
694  // assume iterator size <= matrix size (no check to be done in AssignItr)
695  assert( size <= R::kSize);
696  AssignItr<T,D1,D2,R>::Evaluate(*this,ibegin,ibegin+size,triang,lower,false);
697 }
698 
699 
700 
701 //==============================================================================
702 // SubMatrices and slices of columns and rows
703 //==============================================================================
704 template <class T, unsigned int D1, unsigned int D2, class R>
705 template <class SubVector>
706 SubVector SMatrix<T,D1,D2,R>::SubRow(unsigned int therow, unsigned int col0 ) const {
707 
708  const unsigned int offset = therow*D2 + col0;
709 
710  STATIC_CHECK( SubVector::kSize <= D2,SVector_dimension_too_small);
711  assert(col0 + SubVector::kSize <= D2);
712 
713  SubVector tmp;
714  for(unsigned int i=0; i<SubVector::kSize; ++i) {
715  tmp[i] = fRep[offset+i];
716  }
717  return tmp;
718 }
719 
720 template <class T, unsigned int D1, unsigned int D2, class R>
721 template <class SubVector>
722 SubVector SMatrix<T,D1,D2,R>::SubCol(unsigned int thecol, unsigned int row0 ) const {
723 
724  const unsigned int offset = thecol + row0*D1;
725 
726  STATIC_CHECK( SubVector::kSize <= D1,SVector_dimension_too_small);
727  assert(row0 + SubVector::kSize <= D1);
728 
729  SubVector tmp;
730  for(unsigned int i=0; i<SubVector::kSize; ++i) {
731  tmp[i] = fRep[offset+i*D1];
732  }
733  return tmp;
734 }
735 
736 // submatrix
737 template <class T, unsigned int D1, unsigned int D2, class R>
738 template <class SubMatrix>
739 SubMatrix SMatrix<T,D1,D2,R>::Sub(unsigned int row0, unsigned int col0) const {
740 
741  SubMatrix tmp;
743  (tmp,*this,row0,col0);
744  return tmp;
745 }
746 
747 //diagonal
748 template <class T, unsigned int D1, unsigned int D2, class R>
750 
751  // only for squared matrices
752  STATIC_CHECK( D1 == D2,SMatrix_NOT_square );
753 
754  SVector<T,D1> tmp;
755  for(unsigned int i=0; i<D1; ++i) {
756  tmp[i] = fRep[ i*D2 + i];
757  }
758  return tmp;
759 }
760 
761 //set diagonal
762 template <class T, unsigned int D1, unsigned int D2, class R>
763 template <class Vector>
764 void SMatrix<T,D1,D2,R>::SetDiagonal( const Vector & v) {
765 
766  // check size that size of vector is correct
767  STATIC_CHECK( ( ( D1 <= D2) && Vector::kSize == D1 ) ||
768  ( ( D2 < D1 ) && Vector::kSize == D2 ), SVector_size_NOT_correct );
769 
770 
771  for(unsigned int i=0; i<Vector::kSize; ++i) {
772  fRep[ i*D2 + i] = v[i];
773  }
774 }
775 
776 // matrix trace
777 template <class T, unsigned int D1, unsigned int D2, class R>
779  unsigned int diagSize = D1;
780  if (D2 < D1) diagSize = D2;
781  T trace = 0;
782  for(unsigned int i=0; i< diagSize; ++i) {
783  trace += fRep[ i*D2 + i] ;
784  }
785  return trace;
786 }
787 
788 //upper block
789 template <class T, unsigned int D1, unsigned int D2, class R>
790 #ifndef UNSUPPORTED_TEMPLATE_EXPRESSION
791 SVector<T, D1 * (D2 +1)/2 > SMatrix<T,D1,D2,R>::UpperBlock( ) const {
792 #else
793 template <class SubVector>
794 SubVector SMatrix<T,D1,D2,R>::UpperBlock( ) const {
795 #endif
796  // only for squared matrices
797  STATIC_CHECK( D1 == D2,SMatrix_NOT_square );
798 
799 #ifndef UNSUPPORTED_TEMPLATE_EXPRESSION
800  SVector<T, D1 * (D2 +1)/2 > tmp;
801 #else
802  // N must be equal D1 *(D1 +1)/2
803  STATIC_CHECK( SubVector::kSize == D1*(D1+1)/2,SVector_Wrong_Size );
804  SubVector tmp;
805 #endif
806 
807  int k = 0;
808  for(unsigned int i=0; i<D1; ++i) {
809  for(unsigned int j=i; j<D2; ++j) {
810  tmp[k] = fRep[ i*D2 + j];
811  k++;
812  }
813  }
814  return tmp;
815 }
816 
817 //lower block
818 template <class T, unsigned int D1, unsigned int D2, class R>
819 #ifndef UNSUPPORTED_TEMPLATE_EXPRESSION
820 SVector<T, D1 * (D2 +1)/2 > SMatrix<T,D1,D2,R>::LowerBlock( ) const {
821 #else
822 template <class SubVector>
823 SubVector SMatrix<T,D1,D2,R>::LowerBlock( ) const {
824 #endif
825 
826  // only for squared matrices
827  STATIC_CHECK( D1 == D2,SMatrix_NOT_square );
828 
829 #ifndef UNSUPPORTED_TEMPLATE_EXPRESSION
830  SVector<T, D1 * (D2 +1)/2 > tmp;
831 #else
832  // N must be equal D1 *(D1 +1)/2
833  STATIC_CHECK( SubVector::kSize == D1*(D1+1)/2,SVector_Wrong_Size );
834  SubVector tmp;
835 #endif
836 
837  int k = 0;
838  for(unsigned int i=0; i<D1; ++i) {
839  for(unsigned int j=0; j<=i; ++j) {
840  tmp[k] = fRep[ i*D2 + j];
841  k++;
842  }
843  }
844  return tmp;
845 }
846 
847 /// construct from upper/lower block
848 
849 //lower block
850 template <class T, unsigned int D1, unsigned int D2, class R>
851 #ifndef UNSUPPORTED_TEMPLATE_EXPRESSION
852 SMatrix<T,D1,D2,R>::SMatrix(const SVector<T, D1*(D2+1)/2 > & v, bool lower ) {
853 #else
854 template <unsigned int N>
855 SMatrix<T,D1,D2,R>::SMatrix(const SVector<T, N > & v, bool lower ) {
856 #endif
857 
858  // only for squared matrices
859  STATIC_CHECK( D1 == D2,SMatrix_NOT_square );
860 
861 #ifdef UNSUPPORTED_TEMPLATE_EXPRESSION
862  STATIC_CHECK( N == D1*(D1+1)/2,SVector_Wrong_Size );
863 #endif
864 
865  int k = 0;
866  if (lower) {
867  // case of lower block
868  for(unsigned int i=0; i<D1; ++i) {
869  for(unsigned int j=0; j<=i; ++j) {
870  fRep[ i*D2 + j] = v[k];
871  if ( i != j) fRep[ j*D2 + i] = v[k];
872  k++;
873  }
874  }
875  } else {
876  // case of upper block
877  for(unsigned int i=0; i<D1; ++i) {
878  for(unsigned int j=i; j<D2; ++j) {
879  fRep[ i*D2 + j] = v[k];
880  if ( i != j) fRep[ j*D2 + i] = v[k];
881  k++;
882  }
883  }
884  }
885 }
886 
887 
888 template <class T, unsigned int D1, unsigned int D2, class R>
889 bool SMatrix<T,D1,D2,R>::IsInUse( const T * p) const {
890  return p == fRep.Array();
891 }
892 
893 
894 
895 
896  } // namespace Math
897 
898 } // namespace ROOT
899 
900 
901 
902 #endif
SMatrix< T, D1, D2, R > InverseChol(int &ifail) const
Invert of a symmetric positive defined Matrix using Choleski decomposition.
Definition: SMatrix.icc:446
SMatrix< T, D1, D2, R > & operator+=(const T &rhs)
addition with a scalar
Definition: SMatrix.icc:191
SMatrix< T, D1, D2, R > & operator*=(const T &rhs)
multiplication with a scalar
Definition: SMatrix.icc:252
void SetDiagonal(const Vector &v)
Set the diagonal elements from a Vector Require that vector implements ::kSize since a check (statica...
Definition: SMatrix.icc:764
SMatrix< T, D1, D2, R > & operator/=(const T &rhs)
division with a scalar
Definition: SMatrix.icc:281
bool operator!=(const T &rhs) const
element wise comparison
Definition: SMatrix.icc:321
Namespace for new ROOT classes and functions.
Definition: TFoamSampler.h:19
bool operator<(const T &rhs) const
element wise comparison
Definition: SMatrix.icc:373
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:222
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:624
bool Invert()
Invert a square Matrix ( this method changes the current matrix).
Definition: SMatrix.icc:406
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:569
bool InvertFast()
Fast Invertion of a square Matrix ( this method changes the current matrix).
Definition: SMatrix.icc:423
bool InvertChol()
Invertion of a symmetric positive defined Matrix using Choleski decomposition.
Definition: SMatrix.icc:440
#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:722
SMatrix< T, D1, D2, R > InverseFast(int &ifail) const
Invert a square Matrix and returns a new matrix.
Definition: SMatrix.icc:430
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:706
bool Det(T &det)
determinant of square Matrix via Dfact.
Definition: SMatrix.icc:460
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:633
#define N
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:889
const T & At(unsigned int i, unsigned int j) const
read only access to matrix element, with indices starting from 0.
Definition: SMatrix.icc:647
Bool_t operator==(const TMatrixTBase< Element > &m1, const TMatrixTBase< Element > &m2)
Check to see if two matrices are identical.
R fRep
Matrix Storage Object containing matrix data.
Definition: SMatrix.h:707
SVector< double, 2 > v
Definition: Dict.h:5
std::ostream & Print(std::ostream &os) const
Print: used by operator<<()
Definition: SMatrix.icc:597
iterator begin()
STL iterator interface.
Definition: SMatrix.icc:664
SVector< T, D1 > Diagonal() const
return diagonal elements of a matrix as a Vector.
Definition: SMatrix.icc:749
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:478
T apply(unsigned int i) const
access the parse tree. Index starts from zero
Definition: SVector.icc:533
bool Det2(T &det) const
determinant of square Matrix via Dfact.
Definition: SMatrix.icc:467
SMatrix< T, D1, D2, R > Inverse(int &ifail) const
Invert a square Matrix and returns a new matrix.
Definition: SMatrix.icc:413
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:739
bool operator>(const T &rhs) const
element wise comparison
Definition: SMatrix.icc:341
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:621
T Trace() const
return the trace of a matrix Sum of the diagonal elements
Definition: SMatrix.icc:778
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.
bool operator==(const T &rhs) const
element wise comparison
Definition: SMatrix.icc:293
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:820
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:546
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:512
void SetElements(InputIterator begin, InputIterator end, bool triang=false, bool lower=true)
Set matrix elements with STL iterator interface.
Definition: SMatrix.icc:686
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:791
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:584
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:669