root/Math/HelperOps.h

0001 // @(#)root/smatrix:$Id$
0002 // Authors: J. Palacios    2006
0003
0004 #ifndef ROOT_Math_HelperOps
0005 #define ROOT_Math_HelperOps 1
0006
0007 // Include files
0008
0009 /** @class HelperOps HelperOps.h Math/HelperOps.h
0010  *
0011  *
0012  *  @author Juan PALACIOS
0013  *  @date   2006-01-11
0014  *
0015  *  Specialised helper classes for binary operators =, +=, -=
0016  *  between SMatrices and Expressions with arbitrary representations.
0017  *  Specialisations at the moment only for Symmetric LHS and Generic RHS
0018  *  and used to throw static assert.
0019  */
0020 #include "Math/StaticCheck.h"
0021 #include <algorithm>  // required by std::copy
0022 #include <cassert>
0023
0024 namespace ROOT {
0025
0026 namespace Math {
0027
0028    template <class T, unsigned int D1, unsigned int D2, class R>
0029    class SMatrix;
0030
0031    template <class A, class T, unsigned int D1, unsigned int D2, class R>
0032    class Expr;
0033
0034    template <class T, unsigned int D>
0035    class MatRepSym;
0036
0037    template <class T, unsigned int D1, unsigned int D2>
0038    class MatRepStd;
0039
0040    //=========================================================================
0041    /**
0042        Structure to assign from an expression based to general matrix to general matrix
0043    */
0044    template <class T,
0045              unsigned int D1, unsigned int D2,
0046              class A, class R1, class R2>
0047
0048    struct Assign
0049    {
0050       /**
0051           Evaluate the expression from general to general matrices.
0052           If the matrix to assign the value is in use in the expression,
0053           a temporary object is created to store the value (case A = B * A)
0054       */
0055       static void Evaluate(SMatrix<T,D1,D2,R1>& lhs,  const Expr<A,T,D1,D2,R2>& rhs)
0056       {
0057          if (! rhs.IsInUse(lhs.begin() )  ) {
0058             unsigned int l = 0;
0059             for(unsigned int i=0; i<D1; ++i)
0060                for(unsigned int j=0; j<D2; ++j) {
0061                   lhs.fRep[l] = rhs(i,j);
0062                   l++;
0063                }
0064          }
0065          // lhs is in use in expression, need to create a temporary with the result
0066          else {
0067             // std::cout << "create temp  for " << typeid(rhs).name() << std::endl;
0068             T tmp[D1*D2];
0069             unsigned int l = 0;
0070             for(unsigned int i=0; i<D1; ++i)
0071                for(unsigned int j=0; j<D2; ++j) {
0072                   tmp[l] = rhs(i,j);
0073                   l++;
0074                }
0075             // copy now the temp object
0076             for(unsigned int i=0; i<D1*D2; ++i) lhs.fRep[i] = tmp[i];
0077          }
0078
0079       }
0080
0081    };
0082
0083    /**
0084        Structure to assign from an expression based to symmetric matrix to symmetric matrix
0085    */
0086    template <class T,
0087              unsigned int D1, unsigned int D2,
0088              class A>
0089
0090    struct Assign<T, D1, D2, A, MatRepSym<T,D1>, MatRepSym<T,D1> >
0091    {
0092       /**
0093           Evaluate the expression from  symmetric to symmetric matrices.
0094           If the matrix to assign the value is in use in the expression,
0095           a temporary object is created to store the value (case A = B * A)
0096       */
0097       static void Evaluate(SMatrix<T,D1,D2,MatRepSym<T,D1> >& lhs,
0098                            const Expr<A,T,D1,D2,MatRepSym<T,D1> >& rhs)
0099       {
0100          if (! rhs.IsInUse(lhs.begin() ) ) {
0101             unsigned int l = 0;
0102             for(unsigned int i=0; i<D1; ++i)
0103                // storage of symmetric matrix is in lower block
0104                for(unsigned int j=0; j<=i; ++j) {
0105                   lhs.fRep.Array()[l] = rhs(i,j);
0106                   l++;
0107                }
0108          }
0109          // create a temporary object to store result
0110          else {
0111             T tmp[MatRepSym<T,D1>::kSize];
0112             unsigned int l = 0;
0113             for(unsigned int i=0; i<D1; ++i)
0114                for(unsigned int j=0; j<=i; ++j) {
0115                   tmp[l] = rhs(i,j);
0116                   l++;
0117                }
0118             // copy now the temp object
0119             for(unsigned int i=0; i<MatRepSym<T,D1>::kSize; ++i) lhs.fRep.Array()[i] = tmp[i];
0120          }
0121       }
0122    };
0123
0124
0125
0126    /**
0127        Dummy Structure which flags an error to avoid assignment from expression based on a
0128        general matrix to a symmetric matrix
0129    */
0130    template <class T, unsigned int D1, unsigned int D2, class A>
0131    struct Assign<T, D1, D2, A, MatRepSym<T,D1>, MatRepStd<T,D1,D2> >
0132    {
0133       static void Evaluate(SMatrix<T,D1,D2,MatRepSym<T,D1> >&,
0134                            const Expr<A,T,D1,D2,MatRepStd<T,D1,D2> >&)
0135       {
0136          STATIC_CHECK(0==1, Cannot_assign_general_to_symmetric_matrix);
0137       }
0138
0139    }; // struct Assign
0140
0141
0142    /**
0143        Force Expression evaluation from general to symmetric.
0144        To be used when is known (like in similarity products) that the result
0145        is symmetric
0146        Note this is function used in the simmilarity product: no check for temporary is
0147        done since in that case is not needed
0148    */
0149    struct AssignSym
0150    {
0151       /// assign a symmetric matrix from an expression
0152       template <class T,
0153                 unsigned int D,
0154                 class A,
0155                 class R>
0156       static void Evaluate(SMatrix<T,D,D,MatRepSym<T,D> >& lhs,  const Expr<A,T,D,D,R>& rhs)
0157       {
0158          //for(unsigned int i=0; i<D1*D2; ++i) lhs.fRep[i] = rhs.apply(i);
0159          unsigned int l = 0;
0160          for(unsigned int i=0; i<D; ++i)
0161             // storage of symmetric matrix is in lower block
0162             for(unsigned int j=0; j<=i; ++j) {
0163                lhs.fRep.Array()[l] = rhs(i,j);
0164                l++;
0165             }
0166       }
0167       /// assign the symmetric matric  from a general matrix
0168       template <class T,
0169                 unsigned int D,
0170                 class R>
0171       static void Evaluate(SMatrix<T,D,D,MatRepSym<T,D> >& lhs,  const SMatrix<T,D,D,R>& rhs)
0172       {
0173          //for(unsigned int i=0; i<D1*D2; ++i) lhs.fRep[i] = rhs.apply(i);
0174          unsigned int l = 0;
0175          for(unsigned int i=0; i<D; ++i)
0176             // storage of symmetric matrix is in lower block
0177             for(unsigned int j=0; j<=i; ++j) {
0178                lhs.fRep.Array()[l] = rhs(i,j);
0179                l++;
0180             }
0181       }
0182
0183
0184    }; // struct AssignSym
0185
0186
0187    //=========================================================================
0188    /**
0189        Evaluate the expression performing a += operation
0190        Need to check whether creating a temporary object with the expression result
0191        (like in op:  A += A * B )
0192    */
0193    template <class T, unsigned int D1, unsigned int D2, class A,
0194              class R1, class R2>
0195    struct PlusEquals
0196    {
0197       static void Evaluate(SMatrix<T,D1,D2,R1>& lhs,  const Expr<A,T,D1,D2,R2>& rhs)
0198       {
0199          if (! rhs.IsInUse(lhs.begin() )  ) {
0200             unsigned int l = 0;
0201             for(unsigned int i=0; i<D1; ++i)
0202                for(unsigned int j=0; j<D2; ++j) {
0203                   lhs.fRep[l] += rhs(i,j);
0204                   l++;
0205                }
0206          }
0207          else {
0208             T tmp[D1*D2];
0209             unsigned int l = 0;
0210             for(unsigned int i=0; i<D1; ++i)
0211                for(unsigned int j=0; j<D2; ++j) {
0212                   tmp[l] = rhs(i,j);
0213                   l++;
0214                }
0215             // += now using the temp object
0216             for(unsigned int i=0; i<D1*D2; ++i) lhs.fRep[i] += tmp[i];
0217          }
0218       }
0219    };
0220
0221    /**
0222        Specialization for symmetric matrices
0223        Evaluate the expression performing a += operation for symmetric matrices
0224        Need to have a separate functions to avoid to modify two times the off-diagonal
0225        elements (i.e applying two times the expression)
0226        Need to check whether creating a temporary object with the expression result
0227        (like in op:  A += A * B )
0228    */
0229    template <class T,
0230              unsigned int D1, unsigned int D2,
0231              class A>
0232    struct PlusEquals<T, D1, D2, A, MatRepSym<T,D1>, MatRepSym<T,D1> >
0233    {
0234       static void Evaluate(SMatrix<T,D1,D2,MatRepSym<T,D1> >& lhs,  const Expr<A,T,D1,D2, MatRepSym<T,D1> >& rhs)
0235       {
0236          if (! rhs.IsInUse(lhs.begin() )  ) {
0237             unsigned int l = 0;  // l span storage of sym matrices
0238             for(unsigned int i=0; i<D1; ++i)
0239                for(unsigned int j=0; j<=i; ++j) {
0240                   lhs.fRep.Array()[l] += rhs(i,j);
0241                   l++;
0242                }
0243          }
0244          else {
0245             T tmp[MatRepSym<T,D1>::kSize];
0246             unsigned int l = 0;
0247             for(unsigned int i=0; i<D1; ++i)
0248                for(unsigned int j=0; j<=i; ++j) {
0249                   tmp[l] = rhs(i,j);
0250                   l++;
0251                }
0252             // += now using the temp object
0253             for(unsigned int i=0; i<MatRepSym<T,D1>::kSize; ++i) lhs.fRep.Array()[i] += tmp[i];
0254          }
0255       }
0256    };
0257    /**
0258       Specialization for symmetrix += general : NOT Allowed operation
0259     */
0260    template <class T, unsigned int D1, unsigned int D2, class A>
0261    struct PlusEquals<T, D1, D2, A, MatRepSym<T,D1>, MatRepStd<T,D1,D2> >
0262    {
0263       static void Evaluate(SMatrix<T,D1,D2,MatRepSym<T,D1> >&,
0264                            const Expr<A,T,D1,D2,MatRepStd<T,D1,D2> >&)
0265       {
0266          STATIC_CHECK(0==1, Cannot_plusEqual_general_to_symmetric_matrix);
0267       }
0268    }; // struct PlusEquals
0269
0270    //=========================================================================
0271
0272    /**
0273        Evaluate the expression performing a -= operation
0274        Need to check whether creating a temporary object with the expression result
0275        (like in op:  A -= A * B )
0276    */
0277    template <class T, unsigned int D1, unsigned int D2, class A,
0278              class R1, class R2>
0279    struct MinusEquals
0280    {
0281       static void Evaluate(SMatrix<T,D1,D2,R1>& lhs,  const Expr<A,T,D1,D2,R2>& rhs)
0282       {
0283          if (! rhs.IsInUse(lhs.begin() )  ) {
0284             unsigned int l = 0;
0285             for(unsigned int i=0; i<D1; ++i)
0286                for(unsigned int j=0; j<D2; ++j) {
0287                   lhs.fRep[l] -= rhs(i,j);
0288                   l++;
0289                }
0290          }
0291          else {
0292             T tmp[D1*D2];
0293             unsigned int l = 0;
0294             for(unsigned int i=0; i<D1; ++i)
0295                for(unsigned int j=0; j<D2; ++j) {
0296                   tmp[l] = rhs(i,j);
0297                   l++;
0298                }
0299             // -= now using the temp object
0300             for(unsigned int i=0; i<D1*D2; ++i) lhs.fRep[i] -= tmp[i];
0301          }
0302       }
0303    };
0304    /**
0305        Specialization for symmetric matrices.
0306        Evaluate the expression performing a -= operation for symmetric matrices
0307        Need to have a separate functions to avoid to modify two times the off-diagonal
0308        elements (i.e applying two times the expression)
0309        Need to check whether creating a temporary object with the expression result
0310        (like in op:  A -= A + B )
0311    */
0312    template <class T,
0313              unsigned int D1, unsigned int D2,
0314              class A>
0315    struct MinusEquals<T, D1, D2, A, MatRepSym<T,D1>, MatRepSym<T,D1> >
0316    {
0317       static void Evaluate(SMatrix<T,D1,D2,MatRepSym<T,D1> >& lhs,  const Expr<A,T,D1,D2, MatRepSym<T,D1> >& rhs)
0318       {
0319          if (! rhs.IsInUse(lhs.begin() )  ) {
0320             unsigned int l = 0;  // l span storage of sym matrices
0321             for(unsigned int i=0; i<D1; ++i)
0322                for(unsigned int j=0; j<=i; ++j) {
0323                   lhs.fRep.Array()[l] -= rhs(i,j);
0324                   l++;
0325                }
0326          }
0327          else {
0328             T tmp[MatRepSym<T,D1>::kSize];
0329             unsigned int l = 0;
0330             for(unsigned int i=0; i<D1; ++i)
0331                for(unsigned int j=0; j<=i; ++j) {
0332                   tmp[l] = rhs(i,j);
0333                   l++;
0334                }
0335             // -= now using the temp object
0336             for(unsigned int i=0; i<MatRepSym<T,D1>::kSize; ++i) lhs.fRep.Array()[i] -= tmp[i];
0337          }
0338       }
0339    };
0340
0341    /**
0342       Specialization for symmetrix -= general : NOT Allowed operation
0343     */
0344    template <class T, unsigned int D1, unsigned int D2, class A>
0345    struct MinusEquals<T, D1, D2, A, MatRepSym<T,D1>, MatRepStd<T,D1,D2> >
0346    {
0347       static void Evaluate(SMatrix<T,D1,D2,MatRepSym<T,D1> >&,
0348                            const Expr<A,T,D1,D2,MatRepStd<T,D1,D2> >&)
0349       {
0350          STATIC_CHECK(0==1, Cannot_minusEqual_general_to_symmetric_matrix);
0351       }
0352    }; // struct MinusEquals
0353
0354
0355    /** Structure to deal when a submatrix is placed in a matrix.
0356        We have different cases according to the matrix representation
0357    */
0358    template <class T, unsigned int D1, unsigned int D2,
0359              unsigned int D3, unsigned int D4,
0360              class R1, class R2>
0361    struct PlaceMatrix
0362    {
0363       static void Evaluate(SMatrix<T,D1,D2,R1>& lhs,  const SMatrix<T,D3,D4,R2>& rhs,
0364                            unsigned int row, unsigned int col) {
0365
0366          assert(row+D3 <= D1 && col+D4 <= D2);
0367          const unsigned int offset = row*D2+col;
0368
0369          for(unsigned int i=0; i<D3*D4; ++i) {
0370             lhs.fRep[offset+(i/D4)*D2+i%D4] = rhs.apply(i);
0371          }
0372
0373       }
0374    }; // struct PlaceMatrix
0375
0376    template <class T, unsigned int D1, unsigned int D2,
0377              unsigned int D3, unsigned int D4,
0378              class A, class R1, class R2>
0379    struct PlaceExpr {
0380       static void Evaluate(SMatrix<T,D1,D2,R1>& lhs,  const Expr<A,T,D3,D4,R2>& rhs,
0381                            unsigned int row, unsigned int col) {
0382
0383          assert(row+D3 <= D1 && col+D4 <= D2);
0384          const unsigned int offset = row*D2+col;
0385
0386          for(unsigned int i=0; i<D3*D4; ++i) {
0387             lhs.fRep[offset+(i/D4)*D2+i%D4] = rhs.apply(i);
0388          }
0389       }
0390    };  // struct PlaceExpr
0391
0392    // specialization for general matrix in symmetric matrices
0393    template <class T, unsigned int D1, unsigned int D2,
0394              unsigned int D3, unsigned int D4 >
0395    struct PlaceMatrix<T, D1, D2, D3, D4, MatRepSym<T,D1>, MatRepStd<T,D3,D4> > {
0396       static void Evaluate(SMatrix<T,D1,D2,MatRepSym<T,D1> >& ,
0397                            const SMatrix<T,D3,D4,MatRepStd<T,D3,D4> >& ,
0398                            unsigned int , unsigned int )
0399       {
0400          STATIC_CHECK(0==1, Cannot_Place_Matrix_general_in_symmetric_matrix);
0401       }
0402    }; // struct PlaceMatrix
0403
0404    // specialization for general expression in symmetric matrices
0405    template <class T, unsigned int D1, unsigned int D2,
0406              unsigned int D3, unsigned int D4, class A >
0407    struct PlaceExpr<T, D1, D2, D3, D4, A, MatRepSym<T,D1>, MatRepStd<T,D3,D4> > {
0408       static void Evaluate(SMatrix<T,D1,D2,MatRepSym<T,D1> >& ,
0409                            const Expr<A,T,D3,D4,MatRepStd<T,D3,D4> >& ,
0410                            unsigned int , unsigned int )
0411       {
0412          STATIC_CHECK(0==1, Cannot_Place_Matrix_general_in_symmetric_matrix);
0413       }
0414    }; // struct PlaceExpr
0415
0416    // specialization for symmetric matrix in symmetric matrices
0417
0418    template <class T, unsigned int D1, unsigned int D2,
0419              unsigned int D3, unsigned int D4 >
0420    struct PlaceMatrix<T, D1, D2, D3, D4, MatRepSym<T,D1>, MatRepSym<T,D3> > {
0421       static void Evaluate(SMatrix<T,D1,D2,MatRepSym<T,D1> >& lhs,
0422                            const SMatrix<T,D3,D4,MatRepSym<T,D3> >& rhs,
0423                            unsigned int row, unsigned int col )
0424       {
0425          // can work only if placed on the diagonal
0426          assert(row == col);
0427
0428          for(unsigned int i=0; i<D3; ++i) {
0429             for(unsigned int j=0; j<=i; ++j)
0430                lhs.fRep(row+i,col+j) = rhs(i,j);
0431          }
0432       }
0433    }; // struct PlaceMatrix
0434
0435    // specialization for symmetric expression in symmetric matrices
0436    template <class T, unsigned int D1, unsigned int D2,
0437              unsigned int D3, unsigned int D4, class A >
0438    struct PlaceExpr<T, D1, D2, D3, D4, A, MatRepSym<T,D1>, MatRepSym<T,D3> > {
0439       static void Evaluate(SMatrix<T,D1,D2,MatRepSym<T,D1> >& lhs,
0440                            const Expr<A,T,D3,D4,MatRepSym<T,D3> >& rhs,
0441                            unsigned int row, unsigned int col )
0442       {
0443          // can work only if placed on the diagonal
0444          assert(row == col);
0445
0446          for(unsigned int i=0; i<D3; ++i) {
0447             for(unsigned int j=0; j<=i; ++j)
0448                lhs.fRep(row+i,col+j) = rhs(i,j);
0449          }
0450       }
0451    }; // struct PlaceExpr
0452
0453
0454
0455    /** Structure for getting sub matrices
0456        We have different cases according to the matrix representations
0457    */
0458    template <class T, unsigned int D1, unsigned int D2,
0459              unsigned int D3, unsigned int D4,
0460              class R1, class R2>
0461    struct RetrieveMatrix
0462    {
0463       static void Evaluate(SMatrix<T,D1,D2,R1>& lhs,  const SMatrix<T,D3,D4,R2>& rhs,
0464                            unsigned int row, unsigned int col) {
0465          STATIC_CHECK( D1 <= D3,Smatrix_nrows_too_small);
0466          STATIC_CHECK( D2 <= D4,Smatrix_ncols_too_small);
0467
0468          assert(row + D1 <= D3);
0469          assert(col + D2 <= D4);
0470
0471          for(unsigned int i=0; i<D1; ++i) {
0472             for(unsigned int j=0; j<D2; ++j)
0473                lhs(i,j) = rhs(i+row,j+col);
0474          }
0475       }
0476    };   // struct RetrieveMatrix
0477
0478    // specialization for getting symmetric matrices from  general matrices (MUST fail)
0479    template <class T, unsigned int D1, unsigned int D2,
0480              unsigned int D3, unsigned int D4 >
0481    struct RetrieveMatrix<T, D1, D2, D3, D4, MatRepSym<T,D1>, MatRepStd<T,D3,D4> > {
0482       static void Evaluate(SMatrix<T,D1,D2,MatRepSym<T,D1> >& ,
0483                            const SMatrix<T,D3,D4,MatRepStd<T,D3,D4> >& ,
0484                            unsigned int , unsigned int )
0485       {
0486          STATIC_CHECK(0==1, Cannot_Sub_Matrix_symmetric_in_general_matrix);
0487       }
0488    }; // struct RetrieveMatrix
0489
0490    // specialization for getting symmetric matrices from  symmetric matrices (OK if row == col)
0491    template <class T, unsigned int D1, unsigned int D2,
0492              unsigned int D3, unsigned int D4 >
0493    struct RetrieveMatrix<T, D1, D2, D3, D4, MatRepSym<T,D1>, MatRepSym<T,D3> > {
0494       static void Evaluate(SMatrix<T,D1,D2,MatRepSym<T,D1> >& lhs,
0495                            const SMatrix<T,D3,D4,MatRepSym<T,D3> >& rhs,
0496                            unsigned int row, unsigned int col )
0497       {
0498          STATIC_CHECK(  D1 <= D3,Smatrix_dimension1_too_small);
0499          // can work only if placed on the diagonal
0500          assert(row == col);
0501          assert(row + D1 <= D3);
0502
0503          for(unsigned int i=0; i<D1; ++i) {
0504             for(unsigned int j=0; j<=i; ++j)
0505                lhs(i,j) = rhs(i+row,j+col );
0506          }
0507       }
0508
0509    }; // struct RetrieveMatrix
0510
0511    /**
0512       Structure for assignment to a general matrix from iterator.
0513       Optionally a check is done that iterator size
0514       is not larger than matrix size
0515     */
0516    template <class T, unsigned int D1, unsigned int D2, class R>
0517    struct AssignItr {
0518       template<class Iterator>
0519       static void Evaluate(SMatrix<T,D1,D2,R>& lhs, Iterator begin, Iterator end,
0520                            bool triang, bool lower,bool check=true) {
0521          // require size match exactly (better)
0522
0523          if (triang) {
0524             Iterator itr = begin;
0525             if (lower) {
0526                for (unsigned int i = 0; i < D1; ++i)
0527                   for (unsigned int j =0; j <= i; ++j) {
0528                      // we assume iterator is well bounded within matrix
0529                      lhs.fRep[i*D2+j] = *itr++;
0530                   }
0531
0532             }
0533             else { // upper
0534                for (unsigned int i = 0; i < D1; ++i)
0535                   for (unsigned int j = i; j <D2; ++j) {
0536                      if (itr != end)
0537                         lhs.fRep[i*D2+j] = *itr++;
0538                      else
0539                         return;
0540                   }
0541
0542             }
0543          }
0544          // case of filling the full matrix
0545          else {
0546             if (check) assert( begin + R::kSize == end);
0547             // copy directly the elements
0548             std::copy(begin, end, lhs.fRep.Array() );
0549          }
0550       }
0551
0552    }; // struct AssignItr
0553
0554    /**
0555       Specialized structure for assignment to a symmetrix matrix from iterator.
0556       Optionally a check is done that iterator size
0557       is the same as the matrix size
0558     */
0559    template <class T, unsigned int D1, unsigned int D2>
0560    struct AssignItr<T, D1, D2, MatRepSym<T,D1> >  {
0561       template<class Iterator>
0562       static void Evaluate(SMatrix<T,D1,D2,MatRepSym<T,D1> >& lhs, Iterator begin, Iterator end, bool , bool lower, bool check = true) {
0563
0564          if (lower) {
0565             if (check) {
0566                assert(begin+ static_cast< int>( MatRepSym<T,D1>::kSize) == end);
0567             }
0568             std::copy(begin, end, lhs.fRep.Array() );
0569          }
0570          else {
0571             Iterator itr = begin;
0572             for (unsigned int i = 0; i < D1; ++i)
0573                for (unsigned int j = i; j <D2; ++j) {
0574                   if (itr != end)
0575                      lhs(i,j) = *itr++;
0576                   else
0577                      return;
0578                }
0579          }
0580       }
0581
0582    }; // struct AssignItr
0583
0584
0585 }  // namespace Math
0586
0587 }  // namespace ROOT
0588
0589 #endif // MATH_HELPEROPS_H