level1_impl.h (3892B)
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 int EIGEN_BLAS_FUNC(axpy)(const int *n, const RealScalar *palpha, const RealScalar *px, const int *incx, RealScalar *py, const int *incy) 13 { 14 const Scalar* x = reinterpret_cast<const Scalar*>(px); 15 Scalar* y = reinterpret_cast<Scalar*>(py); 16 Scalar alpha = *reinterpret_cast<const Scalar*>(palpha); 17 18 if(*n<=0) return 0; 19 20 if(*incx==1 && *incy==1) make_vector(y,*n) += alpha * make_vector(x,*n); 21 else if(*incx>0 && *incy>0) make_vector(y,*n,*incy) += alpha * make_vector(x,*n,*incx); 22 else if(*incx>0 && *incy<0) make_vector(y,*n,-*incy).reverse() += alpha * make_vector(x,*n,*incx); 23 else if(*incx<0 && *incy>0) make_vector(y,*n,*incy) += alpha * make_vector(x,*n,-*incx).reverse(); 24 else if(*incx<0 && *incy<0) make_vector(y,*n,-*incy).reverse() += alpha * make_vector(x,*n,-*incx).reverse(); 25 26 return 0; 27 } 28 29 int EIGEN_BLAS_FUNC(copy)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy) 30 { 31 if(*n<=0) return 0; 32 33 Scalar* x = reinterpret_cast<Scalar*>(px); 34 Scalar* y = reinterpret_cast<Scalar*>(py); 35 36 // be careful, *incx==0 is allowed !! 37 if(*incx==1 && *incy==1) 38 make_vector(y,*n) = make_vector(x,*n); 39 else 40 { 41 if(*incx<0) x = x - (*n-1)*(*incx); 42 if(*incy<0) y = y - (*n-1)*(*incy); 43 for(int i=0;i<*n;++i) 44 { 45 *y = *x; 46 x += *incx; 47 y += *incy; 48 } 49 } 50 51 return 0; 52 } 53 54 int EIGEN_BLAS_FUNC(rotg)(RealScalar *pa, RealScalar *pb, RealScalar *pc, RealScalar *ps) 55 { 56 using std::sqrt; 57 using std::abs; 58 59 Scalar& a = *reinterpret_cast<Scalar*>(pa); 60 Scalar& b = *reinterpret_cast<Scalar*>(pb); 61 RealScalar* c = pc; 62 Scalar* s = reinterpret_cast<Scalar*>(ps); 63 64 #if !ISCOMPLEX 65 Scalar r,z; 66 Scalar aa = abs(a); 67 Scalar ab = abs(b); 68 if((aa+ab)==Scalar(0)) 69 { 70 *c = 1; 71 *s = 0; 72 r = 0; 73 z = 0; 74 } 75 else 76 { 77 r = sqrt(a*a + b*b); 78 Scalar amax = aa>ab ? a : b; 79 r = amax>0 ? r : -r; 80 *c = a/r; 81 *s = b/r; 82 z = 1; 83 if (aa > ab) z = *s; 84 if (ab > aa && *c!=RealScalar(0)) 85 z = Scalar(1)/ *c; 86 } 87 *pa = r; 88 *pb = z; 89 #else 90 Scalar alpha; 91 RealScalar norm,scale; 92 if(abs(a)==RealScalar(0)) 93 { 94 *c = RealScalar(0); 95 *s = Scalar(1); 96 a = b; 97 } 98 else 99 { 100 scale = abs(a) + abs(b); 101 norm = scale*sqrt((numext::abs2(a/scale)) + (numext::abs2(b/scale))); 102 alpha = a/abs(a); 103 *c = abs(a)/norm; 104 *s = alpha*numext::conj(b)/norm; 105 a = alpha*norm; 106 } 107 #endif 108 109 // JacobiRotation<Scalar> r; 110 // r.makeGivens(a,b); 111 // *c = r.c(); 112 // *s = r.s(); 113 114 return 0; 115 } 116 117 int EIGEN_BLAS_FUNC(scal)(int *n, RealScalar *palpha, RealScalar *px, int *incx) 118 { 119 if(*n<=0) return 0; 120 121 Scalar* x = reinterpret_cast<Scalar*>(px); 122 Scalar alpha = *reinterpret_cast<Scalar*>(palpha); 123 124 if(*incx==1) make_vector(x,*n) *= alpha; 125 else make_vector(x,*n,std::abs(*incx)) *= alpha; 126 127 return 0; 128 } 129 130 int EIGEN_BLAS_FUNC(swap)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy) 131 { 132 if(*n<=0) return 0; 133 134 Scalar* x = reinterpret_cast<Scalar*>(px); 135 Scalar* y = reinterpret_cast<Scalar*>(py); 136 137 if(*incx==1 && *incy==1) make_vector(y,*n).swap(make_vector(x,*n)); 138 else if(*incx>0 && *incy>0) make_vector(y,*n,*incy).swap(make_vector(x,*n,*incx)); 139 else if(*incx>0 && *incy<0) make_vector(y,*n,-*incy).reverse().swap(make_vector(x,*n,*incx)); 140 else if(*incx<0 && *incy>0) make_vector(y,*n,*incy).swap(make_vector(x,*n,-*incx).reverse()); 141 else if(*incx<0 && *incy<0) make_vector(y,*n,-*incy).reverse().swap(make_vector(x,*n,-*incx).reverse()); 142 143 return 1; 144 }