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 12int 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 29int 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 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 54int 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) make_vector(x,*n).cwiseAbs().maxCoeff(&ret); 61 else make_vector(x,*n,std::abs(*incx)).cwiseAbs().maxCoeff(&ret); 62 return int(ret)+1; 63} 64 65int 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) make_vector(x,*n).cwiseAbs().minCoeff(&ret); 72 else make_vector(x,*n,std::abs(*incx)).cwiseAbs().minCoeff(&ret); 73 return int(ret)+1; 74} 75 76int EIGEN_BLAS_FUNC(rotg)(RealScalar *pa, RealScalar *pb, RealScalar *pc, RealScalar *ps) 77{ 78 using std::sqrt; 79 using std::abs; 80 81 Scalar& a = *reinterpret_cast<Scalar*>(pa); 82 Scalar& b = *reinterpret_cast<Scalar*>(pb); 83 RealScalar* c = pc; 84 Scalar* s = reinterpret_cast<Scalar*>(ps); 85 86 #if !ISCOMPLEX 87 Scalar r,z; 88 Scalar aa = abs(a); 89 Scalar ab = abs(b); 90 if((aa+ab)==Scalar(0)) 91 { 92 *c = 1; 93 *s = 0; 94 r = 0; 95 z = 0; 96 } 97 else 98 { 99 r = sqrt(a*a + b*b); 100 Scalar amax = aa>ab ? a : b; 101 r = amax>0 ? r : -r; 102 *c = a/r; 103 *s = b/r; 104 z = 1; 105 if (aa > ab) z = *s; 106 if (ab > aa && *c!=RealScalar(0)) 107 z = Scalar(1)/ *c; 108 } 109 *pa = r; 110 *pb = z; 111 #else 112 Scalar alpha; 113 RealScalar norm,scale; 114 if(abs(a)==RealScalar(0)) 115 { 116 *c = RealScalar(0); 117 *s = Scalar(1); 118 a = b; 119 } 120 else 121 { 122 scale = abs(a) + abs(b); 123 norm = scale*sqrt((numext::abs2(a/scale)) + (numext::abs2(b/scale))); 124 alpha = a/abs(a); 125 *c = abs(a)/norm; 126 *s = alpha*numext::conj(b)/norm; 127 a = alpha*norm; 128 } 129 #endif 130 131// JacobiRotation<Scalar> r; 132// r.makeGivens(a,b); 133// *c = r.c(); 134// *s = r.s(); 135 136 return 0; 137} 138 139int EIGEN_BLAS_FUNC(scal)(int *n, RealScalar *palpha, RealScalar *px, int *incx) 140{ 141 if(*n<=0) return 0; 142 143 Scalar* x = reinterpret_cast<Scalar*>(px); 144 Scalar alpha = *reinterpret_cast<Scalar*>(palpha); 145 146 if(*incx==1) make_vector(x,*n) *= alpha; 147 else make_vector(x,*n,std::abs(*incx)) *= alpha; 148 149 return 0; 150} 151 152int EIGEN_BLAS_FUNC(swap)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy) 153{ 154 if(*n<=0) return 0; 155 156 Scalar* x = reinterpret_cast<Scalar*>(px); 157 Scalar* y = reinterpret_cast<Scalar*>(py); 158 159 if(*incx==1 && *incy==1) make_vector(y,*n).swap(make_vector(x,*n)); 160 else if(*incx>0 && *incy>0) make_vector(y,*n,*incy).swap(make_vector(x,*n,*incx)); 161 else if(*incx>0 && *incy<0) make_vector(y,*n,-*incy).reverse().swap(make_vector(x,*n,*incx)); 162 else if(*incx<0 && *incy>0) make_vector(y,*n,*incy).swap(make_vector(x,*n,-*incx).reverse()); 163 else if(*incx<0 && *incy<0) make_vector(y,*n,-*incy).reverse().swap(make_vector(x,*n,-*incx).reverse()); 164 165 return 1; 166} 167