1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2009-2010 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#include "common.h"
11
12struct 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 internal::norm1(a); }
16};
17namespace Eigen {
18  namespace internal {
19    template<> struct functor_traits<scalar_norm1_op >
20    {
21      enum { Cost = 3 * NumTraits<Scalar>::AddCost, PacketAccess = 0 };
22    };
23  }
24}
25
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
28RealScalar EIGEN_CAT(EIGEN_CAT(REAL_SCALAR_SUFFIX,SCALAR_SUFFIX),asum_)(int *n, RealScalar *px, int *incx)
29{
30//   std::cerr << "__asum " << *n << " " << *incx << "\n";
31  Complex* x = reinterpret_cast<Complex*>(px);
32
33  if(*n<=0) return 0;
34
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();
37}
38
39// computes a dot product of a conjugated vector with another vector.
40int EIGEN_BLAS_FUNC(dotcw)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar* pres)
41{
42//   std::cerr << "_dotc " << *n << " " << *incx << " " << *incy << "\n";
43
44  if(*n<=0) return 0;
45
46  Scalar* x = reinterpret_cast<Scalar*>(px);
47  Scalar* y = reinterpret_cast<Scalar*>(py);
48  Scalar* res = reinterpret_cast<Scalar*>(pres);
49
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()));
55  return 0;
56}
57
58// computes a vector-vector dot product without complex conjugation.
59int EIGEN_BLAS_FUNC(dotuw)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar* pres)
60{
61//   std::cerr << "_dotu " << *n << " " << *incx << " " << *incy << "\n";
62
63  if(*n<=0) return 0;
64
65  Scalar* x = reinterpret_cast<Scalar*>(px);
66  Scalar* y = reinterpret_cast<Scalar*>(py);
67  Scalar* res = reinterpret_cast<Scalar*>(pres);
68
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();
74  return 0;
75}
76
77RealScalar EIGEN_CAT(EIGEN_CAT(REAL_SCALAR_SUFFIX,SCALAR_SUFFIX),nrm2_)(int *n, RealScalar *px, int *incx)
78{
79//   std::cerr << "__nrm2 " << *n << " " << *incx << "\n";
80  if(*n<=0) return 0;
81
82  Scalar* x = reinterpret_cast<Scalar*>(px);
83
84  if(*incx==1)
85    return vector(x,*n).stableNorm();
86
87  return vector(x,*n,*incx).stableNorm();
88}
89
90int EIGEN_CAT(EIGEN_CAT(SCALAR_SUFFIX,REAL_SCALAR_SUFFIX),rot_)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pc, RealScalar *ps)
91{
92  if(*n<=0) return 0;
93
94  Scalar* x = reinterpret_cast<Scalar*>(px);
95  Scalar* y = reinterpret_cast<Scalar*>(py);
96  RealScalar c = *pc;
97  RealScalar s = *ps;
98
99  StridedVectorType vx(vector(x,*n,std::abs(*incx)));
100  StridedVectorType vy(vector(y,*n,std::abs(*incy)));
101
102  Reverse<StridedVectorType> rvx(vx);
103  Reverse<StridedVectorType> rvy(vy);
104
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));
109
110  return 0;
111}
112
113int EIGEN_CAT(EIGEN_CAT(SCALAR_SUFFIX,REAL_SCALAR_SUFFIX),scal_)(int *n, RealScalar *palpha, RealScalar *px, int *incx)
114{
115  if(*n<=0) return 0;
116
117  Scalar* x = reinterpret_cast<Scalar*>(px);
118  RealScalar alpha = *palpha;
119
120//   std::cerr << "__scal " << *n << " " << alpha << " " << *incx << "\n";
121
122  if(*incx==1)  vector(x,*n) *= alpha;
123  else          vector(x,*n,std::abs(*incx)) *= alpha;
124
125  return 0;
126}
127
128