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
EIGEN_BLAS_FUNC(axpy)12 int EIGEN_BLAS_FUNC(axpy)(const int *n, const RealScalar *palpha, const RealScalar *px, const int *incx, RealScalar *py, const int *incy)
13 {
14 const Scalar* x = reinterpret_cast<const Scalar*>(px);
15 Scalar* y = reinterpret_cast<Scalar*>(py);
16 Scalar alpha = *reinterpret_cast<const Scalar*>(palpha);
17
18 if(*n<=0) return 0;
19
20 if(*incx==1 && *incy==1) make_vector(y,*n) += alpha * make_vector(x,*n);
21 else if(*incx>0 && *incy>0) make_vector(y,*n,*incy) += alpha * make_vector(x,*n,*incx);
22 else if(*incx>0 && *incy<0) make_vector(y,*n,-*incy).reverse() += alpha * make_vector(x,*n,*incx);
23 else if(*incx<0 && *incy>0) make_vector(y,*n,*incy) += alpha * make_vector(x,*n,-*incx).reverse();
24 else if(*incx<0 && *incy<0) make_vector(y,*n,-*incy).reverse() += alpha * make_vector(x,*n,-*incx).reverse();
25
26 return 0;
27 }
28
EIGEN_BLAS_FUNC(copy)29 int EIGEN_BLAS_FUNC(copy)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy)
30 {
31 if(*n<=0) return 0;
32
33 Scalar* x = reinterpret_cast<Scalar*>(px);
34 Scalar* y = reinterpret_cast<Scalar*>(py);
35
36 // be carefull, *incx==0 is allowed !!
37 if(*incx==1 && *incy==1)
38 make_vector(y,*n) = make_vector(x,*n);
39 else
40 {
41 if(*incx<0) x = x - (*n-1)*(*incx);
42 if(*incy<0) y = y - (*n-1)*(*incy);
43 for(int i=0;i<*n;++i)
44 {
45 *y = *x;
46 x += *incx;
47 y += *incy;
48 }
49 }
50
51 return 0;
52 }
53
EIGEN_CAT(EIGEN_CAT (i,SCALAR_SUFFIX),amax_)54 int EIGEN_CAT(EIGEN_CAT(i,SCALAR_SUFFIX),amax_)(int *n, RealScalar *px, int *incx)
55 {
56 if(*n<=0) return 0;
57 Scalar* x = reinterpret_cast<Scalar*>(px);
58
59 DenseIndex ret;
60 if(*incx==1) make_vector(x,*n).cwiseAbs().maxCoeff(&ret);
61 else make_vector(x,*n,std::abs(*incx)).cwiseAbs().maxCoeff(&ret);
62 return int(ret)+1;
63 }
64
EIGEN_CAT(EIGEN_CAT (i,SCALAR_SUFFIX),amin_)65 int EIGEN_CAT(EIGEN_CAT(i,SCALAR_SUFFIX),amin_)(int *n, RealScalar *px, int *incx)
66 {
67 if(*n<=0) return 0;
68 Scalar* x = reinterpret_cast<Scalar*>(px);
69
70 DenseIndex ret;
71 if(*incx==1) make_vector(x,*n).cwiseAbs().minCoeff(&ret);
72 else make_vector(x,*n,std::abs(*incx)).cwiseAbs().minCoeff(&ret);
73 return int(ret)+1;
74 }
75
EIGEN_BLAS_FUNC(rotg)76 int EIGEN_BLAS_FUNC(rotg)(RealScalar *pa, RealScalar *pb, RealScalar *pc, RealScalar *ps)
77 {
78 using std::sqrt;
79 using std::abs;
80
81 Scalar& a = *reinterpret_cast<Scalar*>(pa);
82 Scalar& b = *reinterpret_cast<Scalar*>(pb);
83 RealScalar* c = pc;
84 Scalar* s = reinterpret_cast<Scalar*>(ps);
85
86 #if !ISCOMPLEX
87 Scalar r,z;
88 Scalar aa = abs(a);
89 Scalar ab = abs(b);
90 if((aa+ab)==Scalar(0))
91 {
92 *c = 1;
93 *s = 0;
94 r = 0;
95 z = 0;
96 }
97 else
98 {
99 r = sqrt(a*a + b*b);
100 Scalar amax = aa>ab ? a : b;
101 r = amax>0 ? r : -r;
102 *c = a/r;
103 *s = b/r;
104 z = 1;
105 if (aa > ab) z = *s;
106 if (ab > aa && *c!=RealScalar(0))
107 z = Scalar(1)/ *c;
108 }
109 *pa = r;
110 *pb = z;
111 #else
112 Scalar alpha;
113 RealScalar norm,scale;
114 if(abs(a)==RealScalar(0))
115 {
116 *c = RealScalar(0);
117 *s = Scalar(1);
118 a = b;
119 }
120 else
121 {
122 scale = abs(a) + abs(b);
123 norm = scale*sqrt((numext::abs2(a/scale)) + (numext::abs2(b/scale)));
124 alpha = a/abs(a);
125 *c = abs(a)/norm;
126 *s = alpha*numext::conj(b)/norm;
127 a = alpha*norm;
128 }
129 #endif
130
131 // JacobiRotation<Scalar> r;
132 // r.makeGivens(a,b);
133 // *c = r.c();
134 // *s = r.s();
135
136 return 0;
137 }
138
EIGEN_BLAS_FUNC(scal)139 int EIGEN_BLAS_FUNC(scal)(int *n, RealScalar *palpha, RealScalar *px, int *incx)
140 {
141 if(*n<=0) return 0;
142
143 Scalar* x = reinterpret_cast<Scalar*>(px);
144 Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
145
146 if(*incx==1) make_vector(x,*n) *= alpha;
147 else make_vector(x,*n,std::abs(*incx)) *= alpha;
148
149 return 0;
150 }
151
EIGEN_BLAS_FUNC(swap)152 int EIGEN_BLAS_FUNC(swap)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy)
153 {
154 if(*n<=0) return 0;
155
156 Scalar* x = reinterpret_cast<Scalar*>(px);
157 Scalar* y = reinterpret_cast<Scalar*>(py);
158
159 if(*incx==1 && *incy==1) make_vector(y,*n).swap(make_vector(x,*n));
160 else if(*incx>0 && *incy>0) make_vector(y,*n,*incy).swap(make_vector(x,*n,*incx));
161 else if(*incx>0 && *incy<0) make_vector(y,*n,-*incy).reverse().swap(make_vector(x,*n,*incx));
162 else if(*incx<0 && *incy>0) make_vector(y,*n,*incy).swap(make_vector(x,*n,-*incx).reverse());
163 else if(*incx<0 && *incy<0) make_vector(y,*n,-*incy).reverse().swap(make_vector(x,*n,-*incx).reverse());
164
165 return 1;
166 }
167