Attachment #648 for bug #1152

This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen




  *     y := alpha*A*x + beta*y,
  *
  *  where alpha and beta are scalars, x and y are n element vectors and
  *  A is an n by n hermitian matrix.
  */
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)
{
  typedef void (*functype)(int, const Scalar*, int, const Scalar*, Scalar*, Scalar);
  static const functype func[2] = {
    // array index: UP
    (internal::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Upper,false,false>::run),
    // array index: LO
    (internal::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Lower,false,false>::run),
  };
      func[k] = 0;

    func[UP] = (internal::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Upper,false,false>::run);
    func[LO] = (internal::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Lower,false,false>::run);

    init = true;
  }

  Scalar* a = reinterpret_cast<Scalar*>(pa);
  Scalar* x = reinterpret_cast<Scalar*>(px);
  Scalar* y = reinterpret_cast<Scalar*>(py);
  Scalar alpha  = *reinterpret_cast<Scalar*>(palpha);
  Scalar beta   = *reinterpret_cast<Scalar*>(pbeta);

  // check arguments

  *     A := alpha*x*conjg( x' ) + A,
  *
  *  where alpha is a real scalar, x is an n element vector and A is an
  *  n by n hermitian matrix, supplied in packed form.
  */
int EIGEN_BLAS_FUNC(hpr)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *pap)
{
  typedef void (*functype)(int, Scalar*, const Scalar*, RealScalar);
  static const functype func[2] = {
    // array index: UP
    (internal::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Upper,false,Conj>::run),
    // array index: LO
    (internal::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Lower,false,Conj>::run),
  };
      func[k] = 0;

    func[UP] = (internal::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Upper,false,Conj>::run);
    func[LO] = (internal::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Lower,false,Conj>::run);

    init = true;
  }

  Scalar* x = reinterpret_cast<Scalar*>(px);
  Scalar* ap = reinterpret_cast<Scalar*>(pap);
  RealScalar alpha = *palpha;

  int info = 0;
  if(UPLO(*uplo)==INVALID)                                            info = 1;
  else if(*n<0)                                                       info = 2;

  *     A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A,
  *
  *  where alpha is a scalar, x and y are n element vectors and A is an
  *  n by n hermitian matrix, supplied in packed form.
  */
int EIGEN_BLAS_FUNC(hpr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pap)
{
  typedef void (*functype)(int, Scalar*, const Scalar*, const Scalar*, Scalar);
  static const functype func[2] = {
    // array index: UP
    (internal::packed_rank2_update_selector<Scalar,int,Upper>::run),
    // array index: LO
    (internal::packed_rank2_update_selector<Scalar,int,Lower>::run),
  };
      func[k] = 0;

    func[UP] = (internal::packed_rank2_update_selector<Scalar,int,Upper>::run);
    func[LO] = (internal::packed_rank2_update_selector<Scalar,int,Lower>::run);

    init = true;
  }

  Scalar* x = reinterpret_cast<Scalar*>(px);
  Scalar* y = reinterpret_cast<Scalar*>(py);
  Scalar* ap = reinterpret_cast<Scalar*>(pap);
  Scalar alpha = *reinterpret_cast<Scalar*>(palpha);

  int info = 0;
  if(UPLO(*uplo)==INVALID)                                            info = 1;

  *     A := alpha*x*conjg( x' ) + A,
  *
  *  where alpha is a real scalar, x is an n element vector and A is an
  *  n by n hermitian matrix.
  */
int EIGEN_BLAS_FUNC(her)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *pa, int *lda)
{
  typedef void (*functype)(int, Scalar*, int, const Scalar*, const Scalar*, const Scalar&);
  static const functype func[2] = {
    // array index: UP
    (selfadjoint_rank1_update<Scalar,int,ColMajor,Upper,false,Conj>::run),
    // array index: LO
    (selfadjoint_rank1_update<Scalar,int,ColMajor,Lower,false,Conj>::run),
  };
      func[k] = 0;

    func[UP] = (selfadjoint_rank1_update<Scalar,int,ColMajor,Upper,false,Conj>::run);
    func[LO] = (selfadjoint_rank1_update<Scalar,int,ColMajor,Lower,false,Conj>::run);

    init = true;
  }

  Scalar* x = reinterpret_cast<Scalar*>(px);
  Scalar* a = reinterpret_cast<Scalar*>(pa);
  RealScalar alpha = *reinterpret_cast<RealScalar*>(palpha);

  int info = 0;
  if(UPLO(*uplo)==INVALID)                                            info = 1;
  else if(*n<0)                                                       info = 2;

  *     A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A,
  *
  *  where alpha is a scalar, x and y are n element vectors and A is an n
  *  by n hermitian matrix.
  */
int EIGEN_BLAS_FUNC(her2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pa, int *lda)
{
  typedef void (*functype)(int, Scalar*, int, const Scalar*, const Scalar*, Scalar);
  static const functype func[2] = {
    // array index: UP
    (internal::rank2_update_selector<Scalar,int,Upper>::run),
    // array index: LO
    (internal::rank2_update_selector<Scalar,int,Lower>::run),
  };
      func[k] = 0;

    func[UP] = (internal::rank2_update_selector<Scalar,int,Upper>::run);
    func[LO] = (internal::rank2_update_selector<Scalar,int,Lower>::run);

    init = true;
  }

  Scalar* x = reinterpret_cast<Scalar*>(px);
  Scalar* y = reinterpret_cast<Scalar*>(py);
  Scalar* a = reinterpret_cast<Scalar*>(pa);
  Scalar alpha = *reinterpret_cast<Scalar*>(palpha);

  int info = 0;
  if(UPLO(*uplo)==INVALID)                                            info = 1;




        <Index,Scalar,LhsMapper,StorageOrder,ConjugateLhs,Scalar,RhsMapper,ConjugateRhs>::run(
        rows, cols, LhsMapper(lhs, lhsStride), RhsMapper(rhs, rhsIncr), res, resIncr, alpha);
  }
};

int EIGEN_BLAS_FUNC(gemv)(char *opa, int *m, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *incb, RealScalar *pbeta, RealScalar *pc, int *incc)
{
  typedef void (*functype)(int, int, const Scalar *, int, const Scalar *, int , Scalar *, int, Scalar);
  static const functype func[4] = {
    // array index: NOTR
    (general_matrix_vector_product_wrapper<int,Scalar,ColMajor,false,false>::run),
    // array index: TR  
    (general_matrix_vector_product_wrapper<int,Scalar,RowMajor,false,false>::run),
    // array index: ADJ 
    (general_matrix_vector_product_wrapper<int,Scalar,RowMajor,Conj ,false>::run),
    0
  };
    func[TR  ] = (general_matrix_vector_product_wrapper<int,Scalar,RowMajor,false,false>::run);
    func[ADJ ] = (general_matrix_vector_product_wrapper<int,Scalar,RowMajor,Conj ,false>::run);

    init = true;
  }

  Scalar* a = reinterpret_cast<Scalar*>(pa);
  Scalar* b = reinterpret_cast<Scalar*>(pb);
  Scalar* c = reinterpret_cast<Scalar*>(pc);
  Scalar alpha  = *reinterpret_cast<Scalar*>(palpha);
  Scalar beta   = *reinterpret_cast<Scalar*>(pbeta);

  // check arguments

  if(actual_c!=c) delete[] copy_back(actual_c,c,actual_m,*incc);

  return 1;
}

