You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
88 lines
3.1 KiB
88 lines
3.1 KiB
// This file is part of Eigen, a lightweight C++ template library |
|
// for linear algebra. |
|
// |
|
// Copyright (C) 2012 Chen-Pang He <jdh8@ms63.hinet.net> |
|
// |
|
// This Source Code Form is subject to the terms of the Mozilla |
|
// Public License v. 2.0. If a copy of the MPL was not distributed |
|
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. |
|
|
|
#ifndef EIGEN_PACKED_TRIANGULAR_SOLVER_VECTOR_H |
|
#define EIGEN_PACKED_TRIANGULAR_SOLVER_VECTOR_H |
|
|
|
namespace internal { |
|
|
|
template<typename LhsScalar, typename RhsScalar, typename Index, int Side, int Mode, bool Conjugate, int StorageOrder> |
|
struct packed_triangular_solve_vector; |
|
|
|
// forward and backward substitution, row-major, rhs is a vector |
|
template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate> |
|
struct packed_triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, Mode, Conjugate, RowMajor> |
|
{ |
|
enum { |
|
IsLower = (Mode&Lower)==Lower |
|
}; |
|
static void run(Index size, const LhsScalar* lhs, RhsScalar* rhs) |
|
{ |
|
internal::conj_if<Conjugate> cj; |
|
typedef Map<const Matrix<LhsScalar,Dynamic,1> > LhsMap; |
|
typedef typename conj_expr_if<Conjugate,LhsMap>::type ConjLhsType; |
|
|
|
lhs += IsLower ? 0 : (size*(size+1)>>1)-1; |
|
for(Index pi=0; pi<size; ++pi) |
|
{ |
|
Index i = IsLower ? pi : size-pi-1; |
|
Index s = IsLower ? 0 : 1; |
|
if (pi>0) |
|
rhs[i] -= (ConjLhsType(LhsMap(lhs+s,pi)) |
|
.cwiseProduct(Map<const Matrix<RhsScalar,Dynamic,1> >(rhs+(IsLower ? 0 : i+1),pi))).sum(); |
|
if (!(Mode & UnitDiag)) |
|
rhs[i] /= cj(lhs[IsLower ? i : 0]); |
|
IsLower ? lhs += pi+1 : lhs -= pi+2; |
|
} |
|
} |
|
}; |
|
|
|
// forward and backward substitution, column-major, rhs is a vector |
|
template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate> |
|
struct packed_triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, Mode, Conjugate, ColMajor> |
|
{ |
|
enum { |
|
IsLower = (Mode&Lower)==Lower |
|
}; |
|
static void run(Index size, const LhsScalar* lhs, RhsScalar* rhs) |
|
{ |
|
internal::conj_if<Conjugate> cj; |
|
typedef Map<const Matrix<LhsScalar,Dynamic,1> > LhsMap; |
|
typedef typename conj_expr_if<Conjugate,LhsMap>::type ConjLhsType; |
|
|
|
lhs += IsLower ? 0 : size*(size-1)>>1; |
|
for(Index pi=0; pi<size; ++pi) |
|
{ |
|
Index i = IsLower ? pi : size-pi-1; |
|
Index r = size - pi - 1; |
|
if (!(Mode & UnitDiag)) |
|
rhs[i] /= cj(lhs[IsLower ? 0 : i]); |
|
if (r>0) |
|
Map<Matrix<RhsScalar,Dynamic,1> >(rhs+(IsLower? i+1 : 0),r) -= |
|
rhs[i] * ConjLhsType(LhsMap(lhs+(IsLower? 1 : 0),r)); |
|
IsLower ? lhs += size-pi : lhs -= r; |
|
} |
|
} |
|
}; |
|
|
|
template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate, int StorageOrder> |
|
struct packed_triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheRight, Mode, Conjugate, StorageOrder> |
|
{ |
|
static void run(Index size, const LhsScalar* lhs, RhsScalar* rhs) |
|
{ |
|
packed_triangular_solve_vector<LhsScalar,RhsScalar,Index,OnTheLeft, |
|
((Mode&Upper)==Upper ? Lower : Upper) | (Mode&UnitDiag), |
|
Conjugate,StorageOrder==RowMajor?ColMajor:RowMajor |
|
>::run(size, lhs, rhs); |
|
} |
|
}; |
|
|
|
} // end namespace internal |
|
|
|
#endif // EIGEN_PACKED_TRIANGULAR_SOLVER_VECTOR_H
|
|
|