Main MRPT website > C++ reference
MRPT logo

TriangularMatrixVector.h

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 } // end namespace internal
00156 
00157 template<int Mode, typename Lhs, typename Rhs>
00158 struct TriangularProduct<Mode,true,Lhs,false,Rhs,true>
00159   : public ProductBase<TriangularProduct<Mode,true,Lhs,false,Rhs,true>, Lhs, Rhs >
00160 {
00161   EIGEN_PRODUCT_PUBLIC_INTERFACE(TriangularProduct)
00162 
00163   TriangularProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {}
00164 
00165   template<typename Dest> void scaleAndAddTo(Dest& dst, Scalar alpha) const
00166   {
00167     eigen_assert(dst.rows()==m_lhs.rows() && dst.cols()==m_rhs.cols());
00168 
00169     const ActualLhsType lhs = LhsBlasTraits::extract(m_lhs);
00170     const ActualRhsType rhs = RhsBlasTraits::extract(m_rhs);
00171 
00172     Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(m_lhs)
00173                                * RhsBlasTraits::extractScalarFactor(m_rhs);
00174 
00175     internal::product_triangular_matrix_vector
00176       <Index,Mode,
00177        typename _ActualLhsType::Scalar, LhsBlasTraits::NeedToConjugate,
00178        typename _ActualRhsType::Scalar, RhsBlasTraits::NeedToConjugate,
00179        (int(internal::traits<Lhs>::Flags)&RowMajorBit) ? RowMajor : ColMajor>
00180       ::run(lhs.rows(),lhs.cols(),lhs.data(),lhs.outerStride(),rhs.data(),rhs.innerStride(),
00181                                   dst.data(),dst.innerStride(),actualAlpha);
00182   }
00183 };
00184 
00185 template<int Mode, typename Lhs, typename Rhs>
00186 struct TriangularProduct<Mode,false,Lhs,true,Rhs,false>
00187   : public ProductBase<TriangularProduct<Mode,false,Lhs,true,Rhs,false>, Lhs, Rhs >
00188 {
00189   EIGEN_PRODUCT_PUBLIC_INTERFACE(TriangularProduct)
00190 
00191   TriangularProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {}
00192 
00193   template<typename Dest> void scaleAndAddTo(Dest& dst, Scalar alpha) const
00194   {
00195 
00196     eigen_assert(dst.rows()==m_lhs.rows() && dst.cols()==m_rhs.cols());
00197 
00198     const ActualLhsType lhs = LhsBlasTraits::extract(m_lhs);
00199     const ActualRhsType rhs = RhsBlasTraits::extract(m_rhs);
00200 
00201     Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(m_lhs)
00202                                * RhsBlasTraits::extractScalarFactor(m_rhs);
00203 
00204     internal::product_triangular_matrix_vector
00205       <Index,(Mode & UnitDiag) | (Mode & Lower) ? Upper : Lower,
00206        typename _ActualRhsType::Scalar, RhsBlasTraits::NeedToConjugate,
00207        typename _ActualLhsType::Scalar, LhsBlasTraits::NeedToConjugate,
00208        (int(internal::traits<Rhs>::Flags)&RowMajorBit) ? ColMajor : RowMajor>
00209       ::run(rhs.rows(),rhs.cols(),rhs.data(),rhs.outerStride(),lhs.data(),lhs.innerStride(),
00210                                   dst.data(),dst.innerStride(),actualAlpha);
00211   }
00212 };
00213 
00214 #endif // EIGEN_TRIANGULARMATRIXVECTOR_H



Page generated by Doxygen 1.7.3 for MRPT 0.9.4 SVN: at Sat Mar 26 06:40:17 UTC 2011