root/Math/SMatrix.icc

0001 // @(#)root/smatrix:$Id$
0002 // Authors: T. Glebe, L. Moneta    2005
0003
0004 #ifndef ROOT_Math_SMatrix_icc
0005 #define ROOT_Math_SMatrix_icc
0006 // ********************************************************************
0007 //
0008 // source:
0009 //
0010 // type:      source code
0011 //
0012 // created:   21. Mar 2001
0013 //
0014 // author:    Thorsten Glebe
0015 //            HERA-B Collaboration
0016 //            Max-Planck-Institut fuer Kernphysik
0017 //            Saupfercheckweg 1
0018 //            69117 Heidelberg
0019 //            Germany
0020 //            E-mail: T.Glebe@mpi-hd.mpg.de
0021 //
0022 // Description: SMatrix implementation file
0023 //
0024 // changes:
0025 // 21 Mar 2001 (TG) creation
0026 // 26 Mar 2001 (TG) place_in_row(), place_in_col() added
0027 // 03 Apr 2001 (TG) invert() added
0028 // 07 Apr 2001 (TG) CTOR from SVertex (dyadic product) added
0029 // 09 Apr 2001 (TG) CTOR from array added
0030 // 25 Mai 2001 (TG) row(), col() added
0031 // 11 Jul 2001 (TG) added #include Functions.hh
0032 // 11 Jan 2002 (TG) added operator==(), operator!=()
0033 // 14 Jan 2002 (TG) added more operator==(), operator!=(), operator>(), operator<()
0034 //
0035 // ********************************************************************
0036
0037 #ifndef ROOT_Math_SMatrix
0038 #error "Do not use SMatrix.icc directly. #include \"Math/SMatrix.h\" instead."
0039 #endif // ROOT_Math_SMatrix
0040
0041 #include <iostream>
0042 #include <iomanip>
0043 #include <assert.h>
0044 //#ifndef ROOT_Math_Dsinv
0045 //#include "Math/Dsinv.h"
0046 //#endif
0047 //#include "Math/Dsinv_array.h"
0048 //#include "Math/Dsfact.h"
0049
0050 #include "Math/Dfact.h"
0051 #include "Math/Dinv.h"
0052 #include "Math/Functions.h"
0053 #include "Math/HelperOps.h"
0054 #include "Math/StaticCheck.h"
0055
0056
0057
0058
0059
0060
0061
0062 namespace ROOT {
0063
0064 namespace Math {
0065
0066
0067
0068 //==============================================================================
0069 // Constructors
0070 //==============================================================================
0071 template <class T, unsigned int D1, unsigned int D2, class R>
0072 SMatrix<T,D1,D2,R>::SMatrix() {
0073    // operator=(0);   // if operator=(T ) is defined
0074    for(unsigned int i=0; i<R::kSize; ++i) fRep.Array()[i] = 0;
0075 }
0076
0077 //identity
0078 template <class T, unsigned int D1, unsigned int D2, class R>
0079 SMatrix<T,D1,D2,R>::SMatrix( SMatrixIdentity ) {
0080    for(unsigned int i=0; i<R::kSize; ++i)
0081       fRep.Array()[i] = 0;
0082    if (D1 <= D2) {
0083       for(unsigned int i=0; i<D1; ++i)
0084          fRep[i*D2+i] = 1;
0085    }
0086    else {
0087       for(unsigned int i=0; i<D2; ++i)
0088          fRep[i*D2+i] = 1;
0089    }
0090 }
0091
0092 template <class T, unsigned int D1, unsigned int D2, class R>
0093 SMatrix<T,D1,D2,R>::SMatrix(const SMatrix<T,D1,D2,R>& rhs) {
0094    fRep = rhs.fRep;
0095 }
0096
0097
0098 template <class T, unsigned int D1, unsigned int D2, class R>
0099 template <class R2>
0100 SMatrix<T,D1,D2,R>::SMatrix(const SMatrix<T,D1,D2,R2>& rhs) {
0101    operator=(rhs);
0102 }
0103
0104
0105 template <class T, unsigned int D1, unsigned int D2, class R>
0106 template <class A, class R2>
0107 SMatrix<T,D1,D2,R>::SMatrix(const Expr<A,T,D1,D2,R2>& rhs) {
0108    operator=(rhs);
0109 }
0110
0111
0112 //=============================================================================
0113 // New Constructors from STL interfaces
0114 //=============================================================================
0115 template <class T, unsigned int D1, unsigned int D2, class R>
0116 template <class InputIterator>
0117 SMatrix<T,D1,D2,R>::SMatrix(InputIterator ibegin, InputIterator iend, bool triang, bool lower) {
0118    // assume iterator size == matrix size
0119    for(unsigned int i=0; i<R::kSize; ++i) fRep.Array()[i] = 0;
0120    AssignItr<T,D1,D2,R>::Evaluate(*this,ibegin,iend,triang,lower);
0121 }
0122
0123 template <class T, unsigned int D1, unsigned int D2, class R>
0124 template <class InputIterator>
0125 SMatrix<T,D1,D2,R>::SMatrix(InputIterator ibegin, unsigned int size, bool triang, bool lower) {
0126    // assume iterator size <= matrix size (no check needed in AssignItr)
0127    assert( size <= R::kSize);
0128    for(unsigned int i=0; i<R::kSize; ++i) fRep.Array()[i] = 0;
0129    AssignItr<T,D1,D2,R>::Evaluate(*this,ibegin,ibegin+size,triang,lower,false);
0130 }
0131
0132
0133 //==============================================================================
0134 // Assignment and operator= for scalar types for matrices of size 1
0135 // compiles only for matrices of size 1
0136 //==============================================================================
0137 template <class T, unsigned int D1, unsigned int D2, class R>
0138 SMatrix<T,D1,D2,R>::SMatrix(const T&  rhs) {
0139    STATIC_CHECK( kSize == 1,SMatrix_NOT_of_size_1 );
0140    fRep[0] = rhs;
0141 }
0142
0143 template <class T, unsigned int D1, unsigned int D2, class R>
0144 SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::operator=(const T&  rhs) {
0145    STATIC_CHECK( kSize == 1,SMatrix_NOT_of_size_1 );
0146    fRep[0] = rhs;
0147    return *this;
0148 }
0149
0150 //=============================================================================
0151 //=============================================================================
0152
0153 template <class T, unsigned int D1, unsigned int D2, class R>
0154 template <class M>
0155 SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::operator=(const M&  rhs) {
0156    fRep = rhs.fRep;
0157    return *this;
0158 }
0159
0160 template <class T, unsigned int D1, unsigned int D2, class R>
0161 SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::operator=(const SMatrix<T,D1,D2,R>&  rhs) {
0162    fRep = rhs.fRep;
0163    return *this;
0164 }
0165
0166 template <class T, unsigned int D1, unsigned int D2, class R>
0167 template <class A, class R2>
0168 SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::operator=(const Expr<A,T,D1,D2,R2>&  rhs) {
0169
0170    Assign<T,D1,D2,A,R,R2>::Evaluate(*this, rhs);
0171    return *this;
0172 }
0173
0174 //=============================================================================
0175 // assign from an identity
0176 template <class T, unsigned int D1, unsigned int D2, class R>
0177 SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::operator= ( SMatrixIdentity ) {
0178    for(unsigned int i=0; i<R::kSize; ++i)
0179       fRep.Array()[i] = 0;
0180    if (D1 <= D2) {
0181       for(unsigned int i=0; i<D1; ++i)
0182          fRep[i*D2+i] = 1;
0183    }
0184    else {
0185       for(unsigned int i=0; i<D2; ++i)
0186          fRep[i*D2+i] = 1;
0187    }
0188    return *this;
0189 }
0190
0191
0192
0193 //=============================================================================
0194 // operator+=
0195 //=============================================================================
0196 template <class T, unsigned int D1, unsigned int D2, class R>
0197 SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::operator+=(const T&  rhs) {
0198    // self-addition with a scalar value
0199    for(unsigned int i=0; i<R::kSize; ++i) {
0200       fRep.Array()[i] += rhs;
0201    }
0202    return *this;
0203 }
0204
0205 template <class T, unsigned int D1, unsigned int D2, class R>
0206 template <class R2>
0207 SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::operator+=(const SMatrix<T,D1,D2,R2>&  rhs) {
0208    // self-addition with another matrix (any representation)
0209    // use operator+= of the representation object
0210    fRep += rhs.fRep;
0211    return *this;
0212 }
0213
0214
0215 template <class T, unsigned int D1, unsigned int D2, class R>
0216 template <class A, class R2>
0217 SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::operator+=(const Expr<A,T,D1,D2,R2>&  rhs) {
0218    // self-addition with an expression
0219    PlusEquals<T,D1,D2,A,R,R2>::Evaluate(*this, rhs);
0220    return *this;
0221 }
0222
0223
0224 //==============================================================================
0225 // operator-=
0226 //==============================================================================
0227 template <class T, unsigned int D1, unsigned int D2, class R>
0228 SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::operator-=(const T&  rhs) {
0229    // self-subtraction with a scalar value
0230    for(unsigned int i=0; i<R::kSize; ++i) {
0231       fRep.Array()[i] -= rhs;
0232    }
0233    return *this;
0234 }
0235
0236 template <class T, unsigned int D1, unsigned int D2, class R>
0237 template <class R2>
0238 SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::operator-=(const SMatrix<T,D1,D2,R2>&  rhs) {
0239    // self-subtraction with another matrix (any representation)
0240    // use operator-= of the representation object
0241    fRep -= rhs.fRep;
0242    return *this;
0243 }
0244
0245
0246 template <class T, unsigned int D1, unsigned int D2, class R>
0247 template <class A, class R2>
0248 SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::operator-=(const Expr<A,T,D1,D2,R2>&  rhs) {
0249    // self-subtraction with an expression
0250    MinusEquals<T,D1,D2,A,R,R2>::Evaluate(*this, rhs);
0251    return *this;
0252 }
0253
0254 //==============================================================================
0255 // operator*=
0256 //==============================================================================
0257 template <class T, unsigned int D1, unsigned int D2, class R>
0258 SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::operator*=(const T&  rhs) {
0259    // case of multiplication with a scalar
0260    for(unsigned int i=0; i<R::kSize; ++i) {
0261       fRep.Array()[i] *= rhs;
0262    }
0263    return  *this;
0264 }
0265
0266 template <class T, unsigned int D1, unsigned int D2, class R>
0267 template <class R2>
0268 SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::operator*=(const SMatrix<T,D1,D2,R2>&  rhs) {
0269    // self-multiplication with another matrix (will work only for square matrices)
0270    // a temporary is needed and will be created automatically to store intermediate result
0271    return operator=(*this * rhs);
0272 }
0273
0274 template <class T, unsigned int D1, unsigned int D2, class R>
0275 template <class A, class R2>
0276 SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::operator*=(const Expr<A,T,D1,D2,R2>&  rhs) {
0277    // self-multiplication with an expression (will work only for square matrices)
0278    // a temporary is needed and will be created automatically to store intermediate result
0279    return operator=(*this * rhs);
0280 }
0281
0282
0283 //==============================================================================
0284 // operator/= (only for scalar values)
0285 //==============================================================================
0286 template <class T, unsigned int D1, unsigned int D2, class R>
0287 SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::operator/=(const T&  rhs) {
0288    // division with a scalar
0289    for(unsigned int i=0; i<R::kSize; ++i) {
0290       fRep.Array()[i] /= rhs;
0291    }
0292    return  *this;
0293 }
0294
0295 //==============================================================================
0296 // operator== (element wise comparison)
0297 //==============================================================================
0298 template <class T, unsigned int D1, unsigned int D2, class R>
0299 bool SMatrix<T,D1,D2,R>::operator==(const T& rhs) const {
0300    bool rc = true;
0301    for(unsigned int i=0; i<R::kSize; ++i) {
0302       rc = rc && (fRep.Array()[i] == rhs);
0303    }
0304    return rc;
0305 }
0306
0307 template <class T, unsigned int D1, unsigned int D2, class R>
0308 template <class R2>
0309 bool SMatrix<T,D1,D2,R>::operator==(const SMatrix<T,D1,D2,R2>& rhs) const {
0310    return fRep == rhs.fRep;
0311 }
0312
0313 template <class T, unsigned int D1, unsigned int D2, class R>
0314 template <class A, class R2>
0315 bool SMatrix<T,D1,D2,R>::operator==(const Expr<A,T,D1,D2,R2>& rhs) const {
0316    bool rc = true;
0317    for(unsigned int i=0; i<D1*D2; ++i) {
0318       rc = rc && (fRep[i] == rhs.apply(i));
0319    }
0320    return rc;
0321 }
0322
0323 //==============================================================================
0324 // operator!= (element wise comparison)
0325 //==============================================================================
0326 template <class T, unsigned int D1, unsigned int D2, class R>
0327 inline bool SMatrix<T,D1,D2,R>::operator!=(const T& rhs) const {
0328    return !operator==(rhs);
0329 }
0330
0331 template <class T, unsigned int D1, unsigned int D2, class R>
0332 inline bool SMatrix<T,D1,D2,R>::operator!=(const SMatrix<T,D1,D2,R>& rhs) const {
0333    return !operator==(rhs);
0334 }
0335
0336 template <class T, unsigned int D1, unsigned int D2, class R>
0337 template <class A, class R2>
0338 inline bool SMatrix<T,D1,D2,R>::operator!=(const Expr<A,T,D1,D2,R2>& rhs) const {
0339    return !operator==(rhs);
0340 }
0341
0342
0343 //==============================================================================
0344 // operator> (element wise comparison)
0345 //==============================================================================
0346 template <class T, unsigned int D1, unsigned int D2, class R>
0347 bool SMatrix<T,D1,D2,R>::operator>(const T& rhs) const {
0348    bool rc = true;
0349    for(unsigned int i=0; i<D1*D2; ++i) {
0350       rc = rc && (fRep[i] > rhs);
0351    }
0352    return rc;
0353 }
0354
0355 template <class T, unsigned int D1, unsigned int D2, class R>
0356 template <class R2>
0357 bool SMatrix<T,D1,D2,R>::operator>(const SMatrix<T,D1,D2,R2>& rhs) const {
0358    bool rc = true;
0359    for(unsigned int i=0; i<D1*D2; ++i) {
0360       rc = rc && (fRep[i] > rhs.fRep[i]);
0361    }
0362    return rc;
0363 }
0364
0365 template <class T, unsigned int D1, unsigned int D2, class R>
0366 template <class A, class R2>
0367 bool SMatrix<T,D1,D2,R>::operator>(const Expr<A,T,D1,D2,R2>& rhs) const {
0368    bool rc = true;
0369    for(unsigned int i=0; i<D1*D2; ++i) {
0370       rc = rc && (fRep[i] > rhs.apply(i));
0371    }
0372    return rc;
0373 }
0374
0375 //==============================================================================
0376 // operator< (element wise comparison)
0377 //==============================================================================
0378 template <class T, unsigned int D1, unsigned int D2, class R>
0379 bool SMatrix<T,D1,D2,R>::operator<(const T& rhs) const {
0380    bool rc = true;
0381    for(unsigned int i=0; i<D1*D2; ++i) {
0382       rc = rc && (fRep[i] < rhs);
0383    }
0384    return rc;
0385 }
0386
0387 template <class T, unsigned int D1, unsigned int D2, class R>
0388 template <class R2>
0389 bool SMatrix<T,D1,D2,R>::operator<(const SMatrix<T,D1,D2,R2>& rhs) const {
0390    bool rc = true;
0391    for(unsigned int i=0; i<D1*D2; ++i) {
0392       rc = rc && (fRep[i] < rhs.fRep[i]);
0393    }
0394    return rc;
0395 }
0396
0397 template <class T, unsigned int D1, unsigned int D2, class R>
0398 template <class A, class R2>
0399 bool SMatrix<T,D1,D2,R>::operator<(const Expr<A,T,D1,D2,R2>& rhs) const {
0400    bool rc = true;
0401    for(unsigned int i=0; i<D1*D2; ++i) {
0402       rc = rc && (fRep[i] < rhs.apply(i));
0403    }
0404    return rc;
0405 }
0406
0407
0408 //==============================================================================
0409 // invert
0410 //==============================================================================
0411 template <class T, unsigned int D1, unsigned int D2, class R>
0412 inline bool SMatrix<T,D1,D2,R>::Invert() {
0413    STATIC_CHECK( D1 == D2,SMatrix_not_square);
0414    return Inverter<D1,D2>::Dinv((*this).fRep);
0415 }
0416
0417 // invert returning a new matrix
0418 template <class T, unsigned int D1, unsigned int D2, class R>
0419 inline SMatrix<T,D1,D2,R> SMatrix<T,D1,D2,R>::Inverse(int & ifail) const {
0420    SMatrix<T,D1,D2,R> tmp(*this);
0421    bool ok = tmp.Invert();
0422    ifail = 0;
0423    if (!ok) ifail = 1;
0424    return tmp;
0425 }
0426
0427 // fast inversion
0428 template <class T, unsigned int D1, unsigned int D2, class R>
0429 inline bool SMatrix<T,D1,D2,R>::InvertFast() {
0430    STATIC_CHECK( D1 == D2,SMatrix_not_square);
0431    return FastInverter<D1,D2>::Dinv((*this).fRep);
0432 }
0433
0434 // fast inversion returning a new matrix
0435 template <class T, unsigned int D1, unsigned int D2, class R>
0436 inline SMatrix<T,D1,D2,R> SMatrix<T,D1,D2,R>::InverseFast(int & ifail) const {
0437    SMatrix<T,D1,D2,R> tmp(*this);
0438    bool ok = tmp.InvertFast();
0439    ifail = 0;
0440    if (!ok) ifail = 1;
0441    return tmp;
0442 }
0443
0444 // Choleski inversion for symmetric and positive defined matrices
0445 template <class T, unsigned int D1, unsigned int D2, class R>
0446 inline bool SMatrix<T,D1,D2, R>::InvertChol() {
0447    STATIC_CHECK( D1 == D2,SMatrix_not_square);
0448    return CholInverter<D1>::Dinv((*this).fRep);
0449 }
0450
0451 template <class T, unsigned int D1, unsigned int D2, class R>
0452 inline SMatrix<T,D1,D2, R>  SMatrix<T,D1,D2, R>::InverseChol(int &ifail) const {
0453    SMatrix<T,D1,D2,R > tmp(*this);
0454    bool ok = tmp.InvertChol();
0455    ifail = 0;
0456    if (!ok) ifail = 1;
0457    return tmp;
0458 }
0459
0460
0461
0462 //==============================================================================
0463 // det
0464 //==============================================================================
0465 template <class T, unsigned int D1, unsigned int D2, class R>
0466 inline bool SMatrix<T,D1,D2,R>::Det(T& det) {
0467    STATIC_CHECK( D1 == D2,SMatrix_not_square);
0468    //  return Dfact<SMatrix<T,D1,D1>, D1, D1>(*this,det);
0469    //return Dfact<R, D1, D1>(fRep, det);
0470    return Determinant<D1,D2>::Dfact(fRep, det);
0471 }
0472 template <class T, unsigned int D1, unsigned int D2, class R>
0473 inline bool SMatrix<T,D1,D2,R>::Det2(T& det) const {
0474    SMatrix<T,D1,D2,R> tmp(*this);
0475    return tmp.Det(det);
0476 }
0477
0478
0479 //==============================================================================
0480 // place_in_row
0481 //==============================================================================
0482 template <class T, unsigned int D1, unsigned int D2, class R>
0483 template <unsigned int D>
0484 SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::Place_in_row(const SVector<T,D>& rhs,
0485                                                      unsigned int row,
0486                                                      unsigned int col) {
0487
0488    assert(col+D <= D2);
0489
0490    for(unsigned int i=row*D2+col, j=0; j<D; ++i, ++j) {
0491       fRep[i] = rhs.apply(j);
0492    }
0493    return *this;
0494 }
0495
0496 //==============================================================================
0497 // place_in_row
0498 //==============================================================================
0499 template <class T, unsigned int D1, unsigned int D2, class R>
0500 template <class A, unsigned int D>
0501 SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::Place_in_row(const VecExpr<A,T,D>& rhs,
0502                                                      unsigned int row,
0503                                                      unsigned int col) {
0504
0505    assert(col+D <= D2);
0506
0507    for(unsigned int i=row*D2+col, j=0; j<D; ++i, ++j) {
0508       fRep[i] = rhs.apply(j);
0509    }
0510    return *this;
0511 }
0512
0513 //==============================================================================
0514 // place_in_col
0515 //==============================================================================
0516 template <class T, unsigned int D1, unsigned int D2, class R>
0517 template <unsigned int D>
0518 SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::Place_in_col(const SVector<T,D>& rhs,
0519                                                      unsigned int row,
0520                                                      unsigned int col) {
0521
0522    assert(row+D <= D1);
0523
0524    for(unsigned int i=row*D2+col, j=0; j<D; i+=D2, ++j) {
0525       fRep[i] = rhs.apply(j);
0526    }
0527    return *this;
0528 }
0529
0530 //==============================================================================
0531 // place_in_col
0532 //==============================================================================
0533 template <class T, unsigned int D1, unsigned int D2, class R>
0534 template <class A, unsigned int D>
0535 SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::Place_in_col(const VecExpr<A,T,D>& rhs,
0536                                                      unsigned int row,
0537                                                      unsigned int col) {
0538
0539    assert(row+D <= D1);
0540
0541    for(unsigned int i=row*D2+col, j=0; j<D; i+=D2, ++j) {
0542       fRep[i] = rhs.apply(j);
0543    }
0544    return *this;
0545 }
0546
0547 //==============================================================================
0548 // place_at
0549 //==============================================================================
0550 template <class T, unsigned int D1, unsigned int D2, class R>
0551 template <unsigned int D3, unsigned int D4, class R2>
0552 SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::Place_at(const SMatrix<T,D3,D4,R2>& rhs,
0553                                                  unsigned int row,
0554                                                  unsigned int col) {
0555    PlaceMatrix<T,D1,D2,D3,D4,R,R2>::Evaluate(*this,rhs,row,col);
0556    return *this;
0557 }
0558
0559 //==============================================================================
0560 // place_at
0561 //==============================================================================
0562 template <class T, unsigned int D1, unsigned int D2, class R>
0563 template <class A, unsigned int D3, unsigned int D4, class R2>
0564 SMatrix<T,D1,D2,R>& SMatrix<T,D1,D2,R>::Place_at(const Expr<A,T,D3,D4,R2>& rhs,
0565                                                  unsigned int row,
0566                                                  unsigned int col) {
0567    PlaceExpr<T,D1,D2,D3,D4,A,R,R2>::Evaluate(*this,rhs,row,col);
0568    return *this;
0569 }
0570
0571 //==============================================================================
0572 // row
0573 //==============================================================================
0574 template <class T, unsigned int D1, unsigned int D2, class R>
0575 SVector<T,D2> SMatrix<T,D1,D2,R>::Row(unsigned int therow) const {
0576
0577    const unsigned int offset = therow*D2;
0578
0579    /*static*/  SVector<T,D2> tmp;
0580    for(unsigned int i=0; i<D2; ++i) {
0581       tmp[i] = fRep[offset+i];
0582    }
0583    return tmp;
0584 }
0585
0586 //==============================================================================
0587 // col
0588 //==============================================================================
0589 template <class T, unsigned int D1, unsigned int D2, class R>
0590 SVector<T,D1> SMatrix<T,D1,D2,R>::Col(unsigned int thecol) const {
0591
0592    /*static */ SVector<T,D1> tmp;
0593    for(unsigned int i=0; i<D1; ++i) {
0594       tmp[i] = fRep[thecol+i*D2];
0595    }
0596    return tmp;
0597 }
0598
0599 //==============================================================================
0600 // print
0601 //==============================================================================
0602 template <class T, unsigned int D1, unsigned int D2, class R>
0603 std::ostream& SMatrix<T,D1,D2,R>::Print(std::ostream& os) const {
0604    const std::ios_base::fmtflags prevFmt = os.setf(std::ios::right,std::ios::adjustfield);
0605    //  os.setf(ios::fixed);
0606
0607    os << "[ ";
0608    for (unsigned int i=0; i < D1; ++i) {
0609       for (unsigned int j=0; j < D2; ++j) {
0610          os << std::setw(12) << fRep[i*D2+j];
0611          if ((!((j+1)%12)) && (j < D2-1))
0612             os << std::endl << "         ...";
0613       }
0614       if (i != D1 - 1)
0615          os << std::endl  << "  ";
0616    }
0617    os << " ]";
0618
0619    if (prevFmt != os.flags() ) os.setf(prevFmt, std::ios::adjustfield);
0620    return os;
0621 }
0622
0623 //==============================================================================
0624 // Access functions
0625 //==============================================================================
0626 template <class T, unsigned int D1, unsigned int D2, class R>
0627 inline T SMatrix<T,D1,D2,R>::apply(unsigned int i) const { return fRep[i]; }
0628
0629 template <class T, unsigned int D1, unsigned int D2, class R>
0630 inline const T* SMatrix<T,D1,D2,R>::Array() const { return fRep.Array(); }
0631
0632 template <class T, unsigned int D1, unsigned int D2, class R>
0633 inline T* SMatrix<T,D1,D2,R>::Array() { return fRep.Array(); }
0634
0635 //==============================================================================
0636 // Operators
0637 //==============================================================================
0638 template <class T, unsigned int D1, unsigned int D2, class R>
0639 inline const T& SMatrix<T,D1,D2,R>::operator()(unsigned int i, unsigned int j) const {
0640    return fRep(i,j);
0641 }
0642
0643 template <class T, unsigned int D1, unsigned int D2, class R>
0644 inline T& SMatrix<T,D1,D2,R>::operator()(unsigned int i, unsigned int j) {
0645    return fRep(i,j);
0646 }
0647
0648
0649 //==============================================================================
0650 // Element access with At()
0651 //==============================================================================
0652 template <class T, unsigned int D1, unsigned int D2, class R>
0653 inline const T& SMatrix<T,D1,D2,R>::At(unsigned int i, unsigned int j) const {
0654    assert(i < D1);
0655    assert(j < D2);
0656    return fRep(i,j);
0657 }
0658
0659 template <class T, unsigned int D1, unsigned int D2, class R>
0660 inline T& SMatrix<T,D1,D2,R>::At(unsigned int i, unsigned int j) {
0661    assert(i < D1);
0662    assert(j < D2);
0663    return fRep(i,j);
0664 }
0665
0666 //==============================================================================
0667 // STL interface
0668 //==============================================================================
0669 template <class T, unsigned int D1, unsigned int D2, class R>
0670 inline T * SMatrix<T,D1,D2,R>::begin() {
0671    return fRep.Array();
0672 }
0673
0674 template <class T, unsigned int D1, unsigned int D2, class R>
0675 inline T * SMatrix<T,D1,D2,R>::end() {
0676    return fRep.Array() + R::kSize;
0677 }
0678
0679 template <class T, unsigned int D1, unsigned int D2, class R>
0680 inline const T * SMatrix<T,D1,D2,R>::begin() const {
0681    return fRep.Array();
0682 }
0683
0684 template <class T, unsigned int D1, unsigned int D2, class R>
0685 inline const T * SMatrix<T,D1,D2,R>::end() const {
0686    return fRep.Array() + R::kSize;
0687 }
0688
0689
0690 template <class T, unsigned int D1, unsigned int D2, class R>
0691 template <class InputIterator>
0692 void SMatrix<T,D1,D2,R>::SetElements(InputIterator ibegin, InputIterator iend, bool triang, bool lower) {
0693    // assume iterator size == matrix size when filling full matrix
0694    AssignItr<T,D1,D2,R>::Evaluate(*this,ibegin,iend,triang,lower);
0695 }
0696
0697 template <class T, unsigned int D1, unsigned int D2, class R>
0698 template <class InputIterator>
0699 void SMatrix<T,D1,D2,R>::SetElements(InputIterator ibegin, unsigned int size, bool triang, bool lower) {
0700    // assume iterator size <= matrix size  (no check to be done in AssignItr)
0701    assert( size <= R::kSize);
0702    AssignItr<T,D1,D2,R>::Evaluate(*this,ibegin,ibegin+size,triang,lower,false);
0703 }
0704
0705
0706
0707 //==============================================================================
0708 // SubMatrices and slices of columns and rows
0709 //==============================================================================
0710 template <class T, unsigned int D1, unsigned int D2, class R>
0711 template <class SubVector>
0712 SubVector SMatrix<T,D1,D2,R>::SubRow(unsigned int therow, unsigned int col0 ) const {
0713
0714    const unsigned int offset = therow*D2 + col0;
0715
0716    STATIC_CHECK( SubVector::kSize <= D2,SVector_dimension_too_small);
0717    assert(col0 + SubVector::kSize <= D2);
0718
0719    SubVector tmp;
0720    for(unsigned int i=0; i<SubVector::kSize; ++i) {
0721       tmp[i] = fRep[offset+i];
0722    }
0723    return tmp;
0724 }
0725
0726 template <class T, unsigned int D1, unsigned int D2, class R>
0727 template <class SubVector>
0728 SubVector SMatrix<T,D1,D2,R>::SubCol(unsigned int thecol, unsigned int row0 ) const {
0729
0730    const unsigned int offset = thecol + row0*D1;
0731
0732    STATIC_CHECK( SubVector::kSize <= D1,SVector_dimension_too_small);
0733    assert(row0 + SubVector::kSize <= D1);
0734
0735    SubVector tmp;
0736    for(unsigned int i=0; i<SubVector::kSize; ++i) {
0737       tmp[i] = fRep[offset+i*D1];
0738    }
0739    return tmp;
0740 }
0741
0742 // submatrix
0743 template <class T, unsigned int D1, unsigned int D2, class R>
0744 template <class SubMatrix>
0745 SubMatrix SMatrix<T,D1,D2,R>::Sub(unsigned int row0, unsigned int col0) const {
0746
0747    SubMatrix tmp;
0748    RetrieveMatrix<T,SubMatrix::kRows, SubMatrix::kCols, D1, D2, typename SubMatrix::rep_type, R>::Evaluate
0749       (tmp,*this,row0,col0);
0750    return tmp;
0751 }
0752
0753 //diagonal
0754 template <class T, unsigned int D1, unsigned int D2, class R>
0755 SVector<T,D1> SMatrix<T,D1,D2,R>::Diagonal( ) const {
0756
0757    // only for squared matrices
0758    STATIC_CHECK( D1 == D2,SMatrix_NOT_square );
0759
0760    SVector<T,D1> tmp;
0761    for(unsigned int i=0; i<D1; ++i) {
0762       tmp[i] = fRep[ i*D2 + i];
0763    }
0764    return tmp;
0765 }
0766
0767 //set diagonal
0768 template <class T, unsigned int D1, unsigned int D2, class R>
0769 template <class Vector>
0770 void SMatrix<T,D1,D2,R>::SetDiagonal( const Vector & v) {
0771
0772    // check size that size of vector is correct
0773    STATIC_CHECK( ( ( D1 <= D2) && Vector::kSize == D1 ) ||
0774                  ( ( D2 < D1 ) && Vector::kSize == D2 ), SVector_size_NOT_correct );
0775
0776
0777    for(unsigned int i=0; i<Vector::kSize; ++i) {
0778       fRep[ i*D2 + i] = v[i];
0779    }
0780 }
0781
0782 // matrix trace
0783 template <class T, unsigned int D1, unsigned int D2, class R>
0784 T SMatrix<T,D1,D2,R>::Trace( ) const {
0785    unsigned int diagSize = D1;
0786    if (D2 < D1) diagSize = D2;
0787    T trace = 0;
0788    for(unsigned int i=0; i< diagSize; ++i) {
0789       trace += fRep[ i*D2 + i] ;
0790    }
0791    return trace;
0792 }
0793
0794 //upper block
0795 template <class T, unsigned int D1, unsigned int D2, class R>
0796 #ifndef UNSUPPORTED_TEMPLATE_EXPRESSION
0797 SVector<T, D1 * (D2 +1)/2 >  SMatrix<T,D1,D2,R>::UpperBlock( ) const {
0798 #else
0799 template <class SubVector>
0800 SubVector SMatrix<T,D1,D2,R>::UpperBlock( ) const {
0801 #endif
0802    // only for squared matrices
0803    STATIC_CHECK( D1 == D2,SMatrix_NOT_square );
0804
0805 #ifndef UNSUPPORTED_TEMPLATE_EXPRESSION
0806    SVector<T, D1 * (D2 +1)/2 >  tmp;
0807 #else
0808    // N must be equal D1 *(D1 +1)/2
0809    STATIC_CHECK( SubVector::kSize == D1*(D1+1)/2,SVector_Wrong_Size );
0810    SubVector tmp;
0811 #endif
0812
0813    int k = 0;
0814    for(unsigned int i=0; i<D1; ++i) {
0815       for(unsigned int j=i; j<D2; ++j) {
0816          tmp[k] = fRep[ i*D2 + j];
0817          k++;
0818       }
0819    }
0820    return tmp;
0821 }
0822
0823 //lower block
0824 template <class T, unsigned int D1, unsigned int D2, class R>
0825 #ifndef UNSUPPORTED_TEMPLATE_EXPRESSION
0826 SVector<T, D1 * (D2 +1)/2 >  SMatrix<T,D1,D2,R>::LowerBlock( ) const {
0827 #else
0828 template <class SubVector>
0829 SubVector SMatrix<T,D1,D2,R>::LowerBlock( ) const {
0830 #endif
0831
0832    // only for squared matrices
0833    STATIC_CHECK( D1 == D2,SMatrix_NOT_square );
0834
0835 #ifndef UNSUPPORTED_TEMPLATE_EXPRESSION
0836    SVector<T, D1 * (D2 +1)/2 >  tmp;
0837 #else
0838    // N must be equal D1 *(D1 +1)/2
0839    STATIC_CHECK( SubVector::kSize == D1*(D1+1)/2,SVector_Wrong_Size );
0840    SubVector tmp;
0841 #endif
0842
0843    int k = 0;
0844    for(unsigned int i=0; i<D1; ++i) {
0845       for(unsigned int j=0; j<=i; ++j) {
0846          tmp[k] = fRep[ i*D2 + j];
0847          k++;
0848       }
0849    }
0850    return tmp;
0851 }
0852
0853 /// construct from upper/lower block
0854
0855 //lower block
0856 template <class T, unsigned int D1, unsigned int D2, class R>
0857 #ifndef UNSUPPORTED_TEMPLATE_EXPRESSION
0858 SMatrix<T,D1,D2,R>::SMatrix(const SVector<T, D1*(D2+1)/2 >  & v, bool lower ) {
0859 #else
0860 template <unsigned int N>
0861 SMatrix<T,D1,D2,R>::SMatrix(const SVector<T, N >  & v, bool lower ) {
0862 #endif
0863
0864    // only for squared matrices
0865    STATIC_CHECK( D1 == D2,SMatrix_NOT_square );
0866
0867 #ifdef UNSUPPORTED_TEMPLATE_EXPRESSION
0868    STATIC_CHECK( N == D1*(D1+1)/2,SVector_Wrong_Size );
0869 #endif
0870
0871    int k = 0;
0872    if (lower) {
0873       // case of lower block
0874       for(unsigned int i=0; i<D1; ++i) {
0875          for(unsigned int j=0; j<=i; ++j) {
0876             fRep[ i*D2 + j] = v[k];
0877             if ( i != j) fRep[ j*D2 + i] = v[k];
0878             k++;
0879          }
0880       }
0881    } else {
0882       // case of upper block
0883       for(unsigned int i=0; i<D1; ++i) {
0884          for(unsigned int j=i; j<D2; ++j) {
0885             fRep[ i*D2 + j] = v[k];
0886             if ( i != j) fRep[ j*D2 + i] = v[k];
0887             k++;
0888          }
0889       }
0890    }
0891 }
0892
0893
0894 template <class T, unsigned int D1, unsigned int D2, class R>
0895 bool SMatrix<T,D1,D2,R>::IsInUse( const T * p) const {
0896    return p == fRep.Array();
0897 }
0898
0899
0900
0901
0902   }  // namespace Math
0903
0904 }  // namespace ROOT
0905
0906
0907
0908 #endif