[Pw_forum] ATLAS DSYGVX problem + mailing list question

Riccardo Di Meo dimeo at democritos.it
Wed May 26 19:12:43 CEST 2010


Riccardo Di Meo wrote:
> Paolo Giannozzi wrote:
>   
>> Riccardo Di Meo wrote:
>>
>>   
>>     
>>> It seems therefore that the culprit is the line:
>>>
>>>      nb = ILAENV( 1, 'DSYTRD', 'U', n, -1, -1, -1 )
>>>
>>> which produces two different results with the same input, although I 
>>> don't know why.
>>>     
>>>       
>> ILAENV is a Lapack auxiliary routine. The Atlas version of ILAENV
>> should return the value of nb that is most suited to the Atlas
>> version of Lapack. Apparently it doesn't even return a correct
>> value. Not a QE problem, apparently.
>>
>> P.
>>   
>>     
> Thank you for the explanation: it seems indeed a lapack/atlas problem.
>
> I'll look a bit into it, and post if I'll find something of some use.
>
> Best regards,
> Riccardo
>   
Dear all,

I digged a bit in the LAPACK documentation/source

http://www.netlib.org/lapack/util/ilaenv.f
http://www.netlib.org/lapack/double/dsygvx.f

and wheter this is an Espresso problem or a Lapack one, depends on how 
the said docs are read.

For the way *I* understand it, the ILAENV routine returns the optimal nb 
value for DSYGVX, but it's responsability of the application (QE) to 
check if said number satisfies the constraints set by DSYGVX itself.

This also makes sense, since as far as I understand, the results of 
ILAENV with a specific 'NAME' argument, fits multiple LAPACK routines, 
and therefore couldn't be possibly checked in advance by ILAENV.

This interpretation is also reinforced by the trailing comment to 
DSYGVX, in which it's required:

LWORK >= max(1,8*N)

Not much point in writing it explicitly, if the value returned by ILAENV 
(which should be the optimal one) would produce the correct lwork anyway...

That said, inserting:

IF (lwork < 8 * n) THEN
   lwork = 8 * n
ENDIF

before the ALLOCATE(work(lwork)) line fixes the problem (and doesn't 
hurt anyway).

So the bug is neither in QE or lapack, but in the lapack documentation IMO.

Best Regards,
RDM





More information about the users mailing list