Home | History | Annotate | Download | only in blas 
      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 (at) 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) (char *uplo, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *px, int *incx, RealScalar *pbeta, RealScalar *py, int *incy)
     14 {
     15   Scalar* a = reinterpret_cast<Scalar*>(pa);
     16   Scalar* x = reinterpret_cast<Scalar*>(px);
     17   Scalar* y = reinterpret_cast<Scalar*>(py);
     18   Scalar alpha  = *reinterpret_cast<Scalar*>(palpha);
     19   Scalar beta   = *reinterpret_cast<Scalar*>(pbeta);
     20 
     21   // check arguments
     22   int info = 0;
     23   if(UPLO(*uplo)==INVALID)        info = 1;
     24   else if(*n<0)                   info = 2;
     25   else if(*lda<std::max(1,*n))    info = 5;
     26   else if(*incx==0)               info = 7;
     27   else if(*incy==0)               info = 10;
     28   if(info)
     29     return xerbla_(SCALAR_SUFFIX_UP"SYMV ",&info,6);
     30 
     31   if(*n==0)
     32     return 0;
     33 
     34   Scalar* actual_x = get_compact_vector(x,*n,*incx);
     35   Scalar* actual_y = get_compact_vector(y,*n,*incy);
     36 
     37   if(beta!=Scalar(1))
     38   {
     39     if(beta==Scalar(0)) vector(actual_y, *n).setZero();
     40     else                vector(actual_y, *n) *= beta;
     41   }
     42 
     43   // TODO performs a direct call to the underlying implementation function
     44        if(UPLO(*uplo)==UP) vector(actual_y,*n).noalias() += matrix(a,*n,*n,*lda).selfadjointView<Upper>() * (alpha * vector(actual_x,*n));
     45   else if(UPLO(*uplo)==LO) vector(actual_y,*n).noalias() += matrix(a,*n,*n,*lda).selfadjointView<Lower>() * (alpha * vector(actual_x,*n));
     46 
     47   if(actual_x!=x) delete[] actual_x;
     48   if(actual_y!=y) delete[] copy_back(actual_y,y,*n,*incy);
     49 
     50   return 1;
     51 }
     52 
     53 // C := alpha*x*x' + C
     54 int EIGEN_BLAS_FUNC(syr)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *pc, int *ldc)
     55 {
     56 
     57 //   typedef void (*functype)(int, const Scalar *, int, Scalar *, int, Scalar);
     58 //   static functype func[2];
     59 
     60 //   static bool init = false;
     61 //   if(!init)
     62 //   {
     63 //     for(int k=0; k<2; ++k)
     64 //       func[k] = 0;
     65 //
     66 //     func[UP] = (internal::selfadjoint_product<Scalar,ColMajor,ColMajor,false,UpperTriangular>::run);
     67 //     func[LO] = (internal::selfadjoint_product<Scalar,ColMajor,ColMajor,false,LowerTriangular>::run);
     68 
     69 //     init = true;
     70 //   }
     71 
     72   Scalar* x = reinterpret_cast<Scalar*>(px);
     73   Scalar* c = reinterpret_cast<Scalar*>(pc);
     74   Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
     75 
     76   int info = 0;
     77   if(UPLO(*uplo)==INVALID)                                            info = 1;
     78   else if(*n<0)                                                       info = 2;
     79   else if(*incx==0)                                                   info = 5;
     80   else if(*ldc<std::max(1,*n))                                        info = 7;
     81   if(info)
     82     return xerbla_(SCALAR_SUFFIX_UP"SYR  ",&info,6);
     83 
     84   if(*n==0 || alpha==Scalar(0)) return 1;
     85 
     86   // if the increment is not 1, let's copy it to a temporary vector to enable vectorization
     87   Scalar* x_cpy = get_compact_vector(x,*n,*incx);
     88 
     89   Matrix<Scalar,Dynamic,Dynamic> m2(matrix(c,*n,*n,*ldc));
     90 
     91   // TODO check why this is not accurate enough for lapack tests
     92 //   if(UPLO(*uplo)==LO)       matrix(c,*n,*n,*ldc).selfadjointView<Lower>().rankUpdate(vector(x_cpy,*n), alpha);
     93 //   else if(UPLO(*uplo)==UP)  matrix(c,*n,*n,*ldc).selfadjointView<Upper>().rankUpdate(vector(x_cpy,*n), alpha);
     94 
     95   if(UPLO(*uplo)==LO)
     96     for(int j=0;j<*n;++j)
     97       matrix(c,*n,*n,*ldc).col(j).tail(*n-j) += alpha * x_cpy[j] * vector(x_cpy+j,*n-j);
     98   else
     99     for(int j=0;j<*n;++j)
    100       matrix(c,*n,*n,*ldc).col(j).head(j+1) += alpha * x_cpy[j] * vector(x_cpy,j+1);
    101 
    102   if(x_cpy!=x)  delete[] x_cpy;
    103 
    104   return 1;
    105 }
    106 
    107 // C := alpha*x*y' + alpha*y*x' + C
    108 int EIGEN_BLAS_FUNC(syr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pc, int *ldc)
    109 {
    110 //   typedef void (*functype)(int, const Scalar *, int, const Scalar *, int, Scalar *, int, Scalar);
    111 //   static functype func[2];
    112 //
    113 //   static bool init = false;
    114 //   if(!init)
    115 //   {
    116 //     for(int k=0; k<2; ++k)
    117 //       func[k] = 0;
    118 //
    119 //     func[UP] = (internal::selfadjoint_product<Scalar,ColMajor,ColMajor,false,UpperTriangular>::run);
    120 //     func[LO] = (internal::selfadjoint_product<Scalar,ColMajor,ColMajor,false,LowerTriangular>::run);
    121 //
    122 //     init = true;
    123 //   }
    124 
    125   Scalar* x = reinterpret_cast<Scalar*>(px);
    126   Scalar* y = reinterpret_cast<Scalar*>(py);
    127   Scalar* c = reinterpret_cast<Scalar*>(pc);
    128   Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
    129 
    130   int info = 0;
    131   if(UPLO(*uplo)==INVALID)                                            info = 1;
    132   else if(*n<0)                                                       info = 2;
    133   else if(*incx==0)                                                   info = 5;
    134   else if(*incy==0)                                                   info = 7;
    135   else if(*ldc<std::max(1,*n))                                        info = 9;
    136   if(info)
    137     return xerbla_(SCALAR_SUFFIX_UP"SYR2 ",&info,6);
    138 
    139   if(alpha==Scalar(0))
    140     return 1;
    141 
    142   Scalar* x_cpy = get_compact_vector(x,*n,*incx);
    143   Scalar* y_cpy = get_compact_vector(y,*n,*incy);
    144 
    145   // TODO perform direct calls to underlying implementation
    146   if(UPLO(*uplo)==LO)       matrix(c,*n,*n,*ldc).selfadjointView<Lower>().rankUpdate(vector(x_cpy,*n), vector(y_cpy,*n), alpha);
    147   else if(UPLO(*uplo)==UP)  matrix(c,*n,*n,*ldc).selfadjointView<Upper>().rankUpdate(vector(x_cpy,*n), vector(y_cpy,*n), alpha);
    148 
    149   if(x_cpy!=x)  delete[] x_cpy;
    150   if(y_cpy!=y)  delete[] y_cpy;
    151 
    152 //   int code = UPLO(*uplo);
    153 //   if(code>=2 || func[code]==0)
    154 //     return 0;
    155 
    156 //   func[code](*n, a, *inca, b, *incb, c, *ldc, alpha);
    157   return 1;
    158 }
    159 
    160 /**  DSBMV  performs the matrix-vector  operation
    161   *
    162   *     y := alpha*A*x + beta*y,
    163   *
    164   *  where alpha and beta are scalars, x and y are n element vectors and
    165   *  A is an n by n symmetric band matrix, with k super-diagonals.
    166   */
    167 // int EIGEN_BLAS_FUNC(sbmv)( char *uplo, int *n, int *k, RealScalar *alpha, RealScalar *a, int *lda,
    168 //                            RealScalar *x, int *incx, RealScalar *beta, RealScalar *y, int *incy)
    169 // {
    170 //   return 1;
    171 // }
    172 
    173 
    174 /**  DSPMV  performs the matrix-vector operation
    175   *
    176   *     y := alpha*A*x + beta*y,
    177   *
    178   *  where alpha and beta are scalars, x and y are n element vectors and
    179   *  A is an n by n symmetric matrix, supplied in packed form.
    180   *
    181   */
    182 // int EIGEN_BLAS_FUNC(spmv)(char *uplo, int *n, RealScalar *alpha, RealScalar *ap, RealScalar *x, int *incx, RealScalar *beta, RealScalar *y, int *incy)
    183 // {
    184 //   return 1;
    185 // }
    186 
    187 /**  DSPR    performs the symmetric rank 1 operation
    188   *
    189   *     A := alpha*x*x' + A,
    190   *
    191   *  where alpha is a real scalar, x is an n element vector and A is an
    192   *  n by n symmetric matrix, supplied in packed form.
    193   */
    194 // int EIGEN_BLAS_FUNC(spr)(char *uplo, int *n, Scalar *alpha, Scalar *x, int *incx, Scalar *ap)
    195 // {
    196 //   return 1;
    197 // }
    198 
    199 /**  DSPR2  performs the symmetric rank 2 operation
    200   *
    201   *     A := alpha*x*y' + alpha*y*x' + A,
    202   *
    203   *  where alpha is a scalar, x and y are n element vectors and A is an
    204   *  n by n symmetric matrix, supplied in packed form.
    205   */
    206 // int EIGEN_BLAS_FUNC(spr2)(char *uplo, int *n, RealScalar *alpha, RealScalar *x, int *incx, RealScalar *y, int *incy, RealScalar *ap)
    207 // {
    208 //   return 1;
    209 // }
    210 
    211