level1_cplx_impl.h (5646B)
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 12 struct scalar_norm1_op { 13 typedef RealScalar result_type; 14 EIGEN_EMPTY_STRUCT_CTOR(scalar_norm1_op) 15 inline RealScalar operator() (const Scalar& a) const { return numext::norm1(a); } 16 }; 17 namespace Eigen { 18 namespace internal { 19 template<> struct functor_traits<scalar_norm1_op > 20 { 21 enum { Cost = 3 * NumTraits<Scalar>::AddCost, PacketAccess = 0 }; 22 }; 23 } 24 } 25 26 // computes the sum of magnitudes of all vector elements or, for a complex vector x, the sum 27 // res = |Rex1| + |Imx1| + |Rex2| + |Imx2| + ... + |Rexn| + |Imxn|, where x is a vector of order n 28 RealScalar EIGEN_CAT(REAL_SCALAR_SUFFIX, EIGEN_BLAS_FUNC(asum))(int *n, RealScalar *px, int *incx) 29 { 30 // std::cerr << "__asum " << *n << " " << *incx << "\n"; 31 Complex* x = reinterpret_cast<Complex*>(px); 32 33 if(*n<=0) return 0; 34 35 if(*incx==1) return make_vector(x,*n).unaryExpr<scalar_norm1_op>().sum(); 36 else return make_vector(x,*n,std::abs(*incx)).unaryExpr<scalar_norm1_op>().sum(); 37 } 38 39 int EIGEN_CAT(i, EIGEN_BLAS_FUNC(amax))(int *n, RealScalar *px, int *incx) 40 { 41 if(*n<=0) return 0; 42 Scalar* x = reinterpret_cast<Scalar*>(px); 43 44 DenseIndex ret; 45 if(*incx==1) make_vector(x,*n).unaryExpr<scalar_norm1_op>().maxCoeff(&ret); 46 else make_vector(x,*n,std::abs(*incx)).unaryExpr<scalar_norm1_op>().maxCoeff(&ret); 47 return int(ret)+1; 48 } 49 50 int EIGEN_CAT(i, EIGEN_BLAS_FUNC(amin))(int *n, RealScalar *px, int *incx) 51 { 52 if(*n<=0) return 0; 53 Scalar* x = reinterpret_cast<Scalar*>(px); 54 55 DenseIndex ret; 56 if(*incx==1) make_vector(x,*n).unaryExpr<scalar_norm1_op>().minCoeff(&ret); 57 else make_vector(x,*n,std::abs(*incx)).unaryExpr<scalar_norm1_op>().minCoeff(&ret); 58 return int(ret)+1; 59 } 60 61 // computes a dot product of a conjugated vector with another vector. 62 int EIGEN_BLAS_FUNC(dotcw)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar* pres) 63 { 64 // std::cerr << "_dotc " << *n << " " << *incx << " " << *incy << "\n"; 65 Scalar* res = reinterpret_cast<Scalar*>(pres); 66 67 if(*n<=0) 68 { 69 *res = Scalar(0); 70 return 0; 71 } 72 73 Scalar* x = reinterpret_cast<Scalar*>(px); 74 Scalar* y = reinterpret_cast<Scalar*>(py); 75 76 if(*incx==1 && *incy==1) *res = (make_vector(x,*n).dot(make_vector(y,*n))); 77 else if(*incx>0 && *incy>0) *res = (make_vector(x,*n,*incx).dot(make_vector(y,*n,*incy))); 78 else if(*incx<0 && *incy>0) *res = (make_vector(x,*n,-*incx).reverse().dot(make_vector(y,*n,*incy))); 79 else if(*incx>0 && *incy<0) *res = (make_vector(x,*n,*incx).dot(make_vector(y,*n,-*incy).reverse())); 80 else if(*incx<0 && *incy<0) *res = (make_vector(x,*n,-*incx).reverse().dot(make_vector(y,*n,-*incy).reverse())); 81 return 0; 82 } 83 84 // computes a vector-vector dot product without complex conjugation. 85 int EIGEN_BLAS_FUNC(dotuw)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar* pres) 86 { 87 Scalar* res = reinterpret_cast<Scalar*>(pres); 88 89 if(*n<=0) 90 { 91 *res = Scalar(0); 92 return 0; 93 } 94 95 Scalar* x = reinterpret_cast<Scalar*>(px); 96 Scalar* y = reinterpret_cast<Scalar*>(py); 97 98 if(*incx==1 && *incy==1) *res = (make_vector(x,*n).cwiseProduct(make_vector(y,*n))).sum(); 99 else if(*incx>0 && *incy>0) *res = (make_vector(x,*n,*incx).cwiseProduct(make_vector(y,*n,*incy))).sum(); 100 else if(*incx<0 && *incy>0) *res = (make_vector(x,*n,-*incx).reverse().cwiseProduct(make_vector(y,*n,*incy))).sum(); 101 else if(*incx>0 && *incy<0) *res = (make_vector(x,*n,*incx).cwiseProduct(make_vector(y,*n,-*incy).reverse())).sum(); 102 else if(*incx<0 && *incy<0) *res = (make_vector(x,*n,-*incx).reverse().cwiseProduct(make_vector(y,*n,-*incy).reverse())).sum(); 103 return 0; 104 } 105 106 RealScalar EIGEN_CAT(REAL_SCALAR_SUFFIX, EIGEN_BLAS_FUNC(nrm2))(int *n, RealScalar *px, int *incx) 107 { 108 // std::cerr << "__nrm2 " << *n << " " << *incx << "\n"; 109 if(*n<=0) return 0; 110 111 Scalar* x = reinterpret_cast<Scalar*>(px); 112 113 if(*incx==1) 114 return make_vector(x,*n).stableNorm(); 115 116 return make_vector(x,*n,*incx).stableNorm(); 117 } 118 119 int EIGEN_BLAS_FUNC(EIGEN_CAT(REAL_SCALAR_SUFFIX, rot))(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pc, RealScalar *ps) 120 { 121 if(*n<=0) return 0; 122 123 Scalar* x = reinterpret_cast<Scalar*>(px); 124 Scalar* y = reinterpret_cast<Scalar*>(py); 125 RealScalar c = *pc; 126 RealScalar s = *ps; 127 128 StridedVectorType vx(make_vector(x,*n,std::abs(*incx))); 129 StridedVectorType vy(make_vector(y,*n,std::abs(*incy))); 130 131 Reverse<StridedVectorType> rvx(vx); 132 Reverse<StridedVectorType> rvy(vy); 133 134 // TODO implement mixed real-scalar rotations 135 if(*incx<0 && *incy>0) internal::apply_rotation_in_the_plane(rvx, vy, JacobiRotation<Scalar>(c,s)); 136 else if(*incx>0 && *incy<0) internal::apply_rotation_in_the_plane(vx, rvy, JacobiRotation<Scalar>(c,s)); 137 else internal::apply_rotation_in_the_plane(vx, vy, JacobiRotation<Scalar>(c,s)); 138 139 return 0; 140 } 141 142 int EIGEN_BLAS_FUNC(EIGEN_CAT(REAL_SCALAR_SUFFIX, scal))(int *n, RealScalar *palpha, RealScalar *px, int *incx) 143 { 144 if(*n<=0) return 0; 145 146 Scalar* x = reinterpret_cast<Scalar*>(px); 147 RealScalar alpha = *palpha; 148 149 // std::cerr << "__scal " << *n << " " << alpha << " " << *incx << "\n"; 150 151 if(*incx==1) make_vector(x,*n) *= alpha; 152 else make_vector(x,*n,std::abs(*incx)) *= alpha; 153 154 return 0; 155 }