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