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)(int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy)
13 {
14 Scalar* x = reinterpret_cast<Scalar*>(px);
15 Scalar* y = reinterpret_cast<Scalar*>(py);
16 Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
17
18 if(*n<=0) return 0;
19
20 if(*incx==1 && *incy==1) vector(y,*n) += alpha * vector(x,*n);
21 else if(*incx>0 && *incy>0) vector(y,*n,*incy) += alpha * vector(x,*n,*incx);
22 else if(*incx>0 && *incy<0) vector(y,*n,-*incy).reverse() += alpha * vector(x,*n,*incx);
23 else if(*incx<0 && *incy>0) vector(y,*n,*incy) += alpha * vector(x,*n,-*incx).reverse();
24 else if(*incx<0 && *incy<0) vector(y,*n,-*incy).reverse() += alpha * 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 vector(y,*n) = 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) vector(x,*n).cwiseAbs().maxCoeff(&ret);
61 else vector(x,*n,std::abs(*incx)).cwiseAbs().maxCoeff(&ret);
62 return 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) vector(x,*n).cwiseAbs().minCoeff(&ret);
72 else vector(x,*n,std::abs(*incx)).cwiseAbs().minCoeff(&ret);
73 return ret+1;
74 }
75
EIGEN_BLAS_FUNC(rotg)76 int EIGEN_BLAS_FUNC(rotg)(RealScalar *pa, RealScalar *pb, RealScalar *pc, RealScalar *ps)
77 {
78 Scalar& a = *reinterpret_cast<Scalar*>(pa);
79 Scalar& b = *reinterpret_cast<Scalar*>(pb);
80 RealScalar* c = pc;
81 Scalar* s = reinterpret_cast<Scalar*>(ps);
82
83 #if !ISCOMPLEX
84 Scalar r,z;
85 Scalar aa = internal::abs(a);
86 Scalar ab = internal::abs(b);
87 if((aa+ab)==Scalar(0))
88 {
89 *c = 1;
90 *s = 0;
91 r = 0;
92 z = 0;
93 }
94 else
95 {
96 r = internal::sqrt(a*a + b*b);
97 Scalar amax = aa>ab ? a : b;
98 r = amax>0 ? r : -r;
99 *c = a/r;
100 *s = b/r;
101 z = 1;
102 if (aa > ab) z = *s;
103 if (ab > aa && *c!=RealScalar(0))
104 z = Scalar(1)/ *c;
105 }
106 *pa = r;
107 *pb = z;
108 #else
109 Scalar alpha;
110 RealScalar norm,scale;
111 if(internal::abs(a)==RealScalar(0))
112 {
113 *c = RealScalar(0);
114 *s = Scalar(1);
115 a = b;
116 }
117 else
118 {
119 scale = internal::abs(a) + internal::abs(b);
120 norm = scale*internal::sqrt((internal::abs2(a/scale))+ (internal::abs2(b/scale)));
121 alpha = a/internal::abs(a);
122 *c = internal::abs(a)/norm;
123 *s = alpha*internal::conj(b)/norm;
124 a = alpha*norm;
125 }
126 #endif
127
128 // JacobiRotation<Scalar> r;
129 // r.makeGivens(a,b);
130 // *c = r.c();
131 // *s = r.s();
132
133 return 0;
134 }
135
EIGEN_BLAS_FUNC(scal)136 int EIGEN_BLAS_FUNC(scal)(int *n, RealScalar *palpha, RealScalar *px, int *incx)
137 {
138 if(*n<=0) return 0;
139
140 Scalar* x = reinterpret_cast<Scalar*>(px);
141 Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
142
143 if(*incx==1) vector(x,*n) *= alpha;
144 else vector(x,*n,std::abs(*incx)) *= alpha;
145
146 return 0;
147 }
148
EIGEN_BLAS_FUNC(swap)149 int EIGEN_BLAS_FUNC(swap)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy)
150 {
151 if(*n<=0) return 0;
152
153 Scalar* x = reinterpret_cast<Scalar*>(px);
154 Scalar* y = reinterpret_cast<Scalar*>(py);
155
156 if(*incx==1 && *incy==1) vector(y,*n).swap(vector(x,*n));
157 else if(*incx>0 && *incy>0) vector(y,*n,*incy).swap(vector(x,*n,*incx));
158 else if(*incx>0 && *incy<0) vector(y,*n,-*incy).reverse().swap(vector(x,*n,*incx));
159 else if(*incx<0 && *incy>0) vector(y,*n,*incy).swap(vector(x,*n,-*incx).reverse());
160 else if(*incx<0 && *incy<0) vector(y,*n,-*incy).reverse().swap(vector(x,*n,-*incx).reverse());
161
162 return 1;
163 }
164
165