Dear Fernando,

Thanks for your reply. I have changed the code as below but I still get a crash. I guess still need your help. Thanks

Kon

The compilation:

/opt/intel/bin/ifort collaeu.f90 -L/opt/intel/composerxe-2011.2.137/mkl/lib/intel64 -I/opt/intel/composerxe-2011.2.137/mkl/include/intel64/ilp64 -lmkl_blas95_ilp64 -lmkl_lapack95_ilp64 -Wl,--start-group /opt/intel/composerxe-2011.2.137/mkl/lib/intel64/libmkl_intel_ilp64.a /opt/intel/composerxe-2011.2.137/mkl/lib/intel64/libmkl_intel_thread.a /opt/intel/composerxe-2011.2.137/mkl/lib/intel64/libmkl_core.a -Wl,--end-group -openmp -lpthread -o collaeu

The code:

program Colleau_algorithm

!$ use OMP_LIB

implicit none

INTEGER :: i, j, k, l, m, n, io, irow, icol, numAnim, ierr, jerr, numOfThreads

REAL(8), DIMENSION(:,:), ALLOCATABLE :: A, Ainv, T, Tinv, Asub

REAL(8), DIMENSION(:), ALLOCATABLE :: D, Dinv, Sol, Rhs, r, v, s, w, y, x

INTEGER, DIMENSION(:,:), ALLOCATABLE :: Ped

INTEGER, DIMENSION(:), ALLOCATABLE :: indx, indx1

real(8) :: a1

real :: start_time, end_time

numOfThreads = 8

!CALL OPM_SET_NUM_THREADS(numOfThreads)

call cpu_time(start_time)

OPEN(1,FILE="/group/adhis/Kon/800Kdata/holstein/old_50K/A.hol")

OPEN(2,FILE="/group/adhis/Kon/800Kdata/holstein/old_50K/Asubset.hol")

n = 0; io = 0

do

read(1,*,iostat=io) i, j, k, a1, l

if(io < 0) exit

if (i == 0.and.j==0.and.k==0) cycle

n = n + 1

enddo

numAnim = n

print*, n, " animals in the pedigree "

ALLOCATE (D(n), Sol(n), Rhs(n), stat=ierr )

if (ierr.ne.0) then

print*," Allocation request for D, Sol, Rhs denied! "

stop

endif

ALLOCATE (Ped(n,3), v(n), indx(n), indx1(n), stat=jerr )

if (jerr.ne.0) then

print*," Allocation request for Ped, v, indx denied! "

stop

endif

rewind 1

n = 0; io = 0

do

read(1,*,iostat=io) i, j, k, a1, l

if(io < 0) exit

if (i == 0.and.j==0.and.k==0) cycle

n = n + 1

Ped(n,1) = i; Ped(n,2) = j; Ped(n,3) = k

D(n) = a1; indx(n) = l

enddo

print*, " second reading : ", n, " animals in the pedigree "

PRINT*, " Commencing Colleau calculations for sub-matrix : "

m = 0

do i = 1, n

if(indx(i) > 0) then

m = m + 1

indx1(m) = i

endif

enddo

print*, " Number of animals in the subset matrix of A is : " , m

ALLOCATE (Asub(m,m), stat=jerr )

if (jerr.ne.0) then

print*," Allocation request for Asub denied! "

stop

endif

!$ call OMP_SET_DYNAMIC(.false.)

!$ if (OMP_GET_DYNAMIC()) print *,'Warning: dynamic adjustment of threads has been set'

!$ call OMP_SET_NUM_THREADS(5)

!$OMP PARALLEL DO &

!$OMP PRIVATE(k,i,j,n,m, Rhs, Sol, v, Ped, D, Asub, indx1)

do k = 1, m

Rhs = 0

Rhs(indx1(k)) = 1.0

Sol = 0

!

do i = n, 1, -1

Sol(i) = Sol(i) + Rhs(i)

if (Ped(i,2) /= 0) then

Sol(Ped(i,2)) = Sol(Ped(i,2)) + Sol(i)*0.5

end if

if (Ped(i,3) /= 0) then

Sol(Ped(i,3)) = Sol(Ped(i,3)) + Sol(i)*0.5

end if

enddo

!

v = 0

do j = 1, n

if (Ped(j,2) == 0 .and. Ped(j,3) == 0) then

v(j) = D(j)*Sol(j)

end if

if (Ped(j,2) /= 0 .and. Ped(j,3) /= 0) then

v(j) = D(j)*Sol(j) + (v(Ped(j,2)) + v(Ped(j,3)))*0.5

end if

if (Ped(j,2) /= 0 .and. Ped(j,3) == 0) then

v(j) = D(j)*Sol(j) + (v(Ped(j,2)))*0.5

end if

if (Ped(j,2) == 0 .and. Ped(j,3) /= 0) then

v(j) = D(j)*Sol(j) + (v(Ped(j,3)))*0.5

end if

end do

do j = 1, m

WRITE(2,FMT='(2i10,f15.5)') k, j, v(indx1(j)

Asub(k, j) = v(indx1(j))

enddo

end do

!$OMP END PARALLEL DO

!WRITE(*,FMT='(" Final subset matrix : ")')

do i = 1, m

do j = 1, m

WRITE(2,FMT='(2i10,f20.10)') i, j, Asub(i,j)

end do

end do

call cpu_time(end_time)

print*, " Elapsed time taken to calculate a subset of A is : ", end_time - start_time

end program