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 int EIGEN_BLAS_FUNC(axpy)(int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy) 13 { 14 Scalar* x = reinterpret_cast<Scalar*>(px); 15 Scalar* y = reinterpret_cast<Scalar*>(py); 16 Scalar alpha = *reinterpret_cast<Scalar*>(palpha); 17 18 if(*n<=0) return 0; 19 20 if(*incx==1 && *incy==1) vector(y,*n) += alpha * vector(x,*n); 21 else if(*incx>0 && *incy>0) vector(y,*n,*incy) += alpha * vector(x,*n,*incx); 22 else if(*incx>0 && *incy<0) vector(y,*n,-*incy).reverse() += alpha * vector(x,*n,*incx); 23 else if(*incx<0 && *incy>0) vector(y,*n,*incy) += alpha * vector(x,*n,-*incx).reverse(); 24 else if(*incx<0 && *incy<0) vector(y,*n,-*incy).reverse() += alpha * 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 carefull, *incx==0 is allowed !! 37 if(*incx==1 && *incy==1) 38 vector(y,*n) = 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_CAT(EIGEN_CAT(i,SCALAR_SUFFIX),amax_)(int *n, RealScalar *px, int *incx) 55 { 56 if(*n<=0) return 0; 57 Scalar* x = reinterpret_cast<Scalar*>(px); 58 59 DenseIndex ret; 60 if(*incx==1) vector(x,*n).cwiseAbs().maxCoeff(&ret); 61 else vector(x,*n,std::abs(*incx)).cwiseAbs().maxCoeff(&ret); 62 return ret+1; 63 } 64 65 int EIGEN_CAT(EIGEN_CAT(i,SCALAR_SUFFIX),amin_)(int *n, RealScalar *px, int *incx) 66 { 67 if(*n<=0) return 0; 68 Scalar* x = reinterpret_cast<Scalar*>(px); 69 70 DenseIndex ret; 71 if(*incx==1) vector(x,*n).cwiseAbs().minCoeff(&ret); 72 else vector(x,*n,std::abs(*incx)).cwiseAbs().minCoeff(&ret); 73 return ret+1; 74 } 75 76 int EIGEN_BLAS_FUNC(rotg)(RealScalar *pa, RealScalar *pb, RealScalar *pc, RealScalar *ps) 77 { 78 Scalar& a = *reinterpret_cast<Scalar*>(pa); 79 Scalar& b = *reinterpret_cast<Scalar*>(pb); 80 RealScalar* c = pc; 81 Scalar* s = reinterpret_cast<Scalar*>(ps); 82 83 #if !ISCOMPLEX 84 Scalar r,z; 85 Scalar aa = internal::abs(a); 86 Scalar ab = internal::abs(b); 87 if((aa+ab)==Scalar(0)) 88 { 89 *c = 1; 90 *s = 0; 91 r = 0; 92 z = 0; 93 } 94 else 95 { 96 r = internal::sqrt(a*a + b*b); 97 Scalar amax = aa>ab ? a : b; 98 r = amax>0 ? r : -r; 99 *c = a/r; 100 *s = b/r; 101 z = 1; 102 if (aa > ab) z = *s; 103 if (ab > aa && *c!=RealScalar(0)) 104 z = Scalar(1)/ *c; 105 } 106 *pa = r; 107 *pb = z; 108 #else 109 Scalar alpha; 110 RealScalar norm,scale; 111 if(internal::abs(a)==RealScalar(0)) 112 { 113 *c = RealScalar(0); 114 *s = Scalar(1); 115 a = b; 116 } 117 else 118 { 119 scale = internal::abs(a) + internal::abs(b); 120 norm = scale*internal::sqrt((internal::abs2(a/scale))+ (internal::abs2(b/scale))); 121 alpha = a/internal::abs(a); 122 *c = internal::abs(a)/norm; 123 *s = alpha*internal::conj(b)/norm; 124 a = alpha*norm; 125 } 126 #endif 127 128 // JacobiRotation<Scalar> r; 129 // r.makeGivens(a,b); 130 // *c = r.c(); 131 // *s = r.s(); 132 133 return 0; 134 } 135 136 int EIGEN_BLAS_FUNC(scal)(int *n, RealScalar *palpha, RealScalar *px, int *incx) 137 { 138 if(*n<=0) return 0; 139 140 Scalar* x = reinterpret_cast<Scalar*>(px); 141 Scalar alpha = *reinterpret_cast<Scalar*>(palpha); 142 143 if(*incx==1) vector(x,*n) *= alpha; 144 else vector(x,*n,std::abs(*incx)) *= alpha; 145 146 return 0; 147 } 148 149 int EIGEN_BLAS_FUNC(swap)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy) 150 { 151 if(*n<=0) return 0; 152 153 Scalar* x = reinterpret_cast<Scalar*>(px); 154 Scalar* y = reinterpret_cast<Scalar*>(py); 155 156 if(*incx==1 && *incy==1) vector(y,*n).swap(vector(x,*n)); 157 else if(*incx>0 && *incy>0) vector(y,*n,*incy).swap(vector(x,*n,*incx)); 158 else if(*incx>0 && *incy<0) vector(y,*n,-*incy).reverse().swap(vector(x,*n,*incx)); 159 else if(*incx<0 && *incy>0) vector(y,*n,*incy).swap(vector(x,*n,-*incx).reverse()); 160 else if(*incx<0 && *incy<0) vector(y,*n,-*incy).reverse().swap(vector(x,*n,-*incx).reverse()); 161 162 return 1; 163 } 164 165