int EIGEN_BLAS_FUNC(trsv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pa, int *lda, RealScalar *pb, int *incb)
{
  typedef void (*functype)(int, const Scalar *, int, Scalar *);
  static const functype func[16] = {
    // array index: NOTR  | (UP << 2) | (NUNIT << 3)
    (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0,       false,ColMajor>::run),
    // array index: TR    | (UP << 2) | (NUNIT << 3)
    (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0,       false,RowMajor>::run),
    // array index: ADJ   | (UP << 2) | (NUNIT << 3)
    (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0,       Conj, RowMajor>::run),
    0,
    // array index: NOTR  | (LO << 2) | (NUNIT << 3)
    (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0,       false,ColMajor>::run),
    // array index: TR    | (LO << 2) | (NUNIT << 3)
    (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0,       false,RowMajor>::run),
    // array index: ADJ   | (LO << 2) | (NUNIT << 3)
    (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0,       Conj, RowMajor>::run),
    0,
    // array index: NOTR  | (UP << 2) | (UNIT  << 3)
    (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,ColMajor>::run),
    // array index: TR    | (UP << 2) | (UNIT  << 3)
    (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,RowMajor>::run),
    // array index: ADJ   | (UP << 2) | (UNIT  << 3)
    (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,Conj, RowMajor>::run),
    0,
    // array index: NOTR  | (LO << 2) | (UNIT  << 3)
    (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,ColMajor>::run),
    // array index: TR    | (LO << 2) | (UNIT  << 3)
    (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,RowMajor>::run),
    // array index: ADJ   | (LO << 2) | (UNIT  << 3)
    (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,Conj, RowMajor>::run),
    0
  };

  Scalar* a = reinterpret_cast<Scalar*>(pa);
  Scalar* b = reinterpret_cast<Scalar*>(pb);

  int info = 0;
  if(UPLO(*uplo)==INVALID)                                            info = 1;
  else if(OP(*opa)==INVALID)                                          info = 2;
  else if(DIAG(*diag)==INVALID)                                       info = 3;

  return 0;
}



