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-08 04:20:39

On Fri, Apr 8, 2011 at 9:10 AM, Rutger ter Borg <rutger_at_[hidden]> wrote:
> On 2011-04-07 14:04, sguazt wrote:
>>
>> So maybe we can return to the code I've posted earlier ;)
>>
>
> Hey Marco,
>
> I've integrated your code and I've committed the results to the sandbox
> repository.
>

Thanks Rutger!

Maybe, just to cover all possible cases, we can just return 1 when at
least one the two dimensions is <= 0, that is:

--- [code] ---
if ( std::min(m,n) <= 0 ) {
  return 1;
}
if ( m >= n ) {
    if ( jobu == 'N' ) {
        return 5*n;
    } else {
        return std::max< $INTEGER_TYPE >(3*n+m,5*n);
    }
} else {
    if ( jobvt == 'N' ) {
        return 5*m;
    } else {
        return std::max< $INTEGER_TYPE >(3*m+n,5*m);
    }
}
--- [/code] ---

Now we should fix the complex case also ;)
After a look at the related subroutines:
  http://www.netlib.org/lapack/complex/cgesvd.f
  http://www.netlib.org/lapack/complex16/zgesvd.f

I have obtained this set of conditions:

--- [code] ---
m = ...; // # rows of A
n = ...; // # columns of A
sz = 0; // will hold the min work size
/* k = ::std::min(m,n); */
if (m >= n /* && k > 0 */) {
 if (jobu == 'N') {
   sz = 3*n;
 } else {
   sz = 2*n+m;
 }
} else /* if (k > 0) */ {
 if (jobvt == 'N') {
   sz = 3*m;
 } else {
   sz = 2*m+n;
 }
}
return ::std::max(sz,1);
--- [code] ---

which should translate in something like this:

--- [code] ---
if ( std::min(m,n) <= 0 ) {
  return 1;
}
if ( m >= n ) {
    if ( jobu == 'N' ) {
        return 3*n;
    } else {
        return 2*n+m;
    }
} else {
    if ( jobvt == 'N' ) {
        return 3*m;
    } else {
        return 2*m+n;
    }
}
--- [/code] ---

Cheers,

-- Marco