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 careful, *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_BLAS_FUNC(rotg)54 int EIGEN_BLAS_FUNC(rotg)(RealScalar *pa, RealScalar *pb, RealScalar *pc, RealScalar *ps)
55 {
56 using std::sqrt;
57 using std::abs;
58
59 Scalar& a = *reinterpret_cast<Scalar*>(pa);
60 Scalar& b = *reinterpret_cast<Scalar*>(pb);
61 RealScalar* c = pc;
62 Scalar* s = reinterpret_cast<Scalar*>(ps);
63
64 #if !ISCOMPLEX
65 Scalar r,z;
66 Scalar aa = abs(a);
67 Scalar ab = abs(b);
68 if((aa+ab)==Scalar(0))
69 {
70 *c = 1;
71 *s = 0;
72 r = 0;
73 z = 0;
74 }
75 else
76 {
77 r = sqrt(a*a + b*b);
78 Scalar amax = aa>ab ? a : b;
79 r = amax>0 ? r : -r;
80 *c = a/r;
81 *s = b/r;
82 z = 1;
83 if (aa > ab) z = *s;
84 if (ab > aa && *c!=RealScalar(0))
85 z = Scalar(1)/ *c;
86 }
87 *pa = r;
88 *pb = z;
89 #else
90 Scalar alpha;
91 RealScalar norm,scale;
92 if(abs(a)==RealScalar(0))
93 {
94 *c = RealScalar(0);
95 *s = Scalar(1);
96 a = b;
97 }
98 else
99 {
100 scale = abs(a) + abs(b);
101 norm = scale*sqrt((numext::abs2(a/scale)) + (numext::abs2(b/scale)));
102 alpha = a/abs(a);
103 *c = abs(a)/norm;
104 *s = alpha*numext::conj(b)/norm;
105 a = alpha*norm;
106 }
107 #endif
108
109 // JacobiRotation<Scalar> r;
110 // r.makeGivens(a,b);
111 // *c = r.c();
112 // *s = r.s();
113
114 return 0;
115 }
116
EIGEN_BLAS_FUNC(scal)117 int EIGEN_BLAS_FUNC(scal)(int *n, RealScalar *palpha, RealScalar *px, int *incx)
118 {
119 if(*n<=0) return 0;
120
121 Scalar* x = reinterpret_cast<Scalar*>(px);
122 Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
123
124 if(*incx==1) make_vector(x,*n) *= alpha;
125 else make_vector(x,*n,std::abs(*incx)) *= alpha;
126
127 return 0;
128 }
129
EIGEN_BLAS_FUNC(swap)130 int EIGEN_BLAS_FUNC(swap)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy)
131 {
132 if(*n<=0) return 0;
133
134 Scalar* x = reinterpret_cast<Scalar*>(px);
135 Scalar* y = reinterpret_cast<Scalar*>(py);
136
137 if(*incx==1 && *incy==1) make_vector(y,*n).swap(make_vector(x,*n));
138 else if(*incx>0 && *incy>0) make_vector(y,*n,*incy).swap(make_vector(x,*n,*incx));
139 else if(*incx>0 && *incy<0) make_vector(y,*n,-*incy).reverse().swap(make_vector(x,*n,*incx));
140 else if(*incx<0 && *incy>0) make_vector(y,*n,*incy).swap(make_vector(x,*n,-*incx).reverse());
141 else if(*incx<0 && *incy<0) make_vector(y,*n,-*incy).reverse().swap(make_vector(x,*n,-*incx).reverse());
142
143 return 1;
144 }
145