Squashed 'third_party/eigen/' content from commit 61d72f6
Change-Id: Iccc90fa0b55ab44037f018046d2fcffd90d9d025
git-subtree-dir: third_party/eigen
git-subtree-split: 61d72f6383cfa842868c53e30e087b0258177257
diff --git a/blas/level1_impl.h b/blas/level1_impl.h
new file mode 100644
index 0000000..b08c2f6
--- /dev/null
+++ b/blas/level1_impl.h
@@ -0,0 +1,167 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2009-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include "common.h"
+
+int EIGEN_BLAS_FUNC(axpy)(int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy)
+{
+ Scalar* x = reinterpret_cast<Scalar*>(px);
+ Scalar* y = reinterpret_cast<Scalar*>(py);
+ Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
+
+ if(*n<=0) return 0;
+
+ if(*incx==1 && *incy==1) vector(y,*n) += alpha * vector(x,*n);
+ else if(*incx>0 && *incy>0) vector(y,*n,*incy) += alpha * vector(x,*n,*incx);
+ else if(*incx>0 && *incy<0) vector(y,*n,-*incy).reverse() += alpha * vector(x,*n,*incx);
+ else if(*incx<0 && *incy>0) vector(y,*n,*incy) += alpha * vector(x,*n,-*incx).reverse();
+ else if(*incx<0 && *incy<0) vector(y,*n,-*incy).reverse() += alpha * vector(x,*n,-*incx).reverse();
+
+ return 0;
+}
+
+int EIGEN_BLAS_FUNC(copy)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy)
+{
+ if(*n<=0) return 0;
+
+ Scalar* x = reinterpret_cast<Scalar*>(px);
+ Scalar* y = reinterpret_cast<Scalar*>(py);
+
+ // be carefull, *incx==0 is allowed !!
+ if(*incx==1 && *incy==1)
+ vector(y,*n) = vector(x,*n);
+ else
+ {
+ if(*incx<0) x = x - (*n-1)*(*incx);
+ if(*incy<0) y = y - (*n-1)*(*incy);
+ for(int i=0;i<*n;++i)
+ {
+ *y = *x;
+ x += *incx;
+ y += *incy;
+ }
+ }
+
+ return 0;
+}
+
+int EIGEN_CAT(EIGEN_CAT(i,SCALAR_SUFFIX),amax_)(int *n, RealScalar *px, int *incx)
+{
+ if(*n<=0) return 0;
+ Scalar* x = reinterpret_cast<Scalar*>(px);
+
+ DenseIndex ret;
+ if(*incx==1) vector(x,*n).cwiseAbs().maxCoeff(&ret);
+ else vector(x,*n,std::abs(*incx)).cwiseAbs().maxCoeff(&ret);
+ return ret+1;
+}
+
+int EIGEN_CAT(EIGEN_CAT(i,SCALAR_SUFFIX),amin_)(int *n, RealScalar *px, int *incx)
+{
+ if(*n<=0) return 0;
+ Scalar* x = reinterpret_cast<Scalar*>(px);
+
+ DenseIndex ret;
+ if(*incx==1) vector(x,*n).cwiseAbs().minCoeff(&ret);
+ else vector(x,*n,std::abs(*incx)).cwiseAbs().minCoeff(&ret);
+ return ret+1;
+}
+
+int EIGEN_BLAS_FUNC(rotg)(RealScalar *pa, RealScalar *pb, RealScalar *pc, RealScalar *ps)
+{
+ using std::sqrt;
+ using std::abs;
+
+ Scalar& a = *reinterpret_cast<Scalar*>(pa);
+ Scalar& b = *reinterpret_cast<Scalar*>(pb);
+ RealScalar* c = pc;
+ Scalar* s = reinterpret_cast<Scalar*>(ps);
+
+ #if !ISCOMPLEX
+ Scalar r,z;
+ Scalar aa = abs(a);
+ Scalar ab = abs(b);
+ if((aa+ab)==Scalar(0))
+ {
+ *c = 1;
+ *s = 0;
+ r = 0;
+ z = 0;
+ }
+ else
+ {
+ r = sqrt(a*a + b*b);
+ Scalar amax = aa>ab ? a : b;
+ r = amax>0 ? r : -r;
+ *c = a/r;
+ *s = b/r;
+ z = 1;
+ if (aa > ab) z = *s;
+ if (ab > aa && *c!=RealScalar(0))
+ z = Scalar(1)/ *c;
+ }
+ *pa = r;
+ *pb = z;
+ #else
+ Scalar alpha;
+ RealScalar norm,scale;
+ if(abs(a)==RealScalar(0))
+ {
+ *c = RealScalar(0);
+ *s = Scalar(1);
+ a = b;
+ }
+ else
+ {
+ scale = abs(a) + abs(b);
+ norm = scale*sqrt((numext::abs2(a/scale)) + (numext::abs2(b/scale)));
+ alpha = a/abs(a);
+ *c = abs(a)/norm;
+ *s = alpha*numext::conj(b)/norm;
+ a = alpha*norm;
+ }
+ #endif
+
+// JacobiRotation<Scalar> r;
+// r.makeGivens(a,b);
+// *c = r.c();
+// *s = r.s();
+
+ return 0;
+}
+
+int EIGEN_BLAS_FUNC(scal)(int *n, RealScalar *palpha, RealScalar *px, int *incx)
+{
+ if(*n<=0) return 0;
+
+ Scalar* x = reinterpret_cast<Scalar*>(px);
+ Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
+
+ if(*incx==1) vector(x,*n) *= alpha;
+ else vector(x,*n,std::abs(*incx)) *= alpha;
+
+ return 0;
+}
+
+int EIGEN_BLAS_FUNC(swap)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy)
+{
+ if(*n<=0) return 0;
+
+ Scalar* x = reinterpret_cast<Scalar*>(px);
+ Scalar* y = reinterpret_cast<Scalar*>(py);
+
+ if(*incx==1 && *incy==1) vector(y,*n).swap(vector(x,*n));
+ else if(*incx>0 && *incy>0) vector(y,*n,*incy).swap(vector(x,*n,*incx));
+ else if(*incx>0 && *incy<0) vector(y,*n,-*incy).reverse().swap(vector(x,*n,*incx));
+ else if(*incx<0 && *incy>0) vector(y,*n,*incy).swap(vector(x,*n,-*incx).reverse());
+ else if(*incx<0 && *incy<0) vector(y,*n,-*incy).reverse().swap(vector(x,*n,-*incx).reverse());
+
+ return 1;
+}
+