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 05:01:36

On Thu, Apr 7, 2011 at 10:41 AM, Rutger ter Borg <rutger_at_[hidden]> wrote:
> On 04/07/2011 10:29 AM, sguazt wrote:
>>
>> 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.
>>
>
> Sorry, I meant if you could suggest the optimal conditions/structure of
> something like
>
> if ( ... ) {
>   return max( 3*(m+n) );
> } else if ( ... ) {
>   return ...;
> } else {
>   return ...;
> }

Hmmm is not so trivial.
Maybe we can take the problem from another perspective, something like:

invoke(const char jobu, const char jobvt, MatrixA& a, VectorS& s,
MatrixU& u, MatrixVT& vt, detail::workspace1< WORK > work)
{
  If (invoke is called after computing the optimal workspace size)  {
    do not call min_size_work
  } else {
    do call min_size_work
  }
}

The condition in the "if" statement is for the case when the invoke
function is called by the other overloaded invoke function (e.g., see
line 204) which computes the optimal workspace size.

Instead, the "else" case cover the case the user directly passes to
gesvd the workspace size.
In principle min_size_work could still fail for this case. A possible
workaround is to let min_size_work to query gesvd for the optimal size
instead of doing itself the computation.

What do you think?

-- Marco