1 //
2 // Copyright (c) 2002--2010
3 // Toon Knapen, Karl Meerbergen, Kresimir Fresl,
4 // Thomas Klimpel and Rutger ter Borg
5 //
6 // Distributed under the Boost Software License, Version 1.0.
7 // (See accompanying file LICENSE_1_0.txt or copy at
8 // http://www.boost.org/LICENSE_1_0.txt)
9 //
10 // THIS FILE IS AUTOMATICALLY GENERATED
11 // PLEASE DO NOT EDIT!
12 //
13
14 #ifndef BOOST_NUMERIC_BINDINGS_LAPACK_DRIVER_GESDD_HPP
15 #define BOOST_NUMERIC_BINDINGS_LAPACK_DRIVER_GESDD_HPP
16
17 #include <boost/assert.hpp>
18 #include <boost/numeric/bindings/begin.hpp>
19 #include <boost/numeric/bindings/detail/array.hpp>
20 #include <boost/numeric/bindings/is_column_major.hpp>
21 #include <boost/numeric/bindings/is_complex.hpp>
22 #include <boost/numeric/bindings/is_mutable.hpp>
23 #include <boost/numeric/bindings/is_real.hpp>
24 #include <boost/numeric/bindings/lapack/workspace.hpp>
25 #include <boost/numeric/bindings/remove_imaginary.hpp>
26 #include <boost/numeric/bindings/size.hpp>
27 #include <boost/numeric/bindings/stride.hpp>
28 #include <boost/numeric/bindings/traits/detail/utils.hpp>
29 #include <boost/numeric/bindings/value_type.hpp>
30 #include <boost/static_assert.hpp>
31 #include <boost/type_traits/is_same.hpp>
32 #include <boost/type_traits/remove_const.hpp>
33 #include <boost/utility/enable_if.hpp>
34
35 //
36 // The LAPACK-backend for gesdd is the netlib-compatible backend.
37 //
38 #include <boost/numeric/bindings/lapack/detail/lapack.h>
39 #include <boost/numeric/bindings/lapack/detail/lapack_option.hpp>
40
41 namespace boost {
42 namespace numeric {
43 namespace bindings {
44 namespace lapack {
45
46 //
47 // The detail namespace contains value-type-overloaded functions that
48 // dispatch to the appropriate back-end LAPACK-routine.
49 //
50 namespace detail {
51
52 //
53 // Overloaded function for dispatching to
54 // * netlib-compatible LAPACK backend (the default), and
55 // * float value-type.
56 //
gesdd(const char jobz,const fortran_int_t m,const fortran_int_t n,float * a,const fortran_int_t lda,float * s,float * u,const fortran_int_t ldu,float * vt,const fortran_int_t ldvt,float * work,const fortran_int_t lwork,fortran_int_t * iwork)57 inline std::ptrdiff_t gesdd( const char jobz, const fortran_int_t m,
58 const fortran_int_t n, float* a, const fortran_int_t lda, float* s,
59 float* u, const fortran_int_t ldu, float* vt,
60 const fortran_int_t ldvt, float* work, const fortran_int_t lwork,
61 fortran_int_t* iwork ) {
62 fortran_int_t info(0);
63 LAPACK_SGESDD( &jobz, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work,
64 &lwork, iwork, &info );
65 return info;
66 }
67
68 //
69 // Overloaded function for dispatching to
70 // * netlib-compatible LAPACK backend (the default), and
71 // * double value-type.
72 //
gesdd(const char jobz,const fortran_int_t m,const fortran_int_t n,double * a,const fortran_int_t lda,double * s,double * u,const fortran_int_t ldu,double * vt,const fortran_int_t ldvt,double * work,const fortran_int_t lwork,fortran_int_t * iwork)73 inline std::ptrdiff_t gesdd( const char jobz, const fortran_int_t m,
74 const fortran_int_t n, double* a, const fortran_int_t lda, double* s,
75 double* u, const fortran_int_t ldu, double* vt,
76 const fortran_int_t ldvt, double* work, const fortran_int_t lwork,
77 fortran_int_t* iwork ) {
78 fortran_int_t info(0);
79 LAPACK_DGESDD( &jobz, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work,
80 &lwork, iwork, &info );
81 return info;
82 }
83
84 //
85 // Overloaded function for dispatching to
86 // * netlib-compatible LAPACK backend (the default), and
87 // * complex<float> value-type.
88 //
gesdd(const char jobz,const fortran_int_t m,const fortran_int_t n,std::complex<float> * a,const fortran_int_t lda,float * s,std::complex<float> * u,const fortran_int_t ldu,std::complex<float> * vt,const fortran_int_t ldvt,std::complex<float> * work,const fortran_int_t lwork,float * rwork,fortran_int_t * iwork)89 inline std::ptrdiff_t gesdd( const char jobz, const fortran_int_t m,
90 const fortran_int_t n, std::complex<float>* a,
91 const fortran_int_t lda, float* s, std::complex<float>* u,
92 const fortran_int_t ldu, std::complex<float>* vt,
93 const fortran_int_t ldvt, std::complex<float>* work,
94 const fortran_int_t lwork, float* rwork, fortran_int_t* iwork ) {
95 fortran_int_t info(0);
96 LAPACK_CGESDD( &jobz, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work,
97 &lwork, rwork, iwork, &info );
98 return info;
99 }
100
101 //
102 // Overloaded function for dispatching to
103 // * netlib-compatible LAPACK backend (the default), and
104 // * complex<double> value-type.
105 //
gesdd(const char jobz,const fortran_int_t m,const fortran_int_t n,std::complex<double> * a,const fortran_int_t lda,double * s,std::complex<double> * u,const fortran_int_t ldu,std::complex<double> * vt,const fortran_int_t ldvt,std::complex<double> * work,const fortran_int_t lwork,double * rwork,fortran_int_t * iwork)106 inline std::ptrdiff_t gesdd( const char jobz, const fortran_int_t m,
107 const fortran_int_t n, std::complex<double>* a,
108 const fortran_int_t lda, double* s, std::complex<double>* u,
109 const fortran_int_t ldu, std::complex<double>* vt,
110 const fortran_int_t ldvt, std::complex<double>* work,
111 const fortran_int_t lwork, double* rwork, fortran_int_t* iwork ) {
112 fortran_int_t info(0);
113 LAPACK_ZGESDD( &jobz, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work,
114 &lwork, rwork, iwork, &info );
115 return info;
116 }
117
118 } // namespace detail
119
120 //
121 // Value-type based template class. Use this class if you need a type
122 // for dispatching to gesdd.
123 //
124 template< typename Value, typename Enable = void >
125 struct gesdd_impl {};
126
127 //
128 // This implementation is enabled if Value is a real type.
129 //
130 template< typename Value >
131 struct gesdd_impl< Value, typename boost::enable_if< is_real< Value > >::type > {
132
133 typedef Value value_type;
134 typedef typename remove_imaginary< Value >::type real_type;
135
136 //
137 // Static member function for user-defined workspaces, that
138 // * Deduces the required arguments for dispatching to LAPACK, and
139 // * Asserts that most arguments make sense.
140 //
141 template< typename MatrixA, typename VectorS, typename MatrixU,
142 typename MatrixVT, typename WORK, typename IWORK >
invokeboost::numeric::bindings::lapack::gesdd_impl143 static std::ptrdiff_t invoke( const char jobz, MatrixA& a, VectorS& s,
144 MatrixU& u, MatrixVT& vt, detail::workspace2< WORK,
145 IWORK > work ) {
146 namespace bindings = ::boost::numeric::bindings;
147 BOOST_STATIC_ASSERT( (bindings::is_column_major< MatrixA >::value) );
148 BOOST_STATIC_ASSERT( (bindings::is_column_major< MatrixU >::value) );
149 BOOST_STATIC_ASSERT( (bindings::is_column_major< MatrixVT >::value) );
150 BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
151 typename bindings::value_type< MatrixA >::type >::type,
152 typename remove_const< typename bindings::value_type<
153 VectorS >::type >::type >::value) );
154 BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
155 typename bindings::value_type< MatrixA >::type >::type,
156 typename remove_const< typename bindings::value_type<
157 MatrixU >::type >::type >::value) );
158 BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
159 typename bindings::value_type< MatrixA >::type >::type,
160 typename remove_const< typename bindings::value_type<
161 MatrixVT >::type >::type >::value) );
162 BOOST_STATIC_ASSERT( (bindings::is_mutable< MatrixA >::value) );
163 BOOST_STATIC_ASSERT( (bindings::is_mutable< VectorS >::value) );
164 BOOST_STATIC_ASSERT( (bindings::is_mutable< MatrixU >::value) );
165 BOOST_STATIC_ASSERT( (bindings::is_mutable< MatrixVT >::value) );
166 std::ptrdiff_t minmn = std::min< std::ptrdiff_t >( size_row(a),
167 size_column(a) );
168 BOOST_ASSERT( bindings::size(s) >= std::min<
169 std::ptrdiff_t >(bindings::size_row(a),
170 bindings::size_column(a)) );
171 BOOST_ASSERT( bindings::size(work.select(fortran_int_t())) >=
172 min_size_iwork( minmn ));
173 BOOST_ASSERT( bindings::size(work.select(real_type())) >=
174 min_size_work( bindings::size_row(a),
175 bindings::size_column(a), jobz, minmn ));
176 BOOST_ASSERT( bindings::size_column(a) >= 0 );
177 BOOST_ASSERT( bindings::size_minor(a) == 1 ||
178 bindings::stride_minor(a) == 1 );
179 BOOST_ASSERT( bindings::size_minor(u) == 1 ||
180 bindings::stride_minor(u) == 1 );
181 BOOST_ASSERT( bindings::size_minor(vt) == 1 ||
182 bindings::stride_minor(vt) == 1 );
183 BOOST_ASSERT( bindings::size_row(a) >= 0 );
184 BOOST_ASSERT( bindings::stride_major(a) >= std::max< std::ptrdiff_t >(1,
185 bindings::size_row(a)) );
186 BOOST_ASSERT( jobz == 'A' || jobz == 'S' || jobz == 'O' ||
187 jobz == 'N' );
188 return detail::gesdd( jobz, bindings::size_row(a),
189 bindings::size_column(a), bindings::begin_value(a),
190 bindings::stride_major(a), bindings::begin_value(s),
191 bindings::begin_value(u), bindings::stride_major(u),
192 bindings::begin_value(vt), bindings::stride_major(vt),
193 bindings::begin_value(work.select(real_type())),
194 bindings::size(work.select(real_type())),
195 bindings::begin_value(work.select(fortran_int_t())) );
196 }
197
198 //
199 // Static member function that
200 // * Figures out the minimal workspace requirements, and passes
201 // the results to the user-defined workspace overload of the
202 // invoke static member function
203 // * Enables the unblocked algorithm (BLAS level 2)
204 //
205 template< typename MatrixA, typename VectorS, typename MatrixU,
206 typename MatrixVT >
invokeboost::numeric::bindings::lapack::gesdd_impl207 static std::ptrdiff_t invoke( const char jobz, MatrixA& a, VectorS& s,
208 MatrixU& u, MatrixVT& vt, minimal_workspace ) {
209 namespace bindings = ::boost::numeric::bindings;
210 std::ptrdiff_t minmn = std::min< std::ptrdiff_t >( size_row(a),
211 size_column(a) );
212 bindings::detail::array< real_type > tmp_work( min_size_work(
213 bindings::size_row(a), bindings::size_column(a), jobz,
214 minmn ) );
215 bindings::detail::array< fortran_int_t > tmp_iwork(
216 min_size_iwork( minmn ) );
217 return invoke( jobz, a, s, u, vt, workspace( tmp_work, tmp_iwork ) );
218 }
219
220 //
221 // Static member function that
222 // * Figures out the optimal workspace requirements, and passes
223 // the results to the user-defined workspace overload of the
224 // invoke static member
225 // * Enables the blocked algorithm (BLAS level 3)
226 //
227 template< typename MatrixA, typename VectorS, typename MatrixU,
228 typename MatrixVT >
invokeboost::numeric::bindings::lapack::gesdd_impl229 static std::ptrdiff_t invoke( const char jobz, MatrixA& a, VectorS& s,
230 MatrixU& u, MatrixVT& vt, optimal_workspace ) {
231 namespace bindings = ::boost::numeric::bindings;
232 return invoke( jobz, a, s, u, vt, minimal_workspace() );
233 }
234
235 //
236 // Static member function that returns the minimum size of
237 // workspace-array work.
238 //
min_size_workboost::numeric::bindings::lapack::gesdd_impl239 static std::ptrdiff_t min_size_work( const std::ptrdiff_t m,
240 const std::ptrdiff_t n, const char jobz,
241 const std::ptrdiff_t minmn ) {
242 if ( n == 0 ) return 1;
243 if ( jobz == 'N' ) return 3*minmn + std::max<
244 std::ptrdiff_t >( std::max< std::ptrdiff_t >(m,n), 7*minmn );
245 if ( jobz == 'O' ) return 3*minmn*minmn + std::max<
246 std::ptrdiff_t >( std::max< std::ptrdiff_t >( m,n ),
247 5*minmn*minmn + 4*minmn );
248 return 3*minmn*minmn + std::max< std::ptrdiff_t >( std::max<
249 std::ptrdiff_t >( m,n ), 4*minmn*minmn + 4*minmn );
250 }
251
252 //
253 // Static member function that returns the minimum size of
254 // workspace-array iwork.
255 //
min_size_iworkboost::numeric::bindings::lapack::gesdd_impl256 static std::ptrdiff_t min_size_iwork( const std::ptrdiff_t minmn ) {
257 return 8*minmn;
258 }
259 };
260
261 //
262 // This implementation is enabled if Value is a complex type.
263 //
264 template< typename Value >
265 struct gesdd_impl< Value, typename boost::enable_if< is_complex< Value > >::type > {
266
267 typedef Value value_type;
268 typedef typename remove_imaginary< Value >::type real_type;
269
270 //
271 // Static member function for user-defined workspaces, that
272 // * Deduces the required arguments for dispatching to LAPACK, and
273 // * Asserts that most arguments make sense.
274 //
275 template< typename MatrixA, typename VectorS, typename MatrixU,
276 typename MatrixVT, typename WORK, typename RWORK, typename IWORK >
invokeboost::numeric::bindings::lapack::gesdd_impl277 static std::ptrdiff_t invoke( const char jobz, MatrixA& a, VectorS& s,
278 MatrixU& u, MatrixVT& vt, detail::workspace3< WORK, RWORK,
279 IWORK > work ) {
280 namespace bindings = ::boost::numeric::bindings;
281 BOOST_STATIC_ASSERT( (bindings::is_column_major< MatrixA >::value) );
282 BOOST_STATIC_ASSERT( (bindings::is_column_major< MatrixU >::value) );
283 BOOST_STATIC_ASSERT( (bindings::is_column_major< MatrixVT >::value) );
284 BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
285 typename bindings::value_type< MatrixA >::type >::type,
286 typename remove_const< typename bindings::value_type<
287 MatrixU >::type >::type >::value) );
288 BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
289 typename bindings::value_type< MatrixA >::type >::type,
290 typename remove_const< typename bindings::value_type<
291 MatrixVT >::type >::type >::value) );
292 BOOST_STATIC_ASSERT( (bindings::is_mutable< MatrixA >::value) );
293 BOOST_STATIC_ASSERT( (bindings::is_mutable< VectorS >::value) );
294 BOOST_STATIC_ASSERT( (bindings::is_mutable< MatrixU >::value) );
295 BOOST_STATIC_ASSERT( (bindings::is_mutable< MatrixVT >::value) );
296 std::ptrdiff_t minmn = std::min< std::ptrdiff_t >( size_row(a),
297 size_column(a) );
298 BOOST_ASSERT( bindings::size(s) >= std::min<
299 std::ptrdiff_t >(bindings::size_row(a),
300 bindings::size_column(a)) );
301 BOOST_ASSERT( bindings::size(work.select(fortran_int_t())) >=
302 min_size_iwork( minmn ));
303 BOOST_ASSERT( bindings::size(work.select(real_type())) >=
304 min_size_rwork( minmn, jobz ));
305 BOOST_ASSERT( bindings::size(work.select(value_type())) >=
306 min_size_work( bindings::size_row(a),
307 bindings::size_column(a), jobz, minmn ));
308 BOOST_ASSERT( bindings::size_column(a) >= 0 );
309 BOOST_ASSERT( bindings::size_minor(a) == 1 ||
310 bindings::stride_minor(a) == 1 );
311 BOOST_ASSERT( bindings::size_minor(u) == 1 ||
312 bindings::stride_minor(u) == 1 );
313 BOOST_ASSERT( bindings::size_minor(vt) == 1 ||
314 bindings::stride_minor(vt) == 1 );
315 BOOST_ASSERT( bindings::size_row(a) >= 0 );
316 BOOST_ASSERT( bindings::stride_major(a) >= std::max< std::ptrdiff_t >(1,
317 bindings::size_row(a)) );
318 BOOST_ASSERT( jobz == 'A' || jobz == 'S' || jobz == 'O' ||
319 jobz == 'N' );
320 return detail::gesdd( jobz, bindings::size_row(a),
321 bindings::size_column(a), bindings::begin_value(a),
322 bindings::stride_major(a), bindings::begin_value(s),
323 bindings::begin_value(u), bindings::stride_major(u),
324 bindings::begin_value(vt), bindings::stride_major(vt),
325 bindings::begin_value(work.select(value_type())),
326 bindings::size(work.select(value_type())),
327 bindings::begin_value(work.select(real_type())),
328 bindings::begin_value(work.select(fortran_int_t())) );
329 }
330
331 //
332 // Static member function that
333 // * Figures out the minimal workspace requirements, and passes
334 // the results to the user-defined workspace overload of the
335 // invoke static member function
336 // * Enables the unblocked algorithm (BLAS level 2)
337 //
338 template< typename MatrixA, typename VectorS, typename MatrixU,
339 typename MatrixVT >
invokeboost::numeric::bindings::lapack::gesdd_impl340 static std::ptrdiff_t invoke( const char jobz, MatrixA& a, VectorS& s,
341 MatrixU& u, MatrixVT& vt, minimal_workspace ) {
342 namespace bindings = ::boost::numeric::bindings;
343 std::ptrdiff_t minmn = std::min< std::ptrdiff_t >( size_row(a),
344 size_column(a) );
345 bindings::detail::array< value_type > tmp_work( min_size_work(
346 bindings::size_row(a), bindings::size_column(a), jobz,
347 minmn ) );
348 bindings::detail::array< real_type > tmp_rwork( min_size_rwork( minmn,
349 jobz ) );
350 bindings::detail::array< fortran_int_t > tmp_iwork(
351 min_size_iwork( minmn ) );
352 return invoke( jobz, a, s, u, vt, workspace( tmp_work, tmp_rwork,
353 tmp_iwork ) );
354 }
355
356 //
357 // Static member function that
358 // * Figures out the optimal workspace requirements, and passes
359 // the results to the user-defined workspace overload of the
360 // invoke static member
361 // * Enables the blocked algorithm (BLAS level 3)
362 //
363 template< typename MatrixA, typename VectorS, typename MatrixU,
364 typename MatrixVT >
invokeboost::numeric::bindings::lapack::gesdd_impl365 static std::ptrdiff_t invoke( const char jobz, MatrixA& a, VectorS& s,
366 MatrixU& u, MatrixVT& vt, optimal_workspace ) {
367 namespace bindings = ::boost::numeric::bindings;
368 std::ptrdiff_t minmn = std::min< std::ptrdiff_t >( size_row(a),
369 size_column(a) );
370 value_type opt_size_work;
371 bindings::detail::array< real_type > tmp_rwork( min_size_rwork( minmn,
372 jobz ) );
373 bindings::detail::array< fortran_int_t > tmp_iwork(
374 min_size_iwork( minmn ) );
375 detail::gesdd( jobz, bindings::size_row(a),
376 bindings::size_column(a), bindings::begin_value(a),
377 bindings::stride_major(a), bindings::begin_value(s),
378 bindings::begin_value(u), bindings::stride_major(u),
379 bindings::begin_value(vt), bindings::stride_major(vt),
380 &opt_size_work, -1, bindings::begin_value(tmp_rwork),
381 bindings::begin_value(tmp_iwork) );
382 bindings::detail::array< value_type > tmp_work(
383 traits::detail::to_int( opt_size_work ) );
384 return invoke( jobz, a, s, u, vt, workspace( tmp_work, tmp_rwork,
385 tmp_iwork ) );
386 }
387
388 //
389 // Static member function that returns the minimum size of
390 // workspace-array work.
391 //
min_size_workboost::numeric::bindings::lapack::gesdd_impl392 static std::ptrdiff_t min_size_work( const std::ptrdiff_t m,
393 const std::ptrdiff_t n, const char jobz,
394 const std::ptrdiff_t minmn ) {
395 if ( n == 0 ) return 1;
396 if ( jobz == 'N' ) return 2*minmn + std::max< std::ptrdiff_t >( m,n );
397 if ( jobz == 'O' ) return 2*(minmn*minmn + minmn) + std::max<
398 std::ptrdiff_t >( m, n );
399 return minmn*minmn + 2*minmn + std::max< std::ptrdiff_t >( m, n );
400 }
401
402 //
403 // Static member function that returns the minimum size of
404 // workspace-array rwork.
405 //
min_size_rworkboost::numeric::bindings::lapack::gesdd_impl406 static std::ptrdiff_t min_size_rwork( const std::ptrdiff_t minmn,
407 const char jobz ) {
408 if ( jobz == 'N' ) return 5*minmn;
409 return 5*minmn*minmn + 7*minmn;
410 }
411
412 //
413 // Static member function that returns the minimum size of
414 // workspace-array iwork.
415 //
min_size_iworkboost::numeric::bindings::lapack::gesdd_impl416 static std::ptrdiff_t min_size_iwork( const std::ptrdiff_t minmn ) {
417 return 8*minmn;
418 }
419 };
420
421
422 //
423 // Functions for direct use. These functions are overloaded for temporaries,
424 // so that wrapped types can still be passed and used for write-access. In
425 // addition, if applicable, they are overloaded for user-defined workspaces.
426 // Calls to these functions are passed to the gesdd_impl classes. In the
427 // documentation, most overloads are collapsed to avoid a large number of
428 // prototypes which are very similar.
429 //
430
431 //
432 // Overloaded function for gesdd. Its overload differs for
433 // * User-defined workspace
434 //
435 template< typename MatrixA, typename VectorS, typename MatrixU,
436 typename MatrixVT, typename Workspace >
437 inline typename boost::enable_if< detail::is_workspace< Workspace >,
438 std::ptrdiff_t >::type
gesdd(const char jobz,MatrixA & a,VectorS & s,MatrixU & u,MatrixVT & vt,Workspace work)439 gesdd( const char jobz, MatrixA& a, VectorS& s, MatrixU& u, MatrixVT& vt,
440 Workspace work ) {
441 return gesdd_impl< typename bindings::value_type<
442 MatrixA >::type >::invoke( jobz, a, s, u, vt, work );
443 }
444
445 //
446 // Overloaded function for gesdd. Its overload differs for
447 // * Default workspace-type (optimal)
448 //
449 template< typename MatrixA, typename VectorS, typename MatrixU,
450 typename MatrixVT >
451 inline typename boost::disable_if< detail::is_workspace< MatrixVT >,
452 std::ptrdiff_t >::type
gesdd(const char jobz,MatrixA & a,VectorS & s,MatrixU & u,MatrixVT & vt)453 gesdd( const char jobz, MatrixA& a, VectorS& s, MatrixU& u,
454 MatrixVT& vt ) {
455 return gesdd_impl< typename bindings::value_type<
456 MatrixA >::type >::invoke( jobz, a, s, u, vt,
457 optimal_workspace() );
458 }
459
460 } // namespace lapack
461 } // namespace bindings
462 } // namespace numeric
463 } // namespace boost
464
465 #endif
466