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

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

Doc

File size: 3.5 KB
Line 
1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr>
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_BAND_TRIANGULARSOLVER_H
11#define EIGEN_BAND_TRIANGULARSOLVER_H
12
13namespace internal {
14
15 /* \internal
16 * Solve Ax=b with A a band triangular matrix
17 * TODO: extend it to matrices for x abd b */
18template<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar, int StorageOrder>
19struct band_solve_triangular_selector;
20
21
22template<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar>
23struct band_solve_triangular_selector<Index,Mode,LhsScalar,ConjLhs,RhsScalar,RowMajor>
24{
25 typedef Map<const Matrix<LhsScalar,Dynamic,Dynamic,RowMajor>, 0, OuterStride<> > LhsMap;
26 typedef Map<Matrix<RhsScalar,Dynamic,1> > RhsMap;
27 enum { IsLower = (Mode&Lower) ? 1 : 0 };
28 static void run(Index size, Index k, const LhsScalar* _lhs, Index lhsStride, RhsScalar* _other)
29 {
30 const LhsMap lhs(_lhs,size,k+1,OuterStride<>(lhsStride));
31 RhsMap other(_other,size,1);
32 typename internal::conditional<
33 ConjLhs,
34 const CwiseUnaryOp<typename internal::scalar_conjugate_op<LhsScalar>,LhsMap>,
35 const LhsMap&>
36 ::type cjLhs(lhs);
37
38 for(int col=0 ; col<other.cols() ; ++col)
39 {
40 for(int ii=0; ii<size; ++ii)
41 {
42 int i = IsLower ? ii : size-ii-1;
43 int actual_k = (std::min)(k,ii);
44 int actual_start = IsLower ? k-actual_k : 1;
45
46 if(actual_k>0)
47 other.coeffRef(i,col) -= cjLhs.row(i).segment(actual_start,actual_k).transpose()
48 .cwiseProduct(other.col(col).segment(IsLower ? i-actual_k : i+1,actual_k)).sum();
49
50 if((Mode&UnitDiag)==0)
51 other.coeffRef(i,col) /= cjLhs(i,IsLower ? k : 0);
52 }
53 }
54 }
55
56};
57
58template<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar>
59struct band_solve_triangular_selector<Index,Mode,LhsScalar,ConjLhs,RhsScalar,ColMajor>
60{
61 typedef Map<const Matrix<LhsScalar,Dynamic,Dynamic,ColMajor>, 0, OuterStride<> > LhsMap;
62 typedef Map<Matrix<RhsScalar,Dynamic,1> > RhsMap;
63 enum { IsLower = (Mode&Lower) ? 1 : 0 };
64 static void run(Index size, Index k, const LhsScalar* _lhs, Index lhsStride, RhsScalar* _other)
65 {
66 const LhsMap lhs(_lhs,k+1,size,OuterStride<>(lhsStride));
67 RhsMap other(_other,size,1);
68 typename internal::conditional<
69 ConjLhs,
70 const CwiseUnaryOp<typename internal::scalar_conjugate_op<LhsScalar>,LhsMap>,
71 const LhsMap&>
72 ::type cjLhs(lhs);
73
74 for(int col=0 ; col<other.cols() ; ++col)
75 {
76 for(int ii=0; ii<size; ++ii)
77 {
78 int i = IsLower ? ii : size-ii-1;
79 int actual_k = (std::min)(k,size-ii-1);
80 int actual_start = IsLower ? 1 : k-actual_k;
81
82 if((Mode&UnitDiag)==0)
83 other.coeffRef(i,col) /= cjLhs(IsLower ? 0 : k, i);
84
85 if(actual_k>0)
86 other.col(col).segment(IsLower ? i+1 : i-actual_k, actual_k)
87 -= other.coeff(i,col) * cjLhs.col(i).segment(actual_start,actual_k);
88
89 }
90 }
91 }
92};
93
94
95} // end namespace internal
96
97#endif // EIGEN_BAND_TRIANGULARSOLVER_H
Note: See TracBrowser for help on using the repository browser.