cblas3.hpp

Go to the documentation of this file.
00001 /*
00002  * 
00003  * Copyright (c) Kresimir Fresl 2002 
00004  *
00005  * Distributed under the Boost Software License, Version 1.0.
00006  * (See accompanying file LICENSE_1_0.txt or copy at
00007  * http://www.boost.org/LICENSE_1_0.txt)
00008  *
00009  * Author acknowledges the support of the Faculty of Civil Engineering, 
00010  * University of Zagreb, Croatia.
00011  *
00012  */
00013 
00014 #ifndef BOOST_NUMERIC_BINDINGS_CBLAS_LEVEL_3_HPP
00015 #define BOOST_NUMERIC_BINDINGS_CBLAS_LEVEL_3_HPP
00016 
00017 #include <cassert>
00018 
00019 #include <boost/numeric/bindings/traits/traits.hpp>
00020 #include <boost/numeric/bindings/traits/type_traits.hpp>
00021 #include <boost/numeric/bindings/atlas/cblas3_overloads.hpp>
00022 #include <boost/numeric/bindings/atlas/cblas_enum.hpp>
00023 #include <boost/type_traits/same_traits.hpp>
00024 #include <boost/mpl/if.hpp>
00025 
00026 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00027 #  include <boost/static_assert.hpp>
00028 #endif
00029 
00030 namespace boost { namespace numeric { namespace bindings { 
00031 
00032   namespace atlas {
00033 
00034     // C <- alpha * op (A) * op (B) + beta * C 
00035     // op (A) == A || A^T || A^H
00036     template <typename T, typename MatrA, typename MatrB, typename MatrC>
00037     inline
00038     void gemm (CBLAS_TRANSPOSE const TransA, CBLAS_TRANSPOSE const TransB, 
00039                T const& alpha, MatrA const& a, MatrB const& b, 
00040                T const& beta, MatrC& c
00041                )
00042     {
00043 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00044       BOOST_STATIC_ASSERT((boost::is_same<
00045         typename traits::matrix_traits<MatrA>::matrix_structure, 
00046         traits::general_t
00047       >::value)); 
00048       BOOST_STATIC_ASSERT((boost::is_same<
00049         typename traits::matrix_traits<MatrB>::matrix_structure, 
00050         traits::general_t
00051       >::value)); 
00052       BOOST_STATIC_ASSERT((boost::is_same<
00053         typename traits::matrix_traits<MatrC>::matrix_structure, 
00054         traits::general_t
00055       >::value)); 
00056 
00057       BOOST_STATIC_ASSERT((boost::is_same<
00058         typename traits::matrix_traits<MatrA>::ordering_type,
00059         typename traits::matrix_traits<MatrB>::ordering_type
00060       >::value)); 
00061       BOOST_STATIC_ASSERT((boost::is_same<
00062         typename traits::matrix_traits<MatrA>::ordering_type,
00063         typename traits::matrix_traits<MatrC>::ordering_type
00064       >::value)); 
00065 #endif 
00066 
00067       assert (TransA == CblasNoTrans 
00068               || TransA == CblasTrans 
00069               || TransA == CblasConjTrans); 
00070       assert (TransB == CblasNoTrans 
00071               || TransB == CblasTrans 
00072               || TransB == CblasConjTrans); 
00073 
00074       int const m = TransA == CblasNoTrans
00075         ? traits::matrix_size1 (a)
00076         : traits::matrix_size2 (a);
00077       int const n = TransB == CblasNoTrans
00078         ? traits::matrix_size2 (b)
00079         : traits::matrix_size1 (b);
00080       int const k = TransA == CblasNoTrans
00081         ? traits::matrix_size2 (a)
00082         : traits::matrix_size1 (a); 
00083       assert (m == traits::matrix_size1 (c)); 
00084       assert (n == traits::matrix_size2 (c)); 
00085 #ifndef NDEBUG
00086       int const k1 = TransB == CblasNoTrans
00087         ? traits::matrix_size1 (b)
00088         : traits::matrix_size2 (b);
00089       assert (k == k1); 
00090 #endif
00091       // .. what about AtlasConj? 
00092 
00093       CBLAS_ORDER const stor_ord
00094         = enum_cast<CBLAS_ORDER const>
00095         (storage_order<
00096 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00097            typename traits::matrix_traits<MatrA>::ordering_type
00098 #else
00099            typename MatrA::orientation_category 
00100 #endif 
00101          >::value); 
00102 
00103       detail::gemm (stor_ord, TransA, TransB, m, n, k, alpha, 
00104 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
00105                     traits::matrix_storage (a), 
00106 #else
00107                     traits::matrix_storage_const (a), 
00108 #endif
00109                     traits::leading_dimension (a),
00110 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
00111                     traits::matrix_storage (b), 
00112 #else
00113                     traits::matrix_storage_const (b), 
00114 #endif
00115                     traits::leading_dimension (b),
00116                     beta, 
00117                     traits::matrix_storage (c), 
00118                     traits::leading_dimension (c)); 
00119     }
00120 
00121 
00122     // C <- alpha * A * B + beta * C 
00123     template <typename T, typename MatrA, typename MatrB, typename MatrC>
00124     inline
00125     void gemm (T const& alpha, MatrA const& a, MatrB const& b, 
00126                T const& beta, MatrC& c) 
00127     {
00128       gemm (CblasNoTrans, CblasNoTrans, alpha, a, b, beta, c) ;
00129     }
00130     
00131 
00132     // C <- A * B 
00133     template <typename MatrA, typename MatrB, typename MatrC>
00134     inline
00135     void gemm (MatrA const& a, MatrB const& b, MatrC& c) {
00136 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00137       typedef typename traits::matrix_traits<MatrC>::value_type val_t; 
00138 #else
00139       typedef typename MatrC::value_type val_t; 
00140 #endif 
00141       gemm (CblasNoTrans, CblasNoTrans, (val_t) 1, a, b, (val_t) 0, c);
00142     }
00143 
00144 
00145 
00146     // C <- alpha * A * B + beta * C 
00147     // C <- alpha * B * A + beta * C 
00148     // A == A^T
00149 
00150     namespace detail {
00151 
00152       template <typename T, typename SymmA, typename MatrB, typename MatrC>
00153       inline
00154       void symm (CBLAS_SIDE const side, CBLAS_UPLO const uplo, 
00155                  T const& alpha, SymmA const& a, MatrB const& b, 
00156                  T const& beta, MatrC& c)
00157       {
00158 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00159         BOOST_STATIC_ASSERT((boost::is_same<
00160           typename traits::matrix_traits<MatrB>::matrix_structure, 
00161           traits::general_t
00162         >::value));
00163         BOOST_STATIC_ASSERT((boost::is_same<
00164           typename traits::matrix_traits<MatrC>::matrix_structure, 
00165           traits::general_t
00166         >::value)); 
00167 
00168         BOOST_STATIC_ASSERT((boost::is_same<
00169           typename traits::matrix_traits<SymmA>::ordering_type,
00170           typename traits::matrix_traits<MatrB>::ordering_type
00171         >::value)); 
00172         BOOST_STATIC_ASSERT((boost::is_same<
00173           typename traits::matrix_traits<SymmA>::ordering_type,
00174           typename traits::matrix_traits<MatrC>::ordering_type
00175         >::value)); 
00176 #endif 
00177 
00178         assert (side == CblasLeft || side == CblasRight);
00179         assert (uplo == CblasUpper || uplo == CblasLower); 
00180 
00181         int const m = traits::matrix_size1 (c);
00182         int const n = traits::matrix_size2 (c);
00183 
00184         assert (side == CblasLeft 
00185                 ? m == traits::matrix_size1 (a) 
00186                   && m == traits::matrix_size2 (a)
00187                 : n == traits::matrix_size1 (a) 
00188                   && n == traits::matrix_size2 (a)); 
00189         assert (m == traits::matrix_size1 (b) 
00190                 && n == traits::matrix_size2 (b)); 
00191 
00192         CBLAS_ORDER const stor_ord
00193           = enum_cast<CBLAS_ORDER const>
00194           (storage_order<
00195 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00196             typename traits::matrix_traits<SymmA>::ordering_type
00197 #else
00198             typename SymmA::orientation_category 
00199 #endif 
00200            >::value); 
00201 
00202         symm (stor_ord, side, uplo,  
00203               m, n, alpha, 
00204 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
00205               traits::matrix_storage (a), 
00206 #else
00207               traits::matrix_storage_const (a), 
00208 #endif
00209               traits::leading_dimension (a),
00210 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
00211               traits::matrix_storage (b), 
00212 #else
00213               traits::matrix_storage_const (b), 
00214 #endif
00215               traits::leading_dimension (b),
00216               beta, 
00217               traits::matrix_storage (c), 
00218               traits::leading_dimension (c)); 
00219       }
00220 
00221     } // detail 
00222  
00223     // C <- alpha * A * B + beta * C 
00224     // C <- alpha * B * A + beta * C 
00225     // A == A^T
00226     template <typename T, typename SymmA, typename MatrB, typename MatrC>
00227     inline
00228     void symm (CBLAS_SIDE const side, CBLAS_UPLO const uplo, 
00229                T const& alpha, SymmA const& a, MatrB const& b, 
00230                T const& beta, MatrC& c)
00231     {
00232 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00233       BOOST_STATIC_ASSERT((boost::is_same<
00234         typename traits::matrix_traits<SymmA>::matrix_structure, 
00235         traits::general_t
00236       >::value)); 
00237 #endif 
00238 
00239       detail::symm (side, uplo, alpha, a, b, beta, c); 
00240     }
00241 
00242     template <typename T, typename SymmA, typename MatrB, typename MatrC>
00243     inline
00244     void symm (CBLAS_SIDE const side, 
00245                T const& alpha, SymmA const& a, MatrB const& b, 
00246                T const& beta, MatrC& c)
00247     {
00248 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00249       BOOST_STATIC_ASSERT((boost::is_same<
00250         typename traits::matrix_traits<SymmA>::matrix_structure, 
00251         traits::symmetric_t
00252       >::value)); 
00253 #endif 
00254 
00255       CBLAS_UPLO const uplo
00256         = enum_cast<CBLAS_UPLO const>
00257         (uplo_triang<
00258 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00259           typename traits::matrix_traits<SymmA>::uplo_type
00260 #else
00261           typename SymmA::packed_category 
00262 #endif 
00263          >::value); 
00264 
00265       detail::symm (side, uplo, alpha, a, b, beta, c); 
00266     }
00267 
00268 
00269 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00270 
00271     namespace detail {
00272 
00273       // C <- alpha * A * B + beta * C ;  A == A^T
00274       struct symm_left {
00275         template <typename T, typename SymmA, typename MatrB, typename MatrC>
00276         static void f (T const& alpha, SymmA const& a, MatrB const& b, 
00277                        T const& beta, MatrC& c) 
00278         {
00279 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00280           BOOST_STATIC_ASSERT((boost::is_same<
00281             typename traits::matrix_traits<SymmA>::matrix_structure, 
00282             traits::symmetric_t
00283           >::value)); 
00284           BOOST_STATIC_ASSERT((boost::is_same<
00285             typename traits::matrix_traits<MatrB>::matrix_structure, 
00286             traits::general_t
00287           >::value));
00288 #endif 
00289 
00290           int const m = traits::matrix_size1 (c);
00291           int const n = traits::matrix_size2 (c);
00292 
00293           assert (m == traits::matrix_size1 (a) 
00294                   && m == traits::matrix_size2 (a)); 
00295           assert (m == traits::matrix_size1 (b) 
00296                   && n == traits::matrix_size2 (b)); 
00297 
00298           CBLAS_ORDER const stor_ord
00299             = enum_cast<CBLAS_ORDER const>
00300             (storage_order<
00301               typename traits::matrix_traits<SymmA>::ordering_type
00302              >::value); 
00303 
00304           CBLAS_UPLO const uplo
00305             = enum_cast<CBLAS_UPLO const>
00306             (uplo_triang<
00307               typename traits::matrix_traits<SymmA>::uplo_type
00308              >::value); 
00309 
00310           symm (stor_ord, CblasLeft, uplo,  
00311                 m, n, alpha, 
00312                 traits::matrix_storage (a), traits::leading_dimension (a),
00313                 traits::matrix_storage (b), traits::leading_dimension (b),
00314                 beta, 
00315                 traits::matrix_storage (c), traits::leading_dimension (c)); 
00316         }
00317       }; 
00318 
00319       // C <- alpha * A * B + beta * C ;  B == B^T
00320       struct symm_right {
00321         template <typename T, typename MatrA, typename SymmB, typename MatrC>
00322         static void f (T const& alpha, MatrA const& a, SymmB const& b, 
00323                        T const& beta, MatrC& c) 
00324         {
00325 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00326           BOOST_STATIC_ASSERT((boost::is_same<
00327             typename traits::matrix_traits<MatrA>::matrix_structure, 
00328             traits::general_t
00329           >::value));
00330           BOOST_STATIC_ASSERT((boost::is_same<
00331             typename traits::matrix_traits<SymmB>::matrix_structure, 
00332             traits::symmetric_t
00333           >::value));
00334 #endif 
00335 
00336           int const m = traits::matrix_size1 (c);
00337           int const n = traits::matrix_size2 (c);
00338 
00339           assert (n == traits::matrix_size1 (b) 
00340                   && n == traits::matrix_size2 (b)); 
00341           assert (m == traits::matrix_size1 (a) 
00342                   && n == traits::matrix_size2 (a)); 
00343 
00344           CBLAS_ORDER const stor_ord
00345             = enum_cast<CBLAS_ORDER const>
00346             (storage_order<
00347               typename traits::matrix_traits<SymmB>::ordering_type
00348              >::value); 
00349  
00350           CBLAS_UPLO const uplo
00351             = enum_cast<CBLAS_UPLO const>
00352             (uplo_triang<
00353               typename traits::matrix_traits<SymmB>::uplo_type
00354              >::value); 
00355 
00356           symm (stor_ord, CblasRight, uplo,  
00357                 m, n, alpha, 
00358                 traits::matrix_storage (b), traits::leading_dimension (b),
00359                 traits::matrix_storage (a), traits::leading_dimension (a),
00360                 beta, 
00361                 traits::matrix_storage (c), traits::leading_dimension (c)); 
00362         }
00363       }; 
00364 
00365     } // detail 
00366     
00367     // C <- alpha * A * B + beta * C 
00368     // C <- alpha * B * A + beta * C 
00369     // A == A^T
00370     template <typename T, typename MatrA, typename MatrB, typename MatrC>
00371     inline
00372     void symm (T const& alpha, MatrA const& a, MatrB const& b, 
00373                T const& beta, MatrC& c)
00374     {
00375 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00376       BOOST_STATIC_ASSERT((boost::is_same<
00377         typename traits::matrix_traits<MatrC>::matrix_structure, 
00378         traits::general_t
00379       >::value)); 
00380 
00381       BOOST_STATIC_ASSERT((boost::is_same<
00382         typename traits::matrix_traits<MatrA>::ordering_type,
00383         typename traits::matrix_traits<MatrB>::ordering_type
00384       >::value)); 
00385       BOOST_STATIC_ASSERT((boost::is_same<
00386         typename traits::matrix_traits<MatrA>::ordering_type,
00387         typename traits::matrix_traits<MatrC>::ordering_type
00388       >::value)); 
00389 #endif 
00390 
00391       typedef typename
00392         boost::mpl::if_c<
00393           boost::is_same<
00394             typename traits::matrix_traits<MatrA>::matrix_structure, 
00395             traits::symmetric_t
00396           >::value,
00397           detail::symm_left, 
00398           detail::symm_right
00399         >::type functor; 
00400 
00401       functor::f (alpha, a, b, beta, c); 
00402     }
00403 
00404     // C <- A * B  
00405     // C <- B * A  
00406     template <typename MatrA, typename MatrB, typename MatrC>
00407     inline
00408     void symm (MatrA const& a, MatrB const& b, MatrC& c) {
00409       typedef typename traits::matrix_traits<MatrC>::value_type val_t; 
00410       symm ((val_t) 1, a, b, (val_t) 0, c);
00411     }
00412 
00413 #endif // BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00414 
00415 
00416 
00417     // C <- alpha * A * B + beta * C 
00418     // C <- alpha * B * A + beta * C 
00419     // A == A^H
00420 
00421     namespace detail {
00422 
00423       template <typename T, typename HermA, typename MatrB, typename MatrC>
00424       inline
00425       void hemm (CBLAS_SIDE const side, CBLAS_UPLO const uplo, 
00426                  T const& alpha, HermA const& a, MatrB const& b, 
00427                  T const& beta, MatrC& c)
00428       {
00429 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00430         BOOST_STATIC_ASSERT((boost::is_same<
00431           typename traits::matrix_traits<MatrB>::matrix_structure, 
00432           traits::general_t
00433         >::value));
00434         BOOST_STATIC_ASSERT((boost::is_same<
00435           typename traits::matrix_traits<MatrC>::matrix_structure, 
00436           traits::general_t
00437         >::value)); 
00438 
00439         BOOST_STATIC_ASSERT((boost::is_same<
00440           typename traits::matrix_traits<HermA>::ordering_type,
00441           typename traits::matrix_traits<MatrB>::ordering_type
00442         >::value)); 
00443         BOOST_STATIC_ASSERT((boost::is_same<
00444           typename traits::matrix_traits<HermA>::ordering_type,
00445           typename traits::matrix_traits<MatrC>::ordering_type
00446         >::value)); 
00447 #endif 
00448 
00449         assert (side == CblasLeft || side == CblasRight);
00450         assert (uplo == CblasUpper || uplo == CblasLower); 
00451 
00452         int const m = traits::matrix_size1 (c);
00453         int const n = traits::matrix_size2 (c);
00454 
00455         assert (side == CblasLeft 
00456                 ? m == traits::matrix_size1 (a) 
00457                   && m == traits::matrix_size2 (a)
00458                 : n == traits::matrix_size1 (a) 
00459                   && n == traits::matrix_size2 (a)); 
00460         assert (m == traits::matrix_size1 (b) 
00461                 && n == traits::matrix_size2 (b)); 
00462 
00463         CBLAS_ORDER const stor_ord
00464           = enum_cast<CBLAS_ORDER const>
00465           (storage_order<
00466 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00467             typename traits::matrix_traits<HermA>::ordering_type
00468 #else
00469             typename HermA::orientation_category 
00470 #endif 
00471            >::value); 
00472 
00473         hemm (stor_ord, side, uplo,  
00474               m, n, alpha, 
00475 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
00476               traits::matrix_storage (a), 
00477 #else
00478               traits::matrix_storage_const (a), 
00479 #endif
00480               traits::leading_dimension (a),
00481 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
00482               traits::matrix_storage (b), 
00483 #else
00484               traits::matrix_storage_const (b), 
00485 #endif
00486               traits::leading_dimension (b),
00487               beta, 
00488               traits::matrix_storage (c), 
00489               traits::leading_dimension (c)); 
00490       }
00491 
00492     } // detail 
00493  
00494     // C <- alpha * A * B + beta * C 
00495     // C <- alpha * B * A + beta * C 
00496     // A == A^H
00497     template <typename T, typename HermA, typename MatrB, typename MatrC>
00498     inline
00499     void hemm (CBLAS_SIDE const side, CBLAS_UPLO const uplo, 
00500                T const& alpha, HermA const& a, MatrB const& b, 
00501                T const& beta, MatrC& c)
00502     {
00503 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00504       BOOST_STATIC_ASSERT((boost::is_same<
00505         typename traits::matrix_traits<HermA>::matrix_structure, 
00506         traits::general_t
00507       >::value)); 
00508 #endif 
00509 
00510       detail::hemm (side, uplo, alpha, a, b, beta, c); 
00511     }
00512 
00513     template <typename T, typename HermA, typename MatrB, typename MatrC>
00514     inline
00515     void hemm (CBLAS_SIDE const side, 
00516                T const& alpha, HermA const& a, MatrB const& b, 
00517                T const& beta, MatrC& c)
00518     {
00519 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00520       BOOST_STATIC_ASSERT((boost::is_same<
00521         typename traits::matrix_traits<HermA>::matrix_structure, 
00522         traits::hermitian_t
00523       >::value)); 
00524 #endif 
00525 
00526       CBLAS_UPLO const uplo
00527         = enum_cast<CBLAS_UPLO const>
00528         (uplo_triang<
00529 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00530           typename traits::matrix_traits<HermA>::uplo_type
00531 #else
00532           typename HermA::packed_category 
00533 #endif 
00534          >::value); 
00535 
00536       detail::hemm (side, uplo, alpha, a, b, beta, c); 
00537     }
00538 
00539 
00540 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00541 
00542     namespace detail {
00543 
00544       // C <- alpha * A * B + beta * C ;  A == A^H
00545       struct hemm_left {
00546         template <typename T, typename HermA, typename MatrB, typename MatrC>
00547         static void f (T const& alpha, HermA const& a, MatrB const& b, 
00548                        T const& beta, MatrC& c) 
00549         {
00550 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00551           BOOST_STATIC_ASSERT((boost::is_same<
00552             typename traits::matrix_traits<HermA>::matrix_structure, 
00553             traits::hermitian_t
00554           >::value)); 
00555           BOOST_STATIC_ASSERT((boost::is_same<
00556             typename traits::matrix_traits<MatrB>::matrix_structure, 
00557             traits::general_t
00558           >::value));
00559 #endif 
00560 
00561           int const m = traits::matrix_size1 (c);
00562           int const n = traits::matrix_size2 (c);
00563 
00564           assert (m == traits::matrix_size1 (a) 
00565                   && m == traits::matrix_size2 (a)); 
00566           assert (m == traits::matrix_size1 (b) 
00567                   && n == traits::matrix_size2 (b)); 
00568 
00569           CBLAS_ORDER const stor_ord
00570             = enum_cast<CBLAS_ORDER const>
00571             (storage_order<
00572               typename traits::matrix_traits<HermA>::ordering_type
00573              >::value); 
00574 
00575           CBLAS_UPLO const uplo
00576             = enum_cast<CBLAS_UPLO const>
00577             (uplo_triang<
00578               typename traits::matrix_traits<HermA>::uplo_type
00579              >::value); 
00580 
00581           hemm (stor_ord, CblasLeft, uplo,  
00582                 m, n, alpha, 
00583                 traits::matrix_storage (a), traits::leading_dimension (a),
00584                 traits::matrix_storage (b), traits::leading_dimension (b),
00585                 beta, 
00586                 traits::matrix_storage (c), traits::leading_dimension (c)); 
00587         }
00588       }; 
00589 
00590       // C <- alpha * A * B + beta * C ;  B == B^H
00591       struct hemm_right {
00592         template <typename T, typename MatrA, typename HermB, typename MatrC>
00593         static void f (T const& alpha, MatrA const& a, HermB const& b, 
00594                        T const& beta, MatrC& c) 
00595         {
00596 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00597           BOOST_STATIC_ASSERT((boost::is_same<
00598             typename traits::matrix_traits<MatrA>::matrix_structure, 
00599             traits::general_t
00600           >::value));
00601           BOOST_STATIC_ASSERT((boost::is_same<
00602             typename traits::matrix_traits<HermB>::matrix_structure, 
00603             traits::hermitian_t
00604           >::value));
00605 #endif 
00606 
00607           int const m = traits::matrix_size1 (c);
00608           int const n = traits::matrix_size2 (c);
00609 
00610           assert (n == traits::matrix_size1 (b) 
00611                   && n == traits::matrix_size2 (b)); 
00612           assert (m == traits::matrix_size1 (a) 
00613                   && n == traits::matrix_size2 (a)); 
00614 
00615           CBLAS_ORDER const stor_ord
00616             = enum_cast<CBLAS_ORDER const>
00617             (storage_order<
00618               typename traits::matrix_traits<HermB>::ordering_type
00619              >::value); 
00620  
00621           CBLAS_UPLO const uplo
00622             = enum_cast<CBLAS_UPLO const>
00623             (uplo_triang<
00624               typename traits::matrix_traits<HermB>::uplo_type
00625              >::value); 
00626 
00627           hemm (stor_ord, CblasRight, uplo,  
00628                 m, n, alpha, 
00629                 traits::matrix_storage (b), traits::leading_dimension (b),
00630                 traits::matrix_storage (a), traits::leading_dimension (a),
00631                 beta, 
00632                 traits::matrix_storage (c), traits::leading_dimension (c)); 
00633         }
00634       }; 
00635 
00636     } 
00637     
00638     template <typename T, typename MatrA, typename MatrB, typename MatrC>
00639     inline
00640     void hemm (T const& alpha, MatrA const& a, MatrB const& b, 
00641                T const& beta, MatrC& c)
00642     {
00643 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00644       BOOST_STATIC_ASSERT((boost::is_same<
00645         typename traits::matrix_traits<MatrC>::matrix_structure, 
00646         traits::general_t
00647       >::value)); 
00648 
00649       BOOST_STATIC_ASSERT((boost::is_same<
00650         typename traits::matrix_traits<MatrA>::ordering_type,
00651         typename traits::matrix_traits<MatrB>::ordering_type
00652       >::value)); 
00653       BOOST_STATIC_ASSERT((boost::is_same<
00654         typename traits::matrix_traits<MatrA>::ordering_type,
00655         typename traits::matrix_traits<MatrC>::ordering_type
00656       >::value)); 
00657 #endif 
00658 
00659       typedef typename
00660         boost::mpl::if_c<
00661           boost::is_same<
00662             typename traits::matrix_traits<MatrA>::matrix_structure, 
00663             traits::hermitian_t
00664           >::value,
00665           detail::hemm_left, 
00666           detail::hemm_right
00667         >::type functor; 
00668 
00669       functor::f (alpha, a, b, beta, c); 
00670     }
00671 
00672     // C <- A * B  
00673     // C <- B * A  
00674     template <typename MatrA, typename MatrB, typename MatrC>
00675     inline
00676     void hemm (MatrA const& a, MatrB const& b, MatrC& c) {
00677       typedef typename traits::matrix_traits<MatrC>::value_type val_t; 
00678       hemm ((val_t) 1, a, b, (val_t) 0, c);
00679     }
00680 
00681 #endif // BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00682 
00683     
00684     // C <- alpha * A * A^T + beta * C
00685     // C <- alpha * A^T * A + beta * C
00686     // C == C^T
00687 
00688     namespace detail {
00689 
00690       template <typename T, typename MatrA, typename SymmC>
00691       inline
00692       void syrk (CBLAS_UPLO const uplo, CBLAS_TRANSPOSE trans, 
00693                  T const& alpha, MatrA const& a, 
00694                  T const& beta, SymmC& c)
00695       {
00696 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00697         BOOST_STATIC_ASSERT((boost::is_same<
00698           typename traits::matrix_traits<MatrA>::matrix_structure, 
00699           traits::general_t
00700         >::value));
00701         BOOST_STATIC_ASSERT((boost::is_same<
00702           typename traits::matrix_traits<MatrA>::ordering_type,
00703           typename traits::matrix_traits<SymmC>::ordering_type
00704         >::value)); 
00705 #endif 
00706 
00707         assert (uplo == CblasUpper || uplo == CblasLower); 
00708         assert (trans == CblasNoTrans 
00709                 || trans == CblasTrans 
00710                 || trans == CblasConjTrans); 
00711 
00712         int const n = traits::matrix_size1 (c);
00713         assert (n == traits::matrix_size2 (c)); 
00714         
00715         int const k = trans == CblasNoTrans
00716           ? traits::matrix_size2 (a)
00717           : traits::matrix_size1 (a); 
00718         assert (n == (trans == CblasNoTrans
00719                       ? traits::matrix_size1 (a)
00720                       : traits::matrix_size2 (a))); 
00721 
00722         CBLAS_ORDER const stor_ord
00723           = enum_cast<CBLAS_ORDER const>
00724           (storage_order<
00725 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00726             typename traits::matrix_traits<SymmC>::ordering_type
00727 #else
00728             typename SymmC::orientation_category 
00729 #endif 
00730            >::value); 
00731 
00732         syrk (stor_ord, uplo, trans, 
00733               n, k, alpha, 
00734 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
00735               traits::matrix_storage (a), 
00736 #else
00737               traits::matrix_storage_const (a), 
00738 #endif
00739               traits::leading_dimension (a),
00740               beta, 
00741               traits::matrix_storage (c), 
00742               traits::leading_dimension (c)); 
00743       }
00744 
00745     } // detail 
00746  
00747     // C <- alpha * A * A^T + beta * C
00748     // C <- alpha * A^T * A + beta * C
00749     // C == C^T
00750     template <typename T, typename MatrA, typename SymmC>
00751     inline
00752     void syrk (CBLAS_UPLO const uplo, CBLAS_TRANSPOSE trans, 
00753                T const& alpha, MatrA const& a, 
00754                T const& beta, SymmC& c)
00755     {
00756 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00757       BOOST_STATIC_ASSERT((boost::is_same<
00758         typename traits::matrix_traits<SymmC>::matrix_structure, 
00759         traits::general_t
00760       >::value)); 
00761 #endif 
00762 
00763       detail::syrk (uplo, trans, alpha, a, beta, c); 
00764     }
00765 
00766     template <typename T, typename MatrA, typename SymmC>
00767     inline
00768     void syrk (CBLAS_TRANSPOSE trans, 
00769                T const& alpha, MatrA const& a, 
00770                T const& beta, SymmC& c)
00771     {
00772 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00773       BOOST_STATIC_ASSERT((boost::is_same<
00774         typename traits::matrix_traits<SymmC>::matrix_structure, 
00775         traits::symmetric_t
00776       >::value)); 
00777 #endif 
00778 
00779       CBLAS_UPLO const uplo
00780         = enum_cast<CBLAS_UPLO const>
00781         (uplo_triang<
00782 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00783           typename traits::matrix_traits<SymmC>::uplo_type
00784 #else
00785           typename SymmC::packed_category 
00786 #endif 
00787          >::value); 
00788 
00789       detail::syrk (uplo, trans, alpha, a, beta, c); 
00790     }
00791 
00792     // C <- A * A^T + C
00793     // C <- A^T * A + C
00794     template <typename MatrA, typename SymmC>
00795     inline
00796     void syrk (CBLAS_TRANSPOSE trans, MatrA const& a, SymmC& c) {
00797 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00798       typedef typename traits::matrix_traits<SymmC>::value_type val_t; 
00799 #else
00800       typedef typename SymmC::value_type val_t; 
00801 #endif 
00802       syrk (trans, (val_t) 1, a, (val_t) 0, c);
00803     }
00804 
00805     
00806     // C <- alpha * A * B^T + conj(alpha) * B * A^T + beta * C
00807     // C <- alpha * A^T * B + conj(alpha) * B^T * A + beta * C
00808     // C == C^T
00809 
00810     namespace detail {
00811 
00812       template <typename T, typename MatrA, typename MatrB, typename SymmC>
00813       inline
00814       void syr2k (CBLAS_UPLO const uplo, CBLAS_TRANSPOSE trans, 
00815                   T const& alpha, MatrA const& a, MatrB const& b, 
00816                   T const& beta, SymmC& c)
00817       {
00818 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00819         BOOST_STATIC_ASSERT((boost::is_same<
00820           typename traits::matrix_traits<MatrA>::matrix_structure, 
00821           traits::general_t
00822         >::value));
00823         BOOST_STATIC_ASSERT((boost::is_same<
00824           typename traits::matrix_traits<MatrB>::matrix_structure, 
00825           traits::general_t
00826         >::value));
00827         BOOST_STATIC_ASSERT((boost::is_same<
00828           typename traits::matrix_traits<MatrA>::ordering_type,
00829           typename traits::matrix_traits<SymmC>::ordering_type
00830         >::value)); 
00831         BOOST_STATIC_ASSERT((boost::is_same<
00832           typename traits::matrix_traits<MatrB>::ordering_type,
00833           typename traits::matrix_traits<SymmC>::ordering_type
00834         >::value)); 
00835 #endif 
00836 
00837         assert (uplo == CblasUpper || uplo == CblasLower); 
00838         assert (trans == CblasNoTrans 
00839                 || trans == CblasTrans 
00840                 || trans == CblasConjTrans); 
00841 
00842         int const n = traits::matrix_size1 (c);
00843         assert (n == traits::matrix_size2 (c)); 
00844         
00845         int const k = trans == CblasNoTrans
00846           ? traits::matrix_size2 (a)
00847           : traits::matrix_size1 (a); 
00848         assert (k == (trans == CblasNoTrans
00849                       ? traits::matrix_size2 (b)
00850                       : traits::matrix_size1 (b))); 
00851         assert (n == (trans == CblasNoTrans
00852                       ? traits::matrix_size1 (a)
00853                       : traits::matrix_size2 (a))); 
00854         assert (n == (trans == CblasNoTrans
00855                       ? traits::matrix_size1 (b)
00856                       : traits::matrix_size2 (b))); 
00857 
00858         CBLAS_ORDER const stor_ord
00859           = enum_cast<CBLAS_ORDER const>
00860           (storage_order<
00861 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00862             typename traits::matrix_traits<SymmC>::ordering_type
00863 #else
00864             typename SymmC::orientation_category 
00865 #endif 
00866            >::value); 
00867 
00868         syr2k (stor_ord, uplo, trans, 
00869                n, k, alpha, 
00870 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
00871                traits::matrix_storage (a), 
00872 #else
00873                traits::matrix_storage_const (a), 
00874 #endif
00875                traits::leading_dimension (a),
00876 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
00877                traits::matrix_storage (b), 
00878 #else
00879                traits::matrix_storage_const (b), 
00880 #endif
00881                traits::leading_dimension (b),
00882                beta, 
00883                traits::matrix_storage (c), 
00884                traits::leading_dimension (c)); 
00885       }
00886 
00887     } // detail 
00888  
00889     // C <- alpha * A * B^T + conj(alpha) * B * A^T + beta * C
00890     // C <- alpha * A^T * B + conj(alpha) * B^T * A + beta * C
00891     // C == C^T
00892     template <typename T, typename MatrA, typename MatrB, typename SymmC>
00893     inline
00894     void syr2k (CBLAS_UPLO const uplo, CBLAS_TRANSPOSE trans, 
00895                T const& alpha, MatrA const& a, MatrB const& b, 
00896                T const& beta, SymmC& c)
00897     {
00898 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00899       BOOST_STATIC_ASSERT((boost::is_same<
00900         typename traits::matrix_traits<SymmC>::matrix_structure, 
00901         traits::general_t
00902       >::value)); 
00903 #endif 
00904 
00905       detail::syr2k (uplo, trans, alpha, a, b, beta, c); 
00906     }
00907 
00908     template <typename T, typename MatrA, typename MatrB, typename SymmC>
00909     inline
00910     void syr2k (CBLAS_TRANSPOSE trans, 
00911                T const& alpha, MatrA const& a, MatrB const& b,  
00912                T const& beta, SymmC& c)
00913     {
00914 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00915       BOOST_STATIC_ASSERT((boost::is_same<
00916         typename traits::matrix_traits<SymmC>::matrix_structure, 
00917         traits::symmetric_t
00918       >::value)); 
00919 #endif 
00920 
00921       CBLAS_UPLO const uplo
00922         = enum_cast<CBLAS_UPLO const>
00923         (uplo_triang<
00924 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00925           typename traits::matrix_traits<SymmC>::uplo_type
00926 #else
00927           typename SymmC::packed_category 
00928 #endif 
00929          >::value); 
00930 
00931       detail::syr2k (uplo, trans, alpha, a, b, beta, c); 
00932     }
00933 
00934     // C <- A * B^T + B * A^T + C
00935     // C <- A^T * B + B^T * A + C
00936     template <typename MatrA, typename MatrB, typename SymmC>
00937     inline
00938     void syr2k (CBLAS_TRANSPOSE trans, 
00939                 MatrA const& a, MatrB const& b, SymmC& c) 
00940     {
00941 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00942       typedef typename traits::matrix_traits<SymmC>::value_type val_t; 
00943 #else
00944       typedef typename SymmC::value_type val_t; 
00945 #endif 
00946       syr2k (trans, (val_t) 1, a, b, (val_t) 0, c);
00947     }
00948 
00949     
00950     // C <- alpha * A * A^H + beta * C
00951     // C <- alpha * A^H * A + beta * C
00952     // C == C^H
00953 
00954     namespace detail {
00955 
00956       template <typename T, typename MatrA, typename HermC>
00957       inline
00958       void herk (CBLAS_UPLO const uplo, CBLAS_TRANSPOSE trans, 
00959                  T const& alpha, MatrA const& a, 
00960                  T const& beta, HermC& c)
00961       {
00962 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
00963         BOOST_STATIC_ASSERT((boost::is_same<
00964           typename traits::matrix_traits<MatrA>::matrix_structure, 
00965           traits::general_t
00966         >::value));
00967         BOOST_STATIC_ASSERT((boost::is_same<
00968           typename traits::matrix_traits<MatrA>::ordering_type,
00969           typename traits::matrix_traits<HermC>::ordering_type
00970         >::value)); 
00971 #endif 
00972 
00973         assert (uplo == CblasUpper || uplo == CblasLower); 
00974         assert (trans == CblasNoTrans || trans == CblasConjTrans); 
00975 
00976         int const n = traits::matrix_size1 (c);
00977         assert (n == traits::matrix_size2 (c)); 
00978         
00979         int const k = trans == CblasNoTrans
00980           ? traits::matrix_size2 (a)
00981           : traits::matrix_size1 (a); 
00982         assert (n == (trans == CblasNoTrans
00983                       ? traits::matrix_size1 (a)
00984                       : traits::matrix_size2 (a))); 
00985 
00986         CBLAS_ORDER const stor_ord
00987           = enum_cast<CBLAS_ORDER const>
00988           (storage_order<
00989 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00990             typename traits::matrix_traits<HermC>::ordering_type
00991 #else
00992             typename HermC::orientation_category 
00993 #endif 
00994            >::value); 
00995 
00996         herk (stor_ord, uplo, trans, 
00997               n, k, alpha, 
00998 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
00999               traits::matrix_storage (a), 
01000 #else
01001               traits::matrix_storage_const (a), 
01002 #endif
01003               traits::leading_dimension (a),
01004               beta, 
01005               traits::matrix_storage (c), 
01006               traits::leading_dimension (c)); 
01007       }
01008 
01009     } // detail 
01010  
01011     // C <- alpha * A * A^H + beta * C
01012     // C <- alpha * A^H * A + beta * C
01013     // C == C^H
01014     template <typename T, typename MatrA, typename HermC>
01015     inline
01016     void herk (CBLAS_UPLO const uplo, CBLAS_TRANSPOSE trans, 
01017                T const& alpha, MatrA const& a, 
01018                T const& beta, HermC& c)
01019     {
01020 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
01021       BOOST_STATIC_ASSERT((boost::is_same<
01022         typename traits::matrix_traits<HermC>::matrix_structure, 
01023         traits::general_t
01024       >::value)); 
01025 #endif 
01026 
01027       detail::herk (uplo, trans, alpha, a, beta, c); 
01028     }
01029 
01030     template <typename T, typename MatrA, typename HermC>
01031     inline
01032     void herk (CBLAS_TRANSPOSE trans, 
01033                T const& alpha, MatrA const& a, 
01034                T const& beta, HermC& c)
01035     {
01036 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
01037       BOOST_STATIC_ASSERT((boost::is_same<
01038         typename traits::matrix_traits<HermC>::matrix_structure, 
01039         traits::hermitian_t
01040       >::value)); 
01041 #endif 
01042 
01043       CBLAS_UPLO const uplo
01044         = enum_cast<CBLAS_UPLO const>
01045         (uplo_triang<
01046 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
01047           typename traits::matrix_traits<HermC>::uplo_type
01048 #else
01049           typename HermC::packed_category 
01050 #endif 
01051          >::value); 
01052 
01053       detail::herk (uplo, trans, alpha, a, beta, c); 
01054     }
01055 
01056     // C <- A * A^H + C
01057     // C <- A^H * A + C
01058     template <typename MatrA, typename HermC>
01059     inline
01060     void herk (CBLAS_TRANSPOSE trans, MatrA const& a, HermC& c) {
01061 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
01062       typedef typename traits::matrix_traits<HermC>::value_type val_t; 
01063 #else
01064       typedef typename HermC::value_type val_t; 
01065 #endif 
01066       typedef typename traits::type_traits<val_t>::real_type real_t; 
01067       herk (trans, (real_t) 1, a, (real_t) 0, c);
01068     }
01069 
01070 
01071     // C <- alpha * A * B^H + conj(alpha) * B * A^H + beta * C
01072     // C <- alpha * A^H * B + conj(alpha) * B^H * A + beta * C
01073     // C == C^H
01074 
01075     namespace detail {
01076 
01077       template <typename T1, typename T2, 
01078                 typename MatrA, typename MatrB, typename HermC>
01079       inline
01080       void her2k (CBLAS_UPLO const uplo, CBLAS_TRANSPOSE trans, 
01081                   T1 const& alpha, MatrA const& a, MatrB const& b, 
01082                   T2 const& beta, HermC& c)
01083       {
01084 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
01085         BOOST_STATIC_ASSERT((boost::is_same<
01086           typename traits::matrix_traits<MatrA>::matrix_structure, 
01087           traits::general_t
01088         >::value));
01089         BOOST_STATIC_ASSERT((boost::is_same<
01090           typename traits::matrix_traits<MatrB>::matrix_structure, 
01091           traits::general_t
01092         >::value));
01093         BOOST_STATIC_ASSERT((boost::is_same<
01094           typename traits::matrix_traits<MatrA>::ordering_type,
01095           typename traits::matrix_traits<HermC>::ordering_type
01096         >::value)); 
01097         BOOST_STATIC_ASSERT((boost::is_same<
01098           typename traits::matrix_traits<MatrB>::ordering_type,
01099           typename traits::matrix_traits<HermC>::ordering_type
01100         >::value)); 
01101 #endif 
01102 
01103         assert (uplo == CblasUpper || uplo == CblasLower); 
01104         assert (trans == CblasNoTrans || trans == CblasConjTrans); 
01105 
01106         int const n = traits::matrix_size1 (c);
01107         assert (n == traits::matrix_size2 (c)); 
01108         
01109         int const k = trans == CblasNoTrans
01110           ? traits::matrix_size2 (a)
01111           : traits::matrix_size1 (a); 
01112         assert (k == (trans == CblasNoTrans
01113                       ? traits::matrix_size2 (b)
01114                       : traits::matrix_size1 (b))); 
01115         assert (n == (trans == CblasNoTrans
01116                       ? traits::matrix_size1 (a)
01117                       : traits::matrix_size2 (a))); 
01118         assert (n == (trans == CblasNoTrans
01119                       ? traits::matrix_size1 (b)
01120                       : traits::matrix_size2 (b))); 
01121 
01122         CBLAS_ORDER const stor_ord
01123           = enum_cast<CBLAS_ORDER const>
01124           (storage_order<
01125 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
01126             typename traits::matrix_traits<HermC>::ordering_type
01127 #else
01128             typename HermC::orientation_category 
01129 #endif 
01130            >::value); 
01131 
01132         her2k (stor_ord, uplo, trans, 
01133                n, k, alpha, 
01134 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
01135                traits::matrix_storage (a), 
01136 #else
01137                traits::matrix_storage_const (a), 
01138 #endif
01139                traits::leading_dimension (a),
01140 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
01141                traits::matrix_storage (b), 
01142 #else
01143                traits::matrix_storage_const (b), 
01144 #endif
01145                traits::leading_dimension (b),
01146                beta, 
01147                traits::matrix_storage (c), 
01148                traits::leading_dimension (c)); 
01149       }
01150 
01151     } // detail 
01152  
01153     // C <- alpha * A * B^H + conj(alpha) * B * A^H + beta * C
01154     // C <- alpha * A^H * B + conj(alpha) * B^H * A + beta * C
01155     // C == C^H
01156     template <typename T1, typename T2, 
01157               typename MatrA, typename MatrB, typename HermC>
01158     inline
01159     void her2k (CBLAS_UPLO const uplo, CBLAS_TRANSPOSE trans, 
01160                T1 const& alpha, MatrA const& a, MatrB const& b, 
01161                T2 const& beta, HermC& c)
01162     {
01163 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
01164       BOOST_STATIC_ASSERT((boost::is_same<
01165         typename traits::matrix_traits<HermC>::matrix_structure, 
01166         traits::general_t
01167       >::value)); 
01168 #endif 
01169 
01170       detail::her2k (uplo, trans, alpha, a, b, beta, c); 
01171     }
01172 
01173     template <typename T1, typename T2, 
01174               typename MatrA, typename MatrB, typename HermC>
01175     inline
01176     void her2k (CBLAS_TRANSPOSE trans, 
01177                T1 const& alpha, MatrA const& a, MatrB const& b,  
01178                T2 const& beta, HermC& c)
01179     {
01180 #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
01181       BOOST_STATIC_ASSERT((boost::is_same<
01182         typename traits::matrix_traits<HermC>::matrix_structure, 
01183         traits::hermitian_t
01184       >::value)); 
01185 #endif 
01186 
01187       CBLAS_UPLO const uplo
01188         = enum_cast<CBLAS_UPLO const>
01189         (uplo_triang<
01190 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
01191           typename traits::matrix_traits<HermC>::uplo_type
01192 #else
01193           typename HermC::packed_category 
01194 #endif 
01195          >::value); 
01196 
01197       detail::her2k (uplo, trans, alpha, a, b, beta, c); 
01198     }
01199 
01200     // C <- A * B^H + B * A^H + C
01201     // C <- A^H * B + B^H * A + C
01202     template <typename MatrA, typename MatrB, typename HermC>
01203     inline
01204     void her2k (CBLAS_TRANSPOSE trans, 
01205                 MatrA const& a, MatrB const& b, HermC& c) 
01206     {
01207 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
01208       typedef typename traits::matrix_traits<HermC>::value_type val_t; 
01209 #else
01210       typedef typename HermC::value_type val_t; 
01211 #endif 
01212       typedef typename traits::type_traits<val_t>::real_type real_t; 
01213       her2k (trans, (val_t) 1, a, b, (real_t) 0, c);
01214     }
01215 
01216     
01217   } // namespace atlas
01218 
01219 }}} 
01220 
01221 #endif // BOOST_NUMERIC_BINDINGS_CBLAS_LEVEL_3_HPP