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