如何在 LAPACK 中找到最佳块大小和 LWORK

Question

我正在尝试使用 Fortran 和 lapack 查找 nxn Hermitian 矩阵的逆函数和特征函数。

如何为lda、lwork、liwork 和lrwork 等参数选择最佳值。我浏览了一些示例并找到了这些选择

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

其中nh 是矩阵的维度。我还找到了lwork=16*nh 的另一个例子。如何确定最佳选择？此时，我正在处理 500x500 Hermitian 矩阵（最大）。

我找到了this documentation，这表明

工作

（工作区）REAL 数组，维度（LWORK）

退出时，如果 INFO = 0，则 WORK(1) 返回最优 LWORK。

工作

（输入）整数

数组WORK的维度。 LWORK  最大（1，N）。

为了获得最佳性能 LWORK  N*NB，其中 NB 是 ILAENV 返回的最佳块大小。

对于给定的矩阵维度，是否可以使用WORK 或ILAENV 找出最佳块大小？

我在 mkl 中同时使用 gfortran 和 ifort。

---

编辑

根据@percusse 和@kvantour's answer 的评论，这里有一个示例代码

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

信息：0

工作：255 255

lrwork: 526 526

liwork: 78 78

Answer 1

我不确定你在暗示什么“是否可以使用WORK 或ILAENV 找出特定机器架构的最佳块大小？”。但是，您可以找到特定问题的最佳值。

例如。如果你想找到一个复杂的 Hermitian 矩阵的特征值，使用cheev，你可以让例程返回值：

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 

然后documentation 明确指出（请注意，过去这更容易阅读）：

WORK 是 COMPLEX 数组，维度 (MAX(1,LWORK)) 退出时，如果INFO = 0，WORK(1) 返回最优LWORK。

LWORK 是 INTEGER 数组WORK 的长度。 LWORK >= max(1,2*N-1)。 为了获得最佳效率，LWORK >= (NB+1)*N， 其中NB 是由ILAENV 返回的CHETRD 的块大小。如果LWORK = -1，则假定为工作区查询；例行公事 仅计算 WORK 数组的最佳大小，返回 这个值作为WORK数组的第一个条目，并且没有错误 与LWORK 相关的消息由XERBLA 发出。

所以你需要做的就是

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)