how to use LWORK = -1 in DGEQRF ?

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

how to use LWORK = -1 in DGEQRF ?

Postby Mat » Thu Aug 24, 2006 5:47 pm

Hello,
i have problems to use the LWORK argument within the dgeqrf_ funtction.
When i set the value as recommended within the fortran files:
LWORK >= max(1,N)
i get another output than if i set the value to -1.
Additionally i don't know how to adapt the WORK array if WORK has the value -1. How big has to be the WORK array?

Thanks a lot for help
Mat
 
Posts: 47
Joined: Sat Aug 19, 2006 9:54 am

Postby Julie » Thu Aug 24, 2006 8:32 pm

Hello Mat,

Let start with generality.

Most of the LAPACK routines require the user to give them a workspace (for example DGEQRF). In general this workspace is called WORK and it is followed by the integer variable LWORK in the calling sequence. This is for example the case for DGEQRF.

LWORK represents the size of the workspace you are providing to LAPACK. So WORK is an array of DOUBLE PRECISION number of size LWORK. To be fully precise, WORK is an array of DOUBLE PRECISION number of size MAX(1,LWORK). WORK need to be of size at least 1.

To go fast, LAPACK subroutines work on blocks of columns, this enables us to use Level 3 BLAS subroutines (for example matrix-matrix products). If we work one column at a time, we can only exploit Level 2 BLAS subroutines (for example matrix-vector products) and consequently we have poor performance.
Now there is a trade-off: when we work on block of columns, our workspace gets bigger than the one we would have needed if we worked on one column at a time.
There is a trade-off performance-memory.

There is two basic configuration of users:
(1) those who are short in memory and they want to provide LAPACK with the smallest workspace required to solve the problem, let's call this value LWORK_MIN
(2) those who have memory and want to go as fast as possible, in this case they want the algorithm to use a larger workspace, let's call this value LWORK_FAST

When you do workspace query to an LAPACK routine (i.e. when you call an LAPACK routine with LWORK=-1), then LAPACK returns the value of LWORK_FAST in WORK(1).

If you want to know LWORK_MIN, you need to read this information in the header of the routine (as you have done).

LWORK_MIN and LWORK_FAST are in general different.

Some remarks (that should hopefully makes more sense now):

(1) If you call a routine with LWORK < LWORK_MIN, this is not valid. The routine tells you that the parameter LWORK is wrong.

(2) LAPACK is allowed to use all the workspace from WORK(1) to WORK(LWORK). So make sure that WORK is always at least of size LWORK.

(3) LWORK_MIN only depends on the algorithm so this can be written in hard in a header file. LWORK_FAST depends mostly on the block size your LAPACK implementation is using. This block size is (or should be) tuned for optimal performance on your machine so we expect it to be machine dependent. Obtaining the blocksize used by the algorithm and then computing the associated workspace is complex and machine depedent so this is why the workspace query to determine LWORK_FAST exists.

(4) When you do a workspace query make sure WORK is an array of size at least 1! Because the routine will try to return WORK_FAST in WORK(1)

(5) all this is valid for most LAPACK routines but not all

(6) The workspace is small compared to a dense n-by-n matrix. Most of the users can afford to allocate LWORK_FAST.

Well I hope this is clearer and reexplain what is written the header of DGEQRF.

When i set the value as recommended within the fortran files:
LWORK >= max(1,N)
i get another output than if i set the value to -1.


This is expected. MAX(1,N) is LWORK_MIN. A workspace query returns LWORK_FAST.

Additionally i don't know how to adapt the WORK array if WORK has the value -1. How big has to be the WORK array?


Do you mean: "Additionally i don't know how to adapt the WORK array if LWORK has the value -1. How big has to be the WORK array?"

First:How big has to be the WORK array?
At least LWORK_MIN.

But otherwise, yes, you can be stuck to LWORK_MIN if you do not have memory allocation capacity from your language (for example F77).

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: No registered users and 4 guests