Avoiding Workspace Query

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

Avoiding Workspace Query

Postby jmlopez » Wed Apr 06, 2011 2:25 pm

There is something that caught my eye in one of the posts:

Re: Inverting the following positive semidefinite matrix

Postby Julien Langou » Tue Apr 05, 2011 11:22 pm
Hello

1) DPOTRF + DPOTRI will not work on that matrix. DPOTRF will fail and return with INFO = 1, this is because your matrix is symmetric but it is not symmetric positive definite. (Since A(1,1) = 0, e_1^T A e_1 = 0, so the matrix is not definite. )

2) DGETRF + DGETRI should work. I am glad you got that working.

3) DSYTRF + DSYTRI should work as well, it works just fine on my machine, see quickly written C code below.

4) we should use the C interface to LAPACK when working in C. See:
http://www.netlib.org/lapack/#_standard ... for_lapack
This would avoid the use for writing the prototypes (there is .h file to include ready for use)
This would avoid the use for the worskpace query and workspace allocation.
All in all that would be nicer. I try to send a following post later. Feel free to do it yourself and post the code for others to look at, if you have time on your end.

Julien


Point number 4: How does the C interface to LAPACK avoids the use for the workspace query and workspace allocation?

I looked at the source code for LAPACKE_dgeev since I'm interested in computing the eigenvalues of a matrix. In any case, yes, I can see that the function does not require me to input the workspace. All of this is done in C.

Is this C interface actually optimal? What I mean is the following, LAPACKE_dgeev calls the LAPACKE_dgeev_work twice (I'm assuming LAPACKE_dgeev_work calls the actual fortran routine).

Code: Select all
    /* Query optimal working array(s) size */
    info = LAPACKE_dgeev_work( matrix_order, jobvl, jobvr, n, a, lda, wr, wi,
                               vl, ldvl, vr, ldvr, &work_query, lwork );
    if( info != 0 ) {
        goto exit_level_0;
    }
    lwork = (lapack_int)work_query;
    /* Allocate memory for work arrays */
    work = (double*)LAPACKE_malloc( sizeof(double) * lwork );
    if( work == NULL ) {
        info = LAPACK_WORK_MEMORY_ERROR;
        goto exit_level_0;
    }
    /* Call middle-level interface */
    info = LAPACKE_dgeev_work( matrix_order, jobvl, jobvr, n, a, lda, wr, wi,
                               vl, ldvl, vr, ldvr, work, lwork );


If we look at part of the code in dgeev.f we can see that if LQUERY is true then we exit the routine, otherwise we continue to compute the eigenvalues and eigenvectors. What bothers me is the fact that the fortran routine always computes the optimal workspace size regardless if lquery is true or false.

This may be due to the way fortran works but I think that a better way would be to avoid the computation of the workspace if lwork is not -1. Recall that

Code: Select all
 LQUERY = ( LWORK.EQ.-1 )


So in C we would first look for the optimal size of the workspace by making the call with lwork = -1 and then once we have the optimal size call the routine again with fortran avoiding the executing of the computation of the workspace size.

I'm considering doing this adjustment myself but I would like to know if someone here knows why the programmers decided to compute the workspace size each time the routine was called.

-Manuel
jmlopez
 
Posts: 3
Joined: Wed Apr 06, 2011 1:43 pm
Location: University of Houston, Houston, TX

Re: Avoiding Workspace Query

Postby admin » Wed Apr 06, 2011 11:33 pm

Manuel, you made a good point.
At the moment, yes we are computing the workspace each time in LAPACK. It is the way it has been done for 20 years.
We are working of having a smarter workspace calculation for the next major release of LAPACK but we are still at the discussion stage.
LAPACK contains lots of routines and that change is not as straight-forward as we can imagine.
Nonetheless, your input is greatly appreciated.

The goal of the new C interface is convenience for users without modifying the original LAPACK library. The cost of the LWORK=-1 call is and the Workspace computation is really really minimal. So yes we can say that the C interface is as optimal as LAPACK. And remember that the new C interface has the nowork routine that does not allocate the space for you if you want to save one call.
I encourage you to make some comparison on decent size matrices. (more than 10) and let us know if there is a performance problem.

Thanks for your input!
admin
Site Admin
 
Posts: 486
Joined: Wed Dec 08, 2004 7:07 pm

Re: Avoiding Workspace Query

Postby Julien Langou » Wed Apr 06, 2011 11:50 pm

If we look at part of the code in dgeev.f we can see that if LQUERY is true then we exit the routine, otherwise we continue to compute the eigenvalues and eigenvectors. What bothers me is the fact that the fortran routine always computes the optimal workspace size regardless if lquery is true or false.


You are right. If LQUERY is false, the LAPACK subroutine computes MAXWRK, (the optimal workspace,) which is used nowhere after.

In the case LQUERY is false, we need to compute MINWRK for the following lines of codes:
Code: Select all
         IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
            INFO = -13
         END IF

But we do not need MAXWRK.

So yes the code could be written differently. (More optimally.) If LQUERY is false, we could only compute MINWRK which is either 4*N or 3*N so is fairly easy to compute. No need to go and compute the optimal workspace MAXWRK which is much more complicated. Granted. This would not be too hard to do. Maybe not top priority at this point though.

In the case of a trusted user, like the C-interface, we could have a mechanism which says: "LWORK is good, do not check it, proceed". For this subroutine, since we are speaking of checking against 3*N or 4*N, I am not sure such a mechanism is worth it.

J.
Julien Langou
 
Posts: 733
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Re: Avoiding Workspace Query

Postby jmlopez » Thu Apr 07, 2011 12:55 pm

At the moment, yes we are computing the workspace each time in LAPACK. It is the way it has been done for 20 years.


I find that hard to believe given that the Fortran code has so many if statements yet it does unnecessary work when LWORK is not -1.

The cost of the LWORK=-1 call is and the Workspace computation is really really minimal. So yes we can say that the C interface is as optimal as LAPACK


This might be the case but the thought of my machine doing useless things really gets me thinking why...would anyone do this?

We are working of having a smarter workspace calculation for the next major release of LAPACK but we are still at the discussion stage.
LAPACK contains lots of routines and that change is not as straight-forward as we can imagine.


I would agree that the calculation of the workspace is not straight-forward, I have no idea what the routine does to calculate this but again, the point is
that some (I have only checked the routines to compute eigenvalues and to invert a matrix) routines do work that is not required (even if it is really really minimal). This could avoided by a simple check at the beginning of the routine, either you want to calculate the workspace in which case you can set most
of the output arguments to a NULL pointer or you want to proceed with the computation.

In any case, for now I will keep using LAPACK as it has been done for the last 20 years and hope that in future releases the code is written more optimally. Thank you.

-Manuel
jmlopez
 
Posts: 3
Joined: Wed Apr 06, 2011 1:43 pm
Location: University of Houston, Houston, TX


Return to User Discussion

Who is online

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