
Ublas : 
Subject: Re: [ublas] [bindings] Problem with assertion in GGES
From: Marco Guazzone (marco.guazzone_at_[hidden])
Date: 20100715 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/explorehtml/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