ground/gcs/src/libs/eigen/blas/level2_real_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 // y = alpha*A*x + beta*y
  13 int EIGEN_BLAS_FUNC(symv) (const char *uplo, const int *n, const RealScalar *palpha, const RealScalar *pa, const int *lda,
  14                            const RealScalar *px, const int *incx, const RealScalar *pbeta, RealScalar *py, const int *incy)
  15 {
  16   typedef void (*functype)(int, const Scalar*, int, const Scalar*, Scalar*, Scalar);
  17   static const functype func[2] = {
  18     // array index: UP
  19     (internal::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Upper,false,false>::run),
  20     // array index: LO
  21     (internal::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Lower,false,false>::run),
  22   };
  23
  24   const Scalar* a = reinterpret_cast<const Scalar*>(pa);
  25   const Scalar* x = reinterpret_cast<const Scalar*>(px);
  26   Scalar* y = reinterpret_cast<Scalar*>(py);
  27   Scalar alpha  = *reinterpret_cast<const Scalar*>(palpha);
  28   Scalar beta   = *reinterpret_cast<const Scalar*>(pbeta);
  29
  30   // check arguments
  31   int info = 0;
  32   if(UPLO(*uplo)==INVALID)        info = 1;
  33   else if(*n<0)                   info = 2;
  34   else if(*lda<std::max(1,*n))    info = 5;
  35   else if(*incx==0)               info = 7;
  36   else if(*incy==0)               info = 10;
  37   if(info)
  38     return xerbla_(SCALAR_SUFFIX_UP"SYMV ",&info,6);
  39
  40   if(*n==0)
  41     return 0;
  42
  43   const Scalar* actual_x = get_compact_vector(x,*n,*incx);
  44   Scalar* actual_y = get_compact_vector(y,*n,*incy);
  45
  46   if(beta!=Scalar(1))
  47   {
  48     if(beta==Scalar(0)) make_vector(actual_y, *n).setZero();
  49     else                make_vector(actual_y, *n) *= beta;
  50   }
  51
  52   int code = UPLO(*uplo);
  53   if(code>=2 || func[code]==0)
  54     return 0;
  55
  56   func[code](*n, a, *lda, actual_x, actual_y, alpha);
  57
  58   if(actual_x!=x) delete[] actual_x;
  59   if(actual_y!=y) delete[] copy_back(actual_y,y,*n,*incy);
  60
  61   return 1;
  62 }
  63
  64 // C := alpha*x*x' + C
  65 int EIGEN_BLAS_FUNC(syr)(const char *uplo, const int *n, const RealScalar *palpha, const RealScalar *px, const int *incx, RealScalar *pc, const int *ldc)
  66 {
  67
  68   typedef void (*functype)(int, Scalar*, int, const Scalar*, const Scalar*, const Scalar&);
  69   static const functype func[2] = {
  70     // array index: UP
  71     (selfadjoint_rank1_update<Scalar,int,ColMajor,Upper,false,Conj>::run),
  72     // array index: LO
  73     (selfadjoint_rank1_update<Scalar,int,ColMajor,Lower,false,Conj>::run),
  74   };
  75
  76   const Scalar* x = reinterpret_cast<const Scalar*>(px);
  77   Scalar* c = reinterpret_cast<Scalar*>(pc);
  78   Scalar alpha = *reinterpret_cast<const Scalar*>(palpha);
  79
  80   int info = 0;
  81   if(UPLO(*uplo)==INVALID)                                            info = 1;
  82   else if(*n<0)                                                       info = 2;
  83   else if(*incx==0)                                                   info = 5;
  84   else if(*ldc<std::max(1,*n))                                        info = 7;
  85   if(info)
  86     return xerbla_(SCALAR_SUFFIX_UP"SYR  ",&info,6);
  87
  88   if(*n==0 || alpha==Scalar(0)) return 1;
  89
  90   // if the increment is not 1, let's copy it to a temporary vector to enable vectorization
  91   const Scalar* x_cpy = get_compact_vector(x,*n,*incx);
  92
  93   int code = UPLO(*uplo);
  94   if(code>=2 || func[code]==0)
  95     return 0;
  96
  97   func[code](*n, c, *ldc, x_cpy, x_cpy, alpha);
  98
  99   if(x_cpy!=x)  delete[] x_cpy;
 100
 101   return 1;
 102 }
 103
 104 // C := alpha*x*y' + alpha*y*x' + C
 105 int EIGEN_BLAS_FUNC(syr2)(const char *uplo, const int *n, const RealScalar *palpha, const RealScalar *px, const int *incx, const RealScalar *py, const int *incy, RealScalar *pc, const int *ldc)
 106 {
 107   typedef void (*functype)(int, Scalar*, int, const Scalar*, const Scalar*, Scalar);
 108   static const functype func[2] = {
 109     // array index: UP
 110     (internal::rank2_update_selector<Scalar,int,Upper>::run),
 111     // array index: LO
 112     (internal::rank2_update_selector<Scalar,int,Lower>::run),
 113   };
 114
 115   const Scalar* x = reinterpret_cast<const Scalar*>(px);
 116   const Scalar* y = reinterpret_cast<const Scalar*>(py);
 117   Scalar* c = reinterpret_cast<Scalar*>(pc);
 118   Scalar alpha = *reinterpret_cast<const Scalar*>(palpha);
 119
 120   int info = 0;
 121   if(UPLO(*uplo)==INVALID)                                            info = 1;
 122   else if(*n<0)                                                       info = 2;
 123   else if(*incx==0)                                                   info = 5;
 124   else if(*incy==0)                                                   info = 7;
 125   else if(*ldc<std::max(1,*n))                                        info = 9;
 126   if(info)
 127     return xerbla_(SCALAR_SUFFIX_UP"SYR2 ",&info,6);
 128
 129   if(alpha==Scalar(0))
 130     return 1;
 131
 132   const Scalar* x_cpy = get_compact_vector(x,*n,*incx);
 133   const Scalar* y_cpy = get_compact_vector(y,*n,*incy);
 134
 135   int code = UPLO(*uplo);
 136   if(code>=2 || func[code]==0)
 137     return 0;
 138
 139   func[code](*n, c, *ldc, x_cpy, y_cpy, alpha);
 140
 141   if(x_cpy!=x)  delete[] x_cpy;
 142   if(y_cpy!=y)  delete[] y_cpy;
 143
 144 //   int code = UPLO(*uplo);
 145 //   if(code>=2 || func[code]==0)
 146 //     return 0;
 147
 148 //   func[code](*n, a, *inca, b, *incb, c, *ldc, alpha);
 149   return 1;
 150 }
 151
 152 /**  DSBMV  performs the matrix-vector  operation
 153   *
 154   *     y := alpha*A*x + beta*y,
 155   *
 156   *  where alpha and beta are scalars, x and y are n element vectors and
 157   *  A is an n by n symmetric band matrix, with k super-diagonals.
 158   */
 159 // int EIGEN_BLAS_FUNC(sbmv)( char *uplo, int *n, int *k, RealScalar *alpha, RealScalar *a, int *lda,
 160 //                            RealScalar *x, int *incx, RealScalar *beta, RealScalar *y, int *incy)
 161 // {
 162 //   return 1;
 163 // }
 164
 165
 166 /**  DSPMV  performs the matrix-vector operation
 167   *
 168   *     y := alpha*A*x + beta*y,
 169   *
 170   *  where alpha and beta are scalars, x and y are n element vectors and
 171   *  A is an n by n symmetric matrix, supplied in packed form.
 172   *
 173   */
 174 // int EIGEN_BLAS_FUNC(spmv)(char *uplo, int *n, RealScalar *alpha, RealScalar *ap, RealScalar *x, int *incx, RealScalar *beta, RealScalar *y, int *incy)
 175 // {
 176 //   return 1;
 177 // }
 178
 179 /**  DSPR    performs the symmetric rank 1 operation
 180   *
 181   *     A := alpha*x*x' + A,
 182   *
 183   *  where alpha is a real scalar, x is an n element vector and A is an
 184   *  n by n symmetric matrix, supplied in packed form.
 185   */
 186 int EIGEN_BLAS_FUNC(spr)(char *uplo, int *n, Scalar *palpha, Scalar *px, int *incx, Scalar *pap)
 187 {
 188   typedef void (*functype)(int, Scalar*, const Scalar*, Scalar);
 189   static const functype func[2] = {
 190     // array index: UP
 191     (internal::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Upper,false,false>::run),
 192     // array index: LO
 193     (internal::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Lower,false,false>::run),
 194   };
 195
 196   Scalar* x = reinterpret_cast<Scalar*>(px);
 197   Scalar* ap = reinterpret_cast<Scalar*>(pap);
 198   Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
 199
 200   int info = 0;
 201   if(UPLO(*uplo)==INVALID)                                            info = 1;
 202   else if(*n<0)                                                       info = 2;
 203   else if(*incx==0)                                                   info = 5;
 204   if(info)
 205     return xerbla_(SCALAR_SUFFIX_UP"SPR  ",&info,6);
 206
 207   if(alpha==Scalar(0))
 208     return 1;
 209
 210   Scalar* x_cpy = get_compact_vector(x, *n, *incx);
 211
 212   int code = UPLO(*uplo);
 213   if(code>=2 || func[code]==0)
 214     return 0;
 215
 216   func[code](*n, ap, x_cpy, alpha);
 217
 218   if(x_cpy!=x)  delete[] x_cpy;
 219
 220   return 1;
 221 }
 222
 223 /**  DSPR2  performs the symmetric rank 2 operation
 224   *
 225   *     A := alpha*x*y' + alpha*y*x' + A,
 226   *
 227   *  where alpha is a scalar, x and y are n element vectors and A is an
 228   *  n by n symmetric matrix, supplied in packed form.
 229   */
 230 int EIGEN_BLAS_FUNC(spr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pap)
 231 {
 232   typedef void (*functype)(int, Scalar*, const Scalar*, const Scalar*, Scalar);
 233   static const functype func[2] = {
 234     // array index: UP
 235     (internal::packed_rank2_update_selector<Scalar,int,Upper>::run),
 236     // array index: LO
 237     (internal::packed_rank2_update_selector<Scalar,int,Lower>::run),
 238   };
 239
 240   Scalar* x = reinterpret_cast<Scalar*>(px);
 241   Scalar* y = reinterpret_cast<Scalar*>(py);
 242   Scalar* ap = reinterpret_cast<Scalar*>(pap);
 243   Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
 244
 245   int info = 0;
 246   if(UPLO(*uplo)==INVALID)                                            info = 1;
 247   else if(*n<0)                                                       info = 2;
 248   else if(*incx==0)                                                   info = 5;
 249   else if(*incy==0)                                                   info = 7;
 250   if(info)
 251     return xerbla_(SCALAR_SUFFIX_UP"SPR2 ",&info,6);
 252
 253   if(alpha==Scalar(0))
 254     return 1;
 255
 256   Scalar* x_cpy = get_compact_vector(x, *n, *incx);
 257   Scalar* y_cpy = get_compact_vector(y, *n, *incy);
 258
 259   int code = UPLO(*uplo);
 260   if(code>=2 || func[code]==0)
 261     return 0;
 262
 263   func[code](*n, ap, x_cpy, y_cpy, alpha);
 264
 265   if(x_cpy!=x)  delete[] x_cpy;
 266   if(y_cpy!=y)  delete[] y_cpy;
 267
 268   return 1;
 269 }
 270
 271 /**  DGER   performs the rank 1 operation
 272   *
 273   *     A := alpha*x*y' + A,
 274   *
 275   *  where alpha is a scalar, x is an m element vector, y is an n element
 276   *  vector and A is an m by n matrix.
 277   */
 278 int EIGEN_BLAS_FUNC(ger)(int *m, int *n, Scalar *palpha, Scalar *px, int *incx, Scalar *py, int *incy, Scalar *pa, int *lda)
 279 {
 280   Scalar* x = reinterpret_cast<Scalar*>(px);
 281   Scalar* y = reinterpret_cast<Scalar*>(py);
 282   Scalar* a = reinterpret_cast<Scalar*>(pa);
 283   Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
 284
 285   int info = 0;
 286        if(*m<0)                                                       info = 1;
 287   else if(*n<0)                                                       info = 2;
 288   else if(*incx==0)                                                   info = 5;
 289   else if(*incy==0)                                                   info = 7;
 290   else if(*lda<std::max(1,*m))                                        info = 9;
 291   if(info)
 292     return xerbla_(SCALAR_SUFFIX_UP"GER  ",&info,6);
 293
 294   if(alpha==Scalar(0))
 295     return 1;
 296
 297   Scalar* x_cpy = get_compact_vector(x,*m,*incx);
 298   Scalar* y_cpy = get_compact_vector(y,*n,*incy);
 299
 300   internal::general_rank1_update<Scalar,int,ColMajor,false,false>::run(*m, *n, a, *lda, x_cpy, y_cpy, alpha);
 301
 302   if(x_cpy!=x)  delete[] x_cpy;
 303   if(y_cpy!=y)  delete[] y_cpy;
 304
 305   return 1;
 306 }