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