← Back to branch summary

~ubuntu-branches/ubuntu/maverick/freecad/maverick

« back to all changes in this revision

Viewing changes to src/3rdParty/boost/numeric/bindings/atlas/cblas3.hpp

  • Committer: Bazaar Package Importer
  • Author(s): Teemu Ikonen
  • Date: 2009-07-16 18:37:41 UTC
  • Revision ID: james.westby@ubuntu.com-20090716183741-oww9kcxqrk991i1n
Tags: upstream-0.8.2237
Import upstream version 0.8.2237

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/3rdParty/boost/numeric/bindings/atlas/cblas3.hpp
 
1
/*
 
2
 * 
 
3
 * Copyright (c) Kresimir Fresl 2002 
 
4
 *
 
5
 * Distributed under the Boost Software License, Version 1.0.
 
6
 * (See accompanying file LICENSE_1_0.txt or copy at
 
7
 * http://www.boost.org/LICENSE_1_0.txt)
 
8
 *
 
9
 * Author acknowledges the support of the Faculty of Civil Engineering, 
 
10
 * University of Zagreb, Croatia.
 
11
 *
 
12
 */
 
13
 
 
14
#ifndef BOOST_NUMERIC_BINDINGS_CBLAS_LEVEL_3_HPP
 
15
#define BOOST_NUMERIC_BINDINGS_CBLAS_LEVEL_3_HPP
 
16
 
 
17
#include <cassert>
 
18
 
 
19
#include <boost/numeric/bindings/traits/traits.hpp>
 
20
#include <boost/numeric/bindings/traits/type_traits.hpp>
 
21
#include <boost/numeric/bindings/atlas/cblas3_overloads.hpp>
 
22
#include <boost/numeric/bindings/atlas/cblas_enum.hpp>
 
23
#include <boost/type_traits/same_traits.hpp>
 
24
#include <boost/mpl/if.hpp>
 
25
 
 
26
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
 
27
#  include <boost/static_assert.hpp>
 
28
#endif
 
29
 
 
30
namespace boost { namespace numeric { namespace bindings { 
 
31
 
 
32
  namespace atlas {
 
33
 
 
34
    // C <- alpha * op (A) * op (B) + beta * C 
 
35
    // op (A) == A || A^T || A^H
 
36
    template <typename T, typename MatrA, typename MatrB, typename MatrC>
 
37
    inline
 
38
    void gemm (CBLAS_TRANSPOSE const TransA, CBLAS_TRANSPOSE const TransB, 
 
39
               T const& alpha, MatrA const& a, MatrB const& b, 
 
40
               T const& beta, MatrC& c
 
41
               )
 
42
    {
 
43
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
 
44
      BOOST_STATIC_ASSERT((boost::is_same<
 
45
        typename traits::matrix_traits<MatrA>::matrix_structure, 
 
46
        traits::general_t
 
47
      >::value)); 
 
48
      BOOST_STATIC_ASSERT((boost::is_same<
 
49
        typename traits::matrix_traits<MatrB>::matrix_structure, 
 
50
        traits::general_t
 
51
      >::value)); 
 
52
      BOOST_STATIC_ASSERT((boost::is_same<
 
53
        typename traits::matrix_traits<MatrC>::matrix_structure, 
 
54
        traits::general_t
 
55
      >::value)); 
 
56
 
 
57
      BOOST_STATIC_ASSERT((boost::is_same<
 
58
        typename traits::matrix_traits<MatrA>::ordering_type,
 
59
        typename traits::matrix_traits<MatrB>::ordering_type
 
60
      >::value)); 
 
61
      BOOST_STATIC_ASSERT((boost::is_same<
 
62
        typename traits::matrix_traits<MatrA>::ordering_type,
 
63
        typename traits::matrix_traits<MatrC>::ordering_type
 
64
      >::value)); 
 
65
#endif 
 
66
 
 
67
      assert (TransA == CblasNoTrans 
 
68
              || TransA == CblasTrans 
 
69
              || TransA == CblasConjTrans); 
 
70
      assert (TransB == CblasNoTrans 
 
71
              || TransB == CblasTrans 
 
72
              || TransB == CblasConjTrans); 
 
73
 
 
74
      int const m = TransA == CblasNoTrans
 
75
        ? traits::matrix_size1 (a)
 
76
        : traits::matrix_size2 (a);
 
77
      int const n = TransB == CblasNoTrans
 
78
        ? traits::matrix_size2 (b)
 
79
        : traits::matrix_size1 (b);
 
80
      int const k = TransA == CblasNoTrans
 
81
        ? traits::matrix_size2 (a)
 
82
        : traits::matrix_size1 (a); 
 
83
      assert (m == traits::matrix_size1 (c)); 
 
84
      assert (n == traits::matrix_size2 (c)); 
 
85
#ifndef NDEBUG
 
86
      int const k1 = TransB == CblasNoTrans
 
87
        ? traits::matrix_size1 (b)
 
88
        : traits::matrix_size2 (b);
 
89
      assert (k == k1); 
 
90
#endif
 
91
      // .. what about AtlasConj? 
 
92
 
 
93
      CBLAS_ORDER const stor_ord
 
94
        = enum_cast<CBLAS_ORDER const>
 
95
        (storage_order<
 
96
#ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
 
97
           typename traits::matrix_traits<MatrA>::ordering_type
 
98
#else
 
99
           typename MatrA::orientation_category 
 
100
#endif 
 
101
         >::value); 
 
102
 
 
103
      detail::gemm (stor_ord, TransA, TransB, m, n, k, alpha, 
 
104
#ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
 
105
                    traits::matrix_storage (a), 
 
106
#else
 
107
                    traits::matrix_storage_const (a), 
 
108
#endif
 
109
                    traits::leading_dimension (a),
 
110
#ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
 
111
                    traits::matrix_storage (b), 
 
112
#else
 
113
                    traits::matrix_storage_const (b), 
 
114
#endif
 
115
                    traits::leading_dimension (b),
 
116
                    beta, 
 
117
                    traits::matrix_storage (c), 
 
118
                    traits::leading_dimension (c)); 
 
119
    }
 
120
 
 
121
 
 
122
    // C <- alpha * A * B + beta * C 
 
123
    template <typename T, typename MatrA, typename MatrB, typename MatrC>
 
124
    inline
 
125
    void gemm (T const& alpha, MatrA const& a, MatrB const& b, 
 
126
               T const& beta, MatrC& c) 
 
127
    {
 
128
      gemm (CblasNoTrans, CblasNoTrans, alpha, a, b, beta, c) ;
 
129
    }
 
130
    
 
131
 
 
132
    // C <- A * B 
 
133
    template <typename MatrA, typename MatrB, typename MatrC>
 
134
    inline
 
135
    void gemm (MatrA const& a, MatrB const& b, MatrC& c) {
 
136
#ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
 
137
      typedef typename traits::matrix_traits<MatrC>::value_type val_t; 
 
138
#else
 
139
      typedef typename MatrC::value_type val_t; 
 
140
#endif 
 
141
      gemm (CblasNoTrans, CblasNoTrans, (val_t) 1, a, b, (val_t) 0, c);
 
142
    }
 
143
 
 
144
 
 
145
 
 
146
    // C <- alpha * A * B + beta * C 
 
147
    // C <- alpha * B * A + beta * C 
 
148
    // A == A^T
 
149
 
 
150
    namespace detail {
 
151
 
 
152
      template <typename T, typename SymmA, typename MatrB, typename MatrC>
 
153
      inline
 
154
      void symm (CBLAS_SIDE const side, CBLAS_UPLO const uplo, 
 
155
                 T const& alpha, SymmA const& a, MatrB const& b, 
 
156
                 T const& beta, MatrC& c)
 
157
      {
 
158
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
 
159
        BOOST_STATIC_ASSERT((boost::is_same<
 
160
          typename traits::matrix_traits<MatrB>::matrix_structure, 
 
161
          traits::general_t
 
162
        >::value));
 
163
        BOOST_STATIC_ASSERT((boost::is_same<
 
164
          typename traits::matrix_traits<MatrC>::matrix_structure, 
 
165
          traits::general_t
 
166
        >::value)); 
 
167
 
 
168
        BOOST_STATIC_ASSERT((boost::is_same<
 
169
          typename traits::matrix_traits<SymmA>::ordering_type,
 
170
          typename traits::matrix_traits<MatrB>::ordering_type
 
171
        >::value)); 
 
172
        BOOST_STATIC_ASSERT((boost::is_same<
 
173
          typename traits::matrix_traits<SymmA>::ordering_type,
 
174
          typename traits::matrix_traits<MatrC>::ordering_type
 
175
        >::value)); 
 
176
#endif 
 
177
 
 
178
        assert (side == CblasLeft || side == CblasRight);
 
179
        assert (uplo == CblasUpper || uplo == CblasLower); 
 
180
 
 
181
        int const m = traits::matrix_size1 (c);
 
182
        int const n = traits::matrix_size2 (c);
 
183
 
 
184
        assert (side == CblasLeft 
 
185
                ? m == traits::matrix_size1 (a) 
 
186
                  && m == traits::matrix_size2 (a)
 
187
                : n == traits::matrix_size1 (a) 
 
188
                  && n == traits::matrix_size2 (a)); 
 
189
        assert (m == traits::matrix_size1 (b) 
 
190
                && n == traits::matrix_size2 (b)); 
 
191
 
 
192
        CBLAS_ORDER const stor_ord
 
193
          = enum_cast<CBLAS_ORDER const>
 
194
          (storage_order<
 
195
#ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
 
196
            typename traits::matrix_traits<SymmA>::ordering_type
 
197
#else
 
198
            typename SymmA::orientation_category 
 
199
#endif 
 
200
           >::value); 
 
201
 
 
202
        symm (stor_ord, side, uplo,  
 
203
              m, n, alpha, 
 
204
#ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
 
205
              traits::matrix_storage (a), 
 
206
#else
 
207
              traits::matrix_storage_const (a), 
 
208
#endif
 
209
              traits::leading_dimension (a),
 
210
#ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
 
211
              traits::matrix_storage (b), 
 
212
#else
 
213
              traits::matrix_storage_const (b), 
 
214
#endif
 
215
              traits::leading_dimension (b),
 
216
              beta, 
 
217
              traits::matrix_storage (c), 
 
218
              traits::leading_dimension (c)); 
 
219
      }
 
220
 
 
221
    } // detail 
 
222
 
 
223
    // C <- alpha * A * B + beta * C 
 
224
    // C <- alpha * B * A + beta * C 
 
225
    // A == A^T
 
226
    template <typename T, typename SymmA, typename MatrB, typename MatrC>
 
227
    inline
 
228
    void symm (CBLAS_SIDE const side, CBLAS_UPLO const uplo, 
 
229
               T const& alpha, SymmA const& a, MatrB const& b, 
 
230
               T const& beta, MatrC& c)
 
231
    {
 
232
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
 
233
      BOOST_STATIC_ASSERT((boost::is_same<
 
234
        typename traits::matrix_traits<SymmA>::matrix_structure, 
 
235
        traits::general_t
 
236
      >::value)); 
 
237
#endif 
 
238
 
 
239
      detail::symm (side, uplo, alpha, a, b, beta, c); 
 
240
    }
 
241
 
 
242
    template <typename T, typename SymmA, typename MatrB, typename MatrC>
 
243
    inline
 
244
    void symm (CBLAS_SIDE const side, 
 
245
               T const& alpha, SymmA const& a, MatrB const& b, 
 
246
               T const& beta, MatrC& c)
 
247
    {
 
248
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
 
249
      BOOST_STATIC_ASSERT((boost::is_same<
 
250
        typename traits::matrix_traits<SymmA>::matrix_structure, 
 
251
        traits::symmetric_t
 
252
      >::value)); 
 
253
#endif 
 
254
 
 
255
      CBLAS_UPLO const uplo
 
256
        = enum_cast<CBLAS_UPLO const>
 
257
        (uplo_triang<
 
258
#ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
 
259
          typename traits::matrix_traits<SymmA>::uplo_type
 
260
#else
 
261
          typename SymmA::packed_category 
 
262
#endif 
 
263
         >::value); 
 
264
 
 
265
      detail::symm (side, uplo, alpha, a, b, beta, c); 
 
266
    }
 
267
 
 
268
 
 
269
#ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
 
270
 
 
271
    namespace detail {
 
272
 
 
273
      // C <- alpha * A * B + beta * C ;  A == A^T
 
274
      struct symm_left {
 
275
        template <typename T, typename SymmA, typename MatrB, typename MatrC>
 
276
        static void f (T const& alpha, SymmA const& a, MatrB const& b, 
 
277
                       T const& beta, MatrC& c) 
 
278
        {
 
279
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
 
280
          BOOST_STATIC_ASSERT((boost::is_same<
 
281
            typename traits::matrix_traits<SymmA>::matrix_structure, 
 
282
            traits::symmetric_t
 
283
          >::value)); 
 
284
          BOOST_STATIC_ASSERT((boost::is_same<
 
285
            typename traits::matrix_traits<MatrB>::matrix_structure, 
 
286
            traits::general_t
 
287
          >::value));
 
288
#endif 
 
289
 
 
290
          int const m = traits::matrix_size1 (c);
 
291
          int const n = traits::matrix_size2 (c);
 
292
 
 
293
          assert (m == traits::matrix_size1 (a) 
 
294
                  && m == traits::matrix_size2 (a)); 
 
295
          assert (m == traits::matrix_size1 (b) 
 
296
                  && n == traits::matrix_size2 (b)); 
 
297
 
 
298
          CBLAS_ORDER const stor_ord
 
299
            = enum_cast<CBLAS_ORDER const>
 
300
            (storage_order<
 
301
              typename traits::matrix_traits<SymmA>::ordering_type
 
302
             >::value); 
 
303
 
 
304
          CBLAS_UPLO const uplo
 
305
            = enum_cast<CBLAS_UPLO const>
 
306
            (uplo_triang<
 
307
              typename traits::matrix_traits<SymmA>::uplo_type
 
308
             >::value); 
 
309
 
 
310
          symm (stor_ord, CblasLeft, uplo,  
 
311
                m, n, alpha, 
 
312
                traits::matrix_storage (a), traits::leading_dimension (a),
 
313
                traits::matrix_storage (b), traits::leading_dimension (b),
 
314
                beta, 
 
315
                traits::matrix_storage (c), traits::leading_dimension (c)); 
 
316
        }
 
317
      }; 
 
318
 
 
319
      // C <- alpha * A * B + beta * C ;  B == B^T
 
320
      struct symm_right {
 
321
        template <typename T, typename MatrA, typename SymmB, typename MatrC>
 
322
        static void f (T const& alpha, MatrA const& a, SymmB const& b, 
 
323
                       T const& beta, MatrC& c) 
 
324
        {
 
325
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
 
326
          BOOST_STATIC_ASSERT((boost::is_same<
 
327
            typename traits::matrix_traits<MatrA>::matrix_structure, 
 
328
            traits::general_t
 
329
          >::value));
 
330
          BOOST_STATIC_ASSERT((boost::is_same<
 
331
            typename traits::matrix_traits<SymmB>::matrix_structure, 
 
332
            traits::symmetric_t
 
333
          >::value));
 
334
#endif 
 
335
 
 
336
          int const m = traits::matrix_size1 (c);
 
337
          int const n = traits::matrix_size2 (c);
 
338
 
 
339
          assert (n == traits::matrix_size1 (b) 
 
340
                  && n == traits::matrix_size2 (b)); 
 
341
          assert (m == traits::matrix_size1 (a) 
 
342
                  && n == traits::matrix_size2 (a)); 
 
343
 
 
344
          CBLAS_ORDER const stor_ord
 
345
            = enum_cast<CBLAS_ORDER const>
 
346
            (storage_order<
 
347
              typename traits::matrix_traits<SymmB>::ordering_type
 
348
             >::value); 
 
349
 
 
350
          CBLAS_UPLO const uplo
 
351
            = enum_cast<CBLAS_UPLO const>
 
352
            (uplo_triang<
 
353
              typename traits::matrix_traits<SymmB>::uplo_type
 
354
             >::value); 
 
355
 
 
356
          symm (stor_ord, CblasRight, uplo,  
 
357
                m, n, alpha, 
 
358
                traits::matrix_storage (b), traits::leading_dimension (b),
 
359
                traits::matrix_storage (a), traits::leading_dimension (a),
 
360
                beta, 
 
361
                traits::matrix_storage (c), traits::leading_dimension (c)); 
 
362
        }
 
363
      }; 
 
364
 
 
365
    } // detail 
 
366
    
 
367
    // C <- alpha * A * B + beta * C 
 
368
    // C <- alpha * B * A + beta * C 
 
369
    // A == A^T
 
370
    template <typename T, typename MatrA, typename MatrB, typename MatrC>
 
371
    inline
 
372
    void symm (T const& alpha, MatrA const& a, MatrB const& b, 
 
373
               T const& beta, MatrC& c)
 
374
    {
 
375
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
 
376
      BOOST_STATIC_ASSERT((boost::is_same<
 
377
        typename traits::matrix_traits<MatrC>::matrix_structure, 
 
378
        traits::general_t
 
379
      >::value)); 
 
380
 
 
381
      BOOST_STATIC_ASSERT((boost::is_same<
 
382
        typename traits::matrix_traits<MatrA>::ordering_type,
 
383
        typename traits::matrix_traits<MatrB>::ordering_type
 
384
      >::value)); 
 
385
      BOOST_STATIC_ASSERT((boost::is_same<
 
386
        typename traits::matrix_traits<MatrA>::ordering_type,
 
387
        typename traits::matrix_traits<MatrC>::ordering_type
 
388
      >::value)); 
 
389
#endif 
 
390
 
 
391
      typedef typename
 
392
        boost::mpl::if_c<
 
393
          boost::is_same<
 
394
            typename traits::matrix_traits<MatrA>::matrix_structure, 
 
395
            traits::symmetric_t
 
396
          >::value,
 
397
          detail::symm_left, 
 
398
          detail::symm_right
 
399
        >::type functor; 
 
400
 
 
401
      functor::f (alpha, a, b, beta, c); 
 
402
    }
 
403
 
 
404
    // C <- A * B  
 
405
    // C <- B * A  
 
406
    template <typename MatrA, typename MatrB, typename MatrC>
 
407
    inline
 
408
    void symm (MatrA const& a, MatrB const& b, MatrC& c) {
 
409
      typedef typename traits::matrix_traits<MatrC>::value_type val_t; 
 
410
      symm ((val_t) 1, a, b, (val_t) 0, c);
 
411
    }
 
412
 
 
413
#endif // BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
 
414
 
 
415
 
 
416
 
 
417
    // C <- alpha * A * B + beta * C 
 
418
    // C <- alpha * B * A + beta * C 
 
419
    // A == A^H
 
420
 
 
421
    namespace detail {
 
422
 
 
423
      template <typename T, typename HermA, typename MatrB, typename MatrC>
 
424
      inline
 
425
      void hemm (CBLAS_SIDE const side, CBLAS_UPLO const uplo, 
 
426
                 T const& alpha, HermA const& a, MatrB const& b, 
 
427
                 T const& beta, MatrC& c)
 
428
      {
 
429
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
 
430
        BOOST_STATIC_ASSERT((boost::is_same<
 
431
          typename traits::matrix_traits<MatrB>::matrix_structure, 
 
432
          traits::general_t
 
433
        >::value));
 
434
        BOOST_STATIC_ASSERT((boost::is_same<
 
435
          typename traits::matrix_traits<MatrC>::matrix_structure, 
 
436
          traits::general_t
 
437
        >::value)); 
 
438
 
 
439
        BOOST_STATIC_ASSERT((boost::is_same<
 
440
          typename traits::matrix_traits<HermA>::ordering_type,
 
441
          typename traits::matrix_traits<MatrB>::ordering_type
 
442
        >::value)); 
 
443
        BOOST_STATIC_ASSERT((boost::is_same<
 
444
          typename traits::matrix_traits<HermA>::ordering_type,
 
445
          typename traits::matrix_traits<MatrC>::ordering_type
 
446
        >::value)); 
 
447
#endif 
 
448
 
 
449
        assert (side == CblasLeft || side == CblasRight);
 
450
        assert (uplo == CblasUpper || uplo == CblasLower); 
 
451
 
 
452
        int const m = traits::matrix_size1 (c);
 
453
        int const n = traits::matrix_size2 (c);
 
454
 
 
455
        assert (side == CblasLeft 
 
456
                ? m == traits::matrix_size1 (a) 
 
457
                  && m == traits::matrix_size2 (a)
 
458
                : n == traits::matrix_size1 (a) 
 
459
                  && n == traits::matrix_size2 (a)); 
 
460
        assert (m == traits::matrix_size1 (b) 
 
461
                && n == traits::matrix_size2 (b)); 
 
462
 
 
463
        CBLAS_ORDER const stor_ord
 
464
          = enum_cast<CBLAS_ORDER const>
 
465
          (storage_order<
 
466
#ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
 
467
            typename traits::matrix_traits<HermA>::ordering_type
 
468
#else
 
469
            typename HermA::orientation_category 
 
470
#endif 
 
471
           >::value); 
 
472
 
 
473
        hemm (stor_ord, side, uplo,  
 
474
              m, n, alpha, 
 
475
#ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
 
476
              traits::matrix_storage (a), 
 
477
#else
 
478
              traits::matrix_storage_const (a), 
 
479
#endif
 
480
              traits::leading_dimension (a),
 
481
#ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
 
482
              traits::matrix_storage (b), 
 
483
#else
 
484
              traits::matrix_storage_const (b), 
 
485
#endif
 
486
              traits::leading_dimension (b),
 
487
              beta, 
 
488
              traits::matrix_storage (c), 
 
489
              traits::leading_dimension (c)); 
 
490
      }
 
491
 
 
492
    } // detail 
 
493
 
 
494
    // C <- alpha * A * B + beta * C 
 
495
    // C <- alpha * B * A + beta * C 
 
496
    // A == A^H
 
497
    template <typename T, typename HermA, typename MatrB, typename MatrC>
 
498
    inline
 
499
    void hemm (CBLAS_SIDE const side, CBLAS_UPLO const uplo, 
 
500
               T const& alpha, HermA const& a, MatrB const& b, 
 
501
               T const& beta, MatrC& c)
 
502
    {
 
503
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
 
504
      BOOST_STATIC_ASSERT((boost::is_same<
 
505
        typename traits::matrix_traits<HermA>::matrix_structure, 
 
506
        traits::general_t
 
507
      >::value)); 
 
508
#endif 
 
509
 
 
510
      detail::hemm (side, uplo, alpha, a, b, beta, c); 
 
511
    }
 
512
 
 
513
    template <typename T, typename HermA, typename MatrB, typename MatrC>
 
514
    inline
 
515
    void hemm (CBLAS_SIDE const side, 
 
516
               T const& alpha, HermA const& a, MatrB const& b, 
 
517
               T const& beta, MatrC& c)
 
518
    {
 
519
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
 
520
      BOOST_STATIC_ASSERT((boost::is_same<
 
521
        typename traits::matrix_traits<HermA>::matrix_structure, 
 
522
        traits::hermitian_t
 
523
      >::value)); 
 
524
#endif 
 
525
 
 
526
      CBLAS_UPLO const uplo
 
527
        = enum_cast<CBLAS_UPLO const>
 
528
        (uplo_triang<
 
529
#ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
 
530
          typename traits::matrix_traits<HermA>::uplo_type
 
531
#else
 
532
          typename HermA::packed_category 
 
533
#endif 
 
534
         >::value); 
 
535
 
 
536
      detail::hemm (side, uplo, alpha, a, b, beta, c); 
 
537
    }
 
538
 
 
539
 
 
540
#ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
 
541
 
 
542
    namespace detail {
 
543
 
 
544
      // C <- alpha * A * B + beta * C ;  A == A^H
 
545
      struct hemm_left {
 
546
        template <typename T, typename HermA, typename MatrB, typename MatrC>
 
547
        static void f (T const& alpha, HermA const& a, MatrB const& b, 
 
548
                       T const& beta, MatrC& c) 
 
549
        {
 
550
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
 
551
          BOOST_STATIC_ASSERT((boost::is_same<
 
552
            typename traits::matrix_traits<HermA>::matrix_structure, 
 
553
            traits::hermitian_t
 
554
          >::value)); 
 
555
          BOOST_STATIC_ASSERT((boost::is_same<
 
556
            typename traits::matrix_traits<MatrB>::matrix_structure, 
 
557
            traits::general_t
 
558
          >::value));
 
559
#endif 
 
560
 
 
561
          int const m = traits::matrix_size1 (c);
 
562
          int const n = traits::matrix_size2 (c);
 
563
 
 
564
          assert (m == traits::matrix_size1 (a) 
 
565
                  && m == traits::matrix_size2 (a)); 
 
566
          assert (m == traits::matrix_size1 (b) 
 
567
                  && n == traits::matrix_size2 (b)); 
 
568
 
 
569
          CBLAS_ORDER const stor_ord
 
570
            = enum_cast<CBLAS_ORDER const>
 
571
            (storage_order<
 
572
              typename traits::matrix_traits<HermA>::ordering_type
 
573
             >::value); 
 
574
 
 
575
          CBLAS_UPLO const uplo
 
576
            = enum_cast<CBLAS_UPLO const>
 
577
            (uplo_triang<
 
578
              typename traits::matrix_traits<HermA>::uplo_type
 
579
             >::value); 
 
580
 
 
581
          hemm (stor_ord, CblasLeft, uplo,  
 
582
                m, n, alpha, 
 
583
                traits::matrix_storage (a), traits::leading_dimension (a),
 
584
                traits::matrix_storage (b), traits::leading_dimension (b),
 
585
                beta, 
 
586
                traits::matrix_storage (c), traits::leading_dimension (c)); 
 
587
        }
 
588
      }; 
 
589
 
 
590
      // C <- alpha * A * B + beta * C ;  B == B^H
 
591
      struct hemm_right {
 
592
        template <typename T, typename MatrA, typename HermB, typename MatrC>
 
593
        static void f (T const& alpha, MatrA const& a, HermB const& b, 
 
594
                       T const& beta, MatrC& c) 
 
595
        {
 
596
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
 
597
          BOOST_STATIC_ASSERT((boost::is_same<
 
598
            typename traits::matrix_traits<MatrA>::matrix_structure, 
 
599
            traits::general_t
 
600
          >::value));
 
601
          BOOST_STATIC_ASSERT((boost::is_same<
 
602
            typename traits::matrix_traits<HermB>::matrix_structure, 
 
603
            traits::hermitian_t
 
604
          >::value));
 
605
#endif 
 
606
 
 
607
          int const m = traits::matrix_size1 (c);
 
608
          int const n = traits::matrix_size2 (c);
 
609
 
 
610
          assert (n == traits::matrix_size1 (b) 
 
611
                  && n == traits::matrix_size2 (b)); 
 
612
          assert (m == traits::matrix_size1 (a) 
 
613
                  && n == traits::matrix_size2 (a)); 
 
614
 
 
615
          CBLAS_ORDER const stor_ord
 
616
            = enum_cast<CBLAS_ORDER const>
 
617
            (storage_order<
 
618
              typename traits::matrix_traits<HermB>::ordering_type
 
619
             >::value); 
 
620
 
 
621
          CBLAS_UPLO const uplo
 
622
            = enum_cast<CBLAS_UPLO const>
 
623
            (uplo_triang<
 
624
              typename traits::matrix_traits<HermB>::uplo_type
 
625
             >::value); 
 
626
 
 
627
          hemm (stor_ord, CblasRight, uplo,  
 
628
                m, n, alpha, 
 
629
                traits::matrix_storage (b), traits::leading_dimension (b),
 
630
                traits::matrix_storage (a), traits::leading_dimension (a),
 
631
                beta, 
 
632
                traits::matrix_storage (c), traits::leading_dimension (c)); 
 
633
        }
 
634
      }; 
 
635
 
 
636
    } 
 
637
    
 
638
    template <typename T, typename MatrA, typename MatrB, typename MatrC>
 
639
    inline
 
640
    void hemm (T const& alpha, MatrA const& a, MatrB const& b, 
 
641
               T const& beta, MatrC& c)
 
642
    {
 
643
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
 
644
      BOOST_STATIC_ASSERT((boost::is_same<
 
645
        typename traits::matrix_traits<MatrC>::matrix_structure, 
 
646
        traits::general_t
 
647
      >::value)); 
 
648
 
 
649
      BOOST_STATIC_ASSERT((boost::is_same<
 
650
        typename traits::matrix_traits<MatrA>::ordering_type,
 
651
        typename traits::matrix_traits<MatrB>::ordering_type
 
652
      >::value)); 
 
653
      BOOST_STATIC_ASSERT((boost::is_same<
 
654
        typename traits::matrix_traits<MatrA>::ordering_type,
 
655
        typename traits::matrix_traits<MatrC>::ordering_type
 
656
      >::value)); 
 
657
#endif 
 
658
 
 
659
      typedef typename
 
660
        boost::mpl::if_c<
 
661
          boost::is_same<
 
662
            typename traits::matrix_traits<MatrA>::matrix_structure, 
 
663
            traits::hermitian_t
 
664
          >::value,
 
665
          detail::hemm_left, 
 
666
          detail::hemm_right
 
667
        >::type functor; 
 
668
 
 
669
      functor::f (alpha, a, b, beta, c); 
 
670
    }
 
671
 
 
672
    // C <- A * B  
 
673
    // C <- B * A  
 
674
    template <typename MatrA, typename MatrB, typename MatrC>
 
675
    inline
 
676
    void hemm (MatrA const& a, MatrB const& b, MatrC& c) {
 
677
      typedef typename traits::matrix_traits<MatrC>::value_type val_t; 
 
678
      hemm ((val_t) 1, a, b, (val_t) 0, c);
 
679
    }
 
680
 
 
681
#endif // BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
 
682
 
 
683
    
 
684
    // C <- alpha * A * A^T + beta * C
 
685
    // C <- alpha * A^T * A + beta * C
 
686
    // C == C^T
 
687
 
 
688
    namespace detail {
 
689
 
 
690
      template <typename T, typename MatrA, typename SymmC>
 
691
      inline
 
692
      void syrk (CBLAS_UPLO const uplo, CBLAS_TRANSPOSE trans, 
 
693
                 T const& alpha, MatrA const& a, 
 
694
                 T const& beta, SymmC& c)
 
695
      {
 
696
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
 
697
        BOOST_STATIC_ASSERT((boost::is_same<
 
698
          typename traits::matrix_traits<MatrA>::matrix_structure, 
 
699
          traits::general_t
 
700
        >::value));
 
701
        BOOST_STATIC_ASSERT((boost::is_same<
 
702
          typename traits::matrix_traits<MatrA>::ordering_type,
 
703
          typename traits::matrix_traits<SymmC>::ordering_type
 
704
        >::value)); 
 
705
#endif 
 
706
 
 
707
        assert (uplo == CblasUpper || uplo == CblasLower); 
 
708
        assert (trans == CblasNoTrans 
 
709
                || trans == CblasTrans 
 
710
                || trans == CblasConjTrans); 
 
711
 
 
712
        int const n = traits::matrix_size1 (c);
 
713
        assert (n == traits::matrix_size2 (c)); 
 
714
        
 
715
        int const k = trans == CblasNoTrans
 
716
          ? traits::matrix_size2 (a)
 
717
          : traits::matrix_size1 (a); 
 
718
        assert (n == (trans == CblasNoTrans
 
719
                      ? traits::matrix_size1 (a)
 
720
                      : traits::matrix_size2 (a))); 
 
721
 
 
722
        CBLAS_ORDER const stor_ord
 
723
          = enum_cast<CBLAS_ORDER const>
 
724
          (storage_order<
 
725
#ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
 
726
            typename traits::matrix_traits<SymmC>::ordering_type
 
727
#else
 
728
            typename SymmC::orientation_category 
 
729
#endif 
 
730
           >::value); 
 
731
 
 
732
        syrk (stor_ord, uplo, trans, 
 
733
              n, k, alpha, 
 
734
#ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
 
735
              traits::matrix_storage (a), 
 
736
#else
 
737
              traits::matrix_storage_const (a), 
 
738
#endif
 
739
              traits::leading_dimension (a),
 
740
              beta, 
 
741
              traits::matrix_storage (c), 
 
742
              traits::leading_dimension (c)); 
 
743
      }
 
744
 
 
745
    } // detail 
 
746
 
 
747
    // C <- alpha * A * A^T + beta * C
 
748
    // C <- alpha * A^T * A + beta * C
 
749
    // C == C^T
 
750
    template <typename T, typename MatrA, typename SymmC>
 
751
    inline
 
752
    void syrk (CBLAS_UPLO const uplo, CBLAS_TRANSPOSE trans, 
 
753
               T const& alpha, MatrA const& a, 
 
754
               T const& beta, SymmC& c)
 
755
    {
 
756
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
 
757
      BOOST_STATIC_ASSERT((boost::is_same<
 
758
        typename traits::matrix_traits<SymmC>::matrix_structure, 
 
759
        traits::general_t
 
760
      >::value)); 
 
761
#endif 
 
762
 
 
763
      detail::syrk (uplo, trans, alpha, a, beta, c); 
 
764
    }
 
765
 
 
766
    template <typename T, typename MatrA, typename SymmC>
 
767
    inline
 
768
    void syrk (CBLAS_TRANSPOSE trans, 
 
769
               T const& alpha, MatrA const& a, 
 
770
               T const& beta, SymmC& c)
 
771
    {
 
772
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
 
773
      BOOST_STATIC_ASSERT((boost::is_same<
 
774
        typename traits::matrix_traits<SymmC>::matrix_structure, 
 
775
        traits::symmetric_t
 
776
      >::value)); 
 
777
#endif 
 
778
 
 
779
      CBLAS_UPLO const uplo
 
780
        = enum_cast<CBLAS_UPLO const>
 
781
        (uplo_triang<
 
782
#ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
 
783
          typename traits::matrix_traits<SymmC>::uplo_type
 
784
#else
 
785
          typename SymmC::packed_category 
 
786
#endif 
 
787
         >::value); 
 
788
 
 
789
      detail::syrk (uplo, trans, alpha, a, beta, c); 
 
790
    }
 
791
 
 
792
    // C <- A * A^T + C
 
793
    // C <- A^T * A + C
 
794
    template <typename MatrA, typename SymmC>
 
795
    inline
 
796
    void syrk (CBLAS_TRANSPOSE trans, MatrA const& a, SymmC& c) {
 
797
#ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
 
798
      typedef typename traits::matrix_traits<SymmC>::value_type val_t; 
 
799
#else
 
800
      typedef typename SymmC::value_type val_t; 
 
801
#endif 
 
802
      syrk (trans, (val_t) 1, a, (val_t) 0, c);
 
803
    }
 
804
 
 
805
    
 
806
    // C <- alpha * A * B^T + conj(alpha) * B * A^T + beta * C
 
807
    // C <- alpha * A^T * B + conj(alpha) * B^T * A + beta * C
 
808
    // C == C^T
 
809
 
 
810
    namespace detail {
 
811
 
 
812
      template <typename T, typename MatrA, typename MatrB, typename SymmC>
 
813
      inline
 
814
      void syr2k (CBLAS_UPLO const uplo, CBLAS_TRANSPOSE trans, 
 
815
                  T const& alpha, MatrA const& a, MatrB const& b, 
 
816
                  T const& beta, SymmC& c)
 
817
      {
 
818
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
 
819
        BOOST_STATIC_ASSERT((boost::is_same<
 
820
          typename traits::matrix_traits<MatrA>::matrix_structure, 
 
821
          traits::general_t
 
822
        >::value));
 
823
        BOOST_STATIC_ASSERT((boost::is_same<
 
824
          typename traits::matrix_traits<MatrB>::matrix_structure, 
 
825
          traits::general_t
 
826
        >::value));
 
827
        BOOST_STATIC_ASSERT((boost::is_same<
 
828
          typename traits::matrix_traits<MatrA>::ordering_type,
 
829
          typename traits::matrix_traits<SymmC>::ordering_type
 
830
        >::value)); 
 
831
        BOOST_STATIC_ASSERT((boost::is_same<
 
832
          typename traits::matrix_traits<MatrB>::ordering_type,
 
833
          typename traits::matrix_traits<SymmC>::ordering_type
 
834
        >::value)); 
 
835
#endif 
 
836
 
 
837
        assert (uplo == CblasUpper || uplo == CblasLower); 
 
838
        assert (trans == CblasNoTrans 
 
839
                || trans == CblasTrans 
 
840
                || trans == CblasConjTrans); 
 
841
 
 
842
        int const n = traits::matrix_size1 (c);
 
843
        assert (n == traits::matrix_size2 (c)); 
 
844
        
 
845
        int const k = trans == CblasNoTrans
 
846
          ? traits::matrix_size2 (a)
 
847
          : traits::matrix_size1 (a); 
 
848
        assert (k == (trans == CblasNoTrans
 
849
                      ? traits::matrix_size2 (b)
 
850
                      : traits::matrix_size1 (b))); 
 
851
        assert (n == (trans == CblasNoTrans
 
852
                      ? traits::matrix_size1 (a)
 
853
                      : traits::matrix_size2 (a))); 
 
854
        assert (n == (trans == CblasNoTrans
 
855
                      ? traits::matrix_size1 (b)
 
856
                      : traits::matrix_size2 (b))); 
 
857
 
 
858
        CBLAS_ORDER const stor_ord
 
859
          = enum_cast<CBLAS_ORDER const>
 
860
          (storage_order<
 
861
#ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
 
862
            typename traits::matrix_traits<SymmC>::ordering_type
 
863
#else
 
864
            typename SymmC::orientation_category 
 
865
#endif 
 
866
           >::value); 
 
867
 
 
868
        syr2k (stor_ord, uplo, trans, 
 
869
               n, k, alpha, 
 
870
#ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
 
871
               traits::matrix_storage (a), 
 
872
#else
 
873
               traits::matrix_storage_const (a), 
 
874
#endif
 
875
               traits::leading_dimension (a),
 
876
#ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
 
877
               traits::matrix_storage (b), 
 
878
#else
 
879
               traits::matrix_storage_const (b), 
 
880
#endif
 
881
               traits::leading_dimension (b),
 
882
               beta, 
 
883
               traits::matrix_storage (c), 
 
884
               traits::leading_dimension (c)); 
 
885
      }
 
886
 
 
887
    } // detail 
 
888
 
 
889
    // C <- alpha * A * B^T + conj(alpha) * B * A^T + beta * C
 
890
    // C <- alpha * A^T * B + conj(alpha) * B^T * A + beta * C
 
891
    // C == C^T
 
892
    template <typename T, typename MatrA, typename MatrB, typename SymmC>
 
893
    inline
 
894
    void syr2k (CBLAS_UPLO const uplo, CBLAS_TRANSPOSE trans, 
 
895
               T const& alpha, MatrA const& a, MatrB const& b, 
 
896
               T const& beta, SymmC& c)
 
897
    {
 
898
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
 
899
      BOOST_STATIC_ASSERT((boost::is_same<
 
900
        typename traits::matrix_traits<SymmC>::matrix_structure, 
 
901
        traits::general_t
 
902
      >::value)); 
 
903
#endif 
 
904
 
 
905
      detail::syr2k (uplo, trans, alpha, a, b, beta, c); 
 
906
    }
 
907
 
 
908
    template <typename T, typename MatrA, typename MatrB, typename SymmC>
 
909
    inline
 
910
    void syr2k (CBLAS_TRANSPOSE trans, 
 
911
               T const& alpha, MatrA const& a, MatrB const& b,  
 
912
               T const& beta, SymmC& c)
 
913
    {
 
914
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
 
915
      BOOST_STATIC_ASSERT((boost::is_same<
 
916
        typename traits::matrix_traits<SymmC>::matrix_structure, 
 
917
        traits::symmetric_t
 
918
      >::value)); 
 
919
#endif 
 
920
 
 
921
      CBLAS_UPLO const uplo
 
922
        = enum_cast<CBLAS_UPLO const>
 
923
        (uplo_triang<
 
924
#ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
 
925
          typename traits::matrix_traits<SymmC>::uplo_type
 
926
#else
 
927
          typename SymmC::packed_category 
 
928
#endif 
 
929
         >::value); 
 
930
 
 
931
      detail::syr2k (uplo, trans, alpha, a, b, beta, c); 
 
932
    }
 
933
 
 
934
    // C <- A * B^T + B * A^T + C
 
935
    // C <- A^T * B + B^T * A + C
 
936
    template <typename MatrA, typename MatrB, typename SymmC>
 
937
    inline
 
938
    void syr2k (CBLAS_TRANSPOSE trans, 
 
939
                MatrA const& a, MatrB const& b, SymmC& c) 
 
940
    {
 
941
#ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
 
942
      typedef typename traits::matrix_traits<SymmC>::value_type val_t; 
 
943
#else
 
944
      typedef typename SymmC::value_type val_t; 
 
945
#endif 
 
946
      syr2k (trans, (val_t) 1, a, b, (val_t) 0, c);
 
947
    }
 
948
 
 
949
    
 
950
    // C <- alpha * A * A^H + beta * C
 
951
    // C <- alpha * A^H * A + beta * C
 
952
    // C == C^H
 
953
 
 
954
    namespace detail {
 
955
 
 
956
      template <typename T, typename MatrA, typename HermC>
 
957
      inline
 
958
      void herk (CBLAS_UPLO const uplo, CBLAS_TRANSPOSE trans, 
 
959
                 T const& alpha, MatrA const& a, 
 
960
                 T const& beta, HermC& c)
 
961
      {
 
962
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
 
963
        BOOST_STATIC_ASSERT((boost::is_same<
 
964
          typename traits::matrix_traits<MatrA>::matrix_structure, 
 
965
          traits::general_t
 
966
        >::value));
 
967
        BOOST_STATIC_ASSERT((boost::is_same<
 
968
          typename traits::matrix_traits<MatrA>::ordering_type,
 
969
          typename traits::matrix_traits<HermC>::ordering_type
 
970
        >::value)); 
 
971
#endif 
 
972
 
 
973
        assert (uplo == CblasUpper || uplo == CblasLower); 
 
974
        assert (trans == CblasNoTrans || trans == CblasConjTrans); 
 
975
 
 
976
        int const n = traits::matrix_size1 (c);
 
977
        assert (n == traits::matrix_size2 (c)); 
 
978
        
 
979
        int const k = trans == CblasNoTrans
 
980
          ? traits::matrix_size2 (a)
 
981
          : traits::matrix_size1 (a); 
 
982
        assert (n == (trans == CblasNoTrans
 
983
                      ? traits::matrix_size1 (a)
 
984
                      : traits::matrix_size2 (a))); 
 
985
 
 
986
        CBLAS_ORDER const stor_ord
 
987
          = enum_cast<CBLAS_ORDER const>
 
988
          (storage_order<
 
989
#ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
 
990
            typename traits::matrix_traits<HermC>::ordering_type
 
991
#else
 
992
            typename HermC::orientation_category 
 
993
#endif 
 
994
           >::value); 
 
995
 
 
996
        herk (stor_ord, uplo, trans, 
 
997
              n, k, alpha, 
 
998
#ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
 
999
              traits::matrix_storage (a), 
 
1000
#else
 
1001
              traits::matrix_storage_const (a), 
 
1002
#endif
 
1003
              traits::leading_dimension (a),
 
1004
              beta, 
 
1005
              traits::matrix_storage (c), 
 
1006
              traits::leading_dimension (c)); 
 
1007
      }
 
1008
 
 
1009
    } // detail 
 
1010
 
 
1011
    // C <- alpha * A * A^H + beta * C
 
1012
    // C <- alpha * A^H * A + beta * C
 
1013
    // C == C^H
 
1014
    template <typename T, typename MatrA, typename HermC>
 
1015
    inline
 
1016
    void herk (CBLAS_UPLO const uplo, CBLAS_TRANSPOSE trans, 
 
1017
               T const& alpha, MatrA const& a, 
 
1018
               T const& beta, HermC& c)
 
1019
    {
 
1020
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
 
1021
      BOOST_STATIC_ASSERT((boost::is_same<
 
1022
        typename traits::matrix_traits<HermC>::matrix_structure, 
 
1023
        traits::general_t
 
1024
      >::value)); 
 
1025
#endif 
 
1026
 
 
1027
      detail::herk (uplo, trans, alpha, a, beta, c); 
 
1028
    }
 
1029
 
 
1030
    template <typename T, typename MatrA, typename HermC>
 
1031
    inline
 
1032
    void herk (CBLAS_TRANSPOSE trans, 
 
1033
               T const& alpha, MatrA const& a, 
 
1034
               T const& beta, HermC& c)
 
1035
    {
 
1036
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
 
1037
      BOOST_STATIC_ASSERT((boost::is_same<
 
1038
        typename traits::matrix_traits<HermC>::matrix_structure, 
 
1039
        traits::hermitian_t
 
1040
      >::value)); 
 
1041
#endif 
 
1042
 
 
1043
      CBLAS_UPLO const uplo
 
1044
        = enum_cast<CBLAS_UPLO const>
 
1045
        (uplo_triang<
 
1046
#ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
 
1047
          typename traits::matrix_traits<HermC>::uplo_type
 
1048
#else
 
1049
          typename HermC::packed_category 
 
1050
#endif 
 
1051
         >::value); 
 
1052
 
 
1053
      detail::herk (uplo, trans, alpha, a, beta, c); 
 
1054
    }
 
1055
 
 
1056
    // C <- A * A^H + C
 
1057
    // C <- A^H * A + C
 
1058
    template <typename MatrA, typename HermC>
 
1059
    inline
 
1060
    void herk (CBLAS_TRANSPOSE trans, MatrA const& a, HermC& c) {
 
1061
#ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
 
1062
      typedef typename traits::matrix_traits<HermC>::value_type val_t; 
 
1063
#else
 
1064
      typedef typename HermC::value_type val_t; 
 
1065
#endif 
 
1066
      typedef typename traits::type_traits<val_t>::real_type real_t; 
 
1067
      herk (trans, (real_t) 1, a, (real_t) 0, c);
 
1068
    }
 
1069
 
 
1070
 
 
1071
    // C <- alpha * A * B^H + conj(alpha) * B * A^H + beta * C
 
1072
    // C <- alpha * A^H * B + conj(alpha) * B^H * A + beta * C
 
1073
    // C == C^H
 
1074
 
 
1075
    namespace detail {
 
1076
 
 
1077
      template <typename T1, typename T2, 
 
1078
                typename MatrA, typename MatrB, typename HermC>
 
1079
      inline
 
1080
      void her2k (CBLAS_UPLO const uplo, CBLAS_TRANSPOSE trans, 
 
1081
                  T1 const& alpha, MatrA const& a, MatrB const& b, 
 
1082
                  T2 const& beta, HermC& c)
 
1083
      {
 
1084
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
 
1085
        BOOST_STATIC_ASSERT((boost::is_same<
 
1086
          typename traits::matrix_traits<MatrA>::matrix_structure, 
 
1087
          traits::general_t
 
1088
        >::value));
 
1089
        BOOST_STATIC_ASSERT((boost::is_same<
 
1090
          typename traits::matrix_traits<MatrB>::matrix_structure, 
 
1091
          traits::general_t
 
1092
        >::value));
 
1093
        BOOST_STATIC_ASSERT((boost::is_same<
 
1094
          typename traits::matrix_traits<MatrA>::ordering_type,
 
1095
          typename traits::matrix_traits<HermC>::ordering_type
 
1096
        >::value)); 
 
1097
        BOOST_STATIC_ASSERT((boost::is_same<
 
1098
          typename traits::matrix_traits<MatrB>::ordering_type,
 
1099
          typename traits::matrix_traits<HermC>::ordering_type
 
1100
        >::value)); 
 
1101
#endif 
 
1102
 
 
1103
        assert (uplo == CblasUpper || uplo == CblasLower); 
 
1104
        assert (trans == CblasNoTrans || trans == CblasConjTrans); 
 
1105
 
 
1106
        int const n = traits::matrix_size1 (c);
 
1107
        assert (n == traits::matrix_size2 (c)); 
 
1108
        
 
1109
        int const k = trans == CblasNoTrans
 
1110
          ? traits::matrix_size2 (a)
 
1111
          : traits::matrix_size1 (a); 
 
1112
        assert (k == (trans == CblasNoTrans
 
1113
                      ? traits::matrix_size2 (b)
 
1114
                      : traits::matrix_size1 (b))); 
 
1115
        assert (n == (trans == CblasNoTrans
 
1116
                      ? traits::matrix_size1 (a)
 
1117
                      : traits::matrix_size2 (a))); 
 
1118
        assert (n == (trans == CblasNoTrans
 
1119
                      ? traits::matrix_size1 (b)
 
1120
                      : traits::matrix_size2 (b))); 
 
1121
 
 
1122
        CBLAS_ORDER const stor_ord
 
1123
          = enum_cast<CBLAS_ORDER const>
 
1124
          (storage_order<
 
1125
#ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
 
1126
            typename traits::matrix_traits<HermC>::ordering_type
 
1127
#else
 
1128
            typename HermC::orientation_category 
 
1129
#endif 
 
1130
           >::value); 
 
1131
 
 
1132
        her2k (stor_ord, uplo, trans, 
 
1133
               n, k, alpha, 
 
1134
#ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
 
1135
               traits::matrix_storage (a), 
 
1136
#else
 
1137
               traits::matrix_storage_const (a), 
 
1138
#endif
 
1139
               traits::leading_dimension (a),
 
1140
#ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
 
1141
               traits::matrix_storage (b), 
 
1142
#else
 
1143
               traits::matrix_storage_const (b), 
 
1144
#endif
 
1145
               traits::leading_dimension (b),
 
1146
               beta, 
 
1147
               traits::matrix_storage (c), 
 
1148
               traits::leading_dimension (c)); 
 
1149
      }
 
1150
 
 
1151
    } // detail 
 
1152
 
 
1153
    // C <- alpha * A * B^H + conj(alpha) * B * A^H + beta * C
 
1154
    // C <- alpha * A^H * B + conj(alpha) * B^H * A + beta * C
 
1155
    // C == C^H
 
1156
    template <typename T1, typename T2, 
 
1157
              typename MatrA, typename MatrB, typename HermC>
 
1158
    inline
 
1159
    void her2k (CBLAS_UPLO const uplo, CBLAS_TRANSPOSE trans, 
 
1160
               T1 const& alpha, MatrA const& a, MatrB const& b, 
 
1161
               T2 const& beta, HermC& c)
 
1162
    {
 
1163
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
 
1164
      BOOST_STATIC_ASSERT((boost::is_same<
 
1165
        typename traits::matrix_traits<HermC>::matrix_structure, 
 
1166
        traits::general_t
 
1167
      >::value)); 
 
1168
#endif 
 
1169
 
 
1170
      detail::her2k (uplo, trans, alpha, a, b, beta, c); 
 
1171
    }
 
1172
 
 
1173
    template <typename T1, typename T2, 
 
1174
              typename MatrA, typename MatrB, typename HermC>
 
1175
    inline
 
1176
    void her2k (CBLAS_TRANSPOSE trans, 
 
1177
               T1 const& alpha, MatrA const& a, MatrB const& b,  
 
1178
               T2 const& beta, HermC& c)
 
1179
    {
 
1180
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
 
1181
      BOOST_STATIC_ASSERT((boost::is_same<
 
1182
        typename traits::matrix_traits<HermC>::matrix_structure, 
 
1183
        traits::hermitian_t
 
1184
      >::value)); 
 
1185
#endif 
 
1186
 
 
1187
      CBLAS_UPLO const uplo
 
1188
        = enum_cast<CBLAS_UPLO const>
 
1189
        (uplo_triang<
 
1190
#ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
 
1191
          typename traits::matrix_traits<HermC>::uplo_type
 
1192
#else
 
1193
          typename HermC::packed_category 
 
1194
#endif 
 
1195
         >::value); 
 
1196
 
 
1197
      detail::her2k (uplo, trans, alpha, a, b, beta, c); 
 
1198
    }
 
1199
 
 
1200
    // C <- A * B^H + B * A^H + C
 
1201
    // C <- A^H * B + B^H * A + C
 
1202
    template <typename MatrA, typename MatrB, typename HermC>
 
1203
    inline
 
1204
    void her2k (CBLAS_TRANSPOSE trans, 
 
1205
                MatrA const& a, MatrB const& b, HermC& c) 
 
1206
    {
 
1207
#ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
 
1208
      typedef typename traits::matrix_traits<HermC>::value_type val_t; 
 
1209
#else
 
1210
      typedef typename HermC::value_type val_t; 
 
1211
#endif 
 
1212
      typedef typename traits::type_traits<val_t>::real_type real_t; 
 
1213
      her2k (trans, (val_t) 1, a, b, (real_t) 0, c);
 
1214
    }
 
1215
 
 
1216
    
 
1217
  } // namespace atlas
 
1218
 
 
1219
}}} 
 
1220
 
 
1221
#endif // BOOST_NUMERIC_BINDINGS_CBLAS_LEVEL_3_HPP