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
12int 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
29int 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
54int 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
65int 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
76int 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
139int 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
152int 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