int EIGEN_BLAS_FUNC(trmv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pa, int *lda, RealScalar *pb, int *incb)
{
  typedef void (*functype)(int, int, const Scalar *, int, const Scalar *, int, Scalar *, int, const Scalar&);
  static const functype func[16] = {
    // array index: NOTR  | (UP << 2) | (NUNIT << 3)
    (internal::triangular_matrix_vector_product<int,Upper|0,       Scalar,false,Scalar,false,ColMajor>::run),
    // array index: TR    | (UP << 2) | (NUNIT << 3)
    (internal::triangular_matrix_vector_product<int,Lower|0,       Scalar,false,Scalar,false,RowMajor>::run),
    // array index: ADJ   | (UP << 2) | (NUNIT << 3)
    (internal::triangular_matrix_vector_product<int,Lower|0,       Scalar,Conj, Scalar,false,RowMajor>::run),
    0,
    // array index: NOTR  | (LO << 2) | (NUNIT << 3)
    (internal::triangular_matrix_vector_product<int,Lower|0,       Scalar,false,Scalar,false,ColMajor>::run),
    // array index: TR    | (LO << 2) | (NUNIT << 3)
    (internal::triangular_matrix_vector_product<int,Upper|0,       Scalar,false,Scalar,false,RowMajor>::run),
    // array index: ADJ   | (LO << 2) | (NUNIT << 3)
    (internal::triangular_matrix_vector_product<int,Upper|0,       Scalar,Conj, Scalar,false,RowMajor>::run),
    0,
    // array index: NOTR  | (UP << 2) | (UNIT  << 3)
    (internal::triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,false,Scalar,false,ColMajor>::run),
    // array index: TR    | (UP << 2) | (UNIT  << 3)
    (internal::triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,false,Scalar,false,RowMajor>::run),
    // array index: ADJ   | (UP << 2) | (UNIT  << 3)
    (internal::triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,Conj, Scalar,false,RowMajor>::run),
    0,
    // array index: NOTR  | (LO << 2) | (UNIT  << 3)
    (internal::triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,false,Scalar,false,ColMajor>::run),
    // array index: TR    | (LO << 2) | (UNIT  << 3)
    (internal::triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,false,Scalar,false,RowMajor>::run),
    // array index: ADJ   | (LO << 2) | (UNIT  << 3)
    (internal::triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,Conj, Scalar,false,RowMajor>::run),
    0
  };

  Scalar* a = reinterpret_cast<Scalar*>(pa);
  Scalar* b = reinterpret_cast<Scalar*>(pb);

  int info = 0;
  if(UPLO(*uplo)==INVALID)                                            info = 1;
  else if(OP(*opa)==INVALID)                                          info = 2;
  else if(DIAG(*diag)==INVALID)                                       info = 3;

  *  diagonals.
  *
  *  No test for singularity or near-singularity is included in this
  *  routine. Such tests must be performed before calling this routine.
  */
int EIGEN_BLAS_FUNC(tbsv)(char *uplo, char *op, char *diag, int *n, int *k, RealScalar *pa, int *lda, RealScalar *px, int *incx)
{
  typedef void (*functype)(int, int, const Scalar *, int, Scalar *);
  static const functype func[16] = {
    // array index: NOTR  | (UP << 2) | (NUNIT << 3)
    (internal::band_solve_triangular_selector<int,Upper|0,       Scalar,false,Scalar,ColMajor>::run),
    // array index: TR    | (UP << 2) | (NUNIT << 3)
    (internal::band_solve_triangular_selector<int,Lower|0,       Scalar,false,Scalar,RowMajor>::run),
    // array index: ADJ   | (UP << 2) | (NUNIT << 3)
    (internal::band_solve_triangular_selector<int,Lower|0,       Scalar,Conj, Scalar,RowMajor>::run),
    0,
    // array index: NOTR  | (LO << 2) | (NUNIT << 3)
    (internal::band_solve_triangular_selector<int,Lower|0,       Scalar,false,Scalar,ColMajor>::run),
    // array index: TR    | (LO << 2) | (NUNIT << 3)
    (internal::band_solve_triangular_selector<int,Upper|0,       Scalar,false,Scalar,RowMajor>::run),
    // array index: ADJ   | (LO << 2) | (NUNIT << 3)
    (internal::band_solve_triangular_selector<int,Upper|0,       Scalar,Conj, Scalar,RowMajor>::run),
    0,
    // array index: NOTR  | (UP << 2) | (UNIT  << 3)
    (internal::band_solve_triangular_selector<int,Upper|UnitDiag,Scalar,false,Scalar,ColMajor>::run),
    // array index: TR    | (UP << 2) | (UNIT  << 3)
    (internal::band_solve_triangular_selector<int,Lower|UnitDiag,Scalar,false,Scalar,RowMajor>::run),
    // array index: ADJ   | (UP << 2) | (UNIT  << 3)
    (internal::band_solve_triangular_selector<int,Lower|UnitDiag,Scalar,Conj, Scalar,RowMajor>::run),
    0,
    // array index: NOTR  | (LO << 2) | (UNIT  << 3)
    (internal::band_solve_triangular_selector<int,Lower|UnitDiag,Scalar,false,Scalar,ColMajor>::run),
    // array index: TR    | (LO << 2) | (UNIT  << 3)
    (internal::band_solve_triangular_selector<int,Upper|UnitDiag,Scalar,false,Scalar,RowMajor>::run),
    // array index: ADJ   | (LO << 2) | (UNIT  << 3)
    (internal::band_solve_triangular_selector<int,Upper|UnitDiag,Scalar,Conj, Scalar,RowMajor>::run),
    0,
  };

  Scalar* a = reinterpret_cast<Scalar*>(pa);
  Scalar* x = reinterpret_cast<Scalar*>(px);
  int coeff_rows = *k+1;

  int info = 0;
       if(UPLO(*uplo)==INVALID)                                       info = 1;
  else if(OP(*op)==INVALID)                                           info = 2;

  *     x := A*x,   or   x := A'*x,
  *
  *  where x is an n element vector and  A is an n by n unit, or non-unit,
  *  upper or lower triangular matrix, supplied in packed form.
  */
int EIGEN_BLAS_FUNC(tpmv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pap, RealScalar *px, int *incx)
{
  typedef void (*functype)(int, const Scalar*, const Scalar*, Scalar*, Scalar);
  static const functype func[16] = {
    // array index: NOTR  | (UP << 2) | (NUNIT << 3)
    (internal::packed_triangular_matrix_vector_product<int,Upper|0,       Scalar,false,Scalar,false,ColMajor>::run),
    // array index: TR    | (UP << 2) | (NUNIT << 3)
    (internal::packed_triangular_matrix_vector_product<int,Lower|0,       Scalar,false,Scalar,false,RowMajor>::run),
    // array index: ADJ   | (UP << 2) | (NUNIT << 3)
    (internal::packed_triangular_matrix_vector_product<int,Lower|0,       Scalar,Conj, Scalar,false,RowMajor>::run),
    0,
    // array index: NOTR  | (LO << 2) | (NUNIT << 3)
    (internal::packed_triangular_matrix_vector_product<int,Lower|0,       Scalar,false,Scalar,false,ColMajor>::run),
    // array index: TR    | (LO << 2) | (NUNIT << 3)
    (internal::packed_triangular_matrix_vector_product<int,Upper|0,       Scalar,false,Scalar,false,RowMajor>::run),
    // array index: ADJ   | (LO << 2) | (NUNIT << 3)
    (internal::packed_triangular_matrix_vector_product<int,Upper|0,       Scalar,Conj, Scalar,false,RowMajor>::run),
    0,
    // array index: NOTR  | (UP << 2) | (UNIT  << 3)
    (internal::packed_triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,false,Scalar,false,ColMajor>::run),
    // array index: TR    | (UP << 2) | (UNIT  << 3)
    (internal::packed_triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,false,Scalar,false,RowMajor>::run),
    // array index: ADJ   | (UP << 2) | (UNIT  << 3)
    (internal::packed_triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,Conj, Scalar,false,RowMajor>::run),
    0,
    // array index: NOTR  | (LO << 2) | (UNIT  << 3)
    (internal::packed_triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,false,Scalar,false,ColMajor>::run),
    // array index: TR    | (LO << 2) | (UNIT  << 3)
    (internal::packed_triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,false,Scalar,false,RowMajor>::run),
    // array index: ADJ   | (LO << 2) | (UNIT  << 3)
    (internal::packed_triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,Conj, Scalar,false,RowMajor>::run),
    0
  };

  Scalar* ap = reinterpret_cast<Scalar*>(pap);
  Scalar* x = reinterpret_cast<Scalar*>(px);

  int info = 0;
  if(UPLO(*uplo)==INVALID)                                            info = 1;
  else if(OP(*opa)==INVALID)                                          info = 2;
  else if(DIAG(*diag)==INVALID)                                       info = 3;

  *  non-unit, upper or lower triangular matrix, supplied in packed form.
  *
  *  No test for singularity or near-singularity is included in this
  *  routine. Such tests must be performed before calling this routine.
  */
int EIGEN_BLAS_FUNC(tpsv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pap, RealScalar *px, int *incx)
{
  typedef void (*functype)(int, const Scalar*, Scalar*);
  static const functype func[16] = {
    // array index: NOTR  | (UP << 2) | (NUNIT << 3)
    (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0,       false,ColMajor>::run),
    // array index: TR    | (UP << 2) | (NUNIT << 3)
    (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0,       false,RowMajor>::run),
    // array index: ADJ   | (UP << 2) | (NUNIT << 3)
    (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0,       Conj, RowMajor>::run),
    0,
    // array index: NOTR  | (LO << 2) | (NUNIT << 3)
    (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0,       false,ColMajor>::run),
    // array index: TR    | (LO << 2) | (NUNIT << 3)
    (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0,       false,RowMajor>::run),
    // array index: ADJ   | (LO << 2) | (NUNIT << 3)
    (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0,       Conj, RowMajor>::run),
    0,
    // array index: NOTR  | (UP << 2) | (UNIT  << 3)
    (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,ColMajor>::run),
    // array index: TR    | (UP << 2) | (UNIT  << 3)
    (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,RowMajor>::run),
    // array index: ADJ   | (UP << 2) | (UNIT  << 3)
    (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,Conj, RowMajor>::run),
    0,
    // array index: NOTR  | (LO << 2) | (UNIT  << 3)
    (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,ColMajor>::run),
    // array index: TR    | (LO << 2) | (UNIT  << 3)
    (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,RowMajor>::run),
    // array index: ADJ   | (LO << 2) | (UNIT  << 3)
    (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,Conj, RowMajor>::run),
    0
  };

  Scalar* ap = reinterpret_cast<Scalar*>(pap);
  Scalar* x = reinterpret_cast<Scalar*>(px);

  int info = 0;
  if(UPLO(*uplo)==INVALID)                                            info = 1;
  else if(OP(*opa)==INVALID)                                          info = 2;
  else if(DIAG(*diag)==INVALID)                                       info = 3;




// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#include "common.h"

// y = alpha*A*x + beta*y
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)
{
  typedef void (*functype)(int, const Scalar*, int, const Scalar*, Scalar*, Scalar);
  static const functype func[2] = {
    // array index: UP
    (internal::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Upper,false,false>::run),
    // array index: LO
    (internal::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Lower,false,false>::run),
  };
      func[k] = 0;

    func[UP] = (internal::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Upper,false,false>::run);
    func[LO] = (internal::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Lower,false,false>::run);

    init = true;
  }

  Scalar* a = reinterpret_cast<Scalar*>(pa);
  Scalar* x = reinterpret_cast<Scalar*>(px);
  Scalar* y = reinterpret_cast<Scalar*>(py);
  Scalar alpha  = *reinterpret_cast<Scalar*>(palpha);
  Scalar beta   = *reinterpret_cast<Scalar*>(pbeta);

  // check arguments


  return 1;
}

// C := alpha*x*x' + C
int EIGEN_BLAS_FUNC(syr)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *pc, int *ldc)
{

//   typedef void (*functype)(int, const Scalar *, int, Scalar *, int, Scalar);
//   static functype func[2];

//   static bool init = false;
//   if(!init)
//   {
//     for(int k=0; k<2; ++k)
//       func[k] = 0;
//
//     func[UP] = (internal::selfadjoint_product<Scalar,ColMajor,ColMajor,false,UpperTriangular>::run);
//     func[LO] = (internal::selfadjoint_product<Scalar,ColMajor,ColMajor,false,LowerTriangular>::run);

//     init = true;
//   }
  typedef void (*functype)(int, Scalar*, int, const Scalar*, const Scalar*, const Scalar&);
  static const functype func[2] = {
    // array index: UP
    (selfadjoint_rank1_update<Scalar,int,ColMajor,Upper,false,Conj>::run),
    // array index: LO
    (selfadjoint_rank1_update<Scalar,int,ColMajor,Lower,false,Conj>::run),
  };
      func[k] = 0;

    func[UP] = (selfadjoint_rank1_update<Scalar,int,ColMajor,Upper,false,Conj>::run);
    func[LO] = (selfadjoint_rank1_update<Scalar,int,ColMajor,Lower,false,Conj>::run);

    init = true;
  }

  Scalar* x = reinterpret_cast<Scalar*>(px);
  Scalar* c = reinterpret_cast<Scalar*>(pc);
  Scalar alpha = *reinterpret_cast<Scalar*>(palpha);

  int info = 0;
  if(UPLO(*uplo)==INVALID)                                            info = 1;
  else if(*n<0)                                                       info = 2;

  if(x_cpy!=x)  delete[] x_cpy;

  return 1;
}

// C := alpha*x*y' + alpha*y*x' + C
int EIGEN_BLAS_FUNC(syr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pc, int *ldc)
{
//   typedef void (*functype)(int, const Scalar *, int, const Scalar *, int, Scalar *, int, Scalar);
//   static functype func[2];
//
//   static bool init = false;
//   if(!init)
//   {
//     for(int k=0; k<2; ++k)
//       func[k] = 0;
//
//     func[UP] = (internal::selfadjoint_product<Scalar,ColMajor,ColMajor,false,UpperTriangular>::run);
//     func[LO] = (internal::selfadjoint_product<Scalar,ColMajor,ColMajor,false,LowerTriangular>::run);
//
//     init = true;
//   }
  typedef void (*functype)(int, Scalar*, int, const Scalar*, const Scalar*, Scalar);
  static const functype func[2] = {
    // array index: UP
    (internal::rank2_update_selector<Scalar,int,Upper>::run),
    // array index: LO
    (internal::rank2_update_selector<Scalar,int,Lower>::run),
  };
      func[k] = 0;

    func[UP] = (internal::rank2_update_selector<Scalar,int,Upper>::run);
    func[LO] = (internal::rank2_update_selector<Scalar,int,Lower>::run);

    init = true;
  }

  Scalar* x = reinterpret_cast<Scalar*>(px);
  Scalar* y = reinterpret_cast<Scalar*>(py);
  Scalar* c = reinterpret_cast<Scalar*>(pc);
  Scalar alpha = *reinterpret_cast<Scalar*>(palpha);

  int info = 0;
  if(UPLO(*uplo)==INVALID)                                            info = 1;

  *     A := alpha*x*x' + A,
  *
  *  where alpha is a real scalar, x is an n element vector and A is an
  *  n by n symmetric matrix, supplied in packed form.
  */
int EIGEN_BLAS_FUNC(spr)(char *uplo, int *n, Scalar *palpha, Scalar *px, int *incx, Scalar *pap)
{
  typedef void (*functype)(int, Scalar*, const Scalar*, Scalar);
  static const functype func[2] = {
    // array index: UP
    (internal::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Upper,false,false>::run),
    // array index: LO
    (internal::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Lower,false,false>::run),
  };
      func[k] = 0;

    func[UP] = (internal::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Upper,false,false>::run);
    func[LO] = (internal::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Lower,false,false>::run);

    init = true;
  }

  Scalar* x = reinterpret_cast<Scalar*>(px);
  Scalar* ap = reinterpret_cast<Scalar*>(pap);
  Scalar alpha = *reinterpret_cast<Scalar*>(palpha);

  int info = 0;
  if(UPLO(*uplo)==INVALID)                                            info = 1;
  else if(*n<0)                                                       info = 2;

  *     A := alpha*x*y' + alpha*y*x' + A,
  *
  *  where alpha is a scalar, x and y are n element vectors and A is an
  *  n by n symmetric matrix, supplied in packed form.
  */
int EIGEN_BLAS_FUNC(spr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pap)
{
  typedef void (*functype)(int, Scalar*, const Scalar*, const Scalar*, Scalar);
  static const functype func[2] = {
    // array index: UP
    (internal::packed_rank2_update_selector<Scalar,int,Upper>::run),
    // array index: LO
    (internal::packed_rank2_update_selector<Scalar,int,Lower>::run),
  };
      func[k] = 0;

    func[UP] = (internal::packed_rank2_update_selector<Scalar,int,Upper>::run);
    func[LO] = (internal::packed_rank2_update_selector<Scalar,int,Lower>::run);

    init = true;
  }

  Scalar* x = reinterpret_cast<Scalar*>(px);
  Scalar* y = reinterpret_cast<Scalar*>(py);
  Scalar* ap = reinterpret_cast<Scalar*>(pap);
  Scalar alpha = *reinterpret_cast<Scalar*>(palpha);

  int info = 0;
  if(UPLO(*uplo)==INVALID)                                            info = 1;




// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#include <iostream>
#include "common.h"

int EIGEN_BLAS_FUNC(gemm)(char *opa, char *opb, int *m, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, RealScalar *pbeta, RealScalar *pc, int *ldc)
{
//   std::cerr << "in gemm " << *opa << " " << *opb << " " << *m << " " << *n << " " << *k << " " << *lda << " " << *ldb << " " << *ldc << " " << *palpha << " " << *pbeta << "\n";
  typedef void (*functype)(DenseIndex, DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, Scalar, internal::level3_blocking<Scalar,Scalar>&, Eigen::internal::GemmParallelInfo<DenseIndex>*);
  static const functype func[12] = {
    // array index: NOTR  | (NOTR << 2)
    (internal::general_matrix_matrix_product<DenseIndex,Scalar,ColMajor,false,Scalar,ColMajor,false,ColMajor>::run),
    // array index: TR    | (NOTR << 2)
    (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,false,Scalar,ColMajor,false,ColMajor>::run),
    // array index: ADJ   | (NOTR << 2)
    (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,false,ColMajor>::run),
    0,
    // array index: NOTR  | (TR   << 2)
    (internal::general_matrix_matrix_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,false,ColMajor>::run),
    // array index: TR    | (TR   << 2)
    (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,false,Scalar,RowMajor,false,ColMajor>::run),
    // array index: ADJ   | (TR   << 2)
    (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,RowMajor,false,ColMajor>::run),
    0,
    // array index: NOTR  | (ADJ  << 2)
    (internal::general_matrix_matrix_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,Conj, ColMajor>::run),
    // array index: TR    | (ADJ  << 2)
    (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,false,Scalar,RowMajor,Conj, ColMajor>::run),
    // array index: ADJ   | (ADJ  << 2)
    (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,RowMajor,Conj, ColMajor>::run),
    0
  };

  Scalar* a = reinterpret_cast<Scalar*>(pa);
  Scalar* b = reinterpret_cast<Scalar*>(pb);
  Scalar* c = reinterpret_cast<Scalar*>(pc);
  Scalar alpha  = *reinterpret_cast<Scalar*>(palpha);
  Scalar beta   = *reinterpret_cast<Scalar*>(pbeta);

  int info = 0;

  func[code](*m, *n, *k, a, *lda, b, *ldb, c, *ldc, alpha, blocking, 0);
  return 0;
}

int EIGEN_BLAS_FUNC(trsm)(char *side, char *uplo, char *opa, char *diag, int *m, int *n, RealScalar *palpha,  RealScalar *pa, int *lda, RealScalar *pb, int *ldb)
{
//   std::cerr << "in trsm " << *side << " " << *uplo << " " << *opa << " " << *diag << " " << *m << "," << *n << " " << *palpha << " " << *lda << " " << *ldb<< "\n";
  typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, internal::level3_blocking<Scalar,Scalar>&);
  static const functype func[32] = {
    // array index: NOTR  | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|0,          false,ColMajor,ColMajor>::run),
    // array index: TR    | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|0,          false,RowMajor,ColMajor>::run),
    // array index: ADJ   | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|0,          Conj, RowMajor,ColMajor>::run),\
    0,
    // array index: NOTR  | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|0,          false,ColMajor,ColMajor>::run),
    // array index: TR    | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|0,          false,RowMajor,ColMajor>::run),
    // array index: ADJ   | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|0,          Conj, RowMajor,ColMajor>::run),
    0,
    // array index: NOTR  | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|0,          false,ColMajor,ColMajor>::run),
    // array index: TR    | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|0,          false,RowMajor,ColMajor>::run),
    // array index: ADJ   | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|0,          Conj, RowMajor,ColMajor>::run),
    0,
    // array index: NOTR  | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|0,          false,ColMajor,ColMajor>::run),
    // array index: TR    | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|0,          false,RowMajor,ColMajor>::run),
    // array index: ADJ   | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|0,          Conj, RowMajor,ColMajor>::run),
    0,
    // array index: NOTR  | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|UnitDiag,false,ColMajor,ColMajor>::run),
    // array index: TR    | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|UnitDiag,false,RowMajor,ColMajor>::run),
    // array index: ADJ   | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|UnitDiag,Conj, RowMajor,ColMajor>::run),
    0,
    // array index: NOTR  | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|UnitDiag,false,ColMajor,ColMajor>::run),
    // array index: TR    | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|UnitDiag,false,RowMajor,ColMajor>::run),
    // array index: ADJ   | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|UnitDiag,Conj, RowMajor,ColMajor>::run),
    0,
    // array index: NOTR  | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|UnitDiag,false,ColMajor,ColMajor>::run),
    // array index: TR    | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|UnitDiag,false,RowMajor,ColMajor>::run),
    // array index: ADJ   | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|UnitDiag,Conj, RowMajor,ColMajor>::run),
    0,
    // array index: NOTR  | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|UnitDiag,false,ColMajor,ColMajor>::run),
    // array index: TR    | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|UnitDiag,false,RowMajor,ColMajor>::run),
    // array index: ADJ   | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)
    (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|UnitDiag,Conj, RowMajor,ColMajor>::run),
    0
  };

  Scalar* a = reinterpret_cast<Scalar*>(pa);
  Scalar* b = reinterpret_cast<Scalar*>(pb);
  Scalar  alpha = *reinterpret_cast<Scalar*>(palpha);

  int info = 0;
  if(SIDE(*side)==INVALID)                                            info = 1;
  else if(UPLO(*uplo)==INVALID)                                       info = 2;



// b = alpha*op(a)*b  for side = 'L'or'l'
// b = alpha*b*op(a)  for side = 'R'or'r'
int EIGEN_BLAS_FUNC(trmm)(char *side, char *uplo, char *opa, char *diag, int *m, int *n, RealScalar *palpha,  RealScalar *pa, int *lda, RealScalar *pb, int *ldb)
{
//   std::cerr << "in trmm " << *side << " " << *uplo << " " << *opa << " " << *diag << " " << *m << " " << *n << " " << *lda << " " << *ldb << " " << *palpha << "\n";
  typedef void (*functype)(DenseIndex, DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, const Scalar&, internal::level3_blocking<Scalar,Scalar>&);
  static const functype func[32] = {
    // array index: NOTR  | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0,          true, ColMajor,false,ColMajor,false,ColMajor>::run),
    // array index: TR    | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0,          true, RowMajor,false,ColMajor,false,ColMajor>::run),
    // array index: ADJ   | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0,          true, RowMajor,Conj, ColMajor,false,ColMajor>::run),
    0,
    // array index: NOTR  | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0,          false,ColMajor,false,ColMajor,false,ColMajor>::run),
    // array index: TR    | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0,          false,ColMajor,false,RowMajor,false,ColMajor>::run),
    // array index: ADJ   | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0,          false,ColMajor,false,RowMajor,Conj, ColMajor>::run),
    0,
    // array index: NOTR  | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0,          true, ColMajor,false,ColMajor,false,ColMajor>::run),
    // array index: TR    | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0,          true, RowMajor,false,ColMajor,false,ColMajor>::run),
    // array index: ADJ   | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0,          true, RowMajor,Conj, ColMajor,false,ColMajor>::run),
    0,
    // array index: NOTR  | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0,          false,ColMajor,false,ColMajor,false,ColMajor>::run),
    // array index: TR    | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0,          false,ColMajor,false,RowMajor,false,ColMajor>::run),
    // array index: ADJ   | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0,          false,ColMajor,false,RowMajor,Conj, ColMajor>::run),
    0,
    // array index: NOTR  | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,true, ColMajor,false,ColMajor,false,ColMajor>::run),
    // array index: TR    | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,true, RowMajor,false,ColMajor,false,ColMajor>::run),
    // array index: ADJ   | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,true, RowMajor,Conj, ColMajor,false,ColMajor>::run),
    0,
    // array index: NOTR  | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,false,ColMajor,false,ColMajor,false,ColMajor>::run),
    // array index: TR    | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,false,ColMajor,false,RowMajor,false,ColMajor>::run),
    // array index: ADJ   | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,false,ColMajor,false,RowMajor,Conj, ColMajor>::run),
    0,
    // array index: NOTR  | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,true, ColMajor,false,ColMajor,false,ColMajor>::run),
    // array index: TR    | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,true, RowMajor,false,ColMajor,false,ColMajor>::run),
    // array index: ADJ   | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,true, RowMajor,Conj, ColMajor,false,ColMajor>::run),
    0,
    // array index: NOTR  | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,false,ColMajor,false,ColMajor,false,ColMajor>::run),
    // array index: TR    | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,false,ColMajor,false,RowMajor,false,ColMajor>::run),
    // array index: ADJ   | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)
    (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,false,ColMajor,false,RowMajor,Conj, ColMajor>::run),
    0
  };

  Scalar* a = reinterpret_cast<Scalar*>(pa);
  Scalar* b = reinterpret_cast<Scalar*>(pb);
  Scalar  alpha = *reinterpret_cast<Scalar*>(palpha);

  int info = 0;
  if(SIDE(*side)==INVALID)                                            info = 1;
  else if(UPLO(*uplo)==INVALID)                                       info = 2;


// c = alpha*a*a' + beta*c  for op = 'N'or'n'
// c = alpha*a'*a + beta*c  for op = 'T'or't','C'or'c'
int EIGEN_BLAS_FUNC(syrk)(char *uplo, char *op, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pbeta, RealScalar *pc, int *ldc)
{
//   std::cerr << "in syrk " << *uplo << " " << *op << " " << *n << " " << *k << " " << *palpha << " " << *lda << " " << *pbeta << " " << *ldc << "\n";
  #if !ISCOMPLEX
  typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, const Scalar&, internal::level3_blocking<Scalar,Scalar>&);
  static const functype func[8] = {
    // array index: NOTR  | (UP << 2)
    (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,ColMajor,Conj, Upper>::run),
    // array index: TR    | (UP << 2)
    (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,false,Scalar,ColMajor,ColMajor,Conj, Upper>::run),
    // array index: ADJ   | (UP << 2)
    (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,ColMajor,false,Upper>::run),
    0,
    // array index: NOTR  | (LO << 2)
    (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,ColMajor,Conj, Lower>::run),
    // array index: TR    | (LO << 2)
    (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,false,Scalar,ColMajor,ColMajor,Conj, Lower>::run),
    // array index: ADJ   | (LO << 2)
    (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,ColMajor,false,Lower>::run),
    0
  };
    init = true;
  }
  #endif

  Scalar* a = reinterpret_cast<Scalar*>(pa);
  Scalar* c = reinterpret_cast<Scalar*>(pc);
  Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
  Scalar beta  = *reinterpret_cast<Scalar*>(pbeta);

  int info = 0;


// c = alpha*a*conj(a') + beta*c  for op = 'N'or'n'
// c = alpha*conj(a')*a + beta*c  for op  = 'C'or'c'
int EIGEN_BLAS_FUNC(herk)(char *uplo, char *op, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pbeta, RealScalar *pc, int *ldc)
{
//   std::cerr << "in herk " << *uplo << " " << *op << " " << *n << " " << *k << " " << *palpha << " " << *lda << " " << *pbeta << " " << *ldc << "\n";

  typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, const Scalar&, internal::level3_blocking<Scalar,Scalar>&);
  static const functype func[8] = {
    // array index: NOTR  | (UP << 2)
    (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,Conj, ColMajor,Upper>::run),
    0,
    // array index: ADJ   | (UP << 2)
    (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,false,ColMajor,Upper>::run),
    0,
    // array index: NOTR  | (LO << 2)
    (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,Conj, ColMajor,Lower>::run),
    0,
    // array index: ADJ   | (LO << 2)
    (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,false,ColMajor,Lower>::run),
    0
  };
    init = true;
  }

  Scalar* a = reinterpret_cast<Scalar*>(pa);
  Scalar* c = reinterpret_cast<Scalar*>(pc);
  RealScalar alpha = *palpha;
  RealScalar beta  = *pbeta;

//   std::cerr << "in herk " << *uplo << " " << *op << " " << *n << " " << *k << " " << alpha << " " << *lda << " " << beta << " " << *ldc << "\n";


Return to bug 1152

Lines 14-42 Link Here

(-)a/blas/level2_cplx_impl.h (-65 / +30 lines)
14	* y := alphaAx + beta*y,	14	* y := alphaAx + beta*y,
15	*	15	*
16	* where alpha and beta are scalars, x and y are n element vectors and	16	* where alpha and beta are scalars, x and y are n element vectors and
17	* A is an n by n hermitian matrix.	17	* A is an n by n hermitian matrix.
18	*/	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)	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	{	20	{
21	typedef void (functype)(int, const Scalar, int, const Scalar, Scalar, Scalar);	21	typedef void (functype)(int, const Scalar, int, const Scalar, Scalar, Scalar);
22	static functype func[2];	22	static const functype func[2] = {
23		23	// array index: UP
24	static bool init = false;	24	(internal::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Upper,false,false>::run),
25	if(!init)	25	// array index: LO
26	{	26	(internal::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Lower,false,false>::run),
27	for(int k=0; k<2; ++k)	27	};
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		28
36	Scalar* a = reinterpret_cast<Scalar*>(pa);	29	Scalar* a = reinterpret_cast<Scalar*>(pa);
37	Scalar* x = reinterpret_cast<Scalar*>(px);	30	Scalar* x = reinterpret_cast<Scalar*>(px);
38	Scalar* y = reinterpret_cast<Scalar*>(py);	31	Scalar* y = reinterpret_cast<Scalar*>(py);
39	Scalar alpha = reinterpret_cast<Scalar>(palpha);	32	Scalar alpha = reinterpret_cast<Scalar>(palpha);
40	Scalar beta = reinterpret_cast<Scalar>(pbeta);	33	Scalar beta = reinterpret_cast<Scalar>(pbeta);
41		34
42	// check arguments	35	// check arguments
Lines 106-134 int EIGEN_BLAS_FUNC(hemv)(char *uplo, in Link Here
106	* A := alphaxconjg( x' ) + A,	99	* A := alphaxconjg( x' ) + A,
107	*	100	*
108	* where alpha is a real scalar, x is an n element vector and A is an	101	* 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.	102	* n by n hermitian matrix, supplied in packed form.
110	*/	103	*/
111	int EIGEN_BLAS_FUNC(hpr)(char uplo, int n, RealScalar palpha, RealScalar px, int incx, RealScalar pap)	104	int EIGEN_BLAS_FUNC(hpr)(char uplo, int n, RealScalar palpha, RealScalar px, int incx, RealScalar pap)
112	{	105	{
113	typedef void (functype)(int, Scalar, const Scalar*, RealScalar);	106	typedef void (functype)(int, Scalar, const Scalar*, RealScalar);
114	static functype func[2];	107	static const functype func[2] = {
115		108	// array index: UP
116	static bool init = false;	109	(internal::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Upper,false,Conj>::run),
117	if(!init)	110	// array index: LO
118	{	111	(internal::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Lower,false,Conj>::run),
119	for(int k=0; k<2; ++k)	112	};
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		113
128	Scalar* x = reinterpret_cast<Scalar*>(px);	114	Scalar* x = reinterpret_cast<Scalar*>(px);
129	Scalar* ap = reinterpret_cast<Scalar*>(pap);	115	Scalar* ap = reinterpret_cast<Scalar*>(pap);
130	RealScalar alpha = *palpha;	116	RealScalar alpha = *palpha;
131		117
132	int info = 0;	118	int info = 0;
133	if(UPLO(*uplo)==INVALID) info = 1;	119	if(UPLO(*uplo)==INVALID) info = 1;
134	else if(*n<0) info = 2;	120	else if(*n<0) info = 2;
Lines 157-185 int EIGEN_BLAS_FUNC(hpr)(char *uplo, int Link Here
157	* A := alphaxconjg( y' ) + conjg( alpha )yconjg( x' ) + A,	143	* A := alphaxconjg( y' ) + conjg( alpha )yconjg( x' ) + A,
158	*	144	*
159	* where alpha is a scalar, x and y are n element vectors and A is an	145	* 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.	146	* n by n hermitian matrix, supplied in packed form.
161	*/	147	*/
162	int EIGEN_BLAS_FUNC(hpr2)(char uplo, int n, RealScalar palpha, RealScalar px, int incx, RealScalar py, int incy, RealScalar pap)	148	int EIGEN_BLAS_FUNC(hpr2)(char uplo, int n, RealScalar palpha, RealScalar px, int incx, RealScalar py, int incy, RealScalar pap)
163	{	149	{
164	typedef void (functype)(int, Scalar, const Scalar, const Scalar, Scalar);	150	typedef void (functype)(int, Scalar, const Scalar, const Scalar, Scalar);
165	static functype func[2];	151	static const functype func[2] = {
166		152	// array index: UP
167	static bool init = false;	153	(internal::packed_rank2_update_selector<Scalar,int,Upper>::run),
168	if(!init)	154	// array index: LO
169	{	155	(internal::packed_rank2_update_selector<Scalar,int,Lower>::run),
170	for(int k=0; k<2; ++k)	156	};
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		157
179	Scalar* x = reinterpret_cast<Scalar*>(px);	158	Scalar* x = reinterpret_cast<Scalar*>(px);
180	Scalar* y = reinterpret_cast<Scalar*>(py);	159	Scalar* y = reinterpret_cast<Scalar*>(py);
181	Scalar* ap = reinterpret_cast<Scalar*>(pap);	160	Scalar* ap = reinterpret_cast<Scalar*>(pap);
182	Scalar alpha = reinterpret_cast<Scalar>(palpha);	161	Scalar alpha = reinterpret_cast<Scalar>(palpha);
183		162
184	int info = 0;	163	int info = 0;
185	if(UPLO(*uplo)==INVALID) info = 1;	164	if(UPLO(*uplo)==INVALID) info = 1;
Lines 212-240 int EIGEN_BLAS_FUNC(hpr2)(char *uplo, in Link Here
212	* A := alphaxconjg( x' ) + A,	191	* A := alphaxconjg( x' ) + A,
213	*	192	*
214	* where alpha is a real scalar, x is an n element vector and A is an	193	* where alpha is a real scalar, x is an n element vector and A is an
215	* n by n hermitian matrix.	194	* n by n hermitian matrix.
216	*/	195	*/
217	int EIGEN_BLAS_FUNC(her)(char uplo, int n, RealScalar palpha, RealScalar px, int incx, RealScalar pa, int *lda)	196	int EIGEN_BLAS_FUNC(her)(char uplo, int n, RealScalar palpha, RealScalar px, int incx, RealScalar pa, int *lda)
218	{	197	{
219	typedef void (functype)(int, Scalar, int, const Scalar, const Scalar, const Scalar&);	198	typedef void (functype)(int, Scalar, int, const Scalar, const Scalar, const Scalar&);
220	static functype func[2];	199	static const functype func[2] = {
221		200	// array index: UP
222	static bool init = false;	201	(selfadjoint_rank1_update<Scalar,int,ColMajor,Upper,false,Conj>::run),
223	if(!init)	202	// array index: LO
224	{	203	(selfadjoint_rank1_update<Scalar,int,ColMajor,Lower,false,Conj>::run),
225	for(int k=0; k<2; ++k)	204	};
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		205
234	Scalar* x = reinterpret_cast<Scalar*>(px);	206	Scalar* x = reinterpret_cast<Scalar*>(px);
235	Scalar* a = reinterpret_cast<Scalar*>(pa);	207	Scalar* a = reinterpret_cast<Scalar*>(pa);
236	RealScalar alpha = reinterpret_cast<RealScalar>(palpha);	208	RealScalar alpha = reinterpret_cast<RealScalar>(palpha);
237		209
238	int info = 0;	210	int info = 0;
239	if(UPLO(*uplo)==INVALID) info = 1;	211	if(UPLO(*uplo)==INVALID) info = 1;
240	else if(*n<0) info = 2;	212	else if(*n<0) info = 2;
Lines 266-294 int EIGEN_BLAS_FUNC(her)(char *uplo, int Link Here
266	* A := alphaxconjg( y' ) + conjg( alpha )yconjg( x' ) + A,	238	* A := alphaxconjg( y' ) + conjg( alpha )yconjg( x' ) + A,
267	*	239	*
268	* where alpha is a scalar, x and y are n element vectors and A is an n	240	* where alpha is a scalar, x and y are n element vectors and A is an n
269	* by n hermitian matrix.	241	* by n hermitian matrix.
270	*/	242	*/
271	int EIGEN_BLAS_FUNC(her2)(char uplo, int n, RealScalar palpha, RealScalar px, int incx, RealScalar py, int incy, RealScalar pa, int *lda)	243	int EIGEN_BLAS_FUNC(her2)(char uplo, int n, RealScalar palpha, RealScalar px, int incx, RealScalar py, int incy, RealScalar pa, int *lda)
272	{	244	{
273	typedef void (functype)(int, Scalar, int, const Scalar, const Scalar, Scalar);	245	typedef void (functype)(int, Scalar, int, const Scalar, const Scalar, Scalar);
274	static functype func[2];	246	static const functype func[2] = {
275		247	// array index: UP
276	static bool init = false;	248	(internal::rank2_update_selector<Scalar,int,Upper>::run),
277	if(!init)	249	// array index: LO
278	{	250	(internal::rank2_update_selector<Scalar,int,Lower>::run),
279	for(int k=0; k<2; ++k)	251	};
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		252
288	Scalar* x = reinterpret_cast<Scalar*>(px);	253	Scalar* x = reinterpret_cast<Scalar*>(px);
289	Scalar* y = reinterpret_cast<Scalar*>(py);	254	Scalar* y = reinterpret_cast<Scalar*>(py);
290	Scalar* a = reinterpret_cast<Scalar*>(pa);	255	Scalar* a = reinterpret_cast<Scalar*>(pa);
291	Scalar alpha = reinterpret_cast<Scalar>(palpha);	256	Scalar alpha = reinterpret_cast<Scalar>(palpha);
292		257
293	int info = 0;	258	int info = 0;
294	if(UPLO(*uplo)==INVALID) info = 1;	259	if(UPLO(*uplo)==INVALID) info = 1;

Lines 8-36 Link Here

(-)a/blas/level2_real_impl.h (-93 / +30 lines)
8	// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.	8	// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9		9
10	#include "common.h"	10	#include "common.h"
11		11
12	// y = alphaAx + beta*y	12	// y = alphaAx + 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)	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	{	14	{
15	typedef void (functype)(int, const Scalar, int, const Scalar, Scalar, Scalar);	15	typedef void (functype)(int, const Scalar, int, const Scalar, Scalar, Scalar);
16	static functype func[2];	16	static const functype func[2] = {
17		17	// array index: UP
18	static bool init = false;	18	(internal::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Upper,false,false>::run),
19	if(!init)	19	// array index: LO
20	{	20	(internal::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Lower,false,false>::run),
21	for(int k=0; k<2; ++k)	21	};
22	func[k] = 0;
23
24	func[UP] = (internal::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Upper,false,false>::run);
25	func[LO] = (internal::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Lower,false,false>::run);
26
27	init = true;
28	}
29		22
30	Scalar* a = reinterpret_cast<Scalar*>(pa);	23	Scalar* a = reinterpret_cast<Scalar*>(pa);
31	Scalar* x = reinterpret_cast<Scalar*>(px);	24	Scalar* x = reinterpret_cast<Scalar*>(px);
32	Scalar* y = reinterpret_cast<Scalar*>(py);	25	Scalar* y = reinterpret_cast<Scalar*>(py);
33	Scalar alpha = reinterpret_cast<Scalar>(palpha);	26	Scalar alpha = reinterpret_cast<Scalar>(palpha);
34	Scalar beta = reinterpret_cast<Scalar>(pbeta);	27	Scalar beta = reinterpret_cast<Scalar>(pbeta);
35		28
36	// check arguments	29	// check arguments
Lines 66-109 int EIGEN_BLAS_FUNC(symv) (char *uplo, i Link Here
66		59
67	return 1;	60	return 1;
68	}	61	}
69		62
70	// C := alphaxx' + C	63	// C := alphaxx' + C
71	int EIGEN_BLAS_FUNC(syr)(char uplo, int n, RealScalar palpha, RealScalar px, int incx, RealScalar pc, int *ldc)	64	int EIGEN_BLAS_FUNC(syr)(char uplo, int n, RealScalar palpha, RealScalar px, int incx, RealScalar pc, int *ldc)
72	{	65	{
73		66
74	// typedef void (functype)(int, const Scalar , int, Scalar *, int, Scalar);
75	// static functype func[2];
76
77	// static bool init = false;
78	// if(!init)
79	// {
80	// for(int k=0; k<2; ++k)
81	// func[k] = 0;
82	//
83	// func[UP] = (internal::selfadjoint_product<Scalar,ColMajor,ColMajor,false,UpperTriangular>::run);
84	// func[LO] = (internal::selfadjoint_product<Scalar,ColMajor,ColMajor,false,LowerTriangular>::run);
85
86	// init = true;
87	// }
88	typedef void (functype)(int, Scalar, int, const Scalar, const Scalar, const Scalar&);	67	typedef void (functype)(int, Scalar, int, const Scalar, const Scalar, const Scalar&);
89	static functype func[2];	68	static const functype func[2] = {
90		69	// array index: UP
91	static bool init = false;	70	(selfadjoint_rank1_update<Scalar,int,ColMajor,Upper,false,Conj>::run),
92	if(!init)	71	// array index: LO
93	{	72	(selfadjoint_rank1_update<Scalar,int,ColMajor,Lower,false,Conj>::run),
94	for(int k=0; k<2; ++k)	73	};
95	func[k] = 0;
96
97	func[UP] = (selfadjoint_rank1_update<Scalar,int,ColMajor,Upper,false,Conj>::run);
98	func[LO] = (selfadjoint_rank1_update<Scalar,int,ColMajor,Lower,false,Conj>::run);
99
100	init = true;
101	}
102		74
103	Scalar* x = reinterpret_cast<Scalar*>(px);	75	Scalar* x = reinterpret_cast<Scalar*>(px);
104	Scalar* c = reinterpret_cast<Scalar*>(pc);	76	Scalar* c = reinterpret_cast<Scalar*>(pc);
105	Scalar alpha = reinterpret_cast<Scalar>(palpha);	77	Scalar alpha = reinterpret_cast<Scalar>(palpha);
106		78
107	int info = 0;	79	int info = 0;
108	if(UPLO(*uplo)==INVALID) info = 1;	80	if(UPLO(*uplo)==INVALID) info = 1;
109	else if(*n<0) info = 2;	81	else if(*n<0) info = 2;
Lines 126-169 int EIGEN_BLAS_FUNC(syr)(char *uplo, int Link Here
126	if(x_cpy!=x) delete[] x_cpy;	98	if(x_cpy!=x) delete[] x_cpy;
127		99
128	return 1;	100	return 1;
129	}	101	}
130		102
131	// C := alphaxy' + alphayx' + C	103	// C := alphaxy' + alphayx' + C
132	int EIGEN_BLAS_FUNC(syr2)(char uplo, int n, RealScalar palpha, RealScalar px, int incx, RealScalar py, int incy, RealScalar pc, int *ldc)	104	int EIGEN_BLAS_FUNC(syr2)(char uplo, int n, RealScalar palpha, RealScalar px, int incx, RealScalar py, int incy, RealScalar pc, int *ldc)
133	{	105	{
134	// typedef void (functype)(int, const Scalar , int, const Scalar , int, Scalar , int, Scalar);
135	// static functype func[2];
136	//
137	// static bool init = false;
138	// if(!init)
139	// {
140	// for(int k=0; k<2; ++k)
141	// func[k] = 0;
142	//
143	// func[UP] = (internal::selfadjoint_product<Scalar,ColMajor,ColMajor,false,UpperTriangular>::run);
144	// func[LO] = (internal::selfadjoint_product<Scalar,ColMajor,ColMajor,false,LowerTriangular>::run);
145	//
146	// init = true;
147	// }
148	typedef void (functype)(int, Scalar, int, const Scalar, const Scalar, Scalar);	106	typedef void (functype)(int, Scalar, int, const Scalar, const Scalar, Scalar);
149	static functype func[2];	107	static const functype func[2] = {
150		108	// array index: UP
151	static bool init = false;	109	(internal::rank2_update_selector<Scalar,int,Upper>::run),
152	if(!init)	110	// array index: LO
153	{	111	(internal::rank2_update_selector<Scalar,int,Lower>::run),
154	for(int k=0; k<2; ++k)	112	};
155	func[k] = 0;
156
157	func[UP] = (internal::rank2_update_selector<Scalar,int,Upper>::run);
158	func[LO] = (internal::rank2_update_selector<Scalar,int,Lower>::run);
159
160	init = true;
161	}
162		113
163	Scalar* x = reinterpret_cast<Scalar*>(px);	114	Scalar* x = reinterpret_cast<Scalar*>(px);
164	Scalar* y = reinterpret_cast<Scalar*>(py);	115	Scalar* y = reinterpret_cast<Scalar*>(py);
165	Scalar* c = reinterpret_cast<Scalar*>(pc);	116	Scalar* c = reinterpret_cast<Scalar*>(pc);
166	Scalar alpha = reinterpret_cast<Scalar>(palpha);	117	Scalar alpha = reinterpret_cast<Scalar>(palpha);
167		118
168	int info = 0;	119	int info = 0;
169	if(UPLO(*uplo)==INVALID) info = 1;	120	if(UPLO(*uplo)==INVALID) info = 1;
Lines 229-257 int EIGEN_BLAS_FUNC(syr2)(char *uplo, in Link Here
229	* A := alphaxx' + A,	180	* A := alphaxx' + A,
230	*	181	*
231	* where alpha is a real scalar, x is an n element vector and A is an	182	* where alpha is a real scalar, x is an n element vector and A is an
232	* n by n symmetric matrix, supplied in packed form.	183	* n by n symmetric matrix, supplied in packed form.
233	*/	184	*/
234	int EIGEN_BLAS_FUNC(spr)(char uplo, int n, Scalar palpha, Scalar px, int incx, Scalar pap)	185	int EIGEN_BLAS_FUNC(spr)(char uplo, int n, Scalar palpha, Scalar px, int incx, Scalar pap)
235	{	186	{
236	typedef void (functype)(int, Scalar, const Scalar*, Scalar);	187	typedef void (functype)(int, Scalar, const Scalar*, Scalar);
237	static functype func[2];	188	static const functype func[2] = {
238		189	// array index: UP
239	static bool init = false;	190	(internal::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Upper,false,false>::run),
240	if(!init)	191	// array index: LO
241	{	192	(internal::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Lower,false,false>::run),
242	for(int k=0; k<2; ++k)	193	};
243	func[k] = 0;
244
245	func[UP] = (internal::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Upper,false,false>::run);
246	func[LO] = (internal::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Lower,false,false>::run);
247
248	init = true;
249	}
250		194
251	Scalar* x = reinterpret_cast<Scalar*>(px);	195	Scalar* x = reinterpret_cast<Scalar*>(px);
252	Scalar* ap = reinterpret_cast<Scalar*>(pap);	196	Scalar* ap = reinterpret_cast<Scalar*>(pap);
253	Scalar alpha = reinterpret_cast<Scalar>(palpha);	197	Scalar alpha = reinterpret_cast<Scalar>(palpha);
254		198
255	int info = 0;	199	int info = 0;
256	if(UPLO(*uplo)==INVALID) info = 1;	200	if(UPLO(*uplo)==INVALID) info = 1;
257	else if(*n<0) info = 2;	201	else if(*n<0) info = 2;
Lines 280-308 int EIGEN_BLAS_FUNC(spr)(char *uplo, int Link Here
280	* A := alphaxy' + alphayx' + A,	224	* A := alphaxy' + alphayx' + A,
281	*	225	*
282	* where alpha is a scalar, x and y are n element vectors and A is an	226	* where alpha is a scalar, x and y are n element vectors and A is an
283	* n by n symmetric matrix, supplied in packed form.	227	* n by n symmetric matrix, supplied in packed form.
284	*/	228	*/
285	int EIGEN_BLAS_FUNC(spr2)(char uplo, int n, RealScalar palpha, RealScalar px, int incx, RealScalar py, int incy, RealScalar pap)	229	int EIGEN_BLAS_FUNC(spr2)(char uplo, int n, RealScalar palpha, RealScalar px, int incx, RealScalar py, int incy, RealScalar pap)
286	{	230	{
287	typedef void (functype)(int, Scalar, const Scalar, const Scalar, Scalar);	231	typedef void (functype)(int, Scalar, const Scalar, const Scalar, Scalar);
288	static functype func[2];	232	static const functype func[2] = {
289		233	// array index: UP
290	static bool init = false;	234	(internal::packed_rank2_update_selector<Scalar,int,Upper>::run),
291	if(!init)	235	// array index: LO
292	{	236	(internal::packed_rank2_update_selector<Scalar,int,Lower>::run),
293	for(int k=0; k<2; ++k)	237	};
294	func[k] = 0;
295
296	func[UP] = (internal::packed_rank2_update_selector<Scalar,int,Upper>::run);
297	func[LO] = (internal::packed_rank2_update_selector<Scalar,int,Lower>::run);
298
299	init = true;
300	}
301		238
302	Scalar* x = reinterpret_cast<Scalar*>(px);	239	Scalar* x = reinterpret_cast<Scalar*>(px);
303	Scalar* y = reinterpret_cast<Scalar*>(py);	240	Scalar* y = reinterpret_cast<Scalar*>(py);
304	Scalar* ap = reinterpret_cast<Scalar*>(pap);	241	Scalar* ap = reinterpret_cast<Scalar*>(pap);
305	Scalar alpha = reinterpret_cast<Scalar>(palpha);	242	Scalar alpha = reinterpret_cast<Scalar>(palpha);
306		243
307	int info = 0;	244	int info = 0;
308	if(UPLO(*uplo)==INVALID) info = 1;	245	if(UPLO(*uplo)==INVALID) info = 1;