ground/gcs/src/libs/eigen/blas/BandTriangularSolver.h

   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@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