3
* Copyright (c) Toon Knapen, Karl Meerbergen & Kresimir Fresl 2003
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)
9
* KF acknowledges the support of the Faculty of Civil Engineering,
10
* University of Zagreb, Croatia.
14
#ifndef BOOST_NUMERIC_BINDINGS_LAPACK_HBEV_HPP
15
#define BOOST_NUMERIC_BINDINGS_LAPACK_HBEV_HPP
17
#include <boost/numeric/bindings/traits/type.hpp>
18
#include <boost/numeric/bindings/traits/traits.hpp>
19
#include <boost/numeric/bindings/traits/type_traits.hpp>
20
#include <boost/numeric/bindings/lapack/lapack.h>
21
#include <boost/numeric/bindings/lapack/workspace.hpp>
22
#include <boost/numeric/bindings/traits/detail/array.hpp>
23
#include <boost/numeric/bindings/traits/detail/utils.hpp>
25
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
26
# include <boost/static_assert.hpp>
27
# include <boost/type_traits.hpp>
31
namespace boost { namespace numeric { namespace bindings {
35
///////////////////////////////////////////////////////////////////
37
// Eigendecomposition of a banded Hermitian matrix.
39
///////////////////////////////////////////////////////////////////
42
* hbev() computes the eigenvalues and optionally the associated
43
* eigenvectors of a banded Hermitian matrix A. A matrix is Hermitian
44
* when herm( A ) == A. When A is real, a Hermitian matrix is also
47
* The eigen decomposition is A = U S * herm(U) where U is a
48
* unitary matrix and S is a diagonal matrix. The eigenvalues of A
49
* are on the main diagonal of S. The eigenvalues are real.
51
* Workspace is organized following the arguments in the calling sequence.
52
* optimal_workspace() : for optimizing use of blas 3 kernels
53
* minimal_workspace() : minimum size of workarrays, but does not allow for optimization
55
* workspace( work ) for real matrices where work is a real array with
56
* vector_size( work ) >= 3*matrix_size1( a ) - 2
57
* workspace( work, rwork ) for complex matrices where work is a complex
58
* array with vector_size( work ) >= matrix_size1( a )
59
* and rwork is a real array with
60
* vector_size( rwork ) >= 3 * matrix_size1( a ) - 2.
64
* If uplo=='L' only the lower triangular part is stored.
65
* If uplo=='U' only the upper triangular part is stored.
67
* The matrix is assumed to be stored in LAPACK band format, i.e.
68
* matrices are stored columnwise, in a compressed format so that when e.g. uplo=='U'
69
* the (i,j) element with j>=i is in position (i-j) + j * (KD+1) + KD where KD is the
70
* half bandwidth of the matrix. For a triadiagonal matrix, KD=1, for a diagonal matrix
72
* When uplo=='L', the (i,j) element with j>=i is in position (i-j) + j * (KD+1).
74
* The matrix A is thus a rectangular matrix with KD+1 rows and N columns.
79
void hbev (char const jobz, char const uplo, int const n, int const kd,
80
float* ab, int const ldab, float* w, float* z, int const ldz,
81
float* work, int& info)
83
//for (int i=0; i<n*kd; ++i) std::cout << *(ab+i) << " " ;
85
LAPACK_SSBEV (&jobz, &uplo, &n, &kd, ab, &ldab, w, z, &ldz,
90
void hbev (char const jobz, char const uplo, int const n, int const kd,
91
double* ab, int const ldab, double* w, double* z, int const ldz,
92
double* work, int& info)
94
LAPACK_DSBEV (&jobz, &uplo, &n, &kd, ab, &ldab, w, z, &ldz,
99
void hbev (char const jobz, char const uplo, int const n, int const kd,
100
traits::complex_f* ab, int const ldab, float* w,
101
traits::complex_f* z, int const ldz,
102
traits::complex_f* work, float* rwork, int& info)
104
LAPACK_CHBEV (&jobz, &uplo, &n, &kd, traits::complex_ptr(ab), &ldab,
105
w, traits::complex_ptr(z), &ldz,
106
traits::complex_ptr(work), rwork, &info);
110
void hbev (char const jobz, char const uplo, int const n, int const kd,
111
traits::complex_d* ab, int const ldab, double* w,
112
traits::complex_d* z, int const ldz,
113
traits::complex_d* work, double* rwork, int& info)
115
LAPACK_ZHBEV (&jobz, &uplo, &n, &kd, traits::complex_ptr(ab), &ldab,
116
w, traits::complex_ptr(z), &ldz,
117
traits::complex_ptr(work), rwork, &info);
127
/// Handling of workspace in the case of one workarray.
130
template <typename T, typename R>
131
void operator() (char const jobz, char const uplo, int const n,
132
int const kd, T* ab, int const ldab, R* w, T* z,
133
int const ldz, minimal_workspace , int& info ) const {
134
traits::detail::array<T> work( 3*n-2 );
135
hbev( jobz, uplo, n, kd, ab, ldab, w, z, ldz,
136
traits::vector_storage( work ),
140
template <typename T, typename R>
141
void operator() (char const jobz, char const uplo, int const n,
142
int const kd, T* ab, int const ldab, R* w, T* z,
143
int const ldz, optimal_workspace , int& info ) const {
144
traits::detail::array<T> work( 3*n-2 );
146
hbev( jobz, uplo, n, kd, ab, ldab, w, z, ldz,
147
traits::vector_storage( work ),
151
template <typename T, typename R, typename W>
152
void operator() (char const jobz, char const uplo, int const n,
153
int const kd, T* ab, int const ldab, R* w, T* z,
154
int const ldz, detail::workspace1<W> work,
156
assert( traits::vector_size( work.w_ ) >= 3*n-2 );
158
hbev( jobz, uplo, n, kd, ab, ldab, w, z, ldz,
159
traits::vector_storage( work.w_ ),
165
/// Handling of workspace in the case of two workarrays.
168
template <typename T, typename R>
169
void operator() (char const jobz, char const uplo, int const n,
170
int const kd, T* ab, int const ldab, R* w, T* z,
171
int const ldz, minimal_workspace , int& info ) const {
172
traits::detail::array<T> work( n );
173
traits::detail::array<R> rwork( 3*n-2 );
175
hbev( jobz, uplo, n, kd, ab, ldab, w, z, ldz,
176
traits::vector_storage( work ),
177
traits::vector_storage( rwork ),
181
template <typename T, typename R>
182
void operator() (char const jobz, char const uplo, int const n,
183
int const kd, T* ab, int const ldab, R* w, T* z,
184
int const ldz, optimal_workspace , int& info ) const {
185
traits::detail::array<T> work( n );
186
traits::detail::array<R> rwork( 3*n-2 );
188
hbev( jobz, uplo, n, kd, ab, ldab, w, z, ldz,
189
traits::vector_storage( work ),
190
traits::vector_storage( rwork ),
194
template <typename T, typename R, typename W, typename RW>
195
void operator() (char const jobz, char const uplo, int const n,
196
int const kd, T* ab, int const ldab, R* w, T* z,
197
int const ldz, detail::workspace2<W,RW> work,
199
assert( traits::vector_size( work.wr_ ) >= 3*n-2 );
200
assert( traits::vector_size( work.w_ ) >= n );
202
hbev( jobz, uplo, n, kd, ab, ldab, w, z, ldz,
203
traits::vector_storage( work.w_ ),
204
traits::vector_storage( work.wr_ ),
211
/// Compute eigendecomposition of the banded Hermitian matrix ab.
212
/// if jobz=='N' only the eigenvalues are computed.
213
/// if jobz=='V' compute the eigenvalues a and the eigenvectors.
215
/// Workspace is organized following the arguments in the calling sequence.
216
/// optimal_workspace() : for optimizing use of blas 3 kernels
217
/// minimal_workspace() : minimum size of workarrays, but does not allow for optimization
218
/// of blas 3 kernels
219
/// workspace( work ) for real matrices where work is a real array with
220
/// vector_size( work ) >= 3*matrix_size1( a )-2
221
/// workspace( work, rwork ) for complex matrices where work is a complex
222
/// array with vector_size( work ) >= matrix_size1( a )
223
/// and rwork is a real array with
224
/// vector_size( rwork ) >= 3*matrix_size1( a )-2.
225
template <typename AB, typename Z, typename W, typename Work>
226
int hbev( char const jobz, AB& ab, W& w, Z& z, Work work ) {
227
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
228
BOOST_STATIC_ASSERT((boost::is_same<
229
typename traits::matrix_traits<AB>::matrix_structure,
234
typedef typename AB::value_type value_type ;
236
int const n = traits::matrix_size2 (ab);
237
assert (n == traits::matrix_size1 (z));
238
assert (n == traits::vector_size (w));
239
assert ( jobz=='N' || jobz=='V' );
242
detail::Hbev< n_workspace_args<value_type>::value >() (jobz,
243
traits::matrix_uplo_tag( ab ), n,
244
traits::matrix_upper_bandwidth(ab),
245
traits::matrix_storage (ab),
246
traits::leading_dimension (ab),
247
traits::vector_storage (w),
248
traits::matrix_storage (z),
249
traits::leading_dimension (z),
254
} // namespace detail
257
/// Compute eigendecomposition without eigenvectors
258
template <typename AB, typename W, typename Work>
260
int hbev (AB& ab, W& w, Work work) {
261
return detail::hbev( 'N', ab, w, ab, work );
265
/// Compute eigendecomposition with eigenvectors
266
template <typename AB, typename W, typename Z, typename Work>
268
int hbev (AB& ab, W& w, Z& z, Work work) {
269
BOOST_STATIC_ASSERT((boost::is_same<
270
typename traits::matrix_traits<Z>::matrix_structure,
273
int const n = traits::matrix_size2 (ab);
274
assert (n == traits::matrix_size1 (z));
275
assert (n == traits::matrix_size2 (z));
276
return detail::hbev( 'V', ab, w, z, work );
added
removed
src/3rdParty/boost/numeric/bindings/lapack/hbev.hpp
