level1_impl.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2009-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 //
00006 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #include "common.h"
00026 
00027 int EIGEN_BLAS_FUNC(axpy)(int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy)
00028 {
00029   Scalar* x = reinterpret_cast<Scalar*>(px);
00030   Scalar* y = reinterpret_cast<Scalar*>(py);
00031   Scalar alpha  = *reinterpret_cast<Scalar*>(palpha);
00032 
00033   if(*n<=0) return 0;
00034 
00035   if(*incx==1 && *incy==1)    vector(y,*n) += alpha * vector(x,*n);
00036   else if(*incx>0 && *incy>0) vector(y,*n,*incy) += alpha * vector(x,*n,*incx);
00037   else if(*incx>0 && *incy<0) vector(y,*n,-*incy).reverse() += alpha * vector(x,*n,*incx);
00038   else if(*incx<0 && *incy>0) vector(y,*n,*incy) += alpha * vector(x,*n,-*incx).reverse();
00039   else if(*incx<0 && *incy<0) vector(y,*n,-*incy).reverse() += alpha * vector(x,*n,-*incx).reverse();
00040 
00041   return 0;
00042 }
00043 
00044 int EIGEN_BLAS_FUNC(copy)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy)
00045 {
00046   if(*n<=0) return 0;
00047 
00048   Scalar* x = reinterpret_cast<Scalar*>(px);
00049   Scalar* y = reinterpret_cast<Scalar*>(py);
00050 
00051   // be carefull, *incx==0 is allowed !!
00052   if(*incx==1 && *incy==1)
00053     vector(y,*n) = vector(x,*n);
00054   else
00055   {
00056     if(*incx<0) x = x - (*n-1)*(*incx);
00057     if(*incy<0) y = y - (*n-1)*(*incy);
00058     for(int i=0;i<*n;++i)
00059     {
00060       *y = *x;
00061       x += *incx;
00062       y += *incy;
00063     }
00064   }
00065 
00066   return 0;
00067 }
00068 
00069 int EIGEN_CAT(EIGEN_CAT(i,SCALAR_SUFFIX),amax_)(int *n, RealScalar *px, int *incx)
00070 {
00071   if(*n<=0) return 0;
00072   Scalar* x = reinterpret_cast<Scalar*>(px);
00073 
00074   DenseIndex ret;
00075   if(*incx==1)  vector(x,*n).cwiseAbs().maxCoeff(&ret);
00076   else          vector(x,*n,std::abs(*incx)).cwiseAbs().maxCoeff(&ret);
00077   return ret+1;
00078 }
00079 
00080 int EIGEN_CAT(EIGEN_CAT(i,SCALAR_SUFFIX),amin_)(int *n, RealScalar *px, int *incx)
00081 {
00082   if(*n<=0) return 0;
00083   Scalar* x = reinterpret_cast<Scalar*>(px);
00084   
00085   DenseIndex ret;
00086   if(*incx==1)  vector(x,*n).cwiseAbs().minCoeff(&ret);
00087   else          vector(x,*n,std::abs(*incx)).cwiseAbs().minCoeff(&ret);
00088   return ret+1;
00089 }
00090 
00091 int EIGEN_BLAS_FUNC(rotg)(RealScalar *pa, RealScalar *pb, RealScalar *pc, RealScalar *ps)
00092 {
00093   Scalar& a = *reinterpret_cast<Scalar*>(pa);
00094   Scalar& b = *reinterpret_cast<Scalar*>(pb);
00095   RealScalar* c = pc;
00096   Scalar* s = reinterpret_cast<Scalar*>(ps);
00097 
00098   #if !ISCOMPLEX
00099   Scalar r,z;
00100   Scalar aa = internal::abs(a);
00101   Scalar ab = internal::abs(b);
00102   if((aa+ab)==Scalar(0))
00103   {
00104     *c = 1;
00105     *s = 0;
00106     r = 0;
00107     z = 0;
00108   }
00109   else
00110   {
00111     r = internal::sqrt(a*a + b*b);
00112     Scalar amax = aa>ab ? a : b;
00113     r = amax>0 ? r : -r;
00114     *c = a/r;
00115     *s = b/r;
00116     z = 1;
00117     if (aa > ab) z = *s;
00118     if (ab > aa && *c!=RealScalar(0))
00119       z = Scalar(1)/ *c;
00120   }
00121   *pa = r;
00122   *pb = z;
00123   #else
00124   Scalar alpha;
00125   RealScalar norm,scale;
00126   if(internal::abs(a)==RealScalar(0))
00127   {
00128     *c = RealScalar(0);
00129     *s = Scalar(1);
00130     a = b;
00131   }
00132   else
00133   {
00134     scale = internal::abs(a) + internal::abs(b);
00135     norm = scale*internal::sqrt((internal::abs2(a/scale))+ (internal::abs2(b/scale)));
00136     alpha = a/internal::abs(a);
00137     *c = internal::abs(a)/norm;
00138     *s = alpha*internal::conj(b)/norm;
00139     a = alpha*norm;
00140   }
00141   #endif
00142 
00143 //   JacobiRotation<Scalar> r;
00144 //   r.makeGivens(a,b);
00145 //   *c = r.c();
00146 //   *s = r.s();
00147 
00148   return 0;
00149 }
00150 
00151 int EIGEN_BLAS_FUNC(scal)(int *n, RealScalar *palpha, RealScalar *px, int *incx)
00152 {
00153   if(*n<=0) return 0;
00154 
00155   Scalar* x = reinterpret_cast<Scalar*>(px);
00156   Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
00157 
00158   if(*incx==1)  vector(x,*n) *= alpha;
00159   else          vector(x,*n,std::abs(*incx)) *= alpha;
00160 
00161   return 0;
00162 }
00163 
00164 int EIGEN_BLAS_FUNC(swap)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy)
00165 {
00166   if(*n<=0) return 0;
00167 
00168   Scalar* x = reinterpret_cast<Scalar*>(px);
00169   Scalar* y = reinterpret_cast<Scalar*>(py);
00170 
00171   if(*incx==1 && *incy==1)    vector(y,*n).swap(vector(x,*n));
00172   else if(*incx>0 && *incy>0) vector(y,*n,*incy).swap(vector(x,*n,*incx));
00173   else if(*incx>0 && *incy<0) vector(y,*n,-*incy).reverse().swap(vector(x,*n,*incx));
00174   else if(*incx<0 && *incy>0) vector(y,*n,*incy).swap(vector(x,*n,-*incx).reverse());
00175   else if(*incx<0 && *incy<0) vector(y,*n,-*incy).reverse().swap(vector(x,*n,-*incx).reverse());
00176 
00177   return 1;
00178 }
00179 


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:32:56