Re: [ublas] [bindings] Problem with assertion in GGES

Subject: Re: [ublas] [bindings] Problem with assertion in GGES
From: Marco Guazzone (marco.guazzone_at_[hidden])
Date: 2010-07-15 03:48:53

On Wed, Jul 14, 2010 at 9:48 PM, Rutger ter Borg <rutger_at_[hidden]> wrote:
> Marco Guazzone wrote:
>
>>
>> Looking inside the gges.hpp code, the failing assertion
>>
>>         BOOST_ASSERT( bindings::size(work.select(real_type())) >=
>> min_size_work( bindings::size_column(a) ));
>>
>> fails because for this example:
>>
>>     bindings::size(work.select(real_type())) --> 28
>>     min_size_work(size_column(a)) --> 32
>>
>> I've tried to comment this assertion and all works and the computations
>> are OK.
>>
>> So, where am I wrong?
>>
>> Thank you for your patience!!
>>
>> Cheers,
>>
>
> I am guessing the documentation in GGES concerning the minimal dimensions of
> the real work array are not correct. I.e., the optimal workspace query (=let
> lapack answer the workspace requirements) is returning something else than
> what is described as minimal requirements.
>
> We have seen this kind of error before (for gelsd, IIRC). It requires
> reading the Fortran code for gges, to figure out the true minimal workspace
> requirements. IOW, is it about the formula in min_size_work (max(1,8n+16))
> being incorrect, but I verified that it does match the Fortran description.
>
> Fix is rather trivial once the correct formula has been found in the Fortran
> code.
>

Hi Rutger.

Just an idea (but note that I'm completely unaware of FORTRAN syntax)...

Here

http://www.netlib.org/lapack/explore-html/a00643_source.html

I see this (line 278)

   MINWRK = MAX( 8*N, 6*N + 16 )

In my case (N=2) this will set MINWRK to 28 (=2*2+16) which is good
for the test done in the assertion.

Maybe we could change "min_size_work" in this way:

    static std::ptrdiff_t min_size_work( const std::ptrdiff_t n ) {
        //return std::max< std::ptrdiff_t >( 1, 8*n + 16 );
        return std::max< std::ptrdiff_t >( 8*n, 6*n + 16 );//<-- NEW CHECK
    }

NOTE: just an observation: the old size computation "max(1,8*n+16)"
could have been replaced by return 8*n+16 since (provided n>0) this
quantity is always greater than 1.

Cheers,

-- Marco