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
