HUMIM Tracker: /home/hauberg/Dokumenter/Capture/humim-tracker-0.1/src/ntk/geometry/Eigen/src/Core/products/TriangularMatrixVector.h Source File

Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 //
00006 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #ifndef EIGEN_TRIANGULARMATRIXVECTOR_H
00026 #define EIGEN_TRIANGULARMATRIXVECTOR_H
00027 
00028 namespace internal {
00029 
00030 template<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar, bool ConjRhs, int StorageOrder>
00031 struct product_triangular_matrix_vector;
00032 
00033 template<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar, bool ConjRhs>
00034 struct product_triangular_matrix_vector<Index,Mode,LhsScalar,ConjLhs,RhsScalar,ConjRhs,ColMajor>
00035 {
00036   typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar;
00037   enum {
00038     IsLower = ((Mode&Lower)==Lower),
00039     HasUnitDiag = (Mode & UnitDiag)==UnitDiag
00040   };
00041   static EIGEN_DONT_INLINE  void run(Index rows, Index cols, const LhsScalar* _lhs, Index lhsStride,
00042                                      const RhsScalar* _rhs, Index rhsIncr, ResScalar* _res, Index resIncr, ResScalar alpha)
00043   {
00044     EIGEN_UNUSED_VARIABLE(resIncr);
00045     eigen_assert(resIncr==1);
00046     
00047     static const Index PanelWidth = EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH;
00048 
00049     typedef Map<const Matrix<LhsScalar,Dynamic,Dynamic,ColMajor>, 0, OuterStride<> > LhsMap;
00050     const LhsMap lhs(_lhs,rows,cols,OuterStride<>(lhsStride));
00051     typename conj_expr_if<ConjLhs,LhsMap>::type cjLhs(lhs);
00052     
00053     typedef Map<const Matrix<RhsScalar,Dynamic,1>, 0, InnerStride<> > RhsMap;
00054     const RhsMap rhs(_rhs,cols,InnerStride<>(rhsIncr));
00055     typename conj_expr_if<ConjRhs,RhsMap>::type cjRhs(rhs);
00056 
00057     typedef Map<Matrix<ResScalar,Dynamic,1> > ResMap;
00058     ResMap res(_res,rows);
00059 
00060     for (Index pi=0; pi<cols; pi+=PanelWidth)
00061     {
00062       Index actualPanelWidth = std::min(PanelWidth, cols-pi);
00063       for (Index k=0; k<actualPanelWidth; ++k)
00064       {
00065         Index i = pi + k;
00066         Index s = IsLower ? (HasUnitDiag ? i+1 : i ) : pi;
00067         Index r = IsLower ? actualPanelWidth-k : k+1;
00068         if ((!HasUnitDiag) || (--r)>0)
00069           res.segment(s,r) += (alpha * cjRhs.coeff(i)) * cjLhs.col(i).segment(s,r);
00070         if (HasUnitDiag)
00071           res.coeffRef(i) += alpha * cjRhs.coeff(i);
00072       }
00073       Index r = IsLower ? cols - pi - actualPanelWidth : pi;
00074       if (r>0)
00075       {
00076         Index s = IsLower ? pi+actualPanelWidth : 0;
00077         general_matrix_vector_product<Index,LhsScalar,ColMajor,ConjLhs,RhsScalar,ConjRhs>::run(
00078             r, actualPanelWidth,
00079             &lhs.coeffRef(s,pi), lhsStride,
00080             &rhs.coeffRef(pi), rhsIncr,
00081             &res.coeffRef(s), resIncr, alpha);
00082       }
00083     }
00084   }
00085 };
00086 
00087 template<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar, bool ConjRhs>
00088 struct product_triangular_matrix_vector<Index,Mode,LhsScalar,ConjLhs,RhsScalar,ConjRhs,RowMajor>
00089 {
00090   typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar;
00091   enum {
00092     IsLower = ((Mode&Lower)==Lower),
00093     HasUnitDiag = (Mode & UnitDiag)==UnitDiag
00094   };
00095   static void run(Index rows, Index cols, const LhsScalar* _lhs, Index lhsStride,
00096                   const RhsScalar* _rhs, Index rhsIncr, ResScalar* _res, Index resIncr, ResScalar alpha)
00097   {
00098     eigen_assert(rhsIncr==1);
00099     EIGEN_UNUSED_VARIABLE(rhsIncr);
00100     
00101     static const Index PanelWidth = EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH;
00102 
00103     typedef Map<const Matrix<LhsScalar,Dynamic,Dynamic,RowMajor>, 0, OuterStride<> > LhsMap;
00104     const LhsMap lhs(_lhs,rows,cols,OuterStride<>(lhsStride));
00105     typename conj_expr_if<ConjLhs,LhsMap>::type cjLhs(lhs);
00106 
00107     typedef Map<const Matrix<RhsScalar,Dynamic,1> > RhsMap;
00108     const RhsMap rhs(_rhs,cols);
00109     typename conj_expr_if<ConjRhs,RhsMap>::type cjRhs(rhs);
00110 
00111     typedef Map<Matrix<ResScalar,Dynamic,1>, 0, InnerStride<> > ResMap;
00112     ResMap res(_res,rows,InnerStride<>(resIncr));
00113     
00114     for (Index pi=0; pi<cols; pi+=PanelWidth)
00115     {
00116       Index actualPanelWidth = std::min(PanelWidth, cols-pi);
00117       for (Index k=0; k<actualPanelWidth; ++k)
00118       {
00119         Index i = pi + k;
00120         Index s = IsLower ? pi  : (HasUnitDiag ? i+1 : i);
00121         Index r = IsLower ? k+1 : actualPanelWidth-k;
00122         if ((!HasUnitDiag) || (--r)>0)
00123           res.coeffRef(i) += alpha * (cjLhs.row(i).segment(s,r).cwiseProduct(cjRhs.segment(s,r).transpose())).sum();
00124         if (HasUnitDiag)
00125           res.coeffRef(i) += alpha * cjRhs.coeff(i);
00126       }
00127       Index r = IsLower ? pi : cols - pi - actualPanelWidth;
00128       if (r>0)
00129       {
00130         Index s = IsLower ? 0 : pi + actualPanelWidth;
00131         general_matrix_vector_product<Index,LhsScalar,RowMajor,ConjLhs,RhsScalar,ConjRhs>::run(
00132             actualPanelWidth, r,
00133             &lhs.coeffRef(pi,s), lhsStride,
00134             &rhs.coeffRef(s), rhsIncr,
00135             &res.coeffRef(pi), resIncr, alpha);
00136       }
00137     }
00138   }
00139 };
00140 
00141 /***************************************************************************
00142 * Wrapper to product_triangular_vector
00143 ***************************************************************************/
00144 
00145 template<int Mode, bool LhsIsTriangular, typename Lhs, typename Rhs>
00146 struct traits<TriangularProduct<Mode,LhsIsTriangular,Lhs,false,Rhs,true> >
00147  : traits<ProductBase<TriangularProduct<Mode,LhsIsTriangular,Lhs,false,Rhs,true>, Lhs, Rhs> >
00148 {};
00149 
00150 template<int Mode, bool LhsIsTriangular, typename Lhs, typename Rhs>
00151 struct traits<TriangularProduct<Mode,LhsIsTriangular,Lhs,true,Rhs,false> >
00152  : traits<ProductBase<TriangularProduct<Mode,LhsIsTriangular,Lhs,true,Rhs,false>, Lhs, Rhs> >
00153 {};
00154 
00155 
00156 template<int StorageOrder>
00157 struct trmv_selector;
00158 
00159 } // end namespace internal
00160 
00161 template<int Mode, typename Lhs, typename Rhs>
00162 struct TriangularProduct<Mode,true,Lhs,false,Rhs,true>
00163   : public ProductBase<TriangularProduct<Mode,true,Lhs,false,Rhs,true>, Lhs, Rhs >
00164 {
00165   EIGEN_PRODUCT_PUBLIC_INTERFACE(TriangularProduct)
00166 
00167   TriangularProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {}
00168 
00169   template<typename Dest> void scaleAndAddTo(Dest& dst, Scalar alpha) const
00170   {
00171     eigen_assert(dst.rows()==m_lhs.rows() && dst.cols()==m_rhs.cols());
00172   
00173     internal::trmv_selector<(int(internal::traits<Lhs>::Flags)&RowMajorBit) ? RowMajor : ColMajor>::run(*this, dst, alpha);
00174   }
00175 };
00176 
00177 template<int Mode, typename Lhs, typename Rhs>
00178 struct TriangularProduct<Mode,false,Lhs,true,Rhs,false>
00179   : public ProductBase<TriangularProduct<Mode,false,Lhs,true,Rhs,false>, Lhs, Rhs >
00180 {
00181   EIGEN_PRODUCT_PUBLIC_INTERFACE(TriangularProduct)
00182 
00183   TriangularProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {}
00184 
00185   template<typename Dest> void scaleAndAddTo(Dest& dst, Scalar alpha) const
00186   {
00187     eigen_assert(dst.rows()==m_lhs.rows() && dst.cols()==m_rhs.cols());
00188     
00189     typedef TriangularProduct<(Mode & UnitDiag) | ((Mode & Lower) ? Upper : Lower),true,Transpose<const Rhs>,false,Transpose<const Lhs>,true> TriangularProductTranspose;
00190     Transpose<Dest> dstT(dst);
00191     internal::trmv_selector<(int(internal::traits<Rhs>::Flags)&RowMajorBit) ? ColMajor : RowMajor>::run(
00192       TriangularProductTranspose(m_rhs.transpose(),m_lhs.transpose()), dstT, alpha);
00193   }
00194 };
00195 
00196 namespace internal {
00197 
00198 // TODO: find a way to factorize this piece of code with gemv_selector since the logic is exactly the same.
00199   
00200 template<> struct trmv_selector<ColMajor>
00201 {
00202   template<int Mode, typename Lhs, typename Rhs, typename Dest>
00203   static void run(const TriangularProduct<Mode,true,Lhs,false,Rhs,true>& prod, Dest& dest, typename TriangularProduct<Mode,true,Lhs,false,Rhs,true>::Scalar alpha)
00204   {
00205     typedef TriangularProduct<Mode,true,Lhs,false,Rhs,true> ProductType;
00206     typedef typename ProductType::Index Index;
00207     typedef typename ProductType::LhsScalar   LhsScalar;
00208     typedef typename ProductType::RhsScalar   RhsScalar;
00209     typedef typename ProductType::Scalar      ResScalar;
00210     typedef typename ProductType::RealScalar  RealScalar;
00211     typedef typename ProductType::ActualLhsType ActualLhsType;
00212     typedef typename ProductType::ActualRhsType ActualRhsType;
00213     typedef typename ProductType::LhsBlasTraits LhsBlasTraits;
00214     typedef typename ProductType::RhsBlasTraits RhsBlasTraits;
00215     typedef Map<Matrix<ResScalar,Dynamic,1>, Aligned> MappedDest;
00216 
00217     const ActualLhsType actualLhs = LhsBlasTraits::extract(prod.lhs());
00218     const ActualRhsType actualRhs = RhsBlasTraits::extract(prod.rhs());
00219 
00220     ResScalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(prod.lhs())
00221                                   * RhsBlasTraits::extractScalarFactor(prod.rhs());
00222 
00223     enum {
00224       // FIXME find a way to allow an inner stride on the result if packet_traits<Scalar>::size==1
00225       // on, the other hand it is good for the cache to pack the vector anyways...
00226       EvalToDestAtCompileTime = Dest::InnerStrideAtCompileTime==1,
00227       ComplexByReal = (NumTraits<LhsScalar>::IsComplex) && (!NumTraits<RhsScalar>::IsComplex),
00228       MightCannotUseDest = (Dest::InnerStrideAtCompileTime!=1) || ComplexByReal
00229     };
00230 
00231     gemv_static_vector_if<ResScalar,Dest::SizeAtCompileTime,Dest::MaxSizeAtCompileTime,MightCannotUseDest> static_dest;
00232 
00233     bool alphaIsCompatible = (!ComplexByReal) || (imag(actualAlpha)==RealScalar(0));
00234     bool evalToDest = EvalToDestAtCompileTime && alphaIsCompatible;
00235     
00236     RhsScalar compatibleAlpha = get_factor<ResScalar,RhsScalar>::run(actualAlpha);
00237 
00238     ResScalar* actualDestPtr;
00239     bool freeDestPtr = false;
00240     if (evalToDest)
00241     {
00242       actualDestPtr = dest.data();
00243     }
00244     else
00245     {
00246       #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN
00247       int size = dest.size();
00248       EIGEN_DENSE_STORAGE_CTOR_PLUGIN
00249       #endif
00250       if((actualDestPtr = static_dest.data())==0)
00251       {
00252         freeDestPtr = true;
00253         actualDestPtr = ei_aligned_stack_new(ResScalar,dest.size());
00254       }
00255       if(!alphaIsCompatible)
00256       {
00257         MappedDest(actualDestPtr, dest.size()).setZero();
00258         compatibleAlpha = RhsScalar(1);
00259       }
00260       else
00261         MappedDest(actualDestPtr, dest.size()) = dest;
00262     }
00263     
00264     internal::product_triangular_matrix_vector
00265       <Index,Mode,
00266        LhsScalar, LhsBlasTraits::NeedToConjugate,
00267        RhsScalar, RhsBlasTraits::NeedToConjugate,
00268        ColMajor>
00269       ::run(actualLhs.rows(),actualLhs.cols(),
00270             actualLhs.data(),actualLhs.outerStride(),
00271             actualRhs.data(),actualRhs.innerStride(),
00272             actualDestPtr,1,compatibleAlpha);
00273 
00274     if (!evalToDest)
00275     {
00276       if(!alphaIsCompatible)
00277         dest += actualAlpha * MappedDest(actualDestPtr, dest.size());
00278       else
00279         dest = MappedDest(actualDestPtr, dest.size());
00280       if(freeDestPtr) ei_aligned_stack_delete(ResScalar, actualDestPtr, dest.size());
00281     }
00282   }
00283 };
00284 
00285 template<> struct trmv_selector<RowMajor>
00286 {
00287   template<int Mode, typename Lhs, typename Rhs, typename Dest>
00288   static void run(const TriangularProduct<Mode,true,Lhs,false,Rhs,true>& prod, Dest& dest, typename TriangularProduct<Mode,true,Lhs,false,Rhs,true>::Scalar alpha)
00289   {
00290     typedef TriangularProduct<Mode,true,Lhs,false,Rhs,true> ProductType;
00291     typedef typename ProductType::LhsScalar LhsScalar;
00292     typedef typename ProductType::RhsScalar RhsScalar;
00293     typedef typename ProductType::Scalar    ResScalar;
00294     typedef typename ProductType::Index Index;
00295     typedef typename ProductType::ActualLhsType ActualLhsType;
00296     typedef typename ProductType::ActualRhsType ActualRhsType;
00297     typedef typename ProductType::_ActualRhsType _ActualRhsType;
00298     typedef typename ProductType::LhsBlasTraits LhsBlasTraits;
00299     typedef typename ProductType::RhsBlasTraits RhsBlasTraits;
00300 
00301     typename add_const<ActualLhsType>::type actualLhs = LhsBlasTraits::extract(prod.lhs());
00302     typename add_const<ActualRhsType>::type actualRhs = RhsBlasTraits::extract(prod.rhs());
00303 
00304     ResScalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(prod.lhs())
00305                                   * RhsBlasTraits::extractScalarFactor(prod.rhs());
00306 
00307     enum {
00308       DirectlyUseRhs = _ActualRhsType::InnerStrideAtCompileTime==1
00309     };
00310 
00311     gemv_static_vector_if<RhsScalar,_ActualRhsType::SizeAtCompileTime,_ActualRhsType::MaxSizeAtCompileTime,!DirectlyUseRhs> static_rhs;
00312 
00313     RhsScalar* actualRhsPtr;
00314     bool freeRhsPtr = false;
00315     if (DirectlyUseRhs)
00316     {
00317       actualRhsPtr = const_cast<RhsScalar*>(actualRhs.data());
00318     }
00319     else
00320     {
00321       #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN
00322       int size = actualRhs.size();
00323       EIGEN_DENSE_STORAGE_CTOR_PLUGIN
00324       #endif
00325       if((actualRhsPtr = static_rhs.data())==0)
00326       {
00327         freeRhsPtr = true;
00328         actualRhsPtr = ei_aligned_stack_new(RhsScalar, actualRhs.size());
00329       }
00330       Map<typename _ActualRhsType::PlainObject>(actualRhsPtr, actualRhs.size()) = actualRhs;
00331     }
00332     
00333     internal::product_triangular_matrix_vector
00334       <Index,Mode,
00335        LhsScalar, LhsBlasTraits::NeedToConjugate,
00336        RhsScalar, RhsBlasTraits::NeedToConjugate,
00337        RowMajor>
00338       ::run(actualLhs.rows(),actualLhs.cols(),
00339             actualLhs.data(),actualLhs.outerStride(),
00340             actualRhsPtr,1,
00341             dest.data(),dest.innerStride(),
00342             actualAlpha);
00343 
00344     if((!DirectlyUseRhs) && freeRhsPtr) ei_aligned_stack_delete(RhsScalar, actualRhsPtr, prod.rhs().size());
00345   }
00346 };
00347 
00348 } // end namespace internal
00349 
00350 #endif // EIGEN_TRIANGULARMATRIXVECTOR_H