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
12 // computes the sum of magnitudes of all vector elements or, for a complex vector x, the sum
13 // res = |Rex1| + |Imx1| + |Rex2| + |Imx2| + ... + |Rexn| + |Imxn|, where x is a vector of order n
EIGEN_BLAS_FUNC(asum)14 RealScalar EIGEN_BLAS_FUNC(asum)(int *n, RealScalar *px, int *incx)
15 {
16 // std::cerr << "_asum " << *n << " " << *incx << "\n";
17
18 Scalar* x = reinterpret_cast<Scalar*>(px);
19
20 if(*n<=0) return 0;
21
22 if(*incx==1) return make_vector(x,*n).cwiseAbs().sum();
23 else return make_vector(x,*n,std::abs(*incx)).cwiseAbs().sum();
24 }
25
26 // computes a vector-vector dot product.
EIGEN_BLAS_FUNC(dot)27 Scalar EIGEN_BLAS_FUNC(dot)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy)
28 {
29 // std::cerr << "_dot " << *n << " " << *incx << " " << *incy << "\n";
30
31 if(*n<=0) return 0;
32
33 Scalar* x = reinterpret_cast<Scalar*>(px);
34 Scalar* y = reinterpret_cast<Scalar*>(py);
35
36 if(*incx==1 && *incy==1) return (make_vector(x,*n).cwiseProduct(make_vector(y,*n))).sum();
37 else if(*incx>0 && *incy>0) return (make_vector(x,*n,*incx).cwiseProduct(make_vector(y,*n,*incy))).sum();
38 else if(*incx<0 && *incy>0) return (make_vector(x,*n,-*incx).reverse().cwiseProduct(make_vector(y,*n,*incy))).sum();
39 else if(*incx>0 && *incy<0) return (make_vector(x,*n,*incx).cwiseProduct(make_vector(y,*n,-*incy).reverse())).sum();
40 else if(*incx<0 && *incy<0) return (make_vector(x,*n,-*incx).reverse().cwiseProduct(make_vector(y,*n,-*incy).reverse())).sum();
41 else return 0;
42 }
43
44 // computes the Euclidean norm of a vector.
45 // FIXME
EIGEN_BLAS_FUNC(nrm2)46 Scalar EIGEN_BLAS_FUNC(nrm2)(int *n, RealScalar *px, int *incx)
47 {
48 // std::cerr << "_nrm2 " << *n << " " << *incx << "\n";
49 if(*n<=0) return 0;
50
51 Scalar* x = reinterpret_cast<Scalar*>(px);
52
53 if(*incx==1) return make_vector(x,*n).stableNorm();
54 else return make_vector(x,*n,std::abs(*incx)).stableNorm();
55 }
56
EIGEN_BLAS_FUNC(rot)57 int EIGEN_BLAS_FUNC(rot)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pc, RealScalar *ps)
58 {
59 // std::cerr << "_rot " << *n << " " << *incx << " " << *incy << "\n";
60 if(*n<=0) return 0;
61
62 Scalar* x = reinterpret_cast<Scalar*>(px);
63 Scalar* y = reinterpret_cast<Scalar*>(py);
64 Scalar c = *reinterpret_cast<Scalar*>(pc);
65 Scalar s = *reinterpret_cast<Scalar*>(ps);
66
67 StridedVectorType vx(make_vector(x,*n,std::abs(*incx)));
68 StridedVectorType vy(make_vector(y,*n,std::abs(*incy)));
69
70 Reverse<StridedVectorType> rvx(vx);
71 Reverse<StridedVectorType> rvy(vy);
72
73 if(*incx<0 && *incy>0) internal::apply_rotation_in_the_plane(rvx, vy, JacobiRotation<Scalar>(c,s));
74 else if(*incx>0 && *incy<0) internal::apply_rotation_in_the_plane(vx, rvy, JacobiRotation<Scalar>(c,s));
75 else internal::apply_rotation_in_the_plane(vx, vy, JacobiRotation<Scalar>(c,s));
76
77
78 return 0;
79 }
80
81 /*
82 // performs rotation of points in the modified plane.
83 int EIGEN_BLAS_FUNC(rotm)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *param)
84 {
85 Scalar* x = reinterpret_cast<Scalar*>(px);
86 Scalar* y = reinterpret_cast<Scalar*>(py);
87
88 // TODO
89
90 return 0;
91 }
92
93 // computes the modified parameters for a Givens rotation.
94 int EIGEN_BLAS_FUNC(rotmg)(RealScalar *d1, RealScalar *d2, RealScalar *x1, RealScalar *x2, RealScalar *param)
95 {
96 // TODO
97
98 return 0;
99 }
100 */
101