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