← 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/heev.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/heev.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_HEEV_HPP
 
15
#define BOOST_NUMERIC_BINDINGS_LAPACK_HEEV_HPP
 
16
 
 
17
#include <boost/numeric/bindings/traits/traits.hpp>
 
18
#include <boost/numeric/bindings/traits/type.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
 
 
29
 
 
30
namespace boost { namespace numeric { namespace bindings { 
 
31
 
 
32
  namespace lapack {
 
33
 
 
34
    ///////////////////////////////////////////////////////////////////
 
35
    //
 
36
    // Eigendecomposition of a complex Hermitian matrix A = Q * D * Q'
 
37
    // 
 
38
    ///////////////////////////////////////////////////////////////////
 
39
 
 
40
    /* 
 
41
     * heev() computes the eigendecomposition of a N x N matrix
 
42
     * A = Q * D * Q',  where Q is a N x N unitary matrix and
 
43
     * D is a diagonal matrix. The diagonal element D(i,i) is an
 
44
     * eigenvalue of A and Q(:,i) is a corresponding eigenvector.
 
45
     * The eigenvalues are stored in ascending order.
 
46
     *
 
47
     * On return of heev, A is overwritten by Q and w contains the main
 
48
     * diagonal of D.
 
49
     *
 
50
     * int heev (char jobz, char uplo, A& a, W& w, minimal_workspace ) ;
 
51
     *    jobz : 'V' : compute eigenvectors
 
52
     *           'N' : do not compute eigenvectors
 
53
     *    uplo : 'U' : only the upper triangular part of A is used on input.
 
54
     *           'L' : only the lower triangular part of A is used on input.
 
55
     */ 
 
56
 
 
57
    namespace detail {
 
58
 
 
59
      inline 
 
60
      void heev (char const jobz, char const uplo, int const n,
 
61
                 traits::complex_f* a, int const lda,
 
62
                 float* w, traits::complex_f* work, int const lwork,
 
63
                 float* rwork, int& info) 
 
64
      {
 
65
        LAPACK_CHEEV (&jobz, &uplo, &n,
 
66
                      traits::complex_ptr(a), &lda, w,
 
67
                      traits::complex_ptr(work), &lwork,
 
68
                      rwork, &info);
 
69
      }
 
70
 
 
71
      inline 
 
72
      void heev (char const jobz, char const uplo, int const n,
 
73
                 traits::complex_d* a, int const lda,
 
74
                 double* w, traits::complex_d* work, int const lwork,
 
75
                 double* rwork, int& info) 
 
76
      {
 
77
        LAPACK_ZHEEV (&jobz, &uplo, &n,
 
78
                      traits::complex_ptr(a), &lda, w,
 
79
                      traits::complex_ptr(work), &lwork,
 
80
                      rwork, &info);
 
81
      }
 
82
 
 
83
 
 
84
      template <typename A, typename W, typename Work, typename RWork>
 
85
      inline
 
86
      int heev (char jobz, char uplo, A& a, W& w, Work& work, RWork& rwork) {
 
87
 
 
88
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
 
89
        BOOST_STATIC_ASSERT((boost::is_same<
 
90
          typename traits::matrix_traits<A>::matrix_structure, 
 
91
          traits::general_t
 
92
        >::value)); 
 
93
#endif 
 
94
 
 
95
        int const n = traits::matrix_size1 (a);
 
96
        assert (traits::matrix_size2 (a)==n); 
 
97
        assert (traits::vector_size (w)==n); 
 
98
        assert (2*n-1 <= traits::vector_size (work)); 
 
99
        assert (3*n-2 <= traits::vector_size (rwork)); 
 
100
        assert ( uplo=='U' || uplo=='L' );
 
101
        assert ( jobz=='N' || jobz=='V' );
 
102
 
 
103
        int info; 
 
104
        detail::heev (jobz, uplo, n,
 
105
                     traits::matrix_storage (a), 
 
106
                     traits::leading_dimension (a),
 
107
                     traits::vector_storage (w),  
 
108
                     traits::vector_storage (work),
 
109
                     traits::vector_size (work),
 
110
                     traits::vector_storage (rwork),
 
111
                     info);
 
112
        return info; 
 
113
      }
 
114
    } // namespace detail
 
115
 
 
116
 
 
117
    // Function that allocates temporary arrays
 
118
    template <typename A, typename W>
 
119
    int heev (char jobz, char uplo, A& a, W& w, minimal_workspace ) {
 
120
       typedef typename A::value_type                              value_type ;
 
121
       typedef typename traits::type_traits<value_type>::real_type real_type ;
 
122
 
 
123
       int const n = traits::matrix_size1 (a);
 
124
 
 
125
       traits::detail::array<value_type> work( std::max(1,2*n-1) );
 
126
       traits::detail::array<real_type> rwork( std::max(3*n-1,1) );
 
127
 
 
128
       return detail::heev( jobz, uplo, a, w, work, rwork );
 
129
    }
 
130
 
 
131
 
 
132
    // Function that allocates temporary arrays
 
133
    template <typename A, typename W>
 
134
    int heev (char jobz, char uplo, A& a, W& w, optimal_workspace ) {
 
135
       typedef typename A::value_type                              value_type ;
 
136
       typedef typename traits::type_traits<value_type>::real_type real_type ;
 
137
 
 
138
       int const n = traits::matrix_size1 (a);
 
139
 
 
140
       traits::detail::array<value_type> work( std::max(1,33*n) );
 
141
       traits::detail::array<real_type> rwork( std::max(3*n-1,1) );
 
142
 
 
143
       return detail::heev( jobz, uplo, a, w, work, rwork );
 
144
    }
 
145
 
 
146
 
 
147
    // Function that uses given workarrays
 
148
    template <typename A, typename W, typename WC, typename WR>
 
149
    int heev (char jobz, char uplo, A& a, W& w, detail::workspace2<WC,WR> workspace ) {
 
150
       typedef typename A::value_type                              value_type ;
 
151
       typedef typename traits::type_traits<value_type>::real_type real_type ;
 
152
 
 
153
       return detail::heev( jobz, uplo, a, w, workspace.w_, workspace.wr_ );
 
154
    }
 
155
 
 
156
  }
 
157
 
 
158
}}}
 
159
 
 
160
#endif