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


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:33:49