ground/gcs/src/libs/eigen/blas/level1_impl.h

   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 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
  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)  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
  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)  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
  76 int 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
 139 int 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
 152 int 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 }