program assign6_2
use printmat2
implicit none
real,allocatable :: a(:,:)
open(1,file='ped.a6-2')
call ainv
print*,'Additive relationship matrix'
call printnice(a,'(20f12.6)')
contains
subroutine ainv
integer :: i,k,l,io,p(3),nped=0
real ::w(3),res(4),mendel
w=(/1., -.5, -.5/)
res=(/2., 4/3., 1., 0./)
do
read(1,*,iostat=io) p
if(io /= 0) exit
nped=nped+1
enddo
print*,'# records in pedigree =',nped
allocate(a(nped,nped))
a=0
rewind 1
do
read(1,*,iostat=io) p
if (io /= 0) exit
i=1
if(p(2) == 0) i=i+1
if(p(3) == 0) i=i+1
mendel=res(i)
do k=1,3
do l=1,3
if (p(k) /= 0 .and. p(l) /=0) then
a(p(k),p(l))=a(p(k),p(l))+w(k)*w(l)*mendel
endif
enddo
enddo
enddo
end subroutine
end program assign6_2