1 // This file is part of Eigen, a lightweight C++ template library
4 // Copyright (C) 2009-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
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/.
12 struct scalar_norm1_op
{
13 typedef RealScalar result_type
;
14 EIGEN_EMPTY_STRUCT_CTOR(scalar_norm1_op
)
15 inline RealScalar
operator() (const Scalar
& a
) const { return numext::norm1(a
); }
19 template<> struct functor_traits
<scalar_norm1_op
>
21 enum { Cost
= 3 * NumTraits
<Scalar
>::AddCost
, PacketAccess
= 0 };
26 // computes the sum of magnitudes of all vector elements or, for a complex vector x, the sum
27 // res = |Rex1| + |Imx1| + |Rex2| + |Imx2| + ... + |Rexn| + |Imxn|, where x is a vector of order n
28 RealScalar
EIGEN_CAT(EIGEN_CAT(REAL_SCALAR_SUFFIX
,SCALAR_SUFFIX
),asum_
)(int *n
, RealScalar
*px
, int *incx
)
30 // std::cerr << "__asum " << *n << " " << *incx << "\n";
31 Complex
* x
= reinterpret_cast<Complex
*>(px
);
35 if(*incx
==1) return vector(x
,*n
).unaryExpr
<scalar_norm1_op
>().sum();
36 else return vector(x
,*n
,std::abs(*incx
)).unaryExpr
<scalar_norm1_op
>().sum();
39 // computes a dot product of a conjugated vector with another vector.
40 int EIGEN_BLAS_FUNC(dotcw
)(int *n
, RealScalar
*px
, int *incx
, RealScalar
*py
, int *incy
, RealScalar
* pres
)
42 // std::cerr << "_dotc " << *n << " " << *incx << " " << *incy << "\n";
46 Scalar
* x
= reinterpret_cast<Scalar
*>(px
);
47 Scalar
* y
= reinterpret_cast<Scalar
*>(py
);
48 Scalar
* res
= reinterpret_cast<Scalar
*>(pres
);
50 if(*incx
==1 && *incy
==1) *res
= (vector(x
,*n
).dot(vector(y
,*n
)));
51 else if(*incx
>0 && *incy
>0) *res
= (vector(x
,*n
,*incx
).dot(vector(y
,*n
,*incy
)));
52 else if(*incx
<0 && *incy
>0) *res
= (vector(x
,*n
,-*incx
).reverse().dot(vector(y
,*n
,*incy
)));
53 else if(*incx
>0 && *incy
<0) *res
= (vector(x
,*n
,*incx
).dot(vector(y
,*n
,-*incy
).reverse()));
54 else if(*incx
<0 && *incy
<0) *res
= (vector(x
,*n
,-*incx
).reverse().dot(vector(y
,*n
,-*incy
).reverse()));
58 // computes a vector-vector dot product without complex conjugation.
59 int EIGEN_BLAS_FUNC(dotuw
)(int *n
, RealScalar
*px
, int *incx
, RealScalar
*py
, int *incy
, RealScalar
* pres
)
61 // std::cerr << "_dotu " << *n << " " << *incx << " " << *incy << "\n";
65 Scalar
* x
= reinterpret_cast<Scalar
*>(px
);
66 Scalar
* y
= reinterpret_cast<Scalar
*>(py
);
67 Scalar
* res
= reinterpret_cast<Scalar
*>(pres
);
69 if(*incx
==1 && *incy
==1) *res
= (vector(x
,*n
).cwiseProduct(vector(y
,*n
))).sum();
70 else if(*incx
>0 && *incy
>0) *res
= (vector(x
,*n
,*incx
).cwiseProduct(vector(y
,*n
,*incy
))).sum();
71 else if(*incx
<0 && *incy
>0) *res
= (vector(x
,*n
,-*incx
).reverse().cwiseProduct(vector(y
,*n
,*incy
))).sum();
72 else if(*incx
>0 && *incy
<0) *res
= (vector(x
,*n
,*incx
).cwiseProduct(vector(y
,*n
,-*incy
).reverse())).sum();
73 else if(*incx
<0 && *incy
<0) *res
= (vector(x
,*n
,-*incx
).reverse().cwiseProduct(vector(y
,*n
,-*incy
).reverse())).sum();
77 RealScalar
EIGEN_CAT(EIGEN_CAT(REAL_SCALAR_SUFFIX
,SCALAR_SUFFIX
),nrm2_
)(int *n
, RealScalar
*px
, int *incx
)
79 // std::cerr << "__nrm2 " << *n << " " << *incx << "\n";
82 Scalar
* x
= reinterpret_cast<Scalar
*>(px
);
85 return vector(x
,*n
).stableNorm();
87 return vector(x
,*n
,*incx
).stableNorm();
90 int EIGEN_CAT(EIGEN_CAT(SCALAR_SUFFIX
,REAL_SCALAR_SUFFIX
),rot_
)(int *n
, RealScalar
*px
, int *incx
, RealScalar
*py
, int *incy
, RealScalar
*pc
, RealScalar
*ps
)
94 Scalar
* x
= reinterpret_cast<Scalar
*>(px
);
95 Scalar
* y
= reinterpret_cast<Scalar
*>(py
);
99 StridedVectorType
vx(vector(x
,*n
,std::abs(*incx
)));
100 StridedVectorType
vy(vector(y
,*n
,std::abs(*incy
)));
102 Reverse
<StridedVectorType
> rvx(vx
);
103 Reverse
<StridedVectorType
> rvy(vy
);
105 // TODO implement mixed real-scalar rotations
106 if(*incx
<0 && *incy
>0) internal::apply_rotation_in_the_plane(rvx
, vy
, JacobiRotation
<Scalar
>(c
,s
));
107 else if(*incx
>0 && *incy
<0) internal::apply_rotation_in_the_plane(vx
, rvy
, JacobiRotation
<Scalar
>(c
,s
));
108 else internal::apply_rotation_in_the_plane(vx
, vy
, JacobiRotation
<Scalar
>(c
,s
));
113 int EIGEN_CAT(EIGEN_CAT(SCALAR_SUFFIX
,REAL_SCALAR_SUFFIX
),scal_
)(int *n
, RealScalar
*palpha
, RealScalar
*px
, int *incx
)
117 Scalar
* x
= reinterpret_cast<Scalar
*>(px
);
118 RealScalar alpha
= *palpha
;
120 // std::cerr << "__scal " << *n << " " << alpha << " " << *incx << "\n";
122 if(*incx
==1) vector(x
,*n
) *= alpha
;
123 else vector(x
,*n
,std::abs(*incx
)) *= alpha
;