TriangularSolverVector.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) 2008-2010 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_TRIANGULAR_SOLVER_VECTOR_H
00026 #define EIGEN_TRIANGULAR_SOLVER_VECTOR_H
00027 
00028 namespace Eigen { 
00029 
00030 namespace internal {
00031 
00032 template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate, int StorageOrder>
00033 struct triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheRight, Mode, Conjugate, StorageOrder>
00034 {
00035   static void run(Index size, const LhsScalar* _lhs, Index lhsStride, RhsScalar* rhs)
00036   {
00037     triangular_solve_vector<LhsScalar,RhsScalar,Index,OnTheLeft,
00038         ((Mode&Upper)==Upper ? Lower : Upper) | (Mode&UnitDiag),
00039         Conjugate,StorageOrder==RowMajor?ColMajor:RowMajor
00040       >::run(size, _lhs, lhsStride, rhs);
00041   }
00042 };
00043     
00044 // forward and backward substitution, row-major, rhs is a vector
00045 template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate>
00046 struct triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, Mode, Conjugate, RowMajor>
00047 {
00048   enum {
00049     IsLower = ((Mode&Lower)==Lower)
00050   };
00051   static void run(Index size, const LhsScalar* _lhs, Index lhsStride, RhsScalar* rhs)
00052   {
00053     typedef Map<const Matrix<LhsScalar,Dynamic,Dynamic,RowMajor>, 0, OuterStride<> > LhsMap;
00054     const LhsMap lhs(_lhs,size,size,OuterStride<>(lhsStride));
00055     typename internal::conditional<
00056                           Conjugate,
00057                           const CwiseUnaryOp<typename internal::scalar_conjugate_op<LhsScalar>,LhsMap>,
00058                           const LhsMap&>
00059                         ::type cjLhs(lhs);
00060     static const Index PanelWidth = EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH;
00061     for(Index pi=IsLower ? 0 : size;
00062         IsLower ? pi<size : pi>0;
00063         IsLower ? pi+=PanelWidth : pi-=PanelWidth)
00064     {
00065       Index actualPanelWidth = (std::min)(IsLower ? size - pi : pi, PanelWidth);
00066 
00067       Index r = IsLower ? pi : size - pi; // remaining size
00068       if (r > 0)
00069       {
00070         // let's directly call the low level product function because:
00071         // 1 - it is faster to compile
00072         // 2 - it is slighlty faster at runtime
00073         Index startRow = IsLower ? pi : pi-actualPanelWidth;
00074         Index startCol = IsLower ? 0 : pi;
00075 
00076         general_matrix_vector_product<Index,LhsScalar,RowMajor,Conjugate,RhsScalar,false>::run(
00077           actualPanelWidth, r,
00078           &lhs.coeffRef(startRow,startCol), lhsStride,
00079           rhs + startCol, 1,
00080           rhs + startRow, 1,
00081           RhsScalar(-1));
00082       }
00083 
00084       for(Index k=0; k<actualPanelWidth; ++k)
00085       {
00086         Index i = IsLower ? pi+k : pi-k-1;
00087         Index s = IsLower ? pi   : i+1;
00088         if (k>0)
00089           rhs[i] -= (cjLhs.row(i).segment(s,k).transpose().cwiseProduct(Map<const Matrix<RhsScalar,Dynamic,1> >(rhs+s,k))).sum();
00090         
00091         if(!(Mode & UnitDiag))
00092           rhs[i] /= cjLhs(i,i);
00093       }
00094     }
00095   }
00096 };
00097 
00098 // forward and backward substitution, column-major, rhs is a vector
00099 template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate>
00100 struct triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, Mode, Conjugate, ColMajor>
00101 {
00102   enum {
00103     IsLower = ((Mode&Lower)==Lower)
00104   };
00105   static void run(Index size, const LhsScalar* _lhs, Index lhsStride, RhsScalar* rhs)
00106   {
00107     typedef Map<const Matrix<LhsScalar,Dynamic,Dynamic,ColMajor>, 0, OuterStride<> > LhsMap;
00108     const LhsMap lhs(_lhs,size,size,OuterStride<>(lhsStride));
00109     typename internal::conditional<Conjugate,
00110                                    const CwiseUnaryOp<typename internal::scalar_conjugate_op<LhsScalar>,LhsMap>,
00111                                    const LhsMap&
00112                                   >::type cjLhs(lhs);
00113     static const Index PanelWidth = EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH;
00114 
00115     for(Index pi=IsLower ? 0 : size;
00116         IsLower ? pi<size : pi>0;
00117         IsLower ? pi+=PanelWidth : pi-=PanelWidth)
00118     {
00119       Index actualPanelWidth = (std::min)(IsLower ? size - pi : pi, PanelWidth);
00120       Index startBlock = IsLower ? pi : pi-actualPanelWidth;
00121       Index endBlock = IsLower ? pi + actualPanelWidth : 0;
00122 
00123       for(Index k=0; k<actualPanelWidth; ++k)
00124       {
00125         Index i = IsLower ? pi+k : pi-k-1;
00126         if(!(Mode & UnitDiag))
00127           rhs[i] /= cjLhs.coeff(i,i);
00128 
00129         Index r = actualPanelWidth - k - 1; // remaining size
00130         Index s = IsLower ? i+1 : i-r;
00131         if (r>0)
00132           Map<Matrix<RhsScalar,Dynamic,1> >(rhs+s,r) -= rhs[i] * cjLhs.col(i).segment(s,r);
00133       }
00134       Index r = IsLower ? size - endBlock : startBlock; // remaining size
00135       if (r > 0)
00136       {
00137         // let's directly call the low level product function because:
00138         // 1 - it is faster to compile
00139         // 2 - it is slighlty faster at runtime
00140         general_matrix_vector_product<Index,LhsScalar,ColMajor,Conjugate,RhsScalar,false>::run(
00141             r, actualPanelWidth,
00142             &lhs.coeffRef(endBlock,startBlock), lhsStride,
00143             rhs+startBlock, 1,
00144             rhs+endBlock, 1, RhsScalar(-1));
00145       }
00146     }
00147   }
00148 };
00149 
00150 } // end namespace internal
00151 
00152 } // end namespace Eigen
00153 
00154 #endif // EIGEN_TRIANGULAR_SOLVER_VECTOR_H