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_GGES_HPP
15 #define BOOST_NUMERIC_BINDINGS_LAPACK_DRIVER_GGES_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 gges 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 //
gges(const char jobvsl,const char jobvsr,const char sort,external_fp selctg,const fortran_int_t n,float * a,const fortran_int_t lda,float * b,const fortran_int_t ldb,fortran_int_t & sdim,float * alphar,float * alphai,float * beta,float * vsl,const fortran_int_t ldvsl,float * vsr,const fortran_int_t ldvsr,float * work,const fortran_int_t lwork,fortran_bool_t * bwork)57 inline std::ptrdiff_t gges( const char jobvsl, const char jobvsr,
58 const char sort, external_fp selctg, const fortran_int_t n, float* a,
59 const fortran_int_t lda, float* b, const fortran_int_t ldb,
60 fortran_int_t& sdim, float* alphar, float* alphai, float* beta,
61 float* vsl, const fortran_int_t ldvsl, float* vsr,
62 const fortran_int_t ldvsr, float* work, const fortran_int_t lwork,
63 fortran_bool_t* bwork ) {
64 fortran_int_t info(0);
65 LAPACK_SGGES( &jobvsl, &jobvsr, &sort, selctg, &n, a, &lda, b, &ldb,
66 &sdim, alphar, alphai, beta, vsl, &ldvsl, vsr, &ldvsr, work,
67 &lwork, bwork, &info );
68 return info;
69 }
70
71 //
72 // Overloaded function for dispatching to
73 // * netlib-compatible LAPACK backend (the default), and
74 // * double value-type.
75 //
gges(const char jobvsl,const char jobvsr,const char sort,external_fp selctg,const fortran_int_t n,double * a,const fortran_int_t lda,double * b,const fortran_int_t ldb,fortran_int_t & sdim,double * alphar,double * alphai,double * beta,double * vsl,const fortran_int_t ldvsl,double * vsr,const fortran_int_t ldvsr,double * work,const fortran_int_t lwork,fortran_bool_t * bwork)76 inline std::ptrdiff_t gges( const char jobvsl, const char jobvsr,
77 const char sort, external_fp selctg, const fortran_int_t n, double* a,
78 const fortran_int_t lda, double* b, const fortran_int_t ldb,
79 fortran_int_t& sdim, double* alphar, double* alphai, double* beta,
80 double* vsl, const fortran_int_t ldvsl, double* vsr,
81 const fortran_int_t ldvsr, double* work, const fortran_int_t lwork,
82 fortran_bool_t* bwork ) {
83 fortran_int_t info(0);
84 LAPACK_DGGES( &jobvsl, &jobvsr, &sort, selctg, &n, a, &lda, b, &ldb,
85 &sdim, alphar, alphai, beta, vsl, &ldvsl, vsr, &ldvsr, work,
86 &lwork, bwork, &info );
87 return info;
88 }
89
90 //
91 // Overloaded function for dispatching to
92 // * netlib-compatible LAPACK backend (the default), and
93 // * complex<float> value-type.
94 //
gges(const char jobvsl,const char jobvsr,const char sort,external_fp selctg,const fortran_int_t n,std::complex<float> * a,const fortran_int_t lda,std::complex<float> * b,const fortran_int_t ldb,fortran_int_t & sdim,std::complex<float> * alpha,std::complex<float> * beta,std::complex<float> * vsl,const fortran_int_t ldvsl,std::complex<float> * vsr,const fortran_int_t ldvsr,std::complex<float> * work,const fortran_int_t lwork,float * rwork,fortran_bool_t * bwork)95 inline std::ptrdiff_t gges( const char jobvsl, const char jobvsr,
96 const char sort, external_fp selctg, const fortran_int_t n,
97 std::complex<float>* a, const fortran_int_t lda,
98 std::complex<float>* b, const fortran_int_t ldb, fortran_int_t& sdim,
99 std::complex<float>* alpha, std::complex<float>* beta,
100 std::complex<float>* vsl, const fortran_int_t ldvsl,
101 std::complex<float>* vsr, const fortran_int_t ldvsr,
102 std::complex<float>* work, const fortran_int_t lwork, float* rwork,
103 fortran_bool_t* bwork ) {
104 fortran_int_t info(0);
105 LAPACK_CGGES( &jobvsl, &jobvsr, &sort, selctg, &n, a, &lda, b, &ldb,
106 &sdim, alpha, beta, vsl, &ldvsl, vsr, &ldvsr, work, &lwork, rwork,
107 bwork, &info );
108 return info;
109 }
110
111 //
112 // Overloaded function for dispatching to
113 // * netlib-compatible LAPACK backend (the default), and
114 // * complex<double> value-type.
115 //
gges(const char jobvsl,const char jobvsr,const char sort,external_fp selctg,const fortran_int_t n,std::complex<double> * a,const fortran_int_t lda,std::complex<double> * b,const fortran_int_t ldb,fortran_int_t & sdim,std::complex<double> * alpha,std::complex<double> * beta,std::complex<double> * vsl,const fortran_int_t ldvsl,std::complex<double> * vsr,const fortran_int_t ldvsr,std::complex<double> * work,const fortran_int_t lwork,double * rwork,fortran_bool_t * bwork)116 inline std::ptrdiff_t gges( const char jobvsl, const char jobvsr,
117 const char sort, external_fp selctg, const fortran_int_t n,
118 std::complex<double>* a, const fortran_int_t lda,
119 std::complex<double>* b, const fortran_int_t ldb, fortran_int_t& sdim,
120 std::complex<double>* alpha, std::complex<double>* beta,
121 std::complex<double>* vsl, const fortran_int_t ldvsl,
122 std::complex<double>* vsr, const fortran_int_t ldvsr,
123 std::complex<double>* work, const fortran_int_t lwork, double* rwork,
124 fortran_bool_t* bwork ) {
125 fortran_int_t info(0);
126 LAPACK_ZGGES( &jobvsl, &jobvsr, &sort, selctg, &n, a, &lda, b, &ldb,
127 &sdim, alpha, beta, vsl, &ldvsl, vsr, &ldvsr, work, &lwork, rwork,
128 bwork, &info );
129 return info;
130 }
131
132 } // namespace detail
133
134 //
135 // Value-type based template class. Use this class if you need a type
136 // for dispatching to gges.
137 //
138 template< typename Value, typename Enable = void >
139 struct gges_impl {};
140
141 //
142 // This implementation is enabled if Value is a real type.
143 //
144 template< typename Value >
145 struct gges_impl< Value, typename boost::enable_if< is_real< Value > >::type > {
146
147 typedef Value value_type;
148 typedef typename remove_imaginary< Value >::type real_type;
149
150 //
151 // Static member function for user-defined workspaces, that
152 // * Deduces the required arguments for dispatching to LAPACK, and
153 // * Asserts that most arguments make sense.
154 //
155 template< typename MatrixA, typename MatrixB, typename VectorALPHAR,
156 typename VectorALPHAI, typename VectorBETA, typename MatrixVSL,
157 typename MatrixVSR, typename WORK, typename BWORK >
invokeboost::numeric::bindings::lapack::gges_impl158 static std::ptrdiff_t invoke( const char jobvsl, const char jobvsr,
159 const char sort, external_fp selctg, MatrixA& a, MatrixB& b,
160 fortran_int_t& sdim, VectorALPHAR& alphar,
161 VectorALPHAI& alphai, VectorBETA& beta, MatrixVSL& vsl,
162 MatrixVSR& vsr, detail::workspace2< WORK, BWORK > work ) {
163 namespace bindings = ::boost::numeric::bindings;
164 BOOST_STATIC_ASSERT( (bindings::is_column_major< MatrixA >::value) );
165 BOOST_STATIC_ASSERT( (bindings::is_column_major< MatrixB >::value) );
166 BOOST_STATIC_ASSERT( (bindings::is_column_major< MatrixVSL >::value) );
167 BOOST_STATIC_ASSERT( (bindings::is_column_major< MatrixVSR >::value) );
168 BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
169 typename bindings::value_type< MatrixA >::type >::type,
170 typename remove_const< typename bindings::value_type<
171 MatrixB >::type >::type >::value) );
172 BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
173 typename bindings::value_type< MatrixA >::type >::type,
174 typename remove_const< typename bindings::value_type<
175 VectorALPHAR >::type >::type >::value) );
176 BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
177 typename bindings::value_type< MatrixA >::type >::type,
178 typename remove_const< typename bindings::value_type<
179 VectorALPHAI >::type >::type >::value) );
180 BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
181 typename bindings::value_type< MatrixA >::type >::type,
182 typename remove_const< typename bindings::value_type<
183 VectorBETA >::type >::type >::value) );
184 BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
185 typename bindings::value_type< MatrixA >::type >::type,
186 typename remove_const< typename bindings::value_type<
187 MatrixVSL >::type >::type >::value) );
188 BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
189 typename bindings::value_type< MatrixA >::type >::type,
190 typename remove_const< typename bindings::value_type<
191 MatrixVSR >::type >::type >::value) );
192 BOOST_STATIC_ASSERT( (bindings::is_mutable< MatrixA >::value) );
193 BOOST_STATIC_ASSERT( (bindings::is_mutable< MatrixB >::value) );
194 BOOST_STATIC_ASSERT( (bindings::is_mutable< VectorALPHAR >::value) );
195 BOOST_STATIC_ASSERT( (bindings::is_mutable< VectorALPHAI >::value) );
196 BOOST_STATIC_ASSERT( (bindings::is_mutable< VectorBETA >::value) );
197 BOOST_STATIC_ASSERT( (bindings::is_mutable< MatrixVSL >::value) );
198 BOOST_STATIC_ASSERT( (bindings::is_mutable< MatrixVSR >::value) );
199 BOOST_ASSERT( bindings::size(alphai) >= bindings::size_column(a) );
200 BOOST_ASSERT( bindings::size(alphar) >= bindings::size_column(a) );
201 BOOST_ASSERT( bindings::size(work.select(fortran_bool_t())) >=
202 min_size_bwork( bindings::size_column(a), sort ));
203 BOOST_ASSERT( bindings::size(work.select(real_type())) >=
204 min_size_work( bindings::size_column(a) ));
205 BOOST_ASSERT( bindings::size_column(a) >= 0 );
206 BOOST_ASSERT( bindings::size_minor(a) == 1 ||
207 bindings::stride_minor(a) == 1 );
208 BOOST_ASSERT( bindings::size_minor(b) == 1 ||
209 bindings::stride_minor(b) == 1 );
210 BOOST_ASSERT( bindings::size_minor(vsl) == 1 ||
211 bindings::stride_minor(vsl) == 1 );
212 BOOST_ASSERT( bindings::size_minor(vsr) == 1 ||
213 bindings::stride_minor(vsr) == 1 );
214 BOOST_ASSERT( bindings::stride_major(a) >= std::max< std::ptrdiff_t >(1,
215 bindings::size_column(a)) );
216 BOOST_ASSERT( bindings::stride_major(b) >= std::max< std::ptrdiff_t >(1,
217 bindings::size_column(a)) );
218 BOOST_ASSERT( jobvsl == 'N' || jobvsl == 'V' );
219 BOOST_ASSERT( jobvsr == 'N' || jobvsr == 'V' );
220 BOOST_ASSERT( sort == 'N' || sort == 'S' );
221 return detail::gges( jobvsl, jobvsr, sort, selctg,
222 bindings::size_column(a), bindings::begin_value(a),
223 bindings::stride_major(a), bindings::begin_value(b),
224 bindings::stride_major(b), sdim,
225 bindings::begin_value(alphar), bindings::begin_value(alphai),
226 bindings::begin_value(beta), bindings::begin_value(vsl),
227 bindings::stride_major(vsl), bindings::begin_value(vsr),
228 bindings::stride_major(vsr),
229 bindings::begin_value(work.select(real_type())),
230 bindings::size(work.select(real_type())),
231 bindings::begin_value(work.select(fortran_bool_t())) );
232 }
233
234 //
235 // Static member function that
236 // * Figures out the minimal workspace requirements, and passes
237 // the results to the user-defined workspace overload of the
238 // invoke static member function
239 // * Enables the unblocked algorithm (BLAS level 2)
240 //
241 template< typename MatrixA, typename MatrixB, typename VectorALPHAR,
242 typename VectorALPHAI, typename VectorBETA, typename MatrixVSL,
243 typename MatrixVSR >
invokeboost::numeric::bindings::lapack::gges_impl244 static std::ptrdiff_t invoke( const char jobvsl, const char jobvsr,
245 const char sort, external_fp selctg, MatrixA& a, MatrixB& b,
246 fortran_int_t& sdim, VectorALPHAR& alphar,
247 VectorALPHAI& alphai, VectorBETA& beta, MatrixVSL& vsl,
248 MatrixVSR& vsr, minimal_workspace ) {
249 namespace bindings = ::boost::numeric::bindings;
250 bindings::detail::array< real_type > tmp_work( min_size_work(
251 bindings::size_column(a) ) );
252 bindings::detail::array< fortran_bool_t > tmp_bwork( min_size_bwork(
253 bindings::size_column(a), sort ) );
254 return invoke( jobvsl, jobvsr, sort, selctg, a, b, sdim, alphar,
255 alphai, beta, vsl, vsr, workspace( tmp_work, tmp_bwork ) );
256 }
257
258 //
259 // Static member function that
260 // * Figures out the optimal workspace requirements, and passes
261 // the results to the user-defined workspace overload of the
262 // invoke static member
263 // * Enables the blocked algorithm (BLAS level 3)
264 //
265 template< typename MatrixA, typename MatrixB, typename VectorALPHAR,
266 typename VectorALPHAI, typename VectorBETA, typename MatrixVSL,
267 typename MatrixVSR >
invokeboost::numeric::bindings::lapack::gges_impl268 static std::ptrdiff_t invoke( const char jobvsl, const char jobvsr,
269 const char sort, external_fp selctg, MatrixA& a, MatrixB& b,
270 fortran_int_t& sdim, VectorALPHAR& alphar,
271 VectorALPHAI& alphai, VectorBETA& beta, MatrixVSL& vsl,
272 MatrixVSR& vsr, optimal_workspace ) {
273 namespace bindings = ::boost::numeric::bindings;
274 real_type opt_size_work;
275 bindings::detail::array< fortran_bool_t > tmp_bwork( min_size_bwork(
276 bindings::size_column(a), sort ) );
277 detail::gges( jobvsl, jobvsr, sort, selctg,
278 bindings::size_column(a), bindings::begin_value(a),
279 bindings::stride_major(a), bindings::begin_value(b),
280 bindings::stride_major(b), sdim,
281 bindings::begin_value(alphar), bindings::begin_value(alphai),
282 bindings::begin_value(beta), bindings::begin_value(vsl),
283 bindings::stride_major(vsl), bindings::begin_value(vsr),
284 bindings::stride_major(vsr), &opt_size_work, -1,
285 bindings::begin_value(tmp_bwork) );
286 bindings::detail::array< real_type > tmp_work(
287 traits::detail::to_int( opt_size_work ) );
288 return invoke( jobvsl, jobvsr, sort, selctg, a, b, sdim, alphar,
289 alphai, beta, vsl, vsr, workspace( tmp_work, tmp_bwork ) );
290 }
291
292 //
293 // Static member function that returns the minimum size of
294 // workspace-array work.
295 //
min_size_workboost::numeric::bindings::lapack::gges_impl296 static std::ptrdiff_t min_size_work( const std::ptrdiff_t n ) {
297 if ( n == 0 ) {
298 return 1;
299 } else {
300 return std::max< std::ptrdiff_t >( 8*n, 6*n + 16 );
301 }
302 }
303
304 //
305 // Static member function that returns the minimum size of
306 // workspace-array bwork.
307 //
min_size_bworkboost::numeric::bindings::lapack::gges_impl308 static std::ptrdiff_t min_size_bwork( const std::ptrdiff_t n,
309 const char sort ) {
310 if ( sort == 'N' )
311 return 0;
312 else
313 return n;
314 }
315 };
316
317 //
318 // This implementation is enabled if Value is a complex type.
319 //
320 template< typename Value >
321 struct gges_impl< Value, typename boost::enable_if< is_complex< Value > >::type > {
322
323 typedef Value value_type;
324 typedef typename remove_imaginary< Value >::type real_type;
325
326 //
327 // Static member function for user-defined workspaces, that
328 // * Deduces the required arguments for dispatching to LAPACK, and
329 // * Asserts that most arguments make sense.
330 //
331 template< typename MatrixA, typename MatrixB, typename VectorALPHA,
332 typename VectorBETA, typename MatrixVSL, typename MatrixVSR,
333 typename WORK, typename RWORK, typename BWORK >
invokeboost::numeric::bindings::lapack::gges_impl334 static std::ptrdiff_t invoke( const char jobvsl, const char jobvsr,
335 const char sort, external_fp selctg, MatrixA& a, MatrixB& b,
336 fortran_int_t& sdim, VectorALPHA& alpha, VectorBETA& beta,
337 MatrixVSL& vsl, MatrixVSR& vsr, detail::workspace3< WORK, RWORK,
338 BWORK > work ) {
339 namespace bindings = ::boost::numeric::bindings;
340 BOOST_STATIC_ASSERT( (bindings::is_column_major< MatrixA >::value) );
341 BOOST_STATIC_ASSERT( (bindings::is_column_major< MatrixB >::value) );
342 BOOST_STATIC_ASSERT( (bindings::is_column_major< MatrixVSL >::value) );
343 BOOST_STATIC_ASSERT( (bindings::is_column_major< MatrixVSR >::value) );
344 BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
345 typename bindings::value_type< MatrixA >::type >::type,
346 typename remove_const< typename bindings::value_type<
347 MatrixB >::type >::type >::value) );
348 BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
349 typename bindings::value_type< MatrixA >::type >::type,
350 typename remove_const< typename bindings::value_type<
351 VectorALPHA >::type >::type >::value) );
352 BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
353 typename bindings::value_type< MatrixA >::type >::type,
354 typename remove_const< typename bindings::value_type<
355 VectorBETA >::type >::type >::value) );
356 BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
357 typename bindings::value_type< MatrixA >::type >::type,
358 typename remove_const< typename bindings::value_type<
359 MatrixVSL >::type >::type >::value) );
360 BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
361 typename bindings::value_type< MatrixA >::type >::type,
362 typename remove_const< typename bindings::value_type<
363 MatrixVSR >::type >::type >::value) );
364 BOOST_STATIC_ASSERT( (bindings::is_mutable< MatrixA >::value) );
365 BOOST_STATIC_ASSERT( (bindings::is_mutable< MatrixB >::value) );
366 BOOST_STATIC_ASSERT( (bindings::is_mutable< VectorALPHA >::value) );
367 BOOST_STATIC_ASSERT( (bindings::is_mutable< VectorBETA >::value) );
368 BOOST_STATIC_ASSERT( (bindings::is_mutable< MatrixVSL >::value) );
369 BOOST_STATIC_ASSERT( (bindings::is_mutable< MatrixVSR >::value) );
370 BOOST_ASSERT( bindings::size(alpha) >= bindings::size_column(a) );
371 BOOST_ASSERT( bindings::size(beta) >= bindings::size_column(a) );
372 BOOST_ASSERT( bindings::size(work.select(fortran_bool_t())) >=
373 min_size_bwork( bindings::size_column(a), sort ));
374 BOOST_ASSERT( bindings::size(work.select(real_type())) >=
375 min_size_rwork( bindings::size_column(a) ));
376 BOOST_ASSERT( bindings::size(work.select(value_type())) >=
377 min_size_work( bindings::size_column(a) ));
378 BOOST_ASSERT( bindings::size_column(a) >= 0 );
379 BOOST_ASSERT( bindings::size_minor(a) == 1 ||
380 bindings::stride_minor(a) == 1 );
381 BOOST_ASSERT( bindings::size_minor(b) == 1 ||
382 bindings::stride_minor(b) == 1 );
383 BOOST_ASSERT( bindings::size_minor(vsl) == 1 ||
384 bindings::stride_minor(vsl) == 1 );
385 BOOST_ASSERT( bindings::size_minor(vsr) == 1 ||
386 bindings::stride_minor(vsr) == 1 );
387 BOOST_ASSERT( bindings::stride_major(a) >= std::max< std::ptrdiff_t >(1,
388 bindings::size_column(a)) );
389 BOOST_ASSERT( bindings::stride_major(b) >= std::max< std::ptrdiff_t >(1,
390 bindings::size_column(a)) );
391 BOOST_ASSERT( jobvsl == 'N' || jobvsl == 'V' );
392 BOOST_ASSERT( jobvsr == 'N' || jobvsr == 'V' );
393 BOOST_ASSERT( sort == 'N' || sort == 'S' );
394 return detail::gges( jobvsl, jobvsr, sort, selctg,
395 bindings::size_column(a), bindings::begin_value(a),
396 bindings::stride_major(a), bindings::begin_value(b),
397 bindings::stride_major(b), sdim, bindings::begin_value(alpha),
398 bindings::begin_value(beta), bindings::begin_value(vsl),
399 bindings::stride_major(vsl), bindings::begin_value(vsr),
400 bindings::stride_major(vsr),
401 bindings::begin_value(work.select(value_type())),
402 bindings::size(work.select(value_type())),
403 bindings::begin_value(work.select(real_type())),
404 bindings::begin_value(work.select(fortran_bool_t())) );
405 }
406
407 //
408 // Static member function that
409 // * Figures out the minimal workspace requirements, and passes
410 // the results to the user-defined workspace overload of the
411 // invoke static member function
412 // * Enables the unblocked algorithm (BLAS level 2)
413 //
414 template< typename MatrixA, typename MatrixB, typename VectorALPHA,
415 typename VectorBETA, typename MatrixVSL, typename MatrixVSR >
invokeboost::numeric::bindings::lapack::gges_impl416 static std::ptrdiff_t invoke( const char jobvsl, const char jobvsr,
417 const char sort, external_fp selctg, MatrixA& a, MatrixB& b,
418 fortran_int_t& sdim, VectorALPHA& alpha, VectorBETA& beta,
419 MatrixVSL& vsl, MatrixVSR& vsr, minimal_workspace ) {
420 namespace bindings = ::boost::numeric::bindings;
421 bindings::detail::array< value_type > tmp_work( min_size_work(
422 bindings::size_column(a) ) );
423 bindings::detail::array< real_type > tmp_rwork( min_size_rwork(
424 bindings::size_column(a) ) );
425 bindings::detail::array< fortran_bool_t > tmp_bwork( min_size_bwork(
426 bindings::size_column(a), sort ) );
427 return invoke( jobvsl, jobvsr, sort, selctg, a, b, sdim, alpha, beta,
428 vsl, vsr, workspace( tmp_work, tmp_rwork, tmp_bwork ) );
429 }
430
431 //
432 // Static member function that
433 // * Figures out the optimal workspace requirements, and passes
434 // the results to the user-defined workspace overload of the
435 // invoke static member
436 // * Enables the blocked algorithm (BLAS level 3)
437 //
438 template< typename MatrixA, typename MatrixB, typename VectorALPHA,
439 typename VectorBETA, typename MatrixVSL, typename MatrixVSR >
invokeboost::numeric::bindings::lapack::gges_impl440 static std::ptrdiff_t invoke( const char jobvsl, const char jobvsr,
441 const char sort, external_fp selctg, MatrixA& a, MatrixB& b,
442 fortran_int_t& sdim, VectorALPHA& alpha, VectorBETA& beta,
443 MatrixVSL& vsl, MatrixVSR& vsr, optimal_workspace ) {
444 namespace bindings = ::boost::numeric::bindings;
445 value_type opt_size_work;
446 bindings::detail::array< real_type > tmp_rwork( min_size_rwork(
447 bindings::size_column(a) ) );
448 bindings::detail::array< fortran_bool_t > tmp_bwork( min_size_bwork(
449 bindings::size_column(a), sort ) );
450 detail::gges( jobvsl, jobvsr, sort, selctg,
451 bindings::size_column(a), bindings::begin_value(a),
452 bindings::stride_major(a), bindings::begin_value(b),
453 bindings::stride_major(b), sdim, bindings::begin_value(alpha),
454 bindings::begin_value(beta), bindings::begin_value(vsl),
455 bindings::stride_major(vsl), bindings::begin_value(vsr),
456 bindings::stride_major(vsr), &opt_size_work, -1,
457 bindings::begin_value(tmp_rwork),
458 bindings::begin_value(tmp_bwork) );
459 bindings::detail::array< value_type > tmp_work(
460 traits::detail::to_int( opt_size_work ) );
461 return invoke( jobvsl, jobvsr, sort, selctg, a, b, sdim, alpha, beta,
462 vsl, vsr, workspace( tmp_work, tmp_rwork, tmp_bwork ) );
463 }
464
465 //
466 // Static member function that returns the minimum size of
467 // workspace-array work.
468 //
min_size_workboost::numeric::bindings::lapack::gges_impl469 static std::ptrdiff_t min_size_work( const std::ptrdiff_t n ) {
470 return std::max< std::ptrdiff_t >( 1, 2*n );
471 }
472
473 //
474 // Static member function that returns the minimum size of
475 // workspace-array rwork.
476 //
min_size_rworkboost::numeric::bindings::lapack::gges_impl477 static std::ptrdiff_t min_size_rwork( const std::ptrdiff_t n ) {
478 return 8*n;
479 }
480
481 //
482 // Static member function that returns the minimum size of
483 // workspace-array bwork.
484 //
min_size_bworkboost::numeric::bindings::lapack::gges_impl485 static std::ptrdiff_t min_size_bwork( const std::ptrdiff_t n,
486 const char sort ) {
487 if ( sort == 'N' )
488 return 0;
489 else
490 return n;
491 }
492 };
493
494
495 //
496 // Functions for direct use. These functions are overloaded for temporaries,
497 // so that wrapped types can still be passed and used for write-access. In
498 // addition, if applicable, they are overloaded for user-defined workspaces.
499 // Calls to these functions are passed to the gges_impl classes. In the
500 // documentation, most overloads are collapsed to avoid a large number of
501 // prototypes which are very similar.
502 //
503
504 //
505 // Overloaded function for gges. Its overload differs for
506 // * User-defined workspace
507 //
508 template< typename MatrixA, typename MatrixB, typename VectorALPHAR,
509 typename VectorALPHAI, typename VectorBETA, typename MatrixVSL,
510 typename MatrixVSR, typename Workspace >
511 inline typename boost::enable_if< detail::is_workspace< Workspace >,
512 std::ptrdiff_t >::type
gges(const char jobvsl,const char jobvsr,const char sort,external_fp selctg,MatrixA & a,MatrixB & b,fortran_int_t & sdim,VectorALPHAR & alphar,VectorALPHAI & alphai,VectorBETA & beta,MatrixVSL & vsl,MatrixVSR & vsr,Workspace work)513 gges( const char jobvsl, const char jobvsr, const char sort,
514 external_fp selctg, MatrixA& a, MatrixB& b, fortran_int_t& sdim,
515 VectorALPHAR& alphar, VectorALPHAI& alphai, VectorBETA& beta,
516 MatrixVSL& vsl, MatrixVSR& vsr, Workspace work ) {
517 return gges_impl< typename bindings::value_type<
518 MatrixA >::type >::invoke( jobvsl, jobvsr, sort, selctg, a, b,
519 sdim, alphar, alphai, beta, vsl, vsr, work );
520 }
521
522 //
523 // Overloaded function for gges. Its overload differs for
524 // * Default workspace-type (optimal)
525 //
526 template< typename MatrixA, typename MatrixB, typename VectorALPHAR,
527 typename VectorALPHAI, typename VectorBETA, typename MatrixVSL,
528 typename MatrixVSR >
529 inline typename boost::disable_if< detail::is_workspace< MatrixVSR >,
530 std::ptrdiff_t >::type
gges(const char jobvsl,const char jobvsr,const char sort,external_fp selctg,MatrixA & a,MatrixB & b,fortran_int_t & sdim,VectorALPHAR & alphar,VectorALPHAI & alphai,VectorBETA & beta,MatrixVSL & vsl,MatrixVSR & vsr)531 gges( const char jobvsl, const char jobvsr, const char sort,
532 external_fp selctg, MatrixA& a, MatrixB& b, fortran_int_t& sdim,
533 VectorALPHAR& alphar, VectorALPHAI& alphai, VectorBETA& beta,
534 MatrixVSL& vsl, MatrixVSR& vsr ) {
535 return gges_impl< typename bindings::value_type<
536 MatrixA >::type >::invoke( jobvsl, jobvsr, sort, selctg, a, b,
537 sdim, alphar, alphai, beta, vsl, vsr, optimal_workspace() );
538 }
539
540 //
541 // Overloaded function for gges. Its overload differs for
542 // * User-defined workspace
543 //
544 template< typename MatrixA, typename MatrixB, typename VectorALPHA,
545 typename VectorBETA, typename MatrixVSL, typename MatrixVSR,
546 typename Workspace >
547 inline typename boost::enable_if< detail::is_workspace< Workspace >,
548 std::ptrdiff_t >::type
gges(const char jobvsl,const char jobvsr,const char sort,external_fp selctg,MatrixA & a,MatrixB & b,fortran_int_t & sdim,VectorALPHA & alpha,VectorBETA & beta,MatrixVSL & vsl,MatrixVSR & vsr,Workspace work)549 gges( const char jobvsl, const char jobvsr, const char sort,
550 external_fp selctg, MatrixA& a, MatrixB& b, fortran_int_t& sdim,
551 VectorALPHA& alpha, VectorBETA& beta, MatrixVSL& vsl, MatrixVSR& vsr,
552 Workspace work ) {
553 return gges_impl< typename bindings::value_type<
554 MatrixA >::type >::invoke( jobvsl, jobvsr, sort, selctg, a, b,
555 sdim, alpha, beta, vsl, vsr, work );
556 }
557
558 //
559 // Overloaded function for gges. Its overload differs for
560 // * Default workspace-type (optimal)
561 //
562 template< typename MatrixA, typename MatrixB, typename VectorALPHA,
563 typename VectorBETA, typename MatrixVSL, typename MatrixVSR >
564 inline typename boost::disable_if< detail::is_workspace< MatrixVSR >,
565 std::ptrdiff_t >::type
gges(const char jobvsl,const char jobvsr,const char sort,external_fp selctg,MatrixA & a,MatrixB & b,fortran_int_t & sdim,VectorALPHA & alpha,VectorBETA & beta,MatrixVSL & vsl,MatrixVSR & vsr)566 gges( const char jobvsl, const char jobvsr, const char sort,
567 external_fp selctg, MatrixA& a, MatrixB& b, fortran_int_t& sdim,
568 VectorALPHA& alpha, VectorBETA& beta, MatrixVSL& vsl,
569 MatrixVSR& vsr ) {
570 return gges_impl< typename bindings::value_type<
571 MatrixA >::type >::invoke( jobvsl, jobvsr, sort, selctg, a, b,
572 sdim, alpha, beta, vsl, vsr, optimal_workspace() );
573 }
574
575 } // namespace lapack
576 } // namespace bindings
577 } // namespace numeric
578 } // namespace boost
579
580 #endif
581