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

   1 // This file is part of Eigen, a lightweight C++ template library
   2 // for linear algebra.
   3 //
   4 // Copyright (C) 2012 Chen-Pang He <jdh8@ms63.hinet.net>
   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_PACKED_TRIANGULAR_SOLVER_VECTOR_H
  11 #define EIGEN_PACKED_TRIANGULAR_SOLVER_VECTOR_H
  12
  13 namespace internal {
  14
  15 template<typename LhsScalar, typename RhsScalar, typename Index, int Side, int Mode, bool Conjugate, int StorageOrder>
  16 struct packed_triangular_solve_vector;
  17
  18 // forward and backward substitution, row-major, rhs is a vector
  19 template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate>
  20 struct packed_triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, Mode, Conjugate, RowMajor>
  21 {
  22   enum {
  23     IsLower = (Mode&Lower)==Lower
  24   };
  25   static void run(Index size, const LhsScalar* lhs, RhsScalar* rhs)
  26   {
  27     internal::conj_if<Conjugate> cj;
  28     typedef Map<const Matrix<LhsScalar,Dynamic,1> > LhsMap;
  29     typedef typename conj_expr_if<Conjugate,LhsMap>::type ConjLhsType;
  30
  31     lhs += IsLower ? 0 : (size*(size+1)>>1)-1;
  32     for(Index pi=0; pi<size; ++pi)
  33     {
  34       Index i = IsLower ? pi : size-pi-1;
  35       Index s = IsLower ? 0 : 1;
  36       if (pi>0)
  37         rhs[i] -= (ConjLhsType(LhsMap(lhs+s,pi))
  38             .cwiseProduct(Map<const Matrix<RhsScalar,Dynamic,1> >(rhs+(IsLower ? 0 : i+1),pi))).sum();
  39       if (!(Mode & UnitDiag))
  40         rhs[i] /= cj(lhs[IsLower ? i : 0]);
  41       IsLower ? lhs += pi+1 : lhs -= pi+2;
  42     }
  43   }
  44 };
  45
  46 // forward and backward substitution, column-major, rhs is a vector
  47 template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate>
  48 struct packed_triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, Mode, Conjugate, ColMajor>
  49 {
  50   enum {
  51     IsLower = (Mode&Lower)==Lower
  52   };
  53   static void run(Index size, const LhsScalar* lhs, RhsScalar* rhs)
  54   {
  55     internal::conj_if<Conjugate> cj;
  56     typedef Map<const Matrix<LhsScalar,Dynamic,1> > LhsMap;
  57     typedef typename conj_expr_if<Conjugate,LhsMap>::type ConjLhsType;
  58
  59     lhs += IsLower ? 0 : size*(size-1)>>1;
  60     for(Index pi=0; pi<size; ++pi)
  61     {
  62       Index i = IsLower ? pi : size-pi-1;
  63       Index r = size - pi - 1;
  64       if (!(Mode & UnitDiag))
  65         rhs[i] /= cj(lhs[IsLower ? 0 : i]);
  66       if (r>0)
  67         Map<Matrix<RhsScalar,Dynamic,1> >(rhs+(IsLower? i+1 : 0),r) -=
  68             rhs[i] * ConjLhsType(LhsMap(lhs+(IsLower? 1 : 0),r));
  69       IsLower ? lhs += size-pi : lhs -= r;
  70     }
  71   }
  72 };
  73
  74 template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate, int StorageOrder>
  75 struct packed_triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheRight, Mode, Conjugate, StorageOrder>
  76 {
  77   static void run(Index size, const LhsScalar* lhs, RhsScalar* rhs)
  78   {
  79     packed_triangular_solve_vector<LhsScalar,RhsScalar,Index,OnTheLeft,
  80         ((Mode&Upper)==Upper ? Lower : Upper) | (Mode&UnitDiag),
  81         Conjugate,StorageOrder==RowMajor?ColMajor:RowMajor
  82       >::run(size, lhs, rhs);
  83   }
  84 };
  85
  86 } // end namespace internal
  87
  88 #endif // EIGEN_PACKED_TRIANGULAR_SOLVER_VECTOR_H