Re: Parallel programming with FORALL?

On Jan 11, 9:43 pm, Gordon Sande <Gordon.Sa...@xxxxxxxxx> wrote:
FORALL is an elaborate arraay assignment to cover the cases that simpler
syntax misses. Think about the diagonal of a 2d matrix. Computes the
full righthandside before doing any assignment.

Interesting. I knew it wasn't a looping construct, but you've added
some details I didn't know.


Could be parallel but
not guaranteed. Just like at lot of other array assigments might be.
Oh the joys of "as if"!

Yeah, that was going to be my follow-up question: Whether I can get
array operations to run in parallel.



If you want to use multiple cores with gfrtran you will need to look at OpenMP
and all that flows from that.

:-(

I just took a quick look at OpenMP. So, basically I have to explicitly
create and destroy threads if i want to use multiple cores. It sounds
like it would be harder to implement any algorithm like that.

I just tested hello-world with OpenMP:

program helloMP
use omp_lib
integer:: id, nthreads

id = omp_get_thread_num()
write (*,*) 'Hello World from thread', id

if ( id == 0 ) then
nthreads = omp_get_num_threads()
write (*,*) 'There are', nthreads, 'threads'
end if
end program


On my computer this produces two threads. I wonder how it decides how
many threads to make. Is it always one thread per CPU core?


Anyways, I tried a stab at writing the program with threads, but I
failed miserably (segmentation fault). My code is below, but I
probably need to start by finding an OpenMP tutorial...

program helloMP

use omp_lib
implicit none

integer, parameter :: N = 10000
integer :: i, j, a, b, id, nthreads
real :: vec(N), mat(N,N)

!
! Initialize "A" before breaking into threads.
!
call seed_prng()
call random_number(vec)


!
! Start threads.
!
id = omp_get_thread_num()


! Pick a slice of B based on my thread ID.
if (id == 0) then
a = 1
b = N/2
else
a = N/2+1
b = N
end if

! Operate on that slice.
forall (i = a:b, j = 1:N)
mat(i,j) = gamma(vec(i) + i/j)
mat(i,j) = cosh(mat(j,i)) + tanh(mat(i,j) + j/i)
end forall



if ( id == 0 ) then
nthreads = omp_get_num_threads()
write (*,*) 'There are', nthreads, 'threads'

! Select a random value to print.
i = floor(1 + vec(1)*N)
j = floor(1 + vec(2)*N)
print *,"i = ",i
print *,"j = ",j
print *,"mat(i,j) = ",mat(i,j)
end if

contains

!!
!! Use system time to seed PRNG.
!!
subroutine seed_prng()
implicit none

! Internal variables
integer :: i = 0, n, clock
integer, dimension(:), allocatable :: seed

! Seed the PRNG
call random_seed(size = n)
allocate(seed(n))

call system_clock(count=clock)
seed = clock + 37 * (/ (i-1, i = 1,n) /)

call random_seed(put = seed)
deallocate(seed)
end subroutine
end program

.