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 (at) 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 internal::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(EIGEN_CAT(REAL_SCALAR_SUFFIX,SCALAR_SUFFIX),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 vector(x,*n).unaryExpr<scalar_norm1_op>().sum(); 36 else return vector(x,*n,std::abs(*incx)).unaryExpr<scalar_norm1_op>().sum(); 37 } 38 39 // computes a dot product of a conjugated vector with another vector. 40 int EIGEN_BLAS_FUNC(dotcw)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar* pres) 41 { 42 // std::cerr << "_dotc " << *n << " " << *incx << " " << *incy << "\n"; 43 44 if(*n<=0) return 0; 45 46 Scalar* x = reinterpret_cast<Scalar*>(px); 47 Scalar* y = reinterpret_cast<Scalar*>(py); 48 Scalar* res = reinterpret_cast<Scalar*>(pres); 49 50 if(*incx==1 && *incy==1) *res = (vector(x,*n).dot(vector(y,*n))); 51 else if(*incx>0 && *incy>0) *res = (vector(x,*n,*incx).dot(vector(y,*n,*incy))); 52 else if(*incx<0 && *incy>0) *res = (vector(x,*n,-*incx).reverse().dot(vector(y,*n,*incy))); 53 else if(*incx>0 && *incy<0) *res = (vector(x,*n,*incx).dot(vector(y,*n,-*incy).reverse())); 54 else if(*incx<0 && *incy<0) *res = (vector(x,*n,-*incx).reverse().dot(vector(y,*n,-*incy).reverse())); 55 return 0; 56 } 57 58 // computes a vector-vector dot product without complex conjugation. 59 int EIGEN_BLAS_FUNC(dotuw)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar* pres) 60 { 61 // std::cerr << "_dotu " << *n << " " << *incx << " " << *incy << "\n"; 62 63 if(*n<=0) return 0; 64 65 Scalar* x = reinterpret_cast<Scalar*>(px); 66 Scalar* y = reinterpret_cast<Scalar*>(py); 67 Scalar* res = reinterpret_cast<Scalar*>(pres); 68 69 if(*incx==1 && *incy==1) *res = (vector(x,*n).cwiseProduct(vector(y,*n))).sum(); 70 else if(*incx>0 && *incy>0) *res = (vector(x,*n,*incx).cwiseProduct(vector(y,*n,*incy))).sum(); 71 else if(*incx<0 && *incy>0) *res = (vector(x,*n,-*incx).reverse().cwiseProduct(vector(y,*n,*incy))).sum(); 72 else if(*incx>0 && *incy<0) *res = (vector(x,*n,*incx).cwiseProduct(vector(y,*n,-*incy).reverse())).sum(); 73 else if(*incx<0 && *incy<0) *res = (vector(x,*n,-*incx).reverse().cwiseProduct(vector(y,*n,-*incy).reverse())).sum(); 74 return 0; 75 } 76 77 RealScalar EIGEN_CAT(EIGEN_CAT(REAL_SCALAR_SUFFIX,SCALAR_SUFFIX),nrm2_)(int *n, RealScalar *px, int *incx) 78 { 79 // std::cerr << "__nrm2 " << *n << " " << *incx << "\n"; 80 if(*n<=0) return 0; 81 82 Scalar* x = reinterpret_cast<Scalar*>(px); 83 84 if(*incx==1) 85 return vector(x,*n).stableNorm(); 86 87 return vector(x,*n,*incx).stableNorm(); 88 } 89 90 int EIGEN_CAT(EIGEN_CAT(SCALAR_SUFFIX,REAL_SCALAR_SUFFIX),rot_)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pc, RealScalar *ps) 91 { 92 if(*n<=0) return 0; 93 94 Scalar* x = reinterpret_cast<Scalar*>(px); 95 Scalar* y = reinterpret_cast<Scalar*>(py); 96 RealScalar c = *pc; 97 RealScalar s = *ps; 98 99 StridedVectorType vx(vector(x,*n,std::abs(*incx))); 100 StridedVectorType vy(vector(y,*n,std::abs(*incy))); 101 102 Reverse<StridedVectorType> rvx(vx); 103 Reverse<StridedVectorType> rvy(vy); 104 105 // TODO implement mixed real-scalar rotations 106 if(*incx<0 && *incy>0) internal::apply_rotation_in_the_plane(rvx, vy, JacobiRotation<Scalar>(c,s)); 107 else if(*incx>0 && *incy<0) internal::apply_rotation_in_the_plane(vx, rvy, JacobiRotation<Scalar>(c,s)); 108 else internal::apply_rotation_in_the_plane(vx, vy, JacobiRotation<Scalar>(c,s)); 109 110 return 0; 111 } 112 113 int EIGEN_CAT(EIGEN_CAT(SCALAR_SUFFIX,REAL_SCALAR_SUFFIX),scal_)(int *n, RealScalar *palpha, RealScalar *px, int *incx) 114 { 115 if(*n<=0) return 0; 116 117 Scalar* x = reinterpret_cast<Scalar*>(px); 118 RealScalar alpha = *palpha; 119 120 // std::cerr << "__scal " << *n << " " << alpha << " " << *incx << "\n"; 121 122 if(*incx==1) vector(x,*n) *= alpha; 123 else vector(x,*n,std::abs(*incx)) *= alpha; 124 125 return 0; 126 } 127 128
