[fortran] 分子血縁行列作成プログラム部品
変量効果の推定とBLUP法 単行本 2007/2
佐々木 義之 (著)
ISBN-13: 978-4876987023
のp124にて掲載された内容を実装したFortran90プログラム部品です。自由にご利用ください。
!tm -------------------------------------------------------
subroutine create_a00
! numerator relatonship matrix subroutine by Fortran90
! 出典:
! 変量効果の推定とBLUP法 単行本 2007/2
! 佐々木 義之 (著)
! ISBN-13: 978-4876987023
! p124
!----------------------------
! a s d sex -> a s d renumbered!
! 1 0 0 m 1 0 0
! 2 1 0 m 2 1 0
! x 1 0 f 3 1 0
! y 1 0 f 4 1 0
! 3 1 0 m 5 2 0
! 4 1 0 m 6 1 0
! 5 3 0 m 7 5 0
! 6 0 x m 8 0 3
! 7 3 y m 9 5 4
!----------------------------
integer, parameter :: nanim = 9
integer :: i, j
real (8) :: amat(0:nanim, 0:nanim)
real (8) :: atmp
integer :: sid(nanim), did(nanim)
print *, '--- create_A00_sasaki ---'
!tm -------------------------
sid(1) = 0
sid(2) = 1
sid(3) = 1
sid(4) = 1
sid(5) = 2
sid(6) = 1
sid(7) = 5
sid(8) = 0
sid(9) = 5
!tm -------------------------
did(1) = 0
did(2) = 0
did(3) = 0
did(4) = 0
did(5) = 0
did(6) = 0
did(7) = 0
did(8) = 3
did(9) = 4
!tm -------------------------
do i = 1, nanim
if (sid(i)>0 .and. did(i)>0) then
amat(i, i) = 1.0d0 + 0.5d0*amat(sid(i), did(i))
else
amat(i, i) = 1.0d0
end if
do j = 1, i - 1
atmp = 0.0d0
if (sid(i)>0) atmp = atmp + 0.5d0*amat(j, sid(i))
if (did(i)>0) atmp = atmp + 0.5d0*amat(j, did(i))
amat(j, i) = atmp
amat(i, j) = atmp
end do
end do
write (*, *) '2016 blup sasaki book p124'
write (*, *) '* a s d sex -> a s d revnum'
write (*, *) '* 1 0 0 m 1 0 0'
write (*, *) '* 2 1 0 m 2 1 0'
write (*, *) '* x 1 0 f 3 1 0'
write (*, *) '* y 1 0 f 4 1 0'
write (*, *) '* 3 1 0 m 5 2 0'
write (*, *) '* 4 1 0 m 6 1 0'
write (*, *) '* 5 3 0 m 7 5 0'
write (*, *) '* 6 0 x m 8 0 3'
write (*, *) '* 7 3 y m 9 5 4'
write (*, *) ''
do i = 1, nanim
write (*, '(20f7.4)')(amat(i,j), j=1, nanim)
end do
! 1.0000 0.5000 0.5000 0.5000 0.2500 0.5000 0.1250 0.2500 0.3750
! 0.5000 1.0000 0.2500 0.2500 0.5000 0.2500 0.2500 0.1250 0.3750
! 0.5000 0.2500 1.0000 0.2500 0.1250 0.2500 0.0625 0.5000 0.1875
! 0.5000 0.2500 0.2500 1.0000 0.1250 0.2500 0.0625 0.1250 0.5625
! 0.2500 0.5000 0.1250 0.1250 1.0000 0.1250 0.5000 0.0625 0.5625
! 0.5000 0.2500 0.2500 0.2500 0.1250 1.0000 0.0625 0.1250 0.1875
! 0.1250 0.2500 0.0625 0.0625 0.5000 0.0625 1.0000 0.0313 0.2813
! 0.2500 0.1250 0.5000 0.1250 0.0625 0.1250 0.0313 1.0000 0.0938
! 0.3750 0.3750 0.1875 0.5625 0.5625 0.1875 0.2813 0.0938 1.0625
end subroutine create_a00
佐々木 義之 (著)
ISBN-13: 978-4876987023
のp124にて掲載された内容を実装したFortran90プログラム部品です。自由にご利用ください。
!tm -------------------------------------------------------
subroutine create_a00
! numerator relatonship matrix subroutine by Fortran90
! 出典:
! 変量効果の推定とBLUP法 単行本 2007/2
! 佐々木 義之 (著)
! ISBN-13: 978-4876987023
! p124
!----------------------------
! a s d sex -> a s d renumbered!
! 1 0 0 m 1 0 0
! 2 1 0 m 2 1 0
! x 1 0 f 3 1 0
! y 1 0 f 4 1 0
! 3 1 0 m 5 2 0
! 4 1 0 m 6 1 0
! 5 3 0 m 7 5 0
! 6 0 x m 8 0 3
! 7 3 y m 9 5 4
!----------------------------
integer, parameter :: nanim = 9
integer :: i, j
real (8) :: amat(0:nanim, 0:nanim)
real (8) :: atmp
integer :: sid(nanim), did(nanim)
print *, '--- create_A00_sasaki ---'
!tm -------------------------
sid(1) = 0
sid(2) = 1
sid(3) = 1
sid(4) = 1
sid(5) = 2
sid(6) = 1
sid(7) = 5
sid(8) = 0
sid(9) = 5
!tm -------------------------
did(1) = 0
did(2) = 0
did(3) = 0
did(4) = 0
did(5) = 0
did(6) = 0
did(7) = 0
did(8) = 3
did(9) = 4
!tm -------------------------
do i = 1, nanim
if (sid(i)>0 .and. did(i)>0) then
amat(i, i) = 1.0d0 + 0.5d0*amat(sid(i), did(i))
else
amat(i, i) = 1.0d0
end if
do j = 1, i - 1
atmp = 0.0d0
if (sid(i)>0) atmp = atmp + 0.5d0*amat(j, sid(i))
if (did(i)>0) atmp = atmp + 0.5d0*amat(j, did(i))
amat(j, i) = atmp
amat(i, j) = atmp
end do
end do
write (*, *) '2016 blup sasaki book p124'
write (*, *) '* a s d sex -> a s d revnum'
write (*, *) '* 1 0 0 m 1 0 0'
write (*, *) '* 2 1 0 m 2 1 0'
write (*, *) '* x 1 0 f 3 1 0'
write (*, *) '* y 1 0 f 4 1 0'
write (*, *) '* 3 1 0 m 5 2 0'
write (*, *) '* 4 1 0 m 6 1 0'
write (*, *) '* 5 3 0 m 7 5 0'
write (*, *) '* 6 0 x m 8 0 3'
write (*, *) '* 7 3 y m 9 5 4'
write (*, *) ''
do i = 1, nanim
write (*, '(20f7.4)')(amat(i,j), j=1, nanim)
end do
! 1.0000 0.5000 0.5000 0.5000 0.2500 0.5000 0.1250 0.2500 0.3750
! 0.5000 1.0000 0.2500 0.2500 0.5000 0.2500 0.2500 0.1250 0.3750
! 0.5000 0.2500 1.0000 0.2500 0.1250 0.2500 0.0625 0.5000 0.1875
! 0.5000 0.2500 0.2500 1.0000 0.1250 0.2500 0.0625 0.1250 0.5625
! 0.2500 0.5000 0.1250 0.1250 1.0000 0.1250 0.5000 0.0625 0.5625
! 0.5000 0.2500 0.2500 0.2500 0.1250 1.0000 0.0625 0.1250 0.1875
! 0.1250 0.2500 0.0625 0.0625 0.5000 0.0625 1.0000 0.0313 0.2813
! 0.2500 0.1250 0.5000 0.1250 0.0625 0.1250 0.0313 1.0000 0.0938
! 0.3750 0.3750 0.1875 0.5625 0.5625 0.1875 0.2813 0.0938 1.0625
end subroutine create_a00