Attachment 647 Details for Bug 1152 – Fix blas static initialization race

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

[patch] Fix blas static initialization race

blas (text/plain), 69.04 KB, created by Benoit Jacob on 2016-01-26 16:04:42 UTC

(hide)

Description:

Filename:

MIME Type:

Creator: Benoit Jacob

Created: 2016-01-26 16:04:42 UTC

Size: 69.04 KB

patch

obsolete

># HG changeset patch
># Parent 48252c4c437ca7abace3bb4479eb6f4e15c56f51
>diff --git a/blas/level2_cplx_impl.h b/blas/level2_cplx_impl.h
>--- a/blas/level2_cplx_impl.h
>+++ b/blas/level2_cplx_impl.h
>@@ -14,29 +14,22 @@
>   *     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 functype func[2];
>-
>-  static bool init = false;
>-  if(!init)
>-  {
>-    for(int k=0; k<2; ++k)
>-      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;
>-  }
>+  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),
>+  };
> 
>   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
>@@ -106,29 +99,22 @@ int EIGEN_BLAS_FUNC(hemv)(char *uplo, in
>   *     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 functype func[2];
>-
>-  static bool init = false;
>-  if(!init)
>-  {
>-    for(int k=0; k<2; ++k)
>-      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;
>-  }
>+  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),
>+  };
> 
>   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;
>@@ -157,29 +143,22 @@ int EIGEN_BLAS_FUNC(hpr)(char *uplo, int
>   *     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 functype func[2];
>-
>-  static bool init = false;
>-  if(!init)
>-  {
>-    for(int k=0; k<2; ++k)
>-      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;
>-  }
>+  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),
>+  };
> 
>   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;
>@@ -212,29 +191,22 @@ int EIGEN_BLAS_FUNC(hpr2)(char *uplo, in
>   *     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 functype func[2];
>-
>-  static bool init = false;
>-  if(!init)
>-  {
>-    for(int k=0; k<2; ++k)
>-      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;
>-  }
>+  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),
>+  };
> 
>   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;
>@@ -266,29 +238,22 @@ int EIGEN_BLAS_FUNC(her)(char *uplo, int
>   *     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 functype func[2];
>-
>-  static bool init = false;
>-  if(!init)
>-  {
>-    for(int k=0; k<2; ++k)
>-      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;
>-  }
>+  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),
>+  };
> 
>   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;
>diff --git a/blas/level2_impl.h b/blas/level2_impl.h
>--- a/blas/level2_impl.h
>+++ b/blas/level2_impl.h
>@@ -21,30 +21,25 @@ struct general_matrix_vector_product_wra
>         <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 functype func[4];
>-
>-  static bool init = false;
>-  if(!init)
>-  {
>-    for(int k=0; k<4; ++k)
>-      func[k] = 0;
>-
>-    func[NOTR] = (general_matrix_vector_product_wrapper<int,Scalar,ColMajor,false,false>::run);
>-    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;
>-  }
>+  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
>+  };
> 
>   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
>@@ -85,42 +80,46 @@ int EIGEN_BLAS_FUNC(gemv)(char *opa, int
>   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 functype func[16];
>-
>-  static bool init = false;
>-  if(!init)
>-  {
>-    for(int k=0; k<16; ++k)
>-      func[k] = 0;
>-
>-    func[NOTR  | (UP << 2) | (NUNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0,       false,ColMajor>::run);
>-    func[TR    | (UP << 2) | (NUNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0,       false,RowMajor>::run);
>-    func[ADJ   | (UP << 2) | (NUNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0,       Conj, RowMajor>::run);
>-
>-    func[NOTR  | (LO << 2) | (NUNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0,       false,ColMajor>::run);
>-    func[TR    | (LO << 2) | (NUNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0,       false,RowMajor>::run);
>-    func[ADJ   | (LO << 2) | (NUNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0,       Conj, RowMajor>::run);
>-
>-    func[NOTR  | (UP << 2) | (UNIT  << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,ColMajor>::run);
>-    func[TR    | (UP << 2) | (UNIT  << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,RowMajor>::run);
>-    func[ADJ   | (UP << 2) | (UNIT  << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,Conj, RowMajor>::run);
>-
>-    func[NOTR  | (LO << 2) | (UNIT  << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,ColMajor>::run);
>-    func[TR    | (LO << 2) | (UNIT  << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,RowMajor>::run);
>-    func[ADJ   | (LO << 2) | (UNIT  << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,Conj, RowMajor>::run);
>-
>-    init = true;
>-  }
>+  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;
>@@ -140,42 +139,46 @@ int EIGEN_BLAS_FUNC(trsv)(char *uplo, ch
>   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 functype func[16];
>-
>-  static bool init = false;
>-  if(!init)
>-  {
>-    for(int k=0; k<16; ++k)
>-      func[k] = 0;
>-
>-    func[NOTR  | (UP << 2) | (NUNIT << 3)] = (internal::triangular_matrix_vector_product<int,Upper|0,       Scalar,false,Scalar,false,ColMajor>::run);
>-    func[TR    | (UP << 2) | (NUNIT << 3)] = (internal::triangular_matrix_vector_product<int,Lower|0,       Scalar,false,Scalar,false,RowMajor>::run);
>-    func[ADJ   | (UP << 2) | (NUNIT << 3)] = (internal::triangular_matrix_vector_product<int,Lower|0,       Scalar,Conj, Scalar,false,RowMajor>::run);
>-
>-    func[NOTR  | (LO << 2) | (NUNIT << 3)] = (internal::triangular_matrix_vector_product<int,Lower|0,       Scalar,false,Scalar,false,ColMajor>::run);
>-    func[TR    | (LO << 2) | (NUNIT << 3)] = (internal::triangular_matrix_vector_product<int,Upper|0,       Scalar,false,Scalar,false,RowMajor>::run);
>-    func[ADJ   | (LO << 2) | (NUNIT << 3)] = (internal::triangular_matrix_vector_product<int,Upper|0,       Scalar,Conj, Scalar,false,RowMajor>::run);
>-
>-    func[NOTR  | (UP << 2) | (UNIT  << 3)] = (internal::triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,false,Scalar,false,ColMajor>::run);
>-    func[TR    | (UP << 2) | (UNIT  << 3)] = (internal::triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,false,Scalar,false,RowMajor>::run);
>-    func[ADJ   | (UP << 2) | (UNIT  << 3)] = (internal::triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,Conj, Scalar,false,RowMajor>::run);
>-
>-    func[NOTR  | (LO << 2) | (UNIT  << 3)] = (internal::triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,false,Scalar,false,ColMajor>::run);
>-    func[TR    | (LO << 2) | (UNIT  << 3)] = (internal::triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,false,Scalar,false,RowMajor>::run);
>-    func[ADJ   | (LO << 2) | (UNIT  << 3)] = (internal::triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,Conj, Scalar,false,RowMajor>::run);
>-
>-    init = true;
>-  }
>+  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;
>@@ -341,42 +344,46 @@ int EIGEN_BLAS_FUNC(tbmv)(char *uplo, ch
>   *  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 functype func[16];
>-
>-  static bool init = false;
>-  if(!init)
>-  {
>-    for(int i=0; i<16; ++i)
>-      func[i] = 0;
>-
>-    func[NOTR  | (UP << 2) | (NUNIT << 3)] = (internal::band_solve_triangular_selector<int,Upper|0,       Scalar,false,Scalar,ColMajor>::run);
>-    func[TR    | (UP << 2) | (NUNIT << 3)] = (internal::band_solve_triangular_selector<int,Lower|0,       Scalar,false,Scalar,RowMajor>::run);
>-    func[ADJ   | (UP << 2) | (NUNIT << 3)] = (internal::band_solve_triangular_selector<int,Lower|0,       Scalar,Conj, Scalar,RowMajor>::run);
>-
>-    func[NOTR  | (LO << 2) | (NUNIT << 3)] = (internal::band_solve_triangular_selector<int,Lower|0,       Scalar,false,Scalar,ColMajor>::run);
>-    func[TR    | (LO << 2) | (NUNIT << 3)] = (internal::band_solve_triangular_selector<int,Upper|0,       Scalar,false,Scalar,RowMajor>::run);
>-    func[ADJ   | (LO << 2) | (NUNIT << 3)] = (internal::band_solve_triangular_selector<int,Upper|0,       Scalar,Conj, Scalar,RowMajor>::run);
>-
>-    func[NOTR  | (UP << 2) | (UNIT  << 3)] = (internal::band_solve_triangular_selector<int,Upper|UnitDiag,Scalar,false,Scalar,ColMajor>::run);
>-    func[TR    | (UP << 2) | (UNIT  << 3)] = (internal::band_solve_triangular_selector<int,Lower|UnitDiag,Scalar,false,Scalar,RowMajor>::run);
>-    func[ADJ   | (UP << 2) | (UNIT  << 3)] = (internal::band_solve_triangular_selector<int,Lower|UnitDiag,Scalar,Conj, Scalar,RowMajor>::run);
>-
>-    func[NOTR  | (LO << 2) | (UNIT  << 3)] = (internal::band_solve_triangular_selector<int,Lower|UnitDiag,Scalar,false,Scalar,ColMajor>::run);
>-    func[TR    | (LO << 2) | (UNIT  << 3)] = (internal::band_solve_triangular_selector<int,Upper|UnitDiag,Scalar,false,Scalar,RowMajor>::run);
>-    func[ADJ   | (LO << 2) | (UNIT  << 3)] = (internal::band_solve_triangular_selector<int,Upper|UnitDiag,Scalar,Conj, Scalar,RowMajor>::run);
>-
>-    init = true;
>-  }
>+  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;
>@@ -411,42 +418,46 @@ int EIGEN_BLAS_FUNC(tbsv)(char *uplo, ch
>   *     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 functype func[16];
>-
>-  static bool init = false;
>-  if(!init)
>-  {
>-    for(int k=0; k<16; ++k)
>-      func[k] = 0;
>-
>-    func[NOTR  | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Upper|0,       Scalar,false,Scalar,false,ColMajor>::run);
>-    func[TR    | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Lower|0,       Scalar,false,Scalar,false,RowMajor>::run);
>-    func[ADJ   | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Lower|0,       Scalar,Conj, Scalar,false,RowMajor>::run);
>-
>-    func[NOTR  | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Lower|0,       Scalar,false,Scalar,false,ColMajor>::run);
>-    func[TR    | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Upper|0,       Scalar,false,Scalar,false,RowMajor>::run);
>-    func[ADJ   | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Upper|0,       Scalar,Conj, Scalar,false,RowMajor>::run);
>-
>-    func[NOTR  | (UP << 2) | (UNIT  << 3)] = (internal::packed_triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,false,Scalar,false,ColMajor>::run);
>-    func[TR    | (UP << 2) | (UNIT  << 3)] = (internal::packed_triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,false,Scalar,false,RowMajor>::run);
>-    func[ADJ   | (UP << 2) | (UNIT  << 3)] = (internal::packed_triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,Conj, Scalar,false,RowMajor>::run);
>-
>-    func[NOTR  | (LO << 2) | (UNIT  << 3)] = (internal::packed_triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,false,Scalar,false,ColMajor>::run);
>-    func[TR    | (LO << 2) | (UNIT  << 3)] = (internal::packed_triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,false,Scalar,false,RowMajor>::run);
>-    func[ADJ   | (LO << 2) | (UNIT  << 3)] = (internal::packed_triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,Conj, Scalar,false,RowMajor>::run);
>-
>-    init = true;
>-  }
>+  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;
>@@ -482,42 +493,46 @@ int EIGEN_BLAS_FUNC(tpmv)(char *uplo, ch
>   *  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 functype func[16];
>-
>-  static bool init = false;
>-  if(!init)
>-  {
>-    for(int k=0; k<16; ++k)
>-      func[k] = 0;
>-
>-    func[NOTR  | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0,       false,ColMajor>::run);
>-    func[TR    | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0,       false,RowMajor>::run);
>-    func[ADJ   | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0,       Conj, RowMajor>::run);
>-
>-    func[NOTR  | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0,       false,ColMajor>::run);
>-    func[TR    | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0,       false,RowMajor>::run);
>-    func[ADJ   | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0,       Conj, RowMajor>::run);
>-
>-    func[NOTR  | (UP << 2) | (UNIT  << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,ColMajor>::run);
>-    func[TR    | (UP << 2) | (UNIT  << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,RowMajor>::run);
>-    func[ADJ   | (UP << 2) | (UNIT  << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,Conj, RowMajor>::run);
>-
>-    func[NOTR  | (LO << 2) | (UNIT  << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,ColMajor>::run);
>-    func[TR    | (LO << 2) | (UNIT  << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,RowMajor>::run);
>-    func[ADJ   | (LO << 2) | (UNIT  << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,Conj, RowMajor>::run);
>-
>-    init = true;
>-  }
>+  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;
>diff --git a/blas/level2_real_impl.h b/blas/level2_real_impl.h
>--- a/blas/level2_real_impl.h
>+++ b/blas/level2_real_impl.h
>@@ -8,29 +8,22 @@
> // 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 functype func[2];
>-
>-  static bool init = false;
>-  if(!init)
>-  {
>-    for(int k=0; k<2; ++k)
>-      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;
>-  }
>+  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),
>+  };
> 
>   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
>@@ -66,44 +59,23 @@ int EIGEN_BLAS_FUNC(symv) (char *uplo, i
> 
>   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 functype func[2];
>-
>-  static bool init = false;
>-  if(!init)
>-  {
>-    for(int k=0; k<2; ++k)
>-      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;
>-  }
>+  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),
>+  };
> 
>   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;
>@@ -126,44 +98,23 @@ int EIGEN_BLAS_FUNC(syr)(char *uplo, int
>   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 functype func[2];
>-
>-  static bool init = false;
>-  if(!init)
>-  {
>-    for(int k=0; k<2; ++k)
>-      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;
>-  }
>+  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),
>+  };
> 
>   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;
>@@ -229,29 +180,22 @@ int EIGEN_BLAS_FUNC(syr2)(char *uplo, in
>   *     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 functype func[2];
>-
>-  static bool init = false;
>-  if(!init)
>-  {
>-    for(int k=0; k<2; ++k)
>-      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;
>-  }
>+  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),
>+  };
> 
>   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;
>@@ -280,29 +224,22 @@ int EIGEN_BLAS_FUNC(spr)(char *uplo, int
>   *     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 functype func[2];
>-
>-  static bool init = false;
>-  if(!init)
>-  {
>-    for(int k=0; k<2; ++k)
>-      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;
>-  }
>+  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),
>+  };
> 
>   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;
>diff --git a/blas/level3_impl.h b/blas/level3_impl.h
>--- a/blas/level3_impl.h
>+++ b/blas/level3_impl.h
>@@ -8,34 +8,39 @@
> // 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 functype func[12];
>-
>-  static bool init = false;
>-  if(!init)
>-  {
>-    for(int i=0; i<12; ++i)
>-      func[i] = 0;
>-    func[NOTR  | (NOTR << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,ColMajor,false,Scalar,ColMajor,false,ColMajor>::run);
>-    func[TR    | (NOTR << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,false,Scalar,ColMajor,false,ColMajor>::run);
>-    func[ADJ   | (NOTR << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,false,ColMajor>::run);
>-    func[NOTR  | (TR   << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,false,ColMajor>::run);
>-    func[TR    | (TR   << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,false,Scalar,RowMajor,false,ColMajor>::run);
>-    func[ADJ   | (TR   << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,RowMajor,false,ColMajor>::run);
>-    func[NOTR  | (ADJ  << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,Conj, ColMajor>::run);
>-    func[TR    | (ADJ  << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,false,Scalar,RowMajor,Conj, ColMajor>::run);
>-    func[ADJ   | (ADJ  << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,RowMajor,Conj, ColMajor>::run);
>-    init = true;
>-  }
>+  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;
>@@ -68,59 +73,74 @@ int EIGEN_BLAS_FUNC(gemm)(char *opa, cha
>   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 functype func[32];
>-
>-  static bool init = false;
>-  if(!init)
>-  {
>-    for(int i=0; i<32; ++i)
>-      func[i] = 0;
>-
>-    func[NOTR  | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|0,          false,ColMajor,ColMajor>::run);
>-    func[TR    | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|0,          false,RowMajor,ColMajor>::run);
>-    func[ADJ   | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|0,          Conj, RowMajor,ColMajor>::run);
>-
>-    func[NOTR  | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|0,          false,ColMajor,ColMajor>::run);
>-    func[TR    | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|0,          false,RowMajor,ColMajor>::run);
>-    func[ADJ   | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|0,          Conj, RowMajor,ColMajor>::run);
>-
>-    func[NOTR  | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|0,          false,ColMajor,ColMajor>::run);
>-    func[TR    | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|0,          false,RowMajor,ColMajor>::run);
>-    func[ADJ   | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|0,          Conj, RowMajor,ColMajor>::run);
>-
>-    func[NOTR  | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|0,          false,ColMajor,ColMajor>::run);
>-    func[TR    | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|0,          false,RowMajor,ColMajor>::run);
>-    func[ADJ   | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|0,          Conj, RowMajor,ColMajor>::run);
>-
>-
>-    func[NOTR  | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|UnitDiag,false,ColMajor,ColMajor>::run);
>-    func[TR    | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|UnitDiag,false,RowMajor,ColMajor>::run);
>-    func[ADJ   | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|UnitDiag,Conj, RowMajor,ColMajor>::run);
>-
>-    func[NOTR  | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|UnitDiag,false,ColMajor,ColMajor>::run);
>-    func[TR    | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|UnitDiag,false,RowMajor,ColMajor>::run);
>-    func[ADJ   | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|UnitDiag,Conj, RowMajor,ColMajor>::run);
>-
>-    func[NOTR  | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|UnitDiag,false,ColMajor,ColMajor>::run);
>-    func[TR    | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|UnitDiag,false,RowMajor,ColMajor>::run);
>-    func[ADJ   | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|UnitDiag,Conj, RowMajor,ColMajor>::run);
>-
>-    func[NOTR  | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|UnitDiag,false,ColMajor,ColMajor>::run);
>-    func[TR    | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|UnitDiag,false,RowMajor,ColMajor>::run);
>-    func[ADJ   | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|UnitDiag,Conj, RowMajor,ColMajor>::run);
>-
>-    init = true;
>-  }
>+  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;
>@@ -157,57 +177,74 @@ int EIGEN_BLAS_FUNC(trsm)(char *side, ch
> 
> 
> // 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 functype func[32];
>-  static bool init = false;
>-  if(!init)
>-  {
>-    for(int k=0; k<32; ++k)
>-      func[k] = 0;
>-
>-    func[NOTR  | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0,          true, ColMajor,false,ColMajor,false,ColMajor>::run);
>-    func[TR    | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0,          true, RowMajor,false,ColMajor,false,ColMajor>::run);
>-    func[ADJ   | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0,          true, RowMajor,Conj, ColMajor,false,ColMajor>::run);
>-
>-    func[NOTR  | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0,          false,ColMajor,false,ColMajor,false,ColMajor>::run);
>-    func[TR    | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0,          false,ColMajor,false,RowMajor,false,ColMajor>::run);
>-    func[ADJ   | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0,          false,ColMajor,false,RowMajor,Conj, ColMajor>::run);
>-
>-    func[NOTR  | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0,          true, ColMajor,false,ColMajor,false,ColMajor>::run);
>-    func[TR    | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0,          true, RowMajor,false,ColMajor,false,ColMajor>::run);
>-    func[ADJ   | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0,          true, RowMajor,Conj, ColMajor,false,ColMajor>::run);
>-
>-    func[NOTR  | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0,          false,ColMajor,false,ColMajor,false,ColMajor>::run);
>-    func[TR    | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0,          false,ColMajor,false,RowMajor,false,ColMajor>::run);
>-    func[ADJ   | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0,          false,ColMajor,false,RowMajor,Conj, ColMajor>::run);
>-
>-    func[NOTR  | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,true, ColMajor,false,ColMajor,false,ColMajor>::run);
>-    func[TR    | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,true, RowMajor,false,ColMajor,false,ColMajor>::run);
>-    func[ADJ   | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,true, RowMajor,Conj, ColMajor,false,ColMajor>::run);
>-
>-    func[NOTR  | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,false,ColMajor,false,ColMajor,false,ColMajor>::run);
>-    func[TR    | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,false,ColMajor,false,RowMajor,false,ColMajor>::run);
>-    func[ADJ   | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,false,ColMajor,false,RowMajor,Conj, ColMajor>::run);
>-
>-    func[NOTR  | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,true, ColMajor,false,ColMajor,false,ColMajor>::run);
>-    func[TR    | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,true, RowMajor,false,ColMajor,false,ColMajor>::run);
>-    func[ADJ   | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,true, RowMajor,Conj, ColMajor,false,ColMajor>::run);
>-
>-    func[NOTR  | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,false,ColMajor,false,ColMajor,false,ColMajor>::run);
>-    func[TR    | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,false,ColMajor,false,RowMajor,false,ColMajor>::run);
>-    func[ADJ   | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,false,ColMajor,false,RowMajor,Conj, ColMajor>::run);
>-
>-    init = true;
>-  }
>+  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;
>@@ -312,16 +349,35 @@ int EIGEN_BLAS_FUNC(symm)(char *side, ch
> }
> 
> // 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
>+<<<<<<< local
>+  typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, const 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
>+  };
>+=======
>   typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, const Scalar&, internal::level3_blocking<Scalar,Scalar>&);
>   static functype func[8];
> 
>   static bool init = false;
>   if(!init)
>   {
>     for(int i=0; i<8; ++i)
>       func[i] = 0;
>@@ -331,16 +387,17 @@ int EIGEN_BLAS_FUNC(syrk)(char *uplo, ch
>     func[ADJ   | (UP << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,ColMajor,false,Upper>::run);
> 
>     func[NOTR  | (LO << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,ColMajor,Conj, Lower>::run);
>     func[TR    | (LO << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,false,Scalar,ColMajor,ColMajor,Conj, Lower>::run);
>     func[ADJ   | (LO << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,ColMajor,false,Lower>::run);
> 
>     init = true;
>   }
>+>>>>>>> other
>   #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;
>@@ -518,16 +575,33 @@ int EIGEN_BLAS_FUNC(hemm)(char *side, ch
> }
> 
> // 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";
> 
>+<<<<<<< local
>+  typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, const 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
>+  };
>+=======
>   typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, const Scalar&, internal::level3_blocking<Scalar,Scalar>&);
>   static functype func[8];
> 
>   static bool init = false;
>   if(!init)
>   {
>     for(int i=0; i<8; ++i)
>       func[i] = 0;
>@@ -535,16 +609,17 @@ int EIGEN_BLAS_FUNC(herk)(char *uplo, ch
>     func[NOTR  | (UP << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,Conj, ColMajor,Upper>::run);
>     func[ADJ   | (UP << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,false,ColMajor,Upper>::run);
> 
>     func[NOTR  | (LO << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,Conj, ColMajor,Lower>::run);
>     func[ADJ   | (LO << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,false,ColMajor,Lower>::run);
> 
>     init = true;
>   }
>+>>>>>>> other
> 
>   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";
>

Actions: View | Diff

Attachments on bug 1152: 647 | 648