← 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/geqrf.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/geqrf.hpp
 
1
/*
 
2
 * 
 
3
 * Copyright (c) Toon Knapen, Karl Meerbergen & Kresimir Fresl 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
 * KF acknowledges the support of the Faculty of Civil Engineering, 
 
10
 * University of Zagreb, Croatia.
 
11
 *
 
12
 */
 
13
 
 
14
#ifndef BOOST_NUMERIC_BINDINGS_LAPACK_GEQRF_HPP
 
15
#define BOOST_NUMERIC_BINDINGS_LAPACK_GEQRF_HPP
 
16
 
 
17
#include <complex>
 
18
#include <boost/numeric/bindings/traits/traits.hpp>
 
19
#include <boost/numeric/bindings/lapack/lapack.h>
 
20
#include <boost/numeric/bindings/lapack/workspace.hpp>
 
21
#include <boost/numeric/bindings/traits/detail/array.hpp>
 
22
// #include <boost/numeric/bindings/traits/std_vector.hpp>
 
23
 
 
24
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
 
25
#  include <boost/static_assert.hpp>
 
26
#  include <boost/type_traits.hpp>
 
27
#endif 
 
28
#include <cassert>
 
29
 
 
30
 
 
31
namespace boost { namespace numeric { namespace bindings { 
 
32
 
 
33
  namespace lapack {
 
34
 
 
35
    ///////////////////////////////////////////////////////////////////
 
36
    //
 
37
    // QR factorization of a general m x n matrix  A = Q * R
 
38
    // 
 
39
    ///////////////////////////////////////////////////////////////////
 
40
 
 
41
    /* 
 
42
     * geqrf() computes the QR factorization of a rectangular matrix
 
43
     * A = Q *  R,  where Q is a M x min(M,N) matrix with orthogonal
 
44
     * and normalized column (i.e. herm(Q) $ Q = I) and R is a
 
45
     * min(M,N) x N upper triangular matrix.
 
46
     *
 
47
     * On return of geqrf, the elements on and above the diagonal of
 
48
     * A contain the min(M,N) x N upper trapezoidal matrix R (R is
 
49
     * upper triangular if  m  >=  n);  the elements below the diagonal,
 
50
     * with the array TAU, represent the orthogonal matrix Q as a product
 
51
     * of min(M,N) elementary reflectors.
 
52
     */ 
 
53
 
 
54
    namespace detail {
 
55
 
 
56
      inline 
 
57
      void geqrf (int const m, int const n,
 
58
                 float* a, int const lda,
 
59
                 float* tau, float* work, int const lwork, int& info) 
 
60
      {
 
61
        LAPACK_SGEQRF (&m, &n, a, &lda, tau, work, &lwork, &info);
 
62
      }
 
63
 
 
64
      inline 
 
65
      void geqrf (int const m, int const n,
 
66
                 double* a, int const lda,
 
67
                 double* tau, double* work, int const lwork, int& info) 
 
68
      {
 
69
        LAPACK_DGEQRF (&m, &n, a, &lda, tau, work, &lwork, &info);
 
70
      }
 
71
 
 
72
      inline 
 
73
      void geqrf (int const m, int const n,
 
74
                  traits::complex_f* a, int const lda,
 
75
                  traits::complex_f* tau, traits::complex_f* work,
 
76
                  int const lwork, int& info) 
 
77
      {
 
78
        LAPACK_CGEQRF (&m, &n,
 
79
                      traits::complex_ptr (a), &lda,
 
80
                      traits::complex_ptr (tau),
 
81
                      traits::complex_ptr (work), &lwork, &info );
 
82
      }
 
83
      
 
84
 
 
85
      inline 
 
86
      void geqrf (int const m, int const n,
 
87
                  traits::complex_d* a, int const lda,
 
88
                  traits::complex_d* tau, traits::complex_d* work,
 
89
                  int const lwork, int& info) 
 
90
      {
 
91
        LAPACK_ZGEQRF (&m, &n,
 
92
                      traits::complex_ptr (a), &lda,
 
93
                      traits::complex_ptr (tau),
 
94
                      traits::complex_ptr (work), &lwork, &info );
 
95
      }
 
96
      
 
97
    } 
 
98
 
 
99
    template <typename A, typename Tau, typename Work>
 
100
    inline
 
101
    int geqrf (A& a, Tau& tau, Work& work) {
 
102
 
 
103
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
 
104
      BOOST_STATIC_ASSERT((boost::is_same<
 
105
        typename traits::matrix_traits<A>::matrix_structure, 
 
106
        traits::general_t
 
107
      >::value)); 
 
108
#endif 
 
109
 
 
110
      int const m = traits::matrix_size1 (a);
 
111
      int const n = traits::matrix_size2 (a);
 
112
      assert (std::min(m,n) <= traits::vector_size (tau)); 
 
113
      assert (n <= traits::vector_size (work)); 
 
114
 
 
115
      int info; 
 
116
      detail::geqrf (m, n,
 
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 the QR factorization.
 
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 geqrf (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(1, n*32));
 
135
       return geqrf( a, tau, work );
 
136
    }
 
137
 
 
138
    // Computation of the QR factorization.
 
139
    // Workspace is allocated dynamically to its minimum size.
 
140
    // Blas 3 calls are not optimal.
 
141
    template <typename A, typename Tau>
 
142
    inline
 
143
    int geqrf (A& a, Tau& tau, minimal_workspace ) {
 
144
       typedef typename A::value_type value_type ;
 
145
       const int n = traits::matrix_size2 (a);
 
146
       traits::detail::array<value_type> work(std::max(1, n));
 
147
       return geqrf( a, tau, work );
 
148
    }
 
149
 
 
150
    // Computation of the QR factorization.
 
151
    // Workspace is taken from the array in workspace.
 
152
    // The calling sequence is
 
153
    // geqrf( a, tau, workspace( work ) ) where work is an array with the same value_type
 
154
    // as a.
 
155
    template <typename A, typename Tau, typename Work>
 
156
    inline
 
157
    int geqrf (A& a, Tau& tau, detail::workspace1<Work> workspace ) {
 
158
       typedef typename A::value_type value_type ;
 
159
       const int n = traits::matrix_size2 (a);
 
160
       traits::detail::array<value_type> work(std::max(1, n));
 
161
       return geqrf( a, tau, workspace.w_ );
 
162
    }
 
163
 
 
164
    // Function without workarray as argument
 
165
    template <typename A, typename Tau>
 
166
    inline
 
167
    int geqrf (A& a, Tau& tau) {
 
168
       return geqrf( a, tau, optimal_workspace() );
 
169
    }
 
170
 
 
171
  }
 
172
 
 
173
}}}
 
174
 
 
175
#endif