Как найти оптимальный размер блока и LWORK в LAPACK

Я пытаюсь найти обратные и собственные функции эрмитовых матриц nxn, используя Фортран с лапаком.

Как выбрать оптимальные значения для таких параметров, как 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 (максимум).

Я нашел эту документацию, в которой предлагается

РАБОТА

(рабочая область) REAL массив, размер (LWORK)

При выходе, если INFO = 0, WORK (1) возвращает оптимальный LWORK.

LWORK

(ввод) ЦЕЛОЕ

Размерность массива РАБОТА. LWORK  max (1, N).

Для оптимальной производительности LWORK  N * NB, где NB - оптимальный размер блока, возвращаемый ILAENV.

Можно ли определить оптимальный размер блока, используя WORK или ILAENV для заданного размера матрицы?

Я использую и gfortran, и ifort с mkl.

ИЗМЕНИТЬ

На основе комментария @percusse и ответа @kvantour вот пример кода

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

работа: 78 78


fortran lapack
person Sumit    schedule 21.03.2018    source источник
comment
Сначала вы вызываете функцию с LWORK=-1 и получаете оптимальный размер блока. Затем вы снова вызываете функцию с этими возвращенными значениями   -  person percusse    schedule 21.03.2018
comment
Существует очень случаев, когда действительно нужно вычислить инверсию матрицы 500x500.   -  person knia    schedule 21.03.2018

Ответы (1)


arrow_upward
3
arrow_downward

Я не уверен, что вы имеете в виду, говоря «Можно ли определить оптимальный размер блока, используя WORK или ILAENV для конкретной архитектуры машины?». Однако вы можете найти оптимальные значения для конкретной проблемы.

Например. Если вы хотите найти собственные значения сложной эрмитовой матрицы, используя 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

Затем в документации четко говорится (имейте в виду, в прошлом это было проще читать):

WORK - это COMPLEX массив, измерение (MAX(1,LWORK)) При выходе, если INFO = 0, WORK(1) возвращает оптимальное LWORK.

LWORK равно INTEGER Длина массива WORK. LWORK >= max(1,2*N-1). Для оптимальной эффективности LWORK >= (NB+1)*N, где NB - размер блока для CHETRD, возвращаемый ILAENV. Если LWORK = -1, то предполагается запрос рабочей области; подпрограмма только вычисляет оптимальный размер массива WORK, возвращает это значение как первую запись массива WORK, и XERBLA не выдает сообщений об ошибках, связанных с LWORK.

Итак, все, что вам нужно сделать, это

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)
person kvantour    schedule 21.03.2018
comment
Спасибо, @kvantour. У меня сложилось впечатление, что эти цифры могут меняться в зависимости от типа процессора и оборудования. Итак, я могу запустить тест с lwork=-1 this для случайной эрмитовой матрицы nxn, получить все остальные параметры (например, lwork, rwork и т. Д.), А затем использовать их для любой эрмитовой матрицы nxn. Это так? - person Sumit; 21.03.2018
comment
@Sumit да, вывод не будет зависеть от содержимого A. Но об аргументах, определяющих проблему (UPLO, JOBZ, LDA, N в случае cheev) - person kvantour; 21.03.2018
comment
Я использую zheevd. Я полагаю, что правило такое же, и мне нужно найти iwork из liwork(1) дополнительно. - person Sumit; 21.03.2018
comment
спасибо @kvantour. Я получаю правильное значение для lwork и lrwork, но получаю неправильный результат для liwork. Я изменил свой вопрос образцом кода. - person Sumit; 22.03.2018
comment
@Sumit ваши аргументы неверны. work является сложным, rwork реальным, iwork целым ... И т. Д. Взгляните на ссылку, которую я отправил. В нем четко указаны типы аргументов. - person kvantour; 22.03.2018
comment
Теперь я вижу настоящую проблему. Мне нужно поспать. - person Sumit; 22.03.2018