GeneralMatrixMatrixTriangular.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-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 //
00006 // This Source Code Form is subject to the terms of the Mozilla
00007 // Public License v. 2.0. If a copy of the MPL was not distributed
00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00009 
00010 #ifndef EIGEN_GENERAL_MATRIX_MATRIX_TRIANGULAR_H
00011 #define EIGEN_GENERAL_MATRIX_MATRIX_TRIANGULAR_H
00012 
00013 namespace Eigen { 
00014 
00015 namespace internal {
00016 
00017 /**********************************************************************
00018 * This file implements a general A * B product while
00019 * evaluating only one triangular part of the product.
00020 * This is more general version of self adjoint product (C += A A^T)
00021 * as the level 3 SYRK Blas routine.
00022 **********************************************************************/
00023 
00024 // forward declarations (defined at the end of this file)
00025 template<typename LhsScalar, typename RhsScalar, typename Index, int mr, int nr, bool ConjLhs, bool ConjRhs, int UpLo>
00026 struct tribb_kernel;
00027   
00028 /* Optimized matrix-matrix product evaluating only one triangular half */
00029 template <typename Index,
00030           typename LhsScalar, int LhsStorageOrder, bool ConjugateLhs,
00031           typename RhsScalar, int RhsStorageOrder, bool ConjugateRhs,
00032                               int ResStorageOrder, int  UpLo, int Version = Specialized>
00033 struct general_matrix_matrix_triangular_product;
00034 
00035 // as usual if the result is row major => we transpose the product
00036 template <typename Index, typename LhsScalar, int LhsStorageOrder, bool ConjugateLhs,
00037                           typename RhsScalar, int RhsStorageOrder, bool ConjugateRhs, int  UpLo, int Version>
00038 struct general_matrix_matrix_triangular_product<Index,LhsScalar,LhsStorageOrder,ConjugateLhs,RhsScalar,RhsStorageOrder,ConjugateRhs,RowMajor,UpLo,Version>
00039 {
00040   typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar;
00041   static EIGEN_STRONG_INLINE void run(Index size, Index depth,const LhsScalar* lhs, Index lhsStride,
00042                                       const RhsScalar* rhs, Index rhsStride, ResScalar* res, Index resStride, ResScalar alpha)
00043   {
00044     general_matrix_matrix_triangular_product<Index,
00045         RhsScalar, RhsStorageOrder==RowMajor ? ColMajor : RowMajor, ConjugateRhs,
00046         LhsScalar, LhsStorageOrder==RowMajor ? ColMajor : RowMajor, ConjugateLhs,
00047         ColMajor, UpLo==Lower?Upper:Lower>
00048       ::run(size,depth,rhs,rhsStride,lhs,lhsStride,res,resStride,alpha);
00049   }
00050 };
00051 
00052 template <typename Index, typename LhsScalar, int LhsStorageOrder, bool ConjugateLhs,
00053                           typename RhsScalar, int RhsStorageOrder, bool ConjugateRhs, int  UpLo, int Version>
00054 struct general_matrix_matrix_triangular_product<Index,LhsScalar,LhsStorageOrder,ConjugateLhs,RhsScalar,RhsStorageOrder,ConjugateRhs,ColMajor,UpLo,Version>
00055 {
00056   typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar;
00057   static EIGEN_STRONG_INLINE void run(Index size, Index depth,const LhsScalar* _lhs, Index lhsStride,
00058                                       const RhsScalar* _rhs, Index rhsStride, ResScalar* res, Index resStride, ResScalar alpha)
00059   {
00060     const_blas_data_mapper<LhsScalar, Index, LhsStorageOrder> lhs(_lhs,lhsStride);
00061     const_blas_data_mapper<RhsScalar, Index, RhsStorageOrder> rhs(_rhs,rhsStride);
00062 
00063     typedef gebp_traits<LhsScalar,RhsScalar> Traits;
00064 
00065     Index kc = depth; // cache block size along the K direction
00066     Index mc = size;  // cache block size along the M direction
00067     Index nc = size;  // cache block size along the N direction
00068     computeProductBlockingSizes<LhsScalar,RhsScalar>(kc, mc, nc);
00069     // !!! mc must be a multiple of nr:
00070     if(mc > Traits::nr)
00071       mc = (mc/Traits::nr)*Traits::nr;
00072 
00073     std::size_t sizeW = kc*Traits::WorkSpaceFactor;
00074     std::size_t sizeB = sizeW + kc*size;
00075     ei_declare_aligned_stack_constructed_variable(LhsScalar, blockA, kc*mc, 0);
00076     ei_declare_aligned_stack_constructed_variable(RhsScalar, allocatedBlockB, sizeB, 0);
00077     RhsScalar* blockB = allocatedBlockB + sizeW;
00078     
00079     gemm_pack_lhs<LhsScalar, Index, Traits::mr, Traits::LhsProgress, LhsStorageOrder> pack_lhs;
00080     gemm_pack_rhs<RhsScalar, Index, Traits::nr, RhsStorageOrder> pack_rhs;
00081     gebp_kernel <LhsScalar, RhsScalar, Index, Traits::mr, Traits::nr, ConjugateLhs, ConjugateRhs> gebp;
00082     tribb_kernel<LhsScalar, RhsScalar, Index, Traits::mr, Traits::nr, ConjugateLhs, ConjugateRhs, UpLo> sybb;
00083 
00084     for(Index k2=0; k2<depth; k2+=kc)
00085     {
00086       const Index actual_kc = (std::min)(k2+kc,depth)-k2;
00087 
00088       // note that the actual rhs is the transpose/adjoint of mat
00089       pack_rhs(blockB, &rhs(k2,0), rhsStride, actual_kc, size);
00090 
00091       for(Index i2=0; i2<size; i2+=mc)
00092       {
00093         const Index actual_mc = (std::min)(i2+mc,size)-i2;
00094 
00095         pack_lhs(blockA, &lhs(i2, k2), lhsStride, actual_kc, actual_mc);
00096 
00097         // the selected actual_mc * size panel of res is split into three different part:
00098         //  1 - before the diagonal => processed with gebp or skipped
00099         //  2 - the actual_mc x actual_mc symmetric block => processed with a special kernel
00100         //  3 - after the diagonal => processed with gebp or skipped
00101         if (UpLo==Lower)
00102           gebp(res+i2, resStride, blockA, blockB, actual_mc, actual_kc, (std::min)(size,i2), alpha,
00103                -1, -1, 0, 0, allocatedBlockB);
00104 
00105         sybb(res+resStride*i2 + i2, resStride, blockA, blockB + actual_kc*i2, actual_mc, actual_kc, alpha, allocatedBlockB);
00106 
00107         if (UpLo==Upper)
00108         {
00109           Index j2 = i2+actual_mc;
00110           gebp(res+resStride*j2+i2, resStride, blockA, blockB+actual_kc*j2, actual_mc, actual_kc, (std::max)(Index(0), size-j2), alpha,
00111                -1, -1, 0, 0, allocatedBlockB);
00112         }
00113       }
00114     }
00115   }
00116 };
00117 
00118 // Optimized packed Block * packed Block product kernel evaluating only one given triangular part
00119 // This kernel is built on top of the gebp kernel:
00120 // - the current destination block is processed per panel of actual_mc x BlockSize
00121 //   where BlockSize is set to the minimal value allowing gebp to be as fast as possible
00122 // - then, as usual, each panel is split into three parts along the diagonal,
00123 //   the sub blocks above and below the diagonal are processed as usual,
00124 //   while the triangular block overlapping the diagonal is evaluated into a
00125 //   small temporary buffer which is then accumulated into the result using a
00126 //   triangular traversal.
00127 template<typename LhsScalar, typename RhsScalar, typename Index, int mr, int nr, bool ConjLhs, bool ConjRhs, int UpLo>
00128 struct tribb_kernel
00129 {
00130   typedef gebp_traits<LhsScalar,RhsScalar,ConjLhs,ConjRhs> Traits;
00131   typedef typename Traits::ResScalar ResScalar;
00132   
00133   enum {
00134     BlockSize  = EIGEN_PLAIN_ENUM_MAX(mr,nr)
00135   };
00136   void operator()(ResScalar* res, Index resStride, const LhsScalar* blockA, const RhsScalar* blockB, Index size, Index depth, ResScalar alpha, RhsScalar* workspace)
00137   {
00138     gebp_kernel<LhsScalar, RhsScalar, Index, mr, nr, ConjLhs, ConjRhs> gebp_kernel;
00139     Matrix<ResScalar,BlockSize,BlockSize,ColMajor> buffer;
00140 
00141     // let's process the block per panel of actual_mc x BlockSize,
00142     // again, each is split into three parts, etc.
00143     for (Index j=0; j<size; j+=BlockSize)
00144     {
00145       Index actualBlockSize = std::min<Index>(BlockSize,size - j);
00146       const RhsScalar* actual_b = blockB+j*depth;
00147 
00148       if(UpLo==Upper)
00149         gebp_kernel(res+j*resStride, resStride, blockA, actual_b, j, depth, actualBlockSize, alpha,
00150                     -1, -1, 0, 0, workspace);
00151 
00152       // selfadjoint micro block
00153       {
00154         Index i = j;
00155         buffer.setZero();
00156         // 1 - apply the kernel on the temporary buffer
00157         gebp_kernel(buffer.data(), BlockSize, blockA+depth*i, actual_b, actualBlockSize, depth, actualBlockSize, alpha,
00158                     -1, -1, 0, 0, workspace);
00159         // 2 - triangular accumulation
00160         for(Index j1=0; j1<actualBlockSize; ++j1)
00161         {
00162           ResScalar* r = res + (j+j1)*resStride + i;
00163           for(Index i1=UpLo==Lower ? j1 : 0;
00164               UpLo==Lower ? i1<actualBlockSize : i1<=j1; ++i1)
00165             r[i1] += buffer(i1,j1);
00166         }
00167       }
00168 
00169       if(UpLo==Lower)
00170       {
00171         Index i = j+actualBlockSize;
00172         gebp_kernel(res+j*resStride+i, resStride, blockA+depth*i, actual_b, size-i, depth, actualBlockSize, alpha,
00173                     -1, -1, 0, 0, workspace);
00174       }
00175     }
00176   }
00177 };
00178 
00179 } // end namespace internal
00180 
00181 // high level API
00182 
00183 template<typename MatrixType, unsigned int UpLo>
00184 template<typename ProductDerived, typename _Lhs, typename _Rhs>
00185 TriangularView<MatrixType,UpLo>& TriangularView<MatrixType,UpLo>::assignProduct(const ProductBase<ProductDerived, _Lhs,_Rhs>& prod, const Scalar& alpha)
00186 {
00187   typedef typename internal::remove_all<typename ProductDerived::LhsNested>::type Lhs;
00188   typedef internal::blas_traits<Lhs> LhsBlasTraits;
00189   typedef typename LhsBlasTraits::DirectLinearAccessType ActualLhs;
00190   typedef typename internal::remove_all<ActualLhs>::type _ActualLhs;
00191   typename internal::add_const_on_value_type<ActualLhs>::type actualLhs = LhsBlasTraits::extract(prod.lhs());
00192   
00193   typedef typename internal::remove_all<typename ProductDerived::RhsNested>::type Rhs;
00194   typedef internal::blas_traits<Rhs> RhsBlasTraits;
00195   typedef typename RhsBlasTraits::DirectLinearAccessType ActualRhs;
00196   typedef typename internal::remove_all<ActualRhs>::type _ActualRhs;
00197   typename internal::add_const_on_value_type<ActualRhs>::type actualRhs = RhsBlasTraits::extract(prod.rhs());
00198 
00199   typename ProductDerived::Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(prod.lhs().derived()) * RhsBlasTraits::extractScalarFactor(prod.rhs().derived());
00200 
00201   internal::general_matrix_matrix_triangular_product<Index,
00202     typename Lhs::Scalar, _ActualLhs::Flags&RowMajorBit ? RowMajor : ColMajor, LhsBlasTraits::NeedToConjugate,
00203     typename Rhs::Scalar, _ActualRhs::Flags&RowMajorBit ? RowMajor : ColMajor, RhsBlasTraits::NeedToConjugate,
00204     MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor, UpLo>
00205     ::run(m_matrix.cols(), actualLhs.cols(),
00206           &actualLhs.coeffRef(0,0), actualLhs.outerStride(), &actualRhs.coeffRef(0,0), actualRhs.outerStride(),
00207           const_cast<Scalar*>(m_matrix.data()), m_matrix.outerStride(), actualAlpha);
00208   
00209   return *this;
00210 }
00211 
00212 } // end namespace Eigen
00213 
00214 #endif // EIGEN_GENERAL_MATRIX_MATRIX_TRIANGULAR_H


win_eigen
Author(s): Daniel Stonier
autogenerated on Mon Oct 6 2014 12:24:34