source: pacpussensors/trunk/Vislab/lib3dv/eigen/blas/PackedTriangularSolverVector.h@ 138

Last change on this file since 138 was 136, checked in by ldecherf, 8 years ago

Doc

File size: 3.1 KB
Line 
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
13namespace internal {
14
15template<typename LhsScalar, typename RhsScalar, typename Index, int Side, int Mode, bool Conjugate, int StorageOrder>
16struct packed_triangular_solve_vector;
17
18// forward and backward substitution, row-major, rhs is a vector
19template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate>
20struct 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
47template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate>
48struct 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
74template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate, int StorageOrder>
75struct 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
Note: See TracBrowser for help on using the repository browser.