Re: [ublas] [bindings] LAPACK ?GESVD and minimum workspace size

Subject: Re: [ublas] [bindings] LAPACK ?GESVD and minimum workspace size
From: sguazt (marco.guazzone_at_[hidden])
Date: 2011-04-07 04:29:48

Hi Rutger,

On Thu, Apr 7, 2011 at 9:26 AM, Rutger ter Borg <rutger_at_[hidden]> wrote:
> On 04/06/2011 05:08 PM, sguazt wrote:
>>
[cut]
>
> You are saying that, in some cases, for gesvd, the minimum workspace size
> returned may be larger than the queried optimum workspace size, and the
> bindings are asserting on that error?
>

Exactly. Specifically, what the problem is in the assertion which
involves the call to the min_size_work function.

> Looking at the Fortran code, I see that there are not less than 10 code
> paths for determining workspace size, but that the variation in minimum
> seems to be limited to
>
> max( 4*n, 5*n )
> max( 3*(m+n), 5*n )
> max( 4*m, 5*n )
>
> Could you come up with a clever branch tree? Let's incorporate that in the
> bindings. This is only for the real variants of gesvd, right?

I've just updated my branch. It is at revision 71072.

For what concerns the problem, it affects all variants of GESVD.
For instance, with complex<double> type, I get the following error:

--- [error] ---
Testing zgesvd...
lapack_gesvd: /home/marco/projects/boost-numeric_bindings/boost/numeric/bindings/lapack/driver/gesvd.hpp:295:
static ptrdiff_t boost::numeric::bindings::lapack::gesvd_impl<Value,
typename boost::enable_if<boost::numeric::bindings::is_complex<T>
>::type>::invoke(char, char, MatrixA&, VectorS&, MatrixU&, MatrixVT&,
boost::numeric::bindings::lapack::detail::workspace2<WORK, BWORK>)
[with MatrixA = boost::numeric::ublas::matrix<std::complex<double>,
boost::numeric::ublas::basic_column_major<>,
boost::numeric::ublas::unbounded_array<std::complex<double>,
std::allocator<std::complex<double> > > >, VectorS =
boost::numeric::ublas::vector<double,
boost::numeric::ublas::unbounded_array<double, std::allocator<double>
> >, MatrixU = boost::numeric::ublas::matrix<std::complex<double>,
boost::numeric::ublas::basic_column_major<>,
boost::numeric::ublas::unbounded_array<std::complex<double>,
std::allocator<std::complex<double> > > >, MatrixVT =
boost::numeric::ublas::matrix<std::complex<double>,
boost::numeric::ublas::basic_column_major<>,
boost::numeric::ublas::unbounded_array<std::complex<double>,
std::allocator<std::complex<double> > > >, WORK =
boost::numeric::bindings::detail::array<std::complex<double> >, RWORK
= boost::numeric::bindings::detail::array<double>, Value =
std::complex<double>, ptrdiff_t = long int]: Assertion
`bindings::size(work.select(value_type())) >= min_size_work(
bindings::size_row(a), bindings::size_column(a), minmn )' failed.
Aborted (core dumped)
--- [/error] ---

For your convenience I attach a simple test.
I wasn't able to catch the assertion at runtime, so to try different
variants of GESVD you need to comment the others one (e.g., to try the
CGESVD you need to comment lines from 50 to 54, limits included).

I've compiled it with G++ 4.5.1:
g++ -Wall -ansi -pedantic -o lapack_gesvd lapack_gesvd.cpp -lm -llapack

Cheers,

-- Marco