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