CODEWITHSUNDEEP
fortranlapack
Sumit
Sumit

Reputation: 221

How to find optimal block size and LWORK in LAPACK

I am trying to find inverse and eigenfunctions of nxn Hermitian matrices using Fortran with lapack.

How do I choose the optimal values for parameters like lda, lwork, liwork and lrwork. I browse through some example and find these choices

integer,parameter::lda=nh
integer,parameter::lwork=2*nh+nh*nh
integer,parameter::liwork=3+5*nh
integer,parameter::lrwork=1 + 5*nh + 2*nh*nh

where nh is the dimension of the matrix. I also find another example with lwork=16*nh. How can I determine the best choice? At this point, I am dealing with 500x500 Hermitian matrices (maximum).

I found this documentation which suggests

WORK

(workspace) REAL array, dimension (LWORK)

On exit, if INFO = 0, then WORK(1) returns the optimal LWORK.

LWORK

(input) INTEGER

The dimension of the array WORK. LWORK  max(1,N).

For optimal performance LWORK  N*NB, where NB is the optimal block size returned by ILAENV.

Is it possible to find out the optimal block size using WORK or ILAENV for a given matrix dimension?

I am using both gfortran and ifort with mkl.

EDIT

Based on the comment by @percusse and @kvantour's answer here is a sample code

character,parameter::jobz="v",uplo="u"
integer, parameter::nh=15
complex*16::m(nh,nh),m1(nh,nh)

integer,parameter::lda=nh
integer::ipiv(nh),info

complex*16::work(1)
real*8::rwork(1), w(nh)
integer::iwork(1)
real*8::x1(nh,nh),x2(nh,nh)

call random_seed()
call random_number(x1)
call random_number(x2)

m=cmplx(x1,x2)
m1=conjg(m)
m1=transpose(m1)
m=(m+m1)/2.0

call zheevd(jobz,uplo,nh,m,lda,w,work,-1,rwork,-1,iwork, -1,info)

print*,"info : ", info
print*,"lwork: ", int(work(1))   , 2*nh+nh*nh
print*,"lrwork:", int(rwork(1))  , 1 + 5*nh + 2*nh*nh
print*,"liwork:", int(iwork(1))  , 3+5*nh

end

info : 0

lwork: 255 255

lrwork: 526 526

liwork: 78 78

Upvotes: 3

Views: 819

Answers (1)

kvantour
kvantour

Reputation: 26511

I'm not sure what you are implying with "Is it possible to find out the optimal block size using WORK or ILAENV for a particular machine architecture?". You can however find the optimal values for a particular problem.

Eg. If you want to find the eigenvalues of a complex Hermitian matrix, using cheev, you can ask the routine to return you the value :

subroutine CHEEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, INFO )
  character                , intent(in)    :: JOBZ
  character                , intent(in)    :: UPLO
  integer                  , intent(in)    ::  N
  complex, dimension(lda,*), intent(inout) :: A
  integer                  , intent(in)    :: LDA
  real   , dimension(*)    , intent(out)   :: W
  complex, dimension(*)    , intent(out)   :: WORK
  integer                  , intent(in)    :: LWORK
  real   , dimension(*)    , intent(out)   :: RWORK
  integer                  , intent(out)   :: INFO

Then the documentation clearly states (be advised, in the past this was easier to read):

WORK is COMPLEX array, dimension (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

LWORK is INTEGER The length of the array WORK. LWORK >= max(1,2*N-1). For optimal efficiency, LWORK >= (NB+1)*N, where NB is the blocksize for CHETRD returned by ILAENV. If LWORK = -1, then a workspace query is assumed; the routine only calculates the optimal size of the WORK array, returns this value as the first entry of the WORK array, and no error message related to LWORK is issued by XERBLA.

So all you need to do is 

call cheev(jobz, uplo, n, a, lda, w, work, -1, rwork, info)
lwork=int(work(1))
dallocate(work)
allocate(work(lwork))
call cheev(jobz, uplo, n, a, lda, w, work, lwork, rwork, info)

Upvotes: 3

Related Questions