diff options
Diffstat (limited to 'blas/level1_real_impl.h')
-rw-r--r-- | blas/level1_real_impl.h | 115 |
1 files changed, 115 insertions, 0 deletions
diff --git a/blas/level1_real_impl.h b/blas/level1_real_impl.h new file mode 100644 index 000000000..e18c40f0a --- /dev/null +++ b/blas/level1_real_impl.h @@ -0,0 +1,115 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009-2010 Gael Guennebaud <gael.guennebaud@inria.fr> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#include "common.h" + +// computes the sum of magnitudes of all vector elements or, for a complex vector x, the sum +// res = |Rex1| + |Imx1| + |Rex2| + |Imx2| + ... + |Rexn| + |Imxn|, where x is a vector of order n +RealScalar EIGEN_BLAS_FUNC(asum)(int *n, RealScalar *px, int *incx) +{ +// std::cerr << "_asum " << *n << " " << *incx << "\n"; + + Scalar* x = reinterpret_cast<Scalar*>(px); + + if(*n<=0) return 0; + + if(*incx==1) return vector(x,*n).cwiseAbs().sum(); + else return vector(x,*n,std::abs(*incx)).cwiseAbs().sum(); +} + +// computes a vector-vector dot product. +Scalar EIGEN_BLAS_FUNC(dot)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy) +{ +// std::cerr << "_dot " << *n << " " << *incx << " " << *incy << "\n"; + + if(*n<=0) return 0; + + Scalar* x = reinterpret_cast<Scalar*>(px); + Scalar* y = reinterpret_cast<Scalar*>(py); + + if(*incx==1 && *incy==1) return (vector(x,*n).cwiseProduct(vector(y,*n))).sum(); + else if(*incx>0 && *incy>0) return (vector(x,*n,*incx).cwiseProduct(vector(y,*n,*incy))).sum(); + else if(*incx<0 && *incy>0) return (vector(x,*n,-*incx).reverse().cwiseProduct(vector(y,*n,*incy))).sum(); + else if(*incx>0 && *incy<0) return (vector(x,*n,*incx).cwiseProduct(vector(y,*n,-*incy).reverse())).sum(); + else if(*incx<0 && *incy<0) return (vector(x,*n,-*incx).reverse().cwiseProduct(vector(y,*n,-*incy).reverse())).sum(); + else return 0; +} + +// computes the Euclidean norm of a vector. +// FIXME +Scalar EIGEN_BLAS_FUNC(nrm2)(int *n, RealScalar *px, int *incx) +{ +// std::cerr << "_nrm2 " << *n << " " << *incx << "\n"; + if(*n<=0) return 0; + + Scalar* x = reinterpret_cast<Scalar*>(px); + + if(*incx==1) return vector(x,*n).stableNorm(); + else return vector(x,*n,std::abs(*incx)).stableNorm(); +} + +int EIGEN_BLAS_FUNC(rot)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pc, RealScalar *ps) +{ +// std::cerr << "_rot " << *n << " " << *incx << " " << *incy << "\n"; + if(*n<=0) return 0; + + Scalar* x = reinterpret_cast<Scalar*>(px); + Scalar* y = reinterpret_cast<Scalar*>(py); + Scalar c = *reinterpret_cast<Scalar*>(pc); + Scalar s = *reinterpret_cast<Scalar*>(ps); + + StridedVectorType vx(vector(x,*n,std::abs(*incx))); + StridedVectorType vy(vector(y,*n,std::abs(*incy))); + + Reverse<StridedVectorType> rvx(vx); + Reverse<StridedVectorType> rvy(vy); + + if(*incx<0 && *incy>0) internal::apply_rotation_in_the_plane(rvx, vy, JacobiRotation<Scalar>(c,s)); + else if(*incx>0 && *incy<0) internal::apply_rotation_in_the_plane(vx, rvy, JacobiRotation<Scalar>(c,s)); + else internal::apply_rotation_in_the_plane(vx, vy, JacobiRotation<Scalar>(c,s)); + + + return 0; +} + +/* +// performs rotation of points in the modified plane. +int EIGEN_BLAS_FUNC(rotm)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *param) +{ + Scalar* x = reinterpret_cast<Scalar*>(px); + Scalar* y = reinterpret_cast<Scalar*>(py); + + // TODO + + return 0; +} + +// computes the modified parameters for a Givens rotation. +int EIGEN_BLAS_FUNC(rotmg)(RealScalar *d1, RealScalar *d2, RealScalar *x1, RealScalar *x2, RealScalar *param) +{ + // TODO + + return 0; +} +*/ |