← 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/lapack/orgqr.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/lapack/orgqr.hpp
 
1
//
 
2
// Copyright Fabien Dekeyser, Quoc-Cuong Pham 2008
 
3
//
 
4
// Distributed under the Boost Software License, Version 1.0.
 
5
// (See accompanying file LICENSE_1_0.txt or copy at
 
6
// http://www.boost.org/LICENSE_1_0.txt)
 
7
//
 
8
 
 
9
/**
 
10
 *
 
11
 * \brief Generates an M-by-N real matrix Q with orthonormal columns,
 
12
 * which is defined as the first N columns of a product of K elementary
 
13
 *        reflectors of order M
 
14
 *        Q  =  H(1) H(2) . . . H(k)
 
15
 *
 
16
 * \warning 
 
17
 * \todo 
 
18
 * \date 2005
 
19
 * \author CEA/DRT/DTSI/SARC 
 
20
 * \author Q.C. PHAM
 
21
 *
 
22
 **/
 
23
 
 
24
 
 
25
#ifndef BOOST_NUMERIC_BINDINGS_LAPACK_ORGQR_HPP
 
26
#define BOOST_NUMERIC_BINDINGS_LAPACK_ORGQR_HPP
 
27
 
 
28
#include <complex>
 
29
#include <boost/numeric/bindings/traits/traits.hpp>
 
30
#include <boost/numeric/bindings/lapack/lapack.h>
 
31
#include <boost/numeric/bindings/lapack/workspace.hpp>
 
32
#include <boost/numeric/bindings/traits/detail/array.hpp>
 
33
// #include <boost/numeric/bindings/traits/std_vector.hpp>
 
34
 
 
35
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
 
36
#  include <boost/static_assert.hpp>
 
37
#  include <boost/type_traits.hpp>
 
38
#endif 
 
39
 
 
40
 
 
41
namespace boost { namespace numeric { namespace bindings { 
 
42
 
 
43
  namespace lapack {
 
44
 
 
45
   
 
46
    ///////////////////////////////////////////////////////////////////
 
47
    /* 
 
48
     * Generates an M-by-N real matrix Q with orthonormal columns,
 
49
         * which is defined as the first N columns of a product of K elementary
 
50
         *        reflectors of order M
 
51
         *        Q  =  H(1) H(2) . . . H(k)
 
52
         * as returned by geqrf().
 
53
         * The init value of matrix Q is the matrix A returned by geqrf()
 
54
     */ 
 
55
         ///////////////////////////////////////////////////////////////////
 
56
    namespace detail {
 
57
 
 
58
 
 
59
      inline 
 
60
      void orgqr(int const m, int const n, int const k,
 
61
                 float* a, int const lda,
 
62
                 float* tau, float* work, int const lwork, int& info) 
 
63
      {
 
64
        LAPACK_SORGQR (&m, &n, &k, a, &lda, tau, work, &lwork, &info);
 
65
      }
 
66
 
 
67
      inline 
 
68
      void orgqr(int const m, int const n, int const k,
 
69
                 double* a, int const lda,
 
70
                 double* tau, double* work, int const lwork, int& info) 
 
71
      {
 
72
        LAPACK_DORGQR (&m, &n, &k, a, &lda, tau, work, &lwork, &info);
 
73
      }
 
74
 
 
75
      inline 
 
76
      void orgqr(int const m, int const n, int const k,
 
77
                 traits::complex_f* a, int const lda,
 
78
                 traits::complex_f* tau, traits::complex_f* work, int const lwork, int& info) 
 
79
      {
 
80
        LAPACK_CUNGQR (&m, &n, &k, traits::complex_ptr(a), &lda, traits::complex_ptr(tau),
 
81
                       traits::complex_ptr(work), &lwork, &info);
 
82
      }
 
83
 
 
84
      inline 
 
85
      void orgqr(int const m, int const n, int const k,
 
86
                 traits::complex_d* a, int const lda,
 
87
                 traits::complex_d* tau, traits::complex_d* work, int const lwork, int& info) 
 
88
      {
 
89
        LAPACK_ZUNGQR (&m, &n, &k, traits::complex_ptr(a), &lda, traits::complex_ptr(tau),
 
90
                       traits::complex_ptr(work), &lwork, &info);
 
91
      }
 
92
 
 
93
    } // fin namespace detail
 
94
 
 
95
 
 
96
        //--------------------------------------------
 
97
 
 
98
   template <typename A, typename Tau, typename Work>
 
99
   inline int orgqr(A& a, Tau& tau, Work &work) {
 
100
 
 
101
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
 
102
      BOOST_STATIC_ASSERT((boost::is_same<
 
103
        typename traits::matrix_traits<A>::matrix_structure, 
 
104
        traits::general_t
 
105
      >::value)); 
 
106
#endif 
 
107
 
 
108
          const int m = traits::matrix_size1 (a);
 
109
          const int n = traits::matrix_size2 (a);
 
110
          const int k = n;
 
111
 
 
112
          assert (std::min<int>(m,n) <= traits::vector_size (tau)); 
 
113
      assert (n <= traits::vector_size (work)); 
 
114
         
 
115
          int info; 
 
116
      detail::orgqr (m, n, k,
 
117
                     traits::matrix_storage (a), 
 
118
                     traits::leading_dimension (a),
 
119
                     traits::vector_storage (tau),  
 
120
                     traits::vector_storage (work),
 
121
                     traits::vector_size (work),
 
122
                     info);
 
123
      return info;
 
124
    }
 
125
 
 
126
    // Computation of Q.
 
127
    // Workspace is allocated dynamically so that the optimization of
 
128
    // blas 3 calls is optimal.
 
129
    template <typename A, typename Tau>
 
130
    inline
 
131
    int orgqr (A& a, Tau& tau, optimal_workspace ) {
 
132
       typedef typename A::value_type value_type ;
 
133
           const int n = traits::matrix_size2 (a);
 
134
       traits::detail::array<value_type> work(std::max<int>(1, n*32));
 
135
       return orgqr( a, tau, work );
 
136
           
 
137
    }
 
138
 
 
139
    // Computation of Q
 
140
    // Workspace is allocated dynamically to its minimum size.
 
141
    // Blas 3 calls are not optimal.
 
142
    template <typename A, typename Tau>
 
143
    inline
 
144
    int orgqr (A& a, Tau& tau, minimal_workspace ) {
 
145
       typedef typename A::value_type value_type ;
 
146
           const int n = traits::matrix_size2 (a);
 
147
       traits::detail::array<value_type> work(std::max<int>(1, n));
 
148
       return orgqr( a, tau, work );
 
149
    }
 
150
 
 
151
    // Computation of the Q
 
152
    // Workspace is taken from the array in workspace.
 
153
   
 
154
    template <typename A, typename Tau, typename Work>
 
155
    inline
 
156
    int orgqr (A& a, Tau& tau, detail::workspace1<Work> workspace ) {
 
157
       typedef typename A::value_type value_type ;
 
158
           const int n = traits::matrix_size2 (a);
 
159
       traits::detail::array<value_type> work(std::max<int>(1, n));
 
160
       return orgqr( a, tau, workspace.w_ );
 
161
    }
 
162
 
 
163
    // Function without workarray as argument
 
164
    template <typename A, typename Tau>
 
165
    inline
 
166
    int orgqr (A& a, Tau& tau) {
 
167
       return orgqr( a, tau, optimal_workspace() );
 
168
    }
 
169
 
 
170
  }
 
171
 
 
172
}}}
 
173
 
 
174
#endif