ground/gcs/src/libs/eigen/blas/level2_cplx_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 /**  ZHEMV  performs the matrix-vector  operation
  13   *
  14   *     y := alpha*A*x + beta*y,
  15   *
  16   *  where alpha and beta are scalars, x and y are n element vectors and
  17   *  A is an n by n hermitian matrix.
  18   */
  19 int EIGEN_BLAS_FUNC(hemv)(const char *uplo, const int *n, const RealScalar *palpha, const RealScalar *pa, const int *lda,
  20                           const RealScalar *px, const int *incx, const RealScalar *pbeta, RealScalar *py, const int *incy)
  21 {
  22   typedef void (*functype)(int, const Scalar*, int, const Scalar*, Scalar*, Scalar);
  23   static const functype func[2] = {
  24     // array index: UP
  25     (internal::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Upper,false,false>::run),
  26     // array index: LO
  27     (internal::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Lower,false,false>::run),
  28   };
  29
  30   const Scalar* a = reinterpret_cast<const Scalar*>(pa);
  31   const Scalar* x = reinterpret_cast<const Scalar*>(px);
  32   Scalar* y = reinterpret_cast<Scalar*>(py);
  33   Scalar alpha  = *reinterpret_cast<const Scalar*>(palpha);
  34   Scalar beta   = *reinterpret_cast<const Scalar*>(pbeta);
  35
  36   // check arguments
  37   int info = 0;
  38   if(UPLO(*uplo)==INVALID)        info = 1;
  39   else if(*n<0)                   info = 2;
  40   else if(*lda<std::max(1,*n))    info = 5;
  41   else if(*incx==0)               info = 7;
  42   else if(*incy==0)               info = 10;
  43   if(info)
  44     return xerbla_(SCALAR_SUFFIX_UP"HEMV ",&info,6);
  45
  46   if(*n==0)
  47     return 1;
  48
  49   const Scalar* actual_x = get_compact_vector(x,*n,*incx);
  50   Scalar* actual_y = get_compact_vector(y,*n,*incy);
  51
  52   if(beta!=Scalar(1))
  53   {
  54     if(beta==Scalar(0)) make_vector(actual_y, *n).setZero();
  55     else                make_vector(actual_y, *n) *= beta;
  56   }
  57
  58   if(alpha!=Scalar(0))
  59   {
  60     int code = UPLO(*uplo);
  61     if(code>=2 || func[code]==0)
  62       return 0;
  63
  64     func[code](*n, a, *lda, actual_x, actual_y, alpha);
  65   }
  66
  67   if(actual_x!=x) delete[] actual_x;
  68   if(actual_y!=y) delete[] copy_back(actual_y,y,*n,*incy);
  69
  70   return 1;
  71 }
  72
  73 /**  ZHBMV  performs the matrix-vector  operation
  74   *
  75   *     y := alpha*A*x + beta*y,
  76   *
  77   *  where alpha and beta are scalars, x and y are n element vectors and
  78   *  A is an n by n hermitian band matrix, with k super-diagonals.
  79   */
  80 // int EIGEN_BLAS_FUNC(hbmv)(char *uplo, int *n, int *k, RealScalar *alpha, RealScalar *a, int *lda,
  81 //                           RealScalar *x, int *incx, RealScalar *beta, RealScalar *y, int *incy)
  82 // {
  83 //   return 1;
  84 // }
  85
  86 /**  ZHPMV  performs the matrix-vector operation
  87   *
  88   *     y := alpha*A*x + beta*y,
  89   *
  90   *  where alpha and beta are scalars, x and y are n element vectors and
  91   *  A is an n by n hermitian matrix, supplied in packed form.
  92   */
  93 // int EIGEN_BLAS_FUNC(hpmv)(char *uplo, int *n, RealScalar *alpha, RealScalar *ap, RealScalar *x, int *incx, RealScalar *beta, RealScalar *y, int *incy)
  94 // {
  95 //   return 1;
  96 // }
  97
  98 /**  ZHPR    performs the hermitian rank 1 operation
  99   *
 100   *     A := alpha*x*conjg( x' ) + A,
 101   *
 102   *  where alpha is a real scalar, x is an n element vector and A is an
 103   *  n by n hermitian matrix, supplied in packed form.
 104   */
 105 int EIGEN_BLAS_FUNC(hpr)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *pap)
 106 {
 107   typedef void (*functype)(int, Scalar*, const Scalar*, RealScalar);
 108   static const functype func[2] = {
 109     // array index: UP
 110     (internal::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Upper,false,Conj>::run),
 111     // array index: LO
 112     (internal::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Lower,false,Conj>::run),
 113   };
 114
 115   Scalar* x = reinterpret_cast<Scalar*>(px);
 116   Scalar* ap = reinterpret_cast<Scalar*>(pap);
 117   RealScalar alpha = *palpha;
 118
 119   int info = 0;
 120   if(UPLO(*uplo)==INVALID)                                            info = 1;
 121   else if(*n<0)                                                       info = 2;
 122   else if(*incx==0)                                                   info = 5;
 123   if(info)
 124     return xerbla_(SCALAR_SUFFIX_UP"HPR  ",&info,6);
 125
 126   if(alpha==Scalar(0))
 127     return 1;
 128
 129   Scalar* x_cpy = get_compact_vector(x, *n, *incx);
 130
 131   int code = UPLO(*uplo);
 132   if(code>=2 || func[code]==0)
 133     return 0;
 134
 135   func[code](*n, ap, x_cpy, alpha);
 136
 137   if(x_cpy!=x)  delete[] x_cpy;
 138
 139   return 1;
 140 }
 141
 142 /**  ZHPR2  performs the hermitian rank 2 operation
 143   *
 144   *     A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A,
 145   *
 146   *  where alpha is a scalar, x and y are n element vectors and A is an
 147   *  n by n hermitian matrix, supplied in packed form.
 148   */
 149 int EIGEN_BLAS_FUNC(hpr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pap)
 150 {
 151   typedef void (*functype)(int, Scalar*, const Scalar*, const Scalar*, Scalar);
 152   static const functype func[2] = {
 153     // array index: UP
 154     (internal::packed_rank2_update_selector<Scalar,int,Upper>::run),
 155     // array index: LO
 156     (internal::packed_rank2_update_selector<Scalar,int,Lower>::run),
 157   };
 158
 159   Scalar* x = reinterpret_cast<Scalar*>(px);
 160   Scalar* y = reinterpret_cast<Scalar*>(py);
 161   Scalar* ap = reinterpret_cast<Scalar*>(pap);
 162   Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
 163
 164   int info = 0;
 165   if(UPLO(*uplo)==INVALID)                                            info = 1;
 166   else if(*n<0)                                                       info = 2;
 167   else if(*incx==0)                                                   info = 5;
 168   else if(*incy==0)                                                   info = 7;
 169   if(info)
 170     return xerbla_(SCALAR_SUFFIX_UP"HPR2 ",&info,6);
 171
 172   if(alpha==Scalar(0))
 173     return 1;
 174
 175   Scalar* x_cpy = get_compact_vector(x, *n, *incx);
 176   Scalar* y_cpy = get_compact_vector(y, *n, *incy);
 177
 178   int code = UPLO(*uplo);
 179   if(code>=2 || func[code]==0)
 180     return 0;
 181
 182   func[code](*n, ap, x_cpy, y_cpy, alpha);
 183
 184   if(x_cpy!=x)  delete[] x_cpy;
 185   if(y_cpy!=y)  delete[] y_cpy;
 186
 187   return 1;
 188 }
 189
 190 /**  ZHER   performs the hermitian rank 1 operation
 191   *
 192   *     A := alpha*x*conjg( x' ) + A,
 193   *
 194   *  where alpha is a real scalar, x is an n element vector and A is an
 195   *  n by n hermitian matrix.
 196   */
 197 int EIGEN_BLAS_FUNC(her)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *pa, int *lda)
 198 {
 199   typedef void (*functype)(int, Scalar*, int, const Scalar*, const Scalar*, const Scalar&);
 200   static const functype func[2] = {
 201     // array index: UP
 202     (selfadjoint_rank1_update<Scalar,int,ColMajor,Upper,false,Conj>::run),
 203     // array index: LO
 204     (selfadjoint_rank1_update<Scalar,int,ColMajor,Lower,false,Conj>::run),
 205   };
 206
 207   Scalar* x = reinterpret_cast<Scalar*>(px);
 208   Scalar* a = reinterpret_cast<Scalar*>(pa);
 209   RealScalar alpha = *reinterpret_cast<RealScalar*>(palpha);
 210
 211   int info = 0;
 212   if(UPLO(*uplo)==INVALID)                                            info = 1;
 213   else if(*n<0)                                                       info = 2;
 214   else if(*incx==0)                                                   info = 5;
 215   else if(*lda<std::max(1,*n))                                        info = 7;
 216   if(info)
 217     return xerbla_(SCALAR_SUFFIX_UP"HER  ",&info,6);
 218
 219   if(alpha==RealScalar(0))
 220     return 1;
 221
 222   Scalar* x_cpy = get_compact_vector(x, *n, *incx);
 223
 224   int code = UPLO(*uplo);
 225   if(code>=2 || func[code]==0)
 226     return 0;
 227
 228   func[code](*n, a, *lda, x_cpy, x_cpy, alpha);
 229
 230   matrix(a,*n,*n,*lda).diagonal().imag().setZero();
 231
 232   if(x_cpy!=x)  delete[] x_cpy;
 233
 234   return 1;
 235 }
 236
 237 /**  ZHER2  performs the hermitian rank 2 operation
 238   *
 239   *     A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A,
 240   *
 241   *  where alpha is a scalar, x and y are n element vectors and A is an n
 242   *  by n hermitian matrix.
 243   */
 244 int EIGEN_BLAS_FUNC(her2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pa, int *lda)
 245 {
 246   typedef void (*functype)(int, Scalar*, int, const Scalar*, const Scalar*, Scalar);
 247   static const functype func[2] = {
 248     // array index: UP
 249     (internal::rank2_update_selector<Scalar,int,Upper>::run),
 250     // array index: LO
 251     (internal::rank2_update_selector<Scalar,int,Lower>::run),
 252   };
 253
 254   Scalar* x = reinterpret_cast<Scalar*>(px);
 255   Scalar* y = reinterpret_cast<Scalar*>(py);
 256   Scalar* a = reinterpret_cast<Scalar*>(pa);
 257   Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
 258
 259   int info = 0;
 260   if(UPLO(*uplo)==INVALID)                                            info = 1;
 261   else if(*n<0)                                                       info = 2;
 262   else if(*incx==0)                                                   info = 5;
 263   else if(*incy==0)                                                   info = 7;
 264   else if(*lda<std::max(1,*n))                                        info = 9;
 265   if(info)
 266     return xerbla_(SCALAR_SUFFIX_UP"HER2 ",&info,6);
 267
 268   if(alpha==Scalar(0))
 269     return 1;
 270
 271   Scalar* x_cpy = get_compact_vector(x, *n, *incx);
 272   Scalar* y_cpy = get_compact_vector(y, *n, *incy);
 273
 274   int code = UPLO(*uplo);
 275   if(code>=2 || func[code]==0)
 276     return 0;
 277
 278   func[code](*n, a, *lda, x_cpy, y_cpy, alpha);
 279
 280   matrix(a,*n,*n,*lda).diagonal().imag().setZero();
 281
 282   if(x_cpy!=x)  delete[] x_cpy;
 283   if(y_cpy!=y)  delete[] y_cpy;
 284
 285   return 1;
 286 }
 287
 288 /**  ZGERU  performs the rank 1 operation
 289   *
 290   *     A := alpha*x*y' + A,
 291   *
 292   *  where alpha is a scalar, x is an m element vector, y is an n element
 293   *  vector and A is an m by n matrix.
 294   */
 295 int EIGEN_BLAS_FUNC(geru)(int *m, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pa, int *lda)
 296 {
 297   Scalar* x = reinterpret_cast<Scalar*>(px);
 298   Scalar* y = reinterpret_cast<Scalar*>(py);
 299   Scalar* a = reinterpret_cast<Scalar*>(pa);
 300   Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
 301
 302   int info = 0;
 303        if(*m<0)                                                       info = 1;
 304   else if(*n<0)                                                       info = 2;
 305   else if(*incx==0)                                                   info = 5;
 306   else if(*incy==0)                                                   info = 7;
 307   else if(*lda<std::max(1,*m))                                        info = 9;
 308   if(info)
 309     return xerbla_(SCALAR_SUFFIX_UP"GERU ",&info,6);
 310
 311   if(alpha==Scalar(0))
 312     return 1;
 313
 314   Scalar* x_cpy = get_compact_vector(x,*m,*incx);
 315   Scalar* y_cpy = get_compact_vector(y,*n,*incy);
 316
 317   internal::general_rank1_update<Scalar,int,ColMajor,false,false>::run(*m, *n, a, *lda, x_cpy, y_cpy, alpha);
 318
 319   if(x_cpy!=x)  delete[] x_cpy;
 320   if(y_cpy!=y)  delete[] y_cpy;
 321
 322   return 1;
 323 }
 324
 325 /**  ZGERC  performs the rank 1 operation
 326   *
 327   *     A := alpha*x*conjg( y' ) + A,
 328   *
 329   *  where alpha is a scalar, x is an m element vector, y is an n element
 330   *  vector and A is an m by n matrix.
 331   */
 332 int EIGEN_BLAS_FUNC(gerc)(int *m, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pa, int *lda)
 333 {
 334   Scalar* x = reinterpret_cast<Scalar*>(px);
 335   Scalar* y = reinterpret_cast<Scalar*>(py);
 336   Scalar* a = reinterpret_cast<Scalar*>(pa);
 337   Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
 338
 339   int info = 0;
 340        if(*m<0)                                                       info = 1;
 341   else if(*n<0)                                                       info = 2;
 342   else if(*incx==0)                                                   info = 5;
 343   else if(*incy==0)                                                   info = 7;
 344   else if(*lda<std::max(1,*m))                                        info = 9;
 345   if(info)
 346     return xerbla_(SCALAR_SUFFIX_UP"GERC ",&info,6);
 347
 348   if(alpha==Scalar(0))
 349     return 1;
 350
 351   Scalar* x_cpy = get_compact_vector(x,*m,*incx);
 352   Scalar* y_cpy = get_compact_vector(y,*n,*incy);
 353
 354   internal::general_rank1_update<Scalar,int,ColMajor,false,Conj>::run(*m, *n, a, *lda, x_cpy, y_cpy, alpha);
 355
 356   if(x_cpy!=x)  delete[] x_cpy;
 357   if(y_cpy!=y)  delete[] y_cpy;
 358
 359   return 1;
 360 }