Home | History | Annotate | Download | only in blas 
      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra.
      3 //
      4 // Copyright (C) 2011 Gael Guennebaud <gael.guennebaud (at) inria.fr>
      5 //
      6 // This Source Code Form is subject to the terms of the Mozilla
      7 // Public License v. 2.0. If a copy of the MPL was not distributed
      8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
      9 
     10 #ifndef EIGEN_BAND_TRIANGULARSOLVER_H
     11 #define EIGEN_BAND_TRIANGULARSOLVER_H
     12 
     13 namespace internal {
     14 
     15  /* \internal
     16   * Solve Ax=b with A a band triangular matrix
     17   * TODO: extend it to matrices for x abd b */
     18 template<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar, int StorageOrder>
     19 struct band_solve_triangular_selector;
     20 
     21 
     22 template<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar>
     23 struct band_solve_triangular_selector<Index,Mode,LhsScalar,ConjLhs,RhsScalar,RowMajor>
     24 {
     25   typedef Map<const Matrix<LhsScalar,Dynamic,Dynamic,RowMajor>, 0, OuterStride<> > LhsMap;
     26   typedef Map<Matrix<RhsScalar,Dynamic,1> > RhsMap;
     27   enum { IsLower = (Mode&Lower) ? 1 : 0 };
     28   static void run(Index size, Index k, const LhsScalar* _lhs, Index lhsStride, RhsScalar* _other)
     29   {
     30     const LhsMap lhs(_lhs,size,k+1,OuterStride<>(lhsStride));
     31     RhsMap other(_other,size,1);
     32     typename internal::conditional<
     33                           ConjLhs,
     34                           const CwiseUnaryOp<typename internal::scalar_conjugate_op<LhsScalar>,LhsMap>,
     35                           const LhsMap&>
     36                         ::type cjLhs(lhs);
     37 
     38     for(int col=0 ; col<other.cols() ; ++col)
     39     {
     40       for(int ii=0; ii<size; ++ii)
     41       {
     42         int i = IsLower ? ii : size-ii-1;
     43         int actual_k = (std::min)(k,ii);
     44         int actual_start = IsLower ? k-actual_k : 1;
     45 
     46         if(actual_k>0)
     47           other.coeffRef(i,col) -= cjLhs.row(i).segment(actual_start,actual_k).transpose()
     48                                   .cwiseProduct(other.col(col).segment(IsLower ? i-actual_k : i+1,actual_k)).sum();
     49 
     50         if((Mode&UnitDiag)==0)
     51           other.coeffRef(i,col) /= cjLhs(i,IsLower ? k : 0);
     52       }
     53     }
     54   }
     55 
     56 };
     57 
     58 template<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar>
     59 struct band_solve_triangular_selector<Index,Mode,LhsScalar,ConjLhs,RhsScalar,ColMajor>
     60 {
     61   typedef Map<const Matrix<LhsScalar,Dynamic,Dynamic,ColMajor>, 0, OuterStride<> > LhsMap;
     62   typedef Map<Matrix<RhsScalar,Dynamic,1> > RhsMap;
     63   enum { IsLower = (Mode&Lower) ? 1 : 0 };
     64   static void run(Index size, Index k, const LhsScalar* _lhs, Index lhsStride, RhsScalar* _other)
     65   {
     66     const LhsMap lhs(_lhs,k+1,size,OuterStride<>(lhsStride));
     67     RhsMap other(_other,size,1);
     68     typename internal::conditional<
     69                           ConjLhs,
     70                           const CwiseUnaryOp<typename internal::scalar_conjugate_op<LhsScalar>,LhsMap>,
     71                           const LhsMap&>
     72                         ::type cjLhs(lhs);
     73 
     74     for(int col=0 ; col<other.cols() ; ++col)
     75     {
     76       for(int ii=0; ii<size; ++ii)
     77       {
     78         int i = IsLower ? ii : size-ii-1;
     79         int actual_k = (std::min)(k,size-ii-1);
     80         int actual_start = IsLower ? 1 : k-actual_k;
     81 
     82         if((Mode&UnitDiag)==0)
     83           other.coeffRef(i,col) /= cjLhs(IsLower ? 0 : k, i);
     84 
     85         if(actual_k>0)
     86           other.col(col).segment(IsLower ? i+1 : i-actual_k, actual_k)
     87               -= other.coeff(i,col) * cjLhs.col(i).segment(actual_start,actual_k);
     88 
     89       }
     90     }
     91   }
     92 };
     93 
     94 
     95 } // end namespace internal
     96 
     97 #endif // EIGEN_BAND_TRIANGULARSOLVER_H
     98