SelfadjointMatrixVector.h
Go to the documentation of this file.
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2009 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_SELFADJOINT_MATRIX_VECTOR_H
11 #define EIGEN_SELFADJOINT_MATRIX_VECTOR_H
12 
13 namespace Eigen {
14 
15 namespace internal {
16 
17 /* Optimized selfadjoint matrix * vector product:
18  * This algorithm processes 2 columns at onces that allows to both reduce
19  * the number of load/stores of the result by a factor 2 and to reduce
20  * the instruction dependency.
21  */
22 
23 template<typename Scalar, typename Index, int StorageOrder, int UpLo, bool ConjugateLhs, bool ConjugateRhs, int Version=Specialized>
25 
26 template<typename Scalar, typename Index, int StorageOrder, int UpLo, bool ConjugateLhs, bool ConjugateRhs, int Version>
28 
29 {
30 static EIGEN_DONT_INLINE void run(
31  Index size,
32  const Scalar* lhs, Index lhsStride,
33  const Scalar* rhs,
34  Scalar* res,
35  Scalar alpha);
36 };
37 
38 template<typename Scalar, typename Index, int StorageOrder, int UpLo, bool ConjugateLhs, bool ConjugateRhs, int Version>
40  Index size,
41  const Scalar* lhs, Index lhsStride,
42  const Scalar* rhs,
43  Scalar* res,
44  Scalar alpha)
45 {
46  typedef typename packet_traits<Scalar>::type Packet;
47  typedef typename NumTraits<Scalar>::Real RealScalar;
48  const Index PacketSize = sizeof(Packet)/sizeof(Scalar);
49 
50  enum {
51  IsRowMajor = StorageOrder==RowMajor ? 1 : 0,
52  IsLower = UpLo == Lower ? 1 : 0,
53  FirstTriangular = IsRowMajor == IsLower
54  };
55 
56  conj_helper<Scalar,Scalar,NumTraits<Scalar>::IsComplex && EIGEN_LOGICAL_XOR(ConjugateLhs, IsRowMajor), ConjugateRhs> cj0;
57  conj_helper<Scalar,Scalar,NumTraits<Scalar>::IsComplex && EIGEN_LOGICAL_XOR(ConjugateLhs, !IsRowMajor), ConjugateRhs> cj1;
59 
60  conj_helper<Packet,Packet,NumTraits<Scalar>::IsComplex && EIGEN_LOGICAL_XOR(ConjugateLhs, IsRowMajor), ConjugateRhs> pcj0;
61  conj_helper<Packet,Packet,NumTraits<Scalar>::IsComplex && EIGEN_LOGICAL_XOR(ConjugateLhs, !IsRowMajor), ConjugateRhs> pcj1;
62 
63  Scalar cjAlpha = ConjugateRhs ? numext::conj(alpha) : alpha;
64 
65 
66  Index bound = (std::max)(Index(0),size-8) & 0xfffffffe;
67  if (FirstTriangular)
68  bound = size - bound;
69 
70  for (Index j=FirstTriangular ? bound : 0;
71  j<(FirstTriangular ? size : bound);j+=2)
72  {
73  const Scalar* EIGEN_RESTRICT A0 = lhs + j*lhsStride;
74  const Scalar* EIGEN_RESTRICT A1 = lhs + (j+1)*lhsStride;
75 
76  Scalar t0 = cjAlpha * rhs[j];
77  Packet ptmp0 = pset1<Packet>(t0);
78  Scalar t1 = cjAlpha * rhs[j+1];
79  Packet ptmp1 = pset1<Packet>(t1);
80 
81  Scalar t2(0);
82  Packet ptmp2 = pset1<Packet>(t2);
83  Scalar t3(0);
84  Packet ptmp3 = pset1<Packet>(t3);
85 
86  Index starti = FirstTriangular ? 0 : j+2;
87  Index endi = FirstTriangular ? j : size;
88  Index alignedStart = (starti) + internal::first_default_aligned(&res[starti], endi-starti);
89  Index alignedEnd = alignedStart + ((endi-alignedStart)/(PacketSize))*(PacketSize);
90 
91  res[j] += cjd.pmul(numext::real(A0[j]), t0);
92  res[j+1] += cjd.pmul(numext::real(A1[j+1]), t1);
93  if(FirstTriangular)
94  {
95  res[j] += cj0.pmul(A1[j], t1);
96  t3 += cj1.pmul(A1[j], rhs[j]);
97  }
98  else
99  {
100  res[j+1] += cj0.pmul(A0[j+1],t0);
101  t2 += cj1.pmul(A0[j+1], rhs[j+1]);
102  }
103 
104  for (Index i=starti; i<alignedStart; ++i)
105  {
106  res[i] += cj0.pmul(A0[i], t0) + cj0.pmul(A1[i],t1);
107  t2 += cj1.pmul(A0[i], rhs[i]);
108  t3 += cj1.pmul(A1[i], rhs[i]);
109  }
110  // Yes this an optimization for gcc 4.3 and 4.4 (=> huge speed up)
111  // gcc 4.2 does this optimization automatically.
112  const Scalar* EIGEN_RESTRICT a0It = A0 + alignedStart;
113  const Scalar* EIGEN_RESTRICT a1It = A1 + alignedStart;
114  const Scalar* EIGEN_RESTRICT rhsIt = rhs + alignedStart;
115  Scalar* EIGEN_RESTRICT resIt = res + alignedStart;
116  for (Index i=alignedStart; i<alignedEnd; i+=PacketSize)
117  {
118  Packet A0i = ploadu<Packet>(a0It); a0It += PacketSize;
119  Packet A1i = ploadu<Packet>(a1It); a1It += PacketSize;
120  Packet Bi = ploadu<Packet>(rhsIt); rhsIt += PacketSize; // FIXME should be aligned in most cases
121  Packet Xi = pload <Packet>(resIt);
122 
123  Xi = pcj0.pmadd(A0i,ptmp0, pcj0.pmadd(A1i,ptmp1,Xi));
124  ptmp2 = pcj1.pmadd(A0i, Bi, ptmp2);
125  ptmp3 = pcj1.pmadd(A1i, Bi, ptmp3);
126  pstore(resIt,Xi); resIt += PacketSize;
127  }
128  for (Index i=alignedEnd; i<endi; i++)
129  {
130  res[i] += cj0.pmul(A0[i], t0) + cj0.pmul(A1[i],t1);
131  t2 += cj1.pmul(A0[i], rhs[i]);
132  t3 += cj1.pmul(A1[i], rhs[i]);
133  }
134 
135  res[j] += alpha * (t2 + predux(ptmp2));
136  res[j+1] += alpha * (t3 + predux(ptmp3));
137  }
138  for (Index j=FirstTriangular ? 0 : bound;j<(FirstTriangular ? bound : size);j++)
139  {
140  const Scalar* EIGEN_RESTRICT A0 = lhs + j*lhsStride;
141 
142  Scalar t1 = cjAlpha * rhs[j];
143  Scalar t2(0);
144  res[j] += cjd.pmul(numext::real(A0[j]), t1);
145  for (Index i=FirstTriangular ? 0 : j+1; i<(FirstTriangular ? j : size); i++)
146  {
147  res[i] += cj0.pmul(A0[i], t1);
148  t2 += cj1.pmul(A0[i], rhs[i]);
149  }
150  res[j] += alpha * t2;
151  }
152 }
153 
154 } // end namespace internal
155 
156 /***************************************************************************
157 * Wrapper to product_selfadjoint_vector
158 ***************************************************************************/
159 
160 namespace internal {
161 
162 template<typename Lhs, int LhsMode, typename Rhs>
163 struct selfadjoint_product_impl<Lhs,LhsMode,false,Rhs,0,true>
164 {
166 
170 
174 
175  enum { LhsUpLo = LhsMode&(Upper|Lower) };
176 
177  template<typename Dest>
178  static void run(Dest& dest, const Lhs &a_lhs, const Rhs &a_rhs, const Scalar& alpha)
179  {
180  typedef typename Dest::Scalar ResScalar;
181  typedef typename Rhs::Scalar RhsScalar;
183 
184  eigen_assert(dest.rows()==a_lhs.rows() && dest.cols()==a_rhs.cols());
185 
186  typename internal::add_const_on_value_type<ActualLhsType>::type lhs = LhsBlasTraits::extract(a_lhs);
187  typename internal::add_const_on_value_type<ActualRhsType>::type rhs = RhsBlasTraits::extract(a_rhs);
188 
189  Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(a_lhs)
190  * RhsBlasTraits::extractScalarFactor(a_rhs);
191 
192  enum {
193  EvalToDest = (Dest::InnerStrideAtCompileTime==1),
194  UseRhs = (ActualRhsTypeCleaned::InnerStrideAtCompileTime==1)
195  };
196 
199 
200  ei_declare_aligned_stack_constructed_variable(ResScalar,actualDestPtr,dest.size(),
201  EvalToDest ? dest.data() : static_dest.data());
202 
203  ei_declare_aligned_stack_constructed_variable(RhsScalar,actualRhsPtr,rhs.size(),
204  UseRhs ? const_cast<RhsScalar*>(rhs.data()) : static_rhs.data());
205 
206  if(!EvalToDest)
207  {
208  #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN
209  Index size = dest.size();
210  EIGEN_DENSE_STORAGE_CTOR_PLUGIN
211  #endif
212  MappedDest(actualDestPtr, dest.size()) = dest;
213  }
214 
215  if(!UseRhs)
216  {
217  #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN
218  Index size = rhs.size();
219  EIGEN_DENSE_STORAGE_CTOR_PLUGIN
220  #endif
221  Map<typename ActualRhsTypeCleaned::PlainObject>(actualRhsPtr, rhs.size()) = rhs;
222  }
223 
224 
226  int(LhsUpLo), bool(LhsBlasTraits::NeedToConjugate), bool(RhsBlasTraits::NeedToConjugate)>::run
227  (
228  lhs.rows(), // size
229  &lhs.coeffRef(0,0), lhs.outerStride(), // lhs info
230  actualRhsPtr, // rhs info
231  actualDestPtr, // result info
232  actualAlpha // scale factor
233  );
234 
235  if(!EvalToDest)
236  dest = MappedDest(actualDestPtr, dest.size());
237  }
238 };
239 
240 template<typename Lhs, typename Rhs, int RhsMode>
241 struct selfadjoint_product_impl<Lhs,0,true,Rhs,RhsMode,false>
242 {
244  enum { RhsUpLo = RhsMode&(Upper|Lower) };
245 
246  template<typename Dest>
247  static void run(Dest& dest, const Lhs &a_lhs, const Rhs &a_rhs, const Scalar& alpha)
248  {
249  // let's simply transpose the product
250  Transpose<Dest> destT(dest);
252  Transpose<const Lhs>, 0, true>::run(destT, a_rhs.transpose(), a_lhs.transpose(), alpha);
253  }
254 };
255 
256 } // end namespace internal
257 
258 } // end namespace Eigen
259 
260 #endif // EIGEN_SELFADJOINT_MATRIX_VECTOR_H
static void run(Dest &dest, const Lhs &a_lhs, const Rhs &a_rhs, const Scalar &alpha)
internal::packet_traits< Scalar >::type Packet
SCALAR Scalar
Definition: bench_gemm.cpp:33
#define max(a, b)
Definition: datatypes.h:20
float real
Definition: datatypes.h:10
Expression of the product of two arbitrary matrices or vectors.
Definition: Product.h:71
return int(ret)+1
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
Expression of the transpose of a matrix.
Definition: Transpose.h:52
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:150
EIGEN_DEVICE_FUNC unpacket_traits< Packet >::type predux(const Packet &a)
#define EIGEN_LOGICAL_XOR(a, b)
Definition: Macros.h:900
double bound(double a, double min, double max)
Definition: PoseRTV.cpp:17
const unsigned int RowMajorBit
Definition: Constants.h:61
#define EIGEN_DONT_INLINE
Definition: Macros.h:517
static Index first_default_aligned(const DenseBase< Derived > &m)
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
static void run(Dest &dest, const Lhs &a_lhs, const Rhs &a_rhs, const Scalar &alpha)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
#define eigen_assert(x)
Definition: Macros.h:579
RealScalar alpha
#define EIGEN_RESTRICT
Definition: Macros.h:796
static EIGEN_DONT_INLINE void run(Index size, const Scalar *lhs, Index lhsStride, const Scalar *rhs, Scalar *res, Scalar alpha)
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:34
#define ei_declare_aligned_stack_constructed_variable(TYPE, NAME, SIZE, BUFFER)
Definition: Memory.h:644
EIGEN_DEVICE_FUNC void pstore(Scalar *to, const Packet &from)
EIGEN_STRONG_INLINE Scalar pmadd(const LhsScalar &x, const RhsScalar &y, const Scalar &c) const
Definition: BlasUtil.h:65
#define EIGEN_PLAIN_ENUM_MIN(a, b)
Definition: Macros.h:875
std::ptrdiff_t j
EIGEN_STRONG_INLINE Scalar pmul(const LhsScalar &x, const RhsScalar &y) const
Definition: BlasUtil.h:68
ScalarWithExceptions conj(const ScalarWithExceptions &x)
Definition: exceptions.cpp:74


gtsam
Author(s):
autogenerated on Sat May 8 2021 02:43:55