В настоящее время я пытаюсь получить код fortran FE (конечный элемент) для работы с openmp. У меня есть петля над всеми элементами, ie, что я хочу работать параллельно. Вот упрощенная часть кода, который не работаетopenmp fortran сокращение и критический не работает для массива
!$omp parallel do default(none) shared(nelm,A,res,enod) private(ie,Fe,B,edof)
do ie=1,nelm
call calcB(B,A(:,ie))
call calcFe(Fe,B)
write(*,*) Fe !writes Fe=40d0, this is correct
call getEdof(edof,enod(:,ie))
!$OMP CRITICAL
res(edof)=res(edof)+fe
!$OMP END CRITICAL
enddo
!$omp end parallel do
Цель кода для вычисления силы Fe, а затем добавить его к res в edof. Сила вычисляется с помощью calcFe, и рассчитанная сила верна, но результат res некорректен после цикла.
Если я заменяю calcFe с просто Fe=40d0 затем добавить его в res результат является правильным после цикла
!$omp parallel do default(none) shared(nelm,A,res,enod) private(ie,Fe,B,edof)
do ie=1,nelm
call calcB(B,A(:,ie))
Fe=40d0
call getEdof(edof,enod(:,ie))
!$OMP CRITICAL
res(edof)=res(edof)+fe
!$OMP END CRITICAL
enddo
!$omp end parallel do
Что вызывает эту ошибку? В обоих случаях Fe=40d0 объявлен private, но только один из них дает правильный результат. Вместо использования !$ CRITICAL я мог бы использовать reduction, но он дает ту же ошибку. В программе также используются несколько больших и разреженных матриц, но они пассивные/не используются во время цикла. У моего диспетчера были проблемы с использованием openmp и разреженных матриц раньше и подозревается, что они используют одну и ту же память. Если ошибка не очевидна, какой отладчик лучше всего использовать? Я новичок и для fortran, и для openmp и для программирования.
Im, используя ifort для компиляции, и моя ОС - ubuntu.
EDIT: Добавлен упрощенный код, который вы можете запустить, хотя этот код работает В коде есть две петли, на параллельный и один последовательный порт, чтобы они должны дать тот же результат, res и res2
program main
use omp_lib
implicit none
integer :: ie, nelm,enod(4,50*50),edof(12),i,j,k
double precision ::B(12,8),fe(12),A(12,12,2500),res(2601*3),res2(2601*3),finish,start
!creates enod
i=1
do j=1,50
ie=j
do k=1,50
nelm=k
enod(:,i)=(/ 51*(nelm-1)+1+ie-1, 51*(nelm-1)+1+ie, 51*(nelm)+1+ie-1, 51*(nelm)+1+ie /)
i=i+1
end do
end do
A=1d0
res2=0d0
nelm=2500
start=omp_get_wtime()
!$omp parallel do default(none) shared(nelm,A,enod) private(ie,fe,edof,B) reduction(+:res2)
do ie=1,nelm
call calcB(B,A(:,:,ie))
call calcFe(fe,B) !the calculated fe is always 2304
!can write fe=2304 to get correct result with real code
call getEdof(edof,enod(:,ie))
res2(edof)=res2(edof)+fe
end do
!$omp end parallel do
finish=omp_get_wtime()
write(*,*) 'time: ', finish-start
res=0d0
nelm=2500
start=omp_get_wtime()
do ie=1,nelm
call calcB(B,A(:,:,ie))
call calcFe(fe,B)
call getEdof(edof,enod(:,ie))
res(edof)=res(edof)+fe
end do
finish=omp_get_wtime()
write(*,*) 'time: ', finish-start
write(*,*) 'difference: ',sum(res2-res)
write(*,*) sum(res2)
stop
end program main
subroutine calcB(B,A)
double precision ::B(12,8),A(12,12),C(12)
integer ::gp
C=1d0
do gp=1,8
B(:,gp)=matmul(A,C)
end do
end subroutine calcB
subroutine calcFe(fe,B)
double precision ::fe(12),B(12,8),D(12,12)
integer ::gp
fe=0d0
D=2d0
do gp=1,8
fe=fe+matmul(D,B(:,gp))
end do
end subroutine calcFe
subroutine getEdof(edof,enod)
implicit none
integer,intent(in) :: enod(4)
integer,intent(out):: edof(12)
edof=0
edof(1:3) =(/ enod(1)*3-2, enod(1)*3-1, enod(1)*3 /)
edof(4:6) =(/ enod(2)*3-2, enod(2)*3-1, enod(2)*3 /)
edof(7:9) =(/ enod(3)*3-2, enod(3)*3-1, enod(3)*3 /)
edof(10:12)=(/ enod(4)*3-2, enod(4)*3-1, enod(4)*3 /)
end subroutine getedof
И сделайте файл
FF = ifort -O3 -openmp
OBJ1 = main.f90
ls: $(FORT_OBJS)
$(FF) -o exec $(OBJ1)
к сожалению, этот кусок кода работает, так что я не в состоянии воспроизвести ошибку. res2 и res рассчитываются последовательно и параллельно. В моей реальной программе я поставил все значения в 1d0, чтобы получить постоянную fe. Калорированный fe является правильным, если я добавлю write(*,*) fe после calcFe, я вижу, что значения верны. Затем я добавляю эти значения в res2 и сравниваю их с серийным номером res. Затем они отличаются большим запасом, поэтому нет числовой ошибки округления. Если я просто объявляю fe=2304 в своей основной программе, я получаю правильный ответ, хотя fe уже есть 2304, когда используется запись.
В моей реальной программе все подпрограммы находятся в разных модулях, нужно ли мне проявлять особую осторожность из-за этого? Также в одном из модулей используются некоторые глобальные переменные, они только для чтения, но поскольку они не объявлены в подпрограмме, они не становятся автоматически закрытыми? Это не должно быть проблемой, так как я положил все переменные, используемые для для calulate fe постоянная, глобальные переменные не используется непосредственно для расчета fe
Используйте сокращение здесь, критическое - плохая идея в этом случае. –
Попробуйте создать полностью компилируемый автономный пример, который мы можем скомпилировать и запустить. Объясните, какое значение является неправильным результатом и какой результат вы ожидаете. –
Я боюсь, что мы ничего не можем сделать без кусочка кода, который * воспроизводит * ошибку. Часто вы сами находите проблему, когда пытаетесь изолировать небольшой рабочий пример. * ", так как они не объявлены в подпрограмме, они не становятся автоматически закрытыми?" * ... нет, но если они не изменены, они могут оставаться разделяемыми, что не является проблемой. –