1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2009-2010 Gael Guennebaud <gael.guennebaud@inria.fr> 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include "common.h" 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathint EIGEN_BLAS_FUNC(axpy)(int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy) 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar* x = reinterpret_cast<Scalar*>(px); 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar* y = reinterpret_cast<Scalar*>(py); 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar alpha = *reinterpret_cast<Scalar*>(palpha); 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(*n<=0) return 0; 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(*incx==1 && *incy==1) vector(y,*n) += alpha * vector(x,*n); 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(*incx>0 && *incy>0) vector(y,*n,*incy) += alpha * vector(x,*n,*incx); 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(*incx>0 && *incy<0) vector(y,*n,-*incy).reverse() += alpha * vector(x,*n,*incx); 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(*incx<0 && *incy>0) vector(y,*n,*incy) += alpha * vector(x,*n,-*incx).reverse(); 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(*incx<0 && *incy<0) vector(y,*n,-*incy).reverse() += alpha * vector(x,*n,-*incx).reverse(); 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathint EIGEN_BLAS_FUNC(copy)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy) 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(*n<=0) return 0; 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar* x = reinterpret_cast<Scalar*>(px); 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar* y = reinterpret_cast<Scalar*>(py); 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // be carefull, *incx==0 is allowed !! 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(*incx==1 && *incy==1) 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath vector(y,*n) = vector(x,*n); 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(*incx<0) x = x - (*n-1)*(*incx); 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(*incy<0) y = y - (*n-1)*(*incy); 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0;i<*n;++i) 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath *y = *x; 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x += *incx; 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath y += *incy; 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathint EIGEN_CAT(EIGEN_CAT(i,SCALAR_SUFFIX),amax_)(int *n, RealScalar *px, int *incx) 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(*n<=0) return 0; 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar* x = reinterpret_cast<Scalar*>(px); 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DenseIndex ret; 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(*incx==1) vector(x,*n).cwiseAbs().maxCoeff(&ret); 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else vector(x,*n,std::abs(*incx)).cwiseAbs().maxCoeff(&ret); 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return ret+1; 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathint EIGEN_CAT(EIGEN_CAT(i,SCALAR_SUFFIX),amin_)(int *n, RealScalar *px, int *incx) 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(*n<=0) return 0; 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar* x = reinterpret_cast<Scalar*>(px); 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DenseIndex ret; 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(*incx==1) vector(x,*n).cwiseAbs().minCoeff(&ret); 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else vector(x,*n,std::abs(*incx)).cwiseAbs().minCoeff(&ret); 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return ret+1; 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathint EIGEN_BLAS_FUNC(rotg)(RealScalar *pa, RealScalar *pb, RealScalar *pc, RealScalar *ps) 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::sqrt; 797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::abs; 807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar& a = *reinterpret_cast<Scalar*>(pa); 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar& b = *reinterpret_cast<Scalar*>(pb); 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar* c = pc; 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar* s = reinterpret_cast<Scalar*>(ps); 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath #if !ISCOMPLEX 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar r,z; 887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar aa = abs(a); 897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar ab = abs(b); 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if((aa+ab)==Scalar(0)) 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath *c = 1; 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath *s = 0; 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath r = 0; 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath z = 0; 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez r = sqrt(a*a + b*b); 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar amax = aa>ab ? a : b; 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath r = amax>0 ? r : -r; 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath *c = a/r; 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath *s = b/r; 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath z = 1; 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (aa > ab) z = *s; 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (ab > aa && *c!=RealScalar(0)) 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath z = Scalar(1)/ *c; 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath *pa = r; 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath *pb = z; 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath #else 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar alpha; 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar norm,scale; 1147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(abs(a)==RealScalar(0)) 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath *c = RealScalar(0); 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath *s = Scalar(1); 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath a = b; 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez scale = abs(a) + abs(b); 1237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez norm = scale*sqrt((numext::abs2(a/scale)) + (numext::abs2(b/scale))); 1247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez alpha = a/abs(a); 1257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez *c = abs(a)/norm; 1267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez *s = alpha*numext::conj(b)/norm; 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath a = alpha*norm; 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath #endif 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// JacobiRotation<Scalar> r; 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// r.makeGivens(a,b); 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// *c = r.c(); 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// *s = r.s(); 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathint EIGEN_BLAS_FUNC(scal)(int *n, RealScalar *palpha, RealScalar *px, int *incx) 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(*n<=0) return 0; 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar* x = reinterpret_cast<Scalar*>(px); 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar alpha = *reinterpret_cast<Scalar*>(palpha); 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(*incx==1) vector(x,*n) *= alpha; 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else vector(x,*n,std::abs(*incx)) *= alpha; 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathint EIGEN_BLAS_FUNC(swap)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy) 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(*n<=0) return 0; 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar* x = reinterpret_cast<Scalar*>(px); 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar* y = reinterpret_cast<Scalar*>(py); 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(*incx==1 && *incy==1) vector(y,*n).swap(vector(x,*n)); 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(*incx>0 && *incy>0) vector(y,*n,*incy).swap(vector(x,*n,*incx)); 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(*incx>0 && *incy<0) vector(y,*n,-*incy).reverse().swap(vector(x,*n,*incx)); 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(*incx<0 && *incy>0) vector(y,*n,*incy).swap(vector(x,*n,-*incx).reverse()); 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(*incx<0 && *incy<0) vector(y,*n,-*incy).reverse().swap(vector(x,*n,-*incx).reverse()); 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 1; 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 168