Я пытался написать программу, которая вычисляет миллионы тензоров диполь-дипольного взаимодействия, а также их производные. Поскольку эти тензоры тривиально распараллеливаемы и часто вырождаются, я решил построить справочную таблицу (LUT) и распределить работу. В конечном итоге они будут объединены в большую матрицу и диагонализированы (в конечном итоге я буду использовать для этого scalapack. Пока что диаг подходит для ноды nersc). Чтобы отслеживать все индексы в фортране, я создал пару производных типов данных.
type dtlrdr_lut
sequence
integer p
integer q
integer s
integer i
integer ind
real(dp), dimension(3,3) :: dtlrdr
end type dtlrdr_lut
type dtlrdh_lut
sequence
integer p
integer q
integer ind
real(dp), dimension(3, 3) :: TLR
real(dp), dimension(3, 3, 3, 3) :: dTLRdh
end type dtlrdh_lut
В моей подпрограмме, где я хочу распараллелить все это, у меня есть:
type(dtlrdr_lut), dimension(:), allocatable :: my_dtlrdr, collected_dtlrdr
type(dtlrdh_lut), dimension(:), allocatable :: my_dtlrdh, collected_dtlrdh
integer :: dh_dtype, dr_dtype, dh_types(5), dr_types(6), dh_blocks(5), dr_blocks(6)
INTEGER(KIND=MPI_ADDRESS_KIND) :: dh_offsets(5), dr_offsets(6)
if(.not.allocated(my_dtlrdh)) allocate(my_dtlrdh(my_num_pairs))
if(.not.allocated(my_dtlrdr)) allocate(my_dtlrdr(my_num_pairs*3*nat))
if(me_image.eq.root_image) then
if(.not.allocated(collected_dtlrdh)) allocate(collected_dtlrdh(num_pairs))
if(.not.allocated(collected_dtlrdr)) allocate(collected_dtlrdr(num_pairs*3*nat))
end if
call mpi_get_address(my_dtlrdr(1)%p, dr_offsets(1), ierr)
call mpi_get_address(my_dtlrdr(1)%q, dr_offsets(2), ierr)
call mpi_get_address(my_dtlrdr(1)%s, dr_offsets(3), ierr)
call mpi_get_address(my_dtlrdr(1)%i, dr_offsets(4), ierr)
call mpi_get_address(my_dtlrdr(1)%ind, dr_offsets(5), ierr)
call mpi_get_address(my_dtlrdr(1)%dtlrdr, dr_offsets(6), ierr)
do i = 2, size(dr_offsets)
dr_offsets(i) = dr_offsets(i) - dr_offsets(1)
end do
dr_offsets(1) = 0
dr_types = (/MPI_INTEGER, MPI_INTEGER, MPI_INTEGER, MPI_INTEGER, MPI_INTEGER, MPI_DOUBLE_PRECISION/)
dr_blocks = (/1, 1, 1, 1, 1, 3*3/)
call mpi_type_struct(6, dr_blocks, dr_offsets, dr_types, dr_dtype, ierr)
call mpi_type_commit(dr_dtype, ierr)
call mpi_get_address(my_dtlrdh(1)%p, dh_offsets(1), ierr)
call mpi_get_address(my_dtlrdh(1)%q, dh_offsets(2), ierr)
call mpi_get_address(my_dtlrdh(1)%ind, dh_offsets(3), ierr)
call mpi_get_address(my_dtlrdh(1)%TLR, dh_offsets(4), ierr)
call mpi_get_address(my_dtlrdh(1)%dTLRdh, dh_offsets(5), ierr)
do i = 2, size(dh_offsets)
dh_offsets(i) = dh_offsets(i) - dh_offsets(1)
end do
dh_offsets(1) = 0
dh_types = (/MPI_INTEGER, MPI_INTEGER, MPI_INTEGER, MPI_DOUBLE_PRECISION, MPI_DOUBLE_PRECISION/)
dh_blocks = (/1, 1, 1, 3*3, 3*3*3*3/)
call mpi_type_struct(5, dh_blocks, dh_offsets, dh_types, dh_dtype, ierr)
call mpi_type_commit(dh_dtype, ierr)
call mpi_gather(my_dtlrdh, my_num_pairs, dh_dtype, &
collected_dtlrdh, num_pairs, dh_dtype, &
root_image, intra_image_comm)
call mp_barrier(intra_image_comm)
call mpi_gather(my_dtlrdr, my_num_pairs*3*nat, dr_dtype, &
collected_dtlrdr, num_pairs*3*nat, dr_dtype, &
root_image, intra_image_comm)
Каков результат кода? Итак, корневой процесс собрался и стал барьерным, а затем сбой seg:
Program received signal SIGSEGV: Segmentation fault - invalid memory reference.
Backtrace for this error:
#0 0x10bac04f2
#1 0x10bac0cae
#2 0x7fff8d7c1f19
В этом моделировании процесса 0:
size(my_dtlrdh) = 97
size(collected_dtlrdh) = 194
size(my_dtlrdr) = 873
size(collected_dtlrdr) = 1746
и в процессе 1
size(my_dtlrdh) = 97
size(collected_dtlrdh) = 3
size(my_dtlrdr) = 873
size(collected_dtlrdr) = 1650521
Когда я печатаю смещения, блоки и т. д. для процесса 0, я получаю:
printing dr vars 0
dr_blocks = 1 1 1 1 1 9
dr_offsets = 0 4 8 12 16 24
dr_types = 7 7 7 7 7 17
dr_dtype = 73
printing dh vars 0
dr_blocks = 1 1 1 9 81
dr_offsets = 0 4 8 16 88
dr_types = 7 7 7 17 17
dr_dtype = 74
и для процесса 1 я получаю:
printing dr vars 1
dr_blocks = 1 1 1 1 1 9
dr_offsets = 0 4 8 12 16 24
dr_types = 7 7 7 7 7 17
dr_dtype = 73
printing dh vars 1
dr_blocks = 1 1 1 9 81
dr_offsets = 0 4 8 16 88
dr_types = 7 7 7 17 17
dr_dtype = 74
Однако случайный размер dtlrdr в proc1 не должен иметь значения, потому что на самом деле он ничего не получает. Кажется, я не могу понять, что происходит или почему процесс 1 не может пройти сбор без недопустимой ссылки на память. Любые идеи? И, пожалуйста, дайте мне знать, если вам, ребята, нужна дополнительная информация от меня.
