Question about lwork in dgesdd

Open discussion regarding features, bugs, issues, vendors, etc.

Question about lwork in dgesdd

Postby bencteux » Sat Dec 06, 2008 5:22 am

Hello everybody,
I want to use dgesdd to compute a SVD decomposition. (JOBZ = 'A').
As suggested in documentation, I make a memory query to know the optimal size of "work".
The numbers I get for lwork does not agree with the formula given in the documentation.
For example : with a (322,288) matrix, the code gives optimal lwork = 250848, although I expected from the formula : lwork >= 3 * 288^2 + 4 * 288^2 + 4 * 288 = 581760 .

Is there an error in the documentation, or in my code ?
Thanks for help
Guy Bencteux
bencteux
 
Posts: 13
Joined: Wed Apr 16, 2008 11:52 am

Re: Question about lwork in dgesdd

Postby Julien Langou » Fri Dec 12, 2008 4:13 pm

Bonjour Guy,

I confirm the numbers you are giving (250848 and 581760).

There is consequently a problem on our side. It's not a big deal but still.

The first dimension of LWORK, the one given by the header of the LAPACK routine:
Code: Select all
            LWORK >= 3*min(M,N)*min(M,N) + max(max(M,N),4*min(M,N)*min(M,N)+4*min(M,N)).

is supposed to be the smallest required workspace for the routine to complete the job. Call it LWORK_OPT_MEM. The dimension of LWORK obtained from the WORKSPACE QUERY mechanism (LWORK=-1, etc.) is supposed to be larger and to enable blocking of the operations and so high performance. Call it LWORK_OPT_PERF.

In the case you gave us, LWORK_OPT_PERF is smaller than LWORK_OPT_MEM and that makes no sense whatsoever. I agree. That ought to be fixed. I will add this to the ERRATA in LAPACK. (Webpage to come soon ... )

Best wishes,
Julien.
Julien Langou
 
Posts: 727
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Re: Question about lwork in dgesdd

Postby Julie » Tue Mar 17, 2009 3:23 pm

Guy,
I just corrected the problem in our repository.
I fixed some formulas in the comments that were wrong
Code: Select all
 
*          If JOBZ = 'N',
 *            LWORK >= 3*min(M,N) + max(max(M,N),7*min(M,N)).
 *          If JOBZ = 'O',
-*            LWORK >= 3*min(M,N)*min(M,N) +
+*            LWORK >= 3*min(M,N) +
 *                     max(max(M,N),5*min(M,N)*min(M,N)+4*min(M,N)).
 *          If JOBZ = 'S' or 'A'
-*            LWORK >= 3*min(M,N)*min(M,N) +
+*            LWORK >= 3*min(M,N) +
 *                     max(max(M,N),4*min(M,N)*min(M,N)+4*min(M,N)).
 *          For good performance, LWORK should generally be larger.
 *          If LWORK = -1 but other input arguments are legal, WORK(1)

Merci de ta contribution
Julie
Julie
 
Posts: 299
Joined: Wed Feb 23, 2005 12:32 am
Location: ICL, Denver. Colorado


Return to User Discussion

Who is online

Users browsing this forum: Google [Bot] and 1 guest