same model, different level, multiple trait animal model 3 - 3형질
능력자료
0 0 1 6 1 7 9 0 9570 9670
1 5 1 6 2 7 10 9219 10530 10730
2 6 4 7 3 8 11 10025 11490 11790
2 6 3 7 4 8 12 7230 7600 8100
1 7 0 0 0 0 13 8121 0 0
위 자료를 dairy_3tr.dat
혈통자료
9 14 17
10 15 17
11 16 18
12 14 17
13 15 18
14 0 0
15 0 0
16 0 0
17 0 0
18 0 0
위 자료를 pedi_3tr.dat
파라미터 자료
30000 1150 1000 31000 1150 32000
8000 4000 3500 8300 4000 8500
1 0 0 0 0 0
0 0 0 1 0 0
1 1 0 1 0 0
0 0 0 0 0 1
1 0 1 0 0 1
0 0 0 1 1 1
1 1 1 1 1 1
위 자료를 vcvtrt_3tr.par로 저장
LHS와RHS를 구한는 프로그램
PROGRAM mt_dl_animal_model3_setup
! program name - multiple trait, different fixed level, animal model setup
! same model but different fixed level
! programmer - Park Byoungho
! usage - mt_df_animal_model3_setup
! purpose : read data file and write left-hand side and right-hand side
! Date - 2009.7.7.
! update - 208.7.14. change pedigree lhs making code
USE gi
IMPLICIT NONE
! data dictionary
INTEGER, PARAMETER :: no_of_trait = 3 ! number of trait
INTEGER, PARAMETER :: no_of_fixed = 2 ! number of fixed effects
INTEGER, PARAMETER :: no_of_eq = 18 ! number of equation
CHARACTER(LEN = 256) :: data_filename ! data file name
CHARACTER(LEN = 256) :: pedi_filename ! pedigree file name
CHARACTER(LEN = 256) :: par_filename ! parameter file name
INTEGER :: no_of_ud ! number of upper digonal
INTEGER :: no_of_tc ! number of trait combination
REAL(KIND = 8), ALLOCATABLE :: gvcv(:) ! upper diagonal part of genetic variance-covariance, (no_of_trait)*(no_of_trait + 1) / 2 dimension
REAL(KIND = 8), ALLOCATABLE :: rvcv(:) ! upper diagonal part of residual variance-covariance, (no_of_trait)*(no_of_trait + 1) / 2 dimension
REAL(KIND = 8), ALLOCATABLE :: temp_rvcv(:) ! temporary rvcv
REAL(KIND = 8), ALLOCATABLE :: rvcv_tc(:,:) ! rvcv according to trait combination and then inverse
REAL(KIND = 8), ALLOCATABLE :: rvcv_tc_m(:,:,:) ! rvcv matrix according to trait combination and then inverse
INTEGER, ALLOCATABLE :: trait_combi(:) ! trait combination
INTEGER :: tc_no ! trait combination number
INTEGER, ALLOCATABLE :: temp_effects(:) ! temporary animal and fixed effect
REAL(KIND = 8) :: milk1, milk2, r1, r2 ! data of trait for each line
INTEGER, ALLOCATABLE :: effects(:,:) ! array for fixed and animal effects
REAL(KINd = 8), ALLOCATABLE :: observations(:) ! array for observations
REAL(KINd = 8), ALLOCATABLE :: observations_mbr(:) ! array for observations multiplied by rvcv
INTEGER :: animal, sire, dam ! animal, sire, dam name for pedigree file
REAL(KIND = 8), ALLOCATABLE :: xpy(:,:) ! right-hand side
INTEGER, ALLOCATABLE :: nobs(:,:) ! number of observations
INTEGER :: i, j, k, l ! do loop
REAL :: temp_value ! temporary storate for swap
INTEGER, DIMENSION(99) :: iw ! for inverse
REAL(KIND = 8) :: z ! for inverse
INTEGER :: mr ! for inverse
INTEGER :: status ! I/O status
CHARACTER(LEN = 40) :: error_msg ! error message
INTEGER :: zero = 0
! get data file name
data_filename = 'dairy_3tr.dat'
! get pedigree file name
pedi_filename = 'pedi_3tr.dat'
! get parameter file name
par_filename = 'vcvtrt_3tr.par'
no_of_ud = (no_of_trait)*(no_of_trait + 1) / 2 ! number of upper diagonal of variance-covariance
no_of_tc = 0 ! number of trait combination
DO i = 1, no_of_trait
no_of_tc = no_of_tc + 2 ** (i - 1)
END DO
ALLOCATE(gvcv(no_of_ud)) ! genetic variance-covariance
ALLOCATE(rvcv(no_of_ud)) ! residual variance-covariance
ALLOCATE(temp_rvcv(no_of_ud)) ! temporary residual variance-covariance
ALLOCATE(rvcv_tc(no_of_tc, no_of_ud)) ! residual variance-covariance according to trait combination
ALLOCATE(rvcv_tc_m(no_of_tc, no_of_trait, no_of_trait)) ! residual variance-covariance matrix according to trait combination
ALLOCATE(trait_combi(no_of_ud)) ! trait combination
ALLOCATE(xpy(no_of_eq, no_of_trait)) ! right hand side
ALLOCATE(nobs(no_of_eq, no_of_trait)) ! number of observations
ALLOCATE(temp_effects(no_of_trait * no_of_fixed + 1)) ! temporary array for fixed and animal effects
ALLOCATE(effects(no_of_trait, no_of_fixed + 1)) ! array for processing fixed and animal effects
ALLOCATE(observations(no_of_trait)) ! array for observations
ALLOCATE(observations_mbr(no_of_trait)) ! array for observations multiplied by rvcv
xpy = 0
nobs = 0
! open file for data processing procedure
OPEN(UNIT = 31, FILE = 'ongoing.dat', STATUS = 'REPLACE', ACTION = 'WRITE')
! open parameter file
OPEN(UNIT = 21, FILE = par_filename, STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg)
IF (status /= 0) THEN ! file open failed
WRITE (*,'(1X, A, A)') 'parameter file open failed -- error message : ', error_msg
STOP
END IF
! residual variance-covariance
READ(UNIT = 21, FMT = *) rvcv
WRITE(31,*) 'residual variance-covariance = ', rvcv
! genetic variance-covariance
READ(UNIT = 21, FMT = *) gvcv
WRITE(31,*) 'genetic variance-covariance = ', gvcv
! trait combination
DO k = 1, no_of_tc
READ(UNIT = 21, FMT = *) trait_combi
WRITE(31,*) 'trait combination ', k , '=', trait_combi
temp_rvcv = rvcv * trait_combi
WRITE(31,*) 'trati combination residual vcv', temp_rvcv
z = 0.0
CALL DJNVHF(temp_rvcv,no_of_trait,iw,z,mr) ! inverse of upper diagonal matrix(residual)
WRITE(31,*) 'inverse of trait combination rvcv', temp_rvcv
rvcv_tc(k,:) = temp_rvcv ! sotre the inverse of upper diagonal matrix(residual)
DO i = 1, no_of_trait
DO j = 1, no_of_trait
rvcv_tc_m(k,i,j) = temp_rvcv(IHMSSF(i,j,no_of_trait))
END DO
END DO
END DO
! open data file
OPEN(UNIT = 11, FILE = data_filename, STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg)
IF (status /= 0) THEN ! file open failed
WRITE (*,'(1X, A, A)') 'Data file open failed -- error message : ', error_msg
STOP
END IF
! open data file for writing left-hand side
OPEN(UNIT = 12, FILE = 'lhs.dat', STATUS = 'REPLACE', ACTION = 'WRITE')
! read each line
DO
READ(UNIT = 11, FMT = *, IOSTAT = status) temp_effects, observations
IF (status /= 0) EXIT ! end of file
! copy temp_effects(fixed) to effects
DO i = 1, no_of_trait
effects(i,1:no_of_fixed) = temp_effects((i-1) * no_of_fixed + 1 : i * no_of_fixed)
END DO
! copy temp_effects(animal) to effects
DO i = 1, no_of_trait
effects(i,no_of_fixed + 1) = temp_effects(no_of_trait * no_of_fixed + 1)
END DO
! find trait combination number
tc_no = 0
DO i = 1, no_of_trait
IF (observations(i) /= 0.0) tc_no = tc_no + 2 ** (i - 1)
END DO
observations_mbr = MATMUL(rvcv_tc_m(tc_no,:,:), observations)
! left hand side
DO i = 1, no_of_trait
IF (effects(i,1) /= 0) THEN
DO j = 1, no_of_fixed + 1
IF (j < no_of_fixed + 1) THEN
WRITE(UNIT = 12, FMT = *) effects(i,j), i, ((effects(k,l),l = 1, no_of_fixed + 1), k = 1, no_of_trait), '1', tc_no
ELSE
WRITE(UNIT = 12, FMT = *) effects(i,j), i, ((effects(k,l),l = 1, no_of_fixed + 1), k = 1, no_of_trait), '2', tc_no
END IF
END DO
END IF
END DO
! right hand side
DO i = 1, no_of_trait
IF (effects(i,1) /= 0) THEN
DO j = 1, no_of_fixed + 1
xpy(effects(i,j),i) = xpy(effects(i,j),i) + observations_mbr(i)
nobs(effects(i,j),i) = nobs(effects(i,j),i) + 1
END DO
END IF
END DO
END DO
CLOSE(11) ! close data file
! open pedigree file
OPEN(UNIT = 13, FILE = pedi_filename, STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg)
IF (status /= 0) THEN ! file open failed
WRITE (*,'(1X, A, A)') 'pedigree file open failed -- error message : ', error_msg
STOP
END IF
! read each line
DO
READ(UNIT = 13, FMT = *, IOSTAT = status) animal, sire, dam
IF (status /= 0) EXIT ! end of file
! write animal effect
! DO i = 1, no_of_trait
IF ( sire /= 0 .AND. dam /= 0) THEN
WRITE(UNIT = 12, FMT = *) animal, zero, animal, sire, dam, (zero,k=1,no_of_trait*(no_of_fixed+1)-3), '1001', zero
WRITE(UNIT = 12, FMT = *) sire, zero, animal, dam, sire, (zero,k=1,no_of_trait*(no_of_fixed+1)-3), '1002', zero
WRITE(UNIT = 12, FMT = *) dam, zero, animal, sire, dam, (zero,k=1,no_of_trait*(no_of_fixed+1)-3), '1002', zero
ELSE IF (sire /= 0 .AND. dam == 0) THEN
WRITE(UNIT = 12, FMT = *) animal, zero, animal, sire, dam, (zero,k=1,no_of_trait*(no_of_fixed+1)-3), '1003', zero
WRITE(UNIT = 12, FMT = *) sire, zero, animal, sire, dam, (zero,k=1,no_of_trait*(no_of_fixed+1)-3), '1004', zero
ELSE IF (sire == 0 .AND. dam /= 0) THEN
WRITE(UNIT = 12, FMT = *) animal, zero, animal, dam, sire, (zero,k=1,no_of_trait*(no_of_fixed+1)-3), '1003', zero
WRITE(UNIT = 12, FMT = *) dam, zero, animal, dam, sire, (zero,k=1,no_of_trait*(no_of_fixed+1)-3), '1004', zero
ELSE
WRITE(UNIT = 12, FMT = *) animal, zero, animal, sire, dam, (zero,k=1,no_of_trait*(no_of_fixed+1)-3), '1005', zero
END IF
! END DO
END DO ! end of reading datas
CLOSE(13) ! close animal file
CLOSE(12) ! close file for left-hand side
! open file for right hand side
OPEN(UNIT = 14, FILE = 'rhs.dat', STATUS = 'REPLACE', ACTION = 'WRITE')
! write xpy
DO i = 1, no_of_eq
WRITE(UNIT = 14, FMT = *) xpy(i,:), nobs(i,:)
END DO
CLOSE(14) ! close file for right-hand side
END PROGRAM mt_dl_animal_model3_setup
위 프로그램을 mt_dl_animal_model3_setup.f95로 저장
역행렬 구하는 프로그램
MODULE gi
CONTAINS
SUBROUTINE geninv(x,gix)
! subroutine name - genive, generalized inverse
! programmer - Park Byoungho
! usage - CALL genive(x, gix)
! purpose - to find out a generalized inverse of matrix using ihmssf and djnvhf
! ihmssf = matrix <-> vector using upper diagonal element
! djnvhf = solve generalized inverse of vector
! Date - 2009.7.7.
! last update - 2009.7.10. : precision of x and gix changes single precision to double precision
REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: x ! square matrix to inverse
REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: gix ! square matirx after inverse
INTEGER :: i, j, no_of_eq, size_x, mr
REAL(KIND = 8), ALLOCATABLE :: upperd(:)
INTEGER, DIMENSION(99) :: iw
REAL(KIND = 8) :: z
no_of_eq = SIZE(x,1)
ALLOCATE(upperd((no_of_eq*(no_of_eq +1)/2)))
upperd = 0
DO i = 1, no_of_eq
DO j = 1, no_of_eq
upperd(ihmssf(i,j,no_of_eq))=x(i,j)
END DO
END DO
!WRITE(*,*) 'upperd = ', upperd
z = 0.0
CALL DJNVHF(upperd,no_of_eq,iw,z,mr)
DO i = 1, no_of_eq
DO j = 1, no_of_eq
gix(i,j) = upperd(ihmssf(i,j,no_of_eq))
END DO
END DO
END SUBROUTINE geninv
SUBROUTINE DJNVHF(A,N,IS,Z,IR)
!
! MATRIX INVERSION SUBROUTINE
!
! 'A' IS A HALF-STORED MATRIX TO BE INVERTED. THE INVERTED
! RESULT IS RETURNED IN 'A'.
! 'N' IS THE ORDER OF THE MATRIX TO BE INVERTED.
! 'IS' IS A WORK VECTOR.
! 'Z' IS SET TO ZERO ON CALLS TO THIS SUBROUTINE FROM THE BEEF
! SIRE MONITORING SYSTEM.
! 'IR' RETURNS THE RANK OF THE MATRIX.
!
! THIS SUBROUTINE CALLS THE SUBROUTINE DREARN.
!
DIMENSION A(1),IS(1)
REAL*8 A,BIG,Z,RECIP,R
201 NP=(N*(N+1))/2
IF(N-1) 21,23,96
21 WRITE(6,22)
22 FORMAT(' N.LT.1')
23 IF(A(1)) 25,24,25
24 IR=0
26 RETURN
25 IR=1
A(1)=1.D0/A(1)
RETURN
96 DO 50 L=1,N
BIG=0.D1
DO 2 I=L,N
II=-(I*(I-3))/2+N*(I-1)
IF(DABS(A(II))-DABS(BIG)) 2,2,1
1 IS(L)=I
BIG=A(II)
2 CONTINUE
IF(DABS(BIG)-Z) 63,63,62
63 IR=L-1
IF(IR) 26,26,98
98 DO 64 I=1,L
LI=-(I*(I-1))/2+N*(I-1)
DO 64 J=L,N
LIJ=LI+J
64 A(LIJ)=-0.D1
IF(L-N) 66,69,67
67 STOP 67
66 LP1=L+1
DO 68 I=LP1,N
LI=-(I*(I-1))/2+N*(I-1)
DO 68 J=I,N
LIJ=LI+J
68 A(LIJ)=-0.D1
69 IF(IR-2) 56,27,27
56 CALL DREARN(A,N,IS,1)
DO 57 I=1,NP
57 A(I)=-A(I)
RETURN
27 LP1=L+1
DO 65 I=2,L
J=LP1-I
65 CALL DREARN(A,N,IS,J)
DO 101 I=1,NP
101 A(I)=-A(I)
RETURN
62 CALL DREARN(A,N,IS,L)
28 K3=N*(L-1)-(L*(L-3))/2
RECIP=1./A(K3)
A(K3)=-RECIP
DO 50 I=1,N
IF(I-L) 6,50,7
6 K11=N*(I-1)-(I*(I-3)/2)
K1=K11+L-I
GO TO 8
7 K11=N*(I-1)-(I*(I-3))/2
K1=K3+I-L
8 R=A(K1)*RECIP
301 DO 12 J=I,N
K4=K11+J-I
IF(J-L) 10,12,11
10 K5=N*(J-1)-(J*(J-3)/2)+L-J
GO TO 14
11 K5=K3+J-L
14 A(K4)=A(K4)-R*A(K5)
12 CONTINUE
A(K1)=R
50 CONTINUE
NP=(N*(N+1))/2
DO 100 I=1,NP
100 A(I)=-A(I)
NP1=N+1
DO 61 I=2,N
L=NP1-I
61 CALL DREARN(A,N,IS,L)
IR=N
RETURN
END SUBROUTINE
SUBROUTINE DREARN(A,N,IS,L)
!
! THIS IS A SUPPORT SUBROUTINE CALLED BY THE MATRIX INVERSION
! SUBROUTINE DJNVHF.
!
DIMENSION A(1),IS(1)
DOUBLE PRECISION A,SAVE
ISL=IS(L)
IF(ISL-L) 3,28,4
3 STOP 3
4 LM1=L-1
IF(LM1) 22,22,5
5 DO 21 I=1,LM1
IL=-(I*(I-1))/2+N*(I-1)+L
IISL=IL-L+ISL
SAVE=A(IL)
A(IL)=A(IISL)
21 A(IISL)=SAVE
22 LP1=L+1
ISLM1=ISL-1
IF(LP1-ISLM1) 23,23,25
23 DO 24 I=LP1,ISLM1
LI=-(L*(L-1))/2+N*(L-1)+I
IISL=-(I*(I-1))/2+N*(I-1)+ISL
SAVE=A(LI)
A(LI)=A(IISL)
24 A(IISL)=SAVE
25 ISLP1=ISL+1
IF(ISLP1-N) 26,26,38
26 DO 27 I=ISLP1,N
ISLI=-(ISL*ISLM1)/2+N*ISLM1+I
LI=-(L*LM1)/2+N*LM1+I
SAVE=A(ISLI)
A(ISLI)=A(LI)
27 A(LI)=SAVE
38 LL=-(L*(L-3))/2+N*(L-1)
ISLISL=-(ISL*(ISL-3))/2+ISLM1*N
SAVE=A(LL)
A(LL)=A(ISLISL)
A(ISLISL)=SAVE
28 RETURN
END SUBROUTINE
FUNCTION IHMSSF(IR,IC,N)
JR=IR
JC=IC
IF(IR-IC)40,40,41
41 JC=IR
JR=IC
40 NR=JR-1
KK=(JR*NR)/2
KK=KK+NR*(N-NR)
KK=KK+JC-NR
IHMSSF=KK
RETURN
END FUNCTION
END MODULE gi
위 프로그램을 inverse.f95로 저장
컴파일 gfortran inverse.f95 mt_dl_animal_model3_setup.f95 -o mt_dl_animal_model3_setup
위 프로그램으로 rhs.dat
0.55533141377722850 0.61532929747454823 0.29148183944825640 2 2 1
0.53367565173619258 0.0000000000000000 0.31501003224321361 2 0 1
0.0000000000000000 0.22801333236648663 0.34632446883858659 0 1 1
0.0000000000000000 0.34632187057857344 0.23792049568881923 0 1 1
0.28463141377722845 0.0000000000000000 0.0000000000000000 1 0 0
0.53367565173619258 0.61532929747454823 0.0000000000000000 2 2 0
0.27070000000000000 0.57433520294506013 0.60649187169147001 1 2 2
0.0000000000000000 0.0000000000000000 0.58424496452740582 0 0 2
0.0000000000000000 0.29789664143982270 0.29148183944825640 0 1 1
0.28463141377722845 0.31743265603472554 0.31501003224321361 1 1 1
0.30934684599986850 0.34632187057857344 0.34632446883858659 1 1 1
0.22432880573632405 0.22801333236648663 0.23792049568881923 1 1 1
0.27070000000000000 0.0000000000000000 0.0000000000000000 1 0 0
0.0000000000000000 0.0000000000000000 0.0000000000000000 0 0 0
0.0000000000000000 0.0000000000000000 0.0000000000000000 0 0 0
0.0000000000000000 0.0000000000000000 0.0000000000000000 0 0 0
0.0000000000000000 0.0000000000000000 0.0000000000000000 0 0 0
0.0000000000000000 0.0000000000000000 0.0000000000000000 0 0 0
위 프로그램으로 구한 lhs.dat
1 2 0 0 9 1 6 9 1 7 9 1 6
6 2 0 0 9 1 6 9 1 7 9 1 6
9 2 0 0 9 1 6 9 1 7 9 2 6
1 3 0 0 9 1 6 9 1 7 9 1 6
7 3 0 0 9 1 6 9 1 7 9 1 6
9 3 0 0 9 1 6 9 1 7 9 2 6
1 1 1 5 10 1 6 10 2 7 10 1 7
5 1 1 5 10 1 6 10 2 7 10 1 7
10 1 1 5 10 1 6 10 2 7 10 2 7
1 2 1 5 10 1 6 10 2 7 10 1 7
6 2 1 5 10 1 6 10 2 7 10 1 7
10 2 1 5 10 1 6 10 2 7 10 2 7
2 3 1 5 10 1 6 10 2 7 10 1 7
7 3 1 5 10 1 6 10 2 7 10 1 7
10 3 1 5 10 1 6 10 2 7 10 2 7
2 1 2 6 11 4 7 11 3 8 11 1 7
6 1 2 6 11 4 7 11 3 8 11 1 7
11 1 2 6 11 4 7 11 3 8 11 2 7
4 2 2 6 11 4 7 11 3 8 11 1 7
7 2 2 6 11 4 7 11 3 8 11 1 7
11 2 2 6 11 4 7 11 3 8 11 2 7
3 3 2 6 11 4 7 11 3 8 11 1 7
8 3 2 6 11 4 7 11 3 8 11 1 7
11 3 2 6 11 4 7 11 3 8 11 2 7
2 1 2 6 12 3 7 12 4 8 12 1 7
6 1 2 6 12 3 7 12 4 8 12 1 7
12 1 2 6 12 3 7 12 4 8 12 2 7
3 2 2 6 12 3 7 12 4 8 12 1 7
7 2 2 6 12 3 7 12 4 8 12 1 7
12 2 2 6 12 3 7 12 4 8 12 2 7
4 3 2 6 12 3 7 12 4 8 12 1 7
8 3 2 6 12 3 7 12 4 8 12 1 7
12 3 2 6 12 3 7 12 4 8 12 2 7
1 1 1 7 13 0 0 13 0 0 13 1 1
7 1 1 7 13 0 0 13 0 0 13 1 1
13 1 1 7 13 0 0 13 0 0 13 2 1
9 0 9 14 17 0 0 0 0 0 0 1001 0
14 0 9 17 14 0 0 0 0 0 0 1002 0
17 0 9 14 17 0 0 0 0 0 0 1002 0
10 0 10 15 17 0 0 0 0 0 0 1001 0
15 0 10 17 15 0 0 0 0 0 0 1002 0
17 0 10 15 17 0 0 0 0 0 0 1002 0
11 0 11 16 18 0 0 0 0 0 0 1001 0
16 0 11 18 16 0 0 0 0 0 0 1002 0
18 0 11 16 18 0 0 0 0 0 0 1002 0
12 0 12 14 17 0 0 0 0 0 0 1001 0
14 0 12 17 14 0 0 0 0 0 0 1002 0
17 0 12 14 17 0 0 0 0 0 0 1002 0
13 0 13 15 18 0 0 0 0 0 0 1001 0
15 0 13 18 15 0 0 0 0 0 0 1002 0
18 0 13 15 18 0 0 0 0 0 0 1002 0
14 0 14 0 0 0 0 0 0 0 0 1005 0
15 0 15 0 0 0 0 0 0 0 0 1005 0
16 0 16 0 0 0 0 0 0 0 0 1005 0
17 0 17 0 0 0 0 0 0 0 0 1005 0
18 0 18 0 0 0 0 0 0 0 0 1005 0
lhs.dat의 정렬 : sort -n -o sorted_lhs.dat lhs.dat
정렬한 lhs.dat : sorted_lhs.dat
1 1 1 5 10 1 6 10 2 7 10 1 7
1 1 1 7 13 0 0 13 0 0 13 1 1
1 2 0 0 9 1 6 9 1 7 9 1 6
1 2 1 5 10 1 6 10 2 7 10 1 7
1 3 0 0 9 1 6 9 1 7 9 1 6
2 1 2 6 11 4 7 11 3 8 11 1 7
2 1 2 6 12 3 7 12 4 8 12 1 7
2 3 1 5 10 1 6 10 2 7 10 1 7
3 2 2 6 12 3 7 12 4 8 12 1 7
3 3 2 6 11 4 7 11 3 8 11 1 7
4 2 2 6 11 4 7 11 3 8 11 1 7
4 3 2 6 12 3 7 12 4 8 12 1 7
5 1 1 5 10 1 6 10 2 7 10 1 7
6 1 2 6 11 4 7 11 3 8 11 1 7
6 1 2 6 12 3 7 12 4 8 12 1 7
6 2 0 0 9 1 6 9 1 7 9 1 6
6 2 1 5 10 1 6 10 2 7 10 1 7
7 1 1 7 13 0 0 13 0 0 13 1 1
7 2 2 6 11 4 7 11 3 8 11 1 7
7 2 2 6 12 3 7 12 4 8 12 1 7
7 3 0 0 9 1 6 9 1 7 9 1 6
7 3 1 5 10 1 6 10 2 7 10 1 7
8 3 2 6 11 4 7 11 3 8 11 1 7
8 3 2 6 12 3 7 12 4 8 12 1 7
9 0 9 14 17 0 0 0 0 0 0 1001 0
9 2 0 0 9 1 6 9 1 7 9 2 6
9 3 0 0 9 1 6 9 1 7 9 2 6
10 0 10 15 17 0 0 0 0 0 0 1001 0
10 1 1 5 10 1 6 10 2 7 10 2 7
10 2 1 5 10 1 6 10 2 7 10 2 7
10 3 1 5 10 1 6 10 2 7 10 2 7
11 0 11 16 18 0 0 0 0 0 0 1001 0
11 1 2 6 11 4 7 11 3 8 11 2 7
11 2 2 6 11 4 7 11 3 8 11 2 7
11 3 2 6 11 4 7 11 3 8 11 2 7
12 0 12 14 17 0 0 0 0 0 0 1001 0
12 1 2 6 12 3 7 12 4 8 12 2 7
12 2 2 6 12 3 7 12 4 8 12 2 7
12 3 2 6 12 3 7 12 4 8 12 2 7
13 0 13 15 18 0 0 0 0 0 0 1001 0
13 1 1 7 13 0 0 13 0 0 13 2 1
14 0 9 17 14 0 0 0 0 0 0 1002 0
14 0 12 17 14 0 0 0 0 0 0 1002 0
14 0 14 0 0 0 0 0 0 0 0 1005 0
15 0 10 17 15 0 0 0 0 0 0 1002 0
15 0 13 18 15 0 0 0 0 0 0 1002 0
15 0 15 0 0 0 0 0 0 0 0 1005 0
16 0 11 18 16 0 0 0 0 0 0 1002 0
16 0 16 0 0 0 0 0 0 0 0 1005 0
17 0 9 14 17 0 0 0 0 0 0 1002 0
17 0 10 15 17 0 0 0 0 0 0 1002 0
17 0 12 14 17 0 0 0 0 0 0 1002 0
17 0 17 0 0 0 0 0 0 0 0 1005 0
18 0 11 16 18 0 0 0 0 0 0 1002 0
18 0 13 15 18 0 0 0 0 0 0 1002 0
18 0 18 0 0 0 0 0 0 0 0 1005 0
해를 구하는 프로그램
PROGRAM mt_dl_animal_model3_solve ! program name - multiple trait, different fixed level, animal_model_solve
! programmer - Park Byoungho
! usage - mt_dl_animal_model3_solve
! purpose : read sorted_lhs.dat and rhs.dat which are made in animal_model_setup.exe and are sorted
! find a solution using Gauss-Seidel iteration
! Date - 2009.7.7
! update - 209.7.14. inv_diag_ele(2,2) -> allocate(inv_diag_ele(no_of_trait,no_of_trait))
! data_type(i,9) -> data_type = lhs(i,2+no_of_trait*(no_of_fixed+1)+1)
! tc_no(i,10) -> tc_no = lhs(i,2+no_of_trait*(no_of_fixed+1)+2) USE gi IMPLICIT NONE
! data dictionary
INTEGER, PARAMETER :: no_of_trait = 3 ! number of trait
INTEGER, PARAMETER :: no_of_fixed = 2 ! number of fixed effects
INTEGER, PARAMETER :: no_of_eq = 18 ! number of equation INTEGER :: no_of_ud ! number of upper digonal
INTEGER :: no_of_tc ! number of trait combination REAL(KIND = 8), ALLOCATABLE :: gvcv(:) ! upper diagonal part of genetic variance-covariance, (no_of_trait)*(no_of_trait + 1) / 2 dimension
REAL(KIND = 8), ALLOCATABLE :: gvcv_m(:,:)
REAL(KIND = 8), ALLOCATABLE :: rvcv(:) ! upper diagonal part of residual variance-covariance, (no_of_trait)*(no_of_trait + 1) / 2 dimension
REAL(KIND = 8), ALLOCATABLE :: temp_rvcv(:) ! temporary rvcv
REAL(KIND = 8), ALLOCATABLE :: rvcv_tc(:,:) ! rvcv according to trait combination and then inverse
REAL(KIND = 8), ALLOCATABLE :: rvcv_tc_m(:,:,:) ! rvcv matrix according to trait combination and then inverse
INTEGER, ALLOCATABLE :: trait_combi(:) ! trait combination
INTEGER :: tc_no ! trait combination number INTEGER, DIMENSION(99) :: iw ! for inverse
REAL(KIND = 8) :: z ! for inverse
INTEGER :: mr ! for inverse
INTEGER :: pre_eq_no ! previous equation number
INTEGER :: cur_eq_no ! current equation number
INTEGER :: pre_trait_no ! previous trait number
INTEGER :: cur_trait_no ! current trati number
INTEGER, ALLOCATABLE :: loc_of_nz(:,:) ! location of none-zero element
INTEGER, ALLOCATABLE :: loc_of_nz_a(:,:) ! location of none-zero element
INTEGER :: data_type ! data type : 1:fixed, 2-6:animal REAL(KIND = 8), ALLOCATABLE :: diag_ele(:,:) ! diagonal element of lhs
REAL(KIND = 8), ALLOCATABLE :: inv_diag_ele(:,:) ! inverse of diagonal INTEGER :: no_of_lhs ! number of lhs lines
INTEGER, ALLOCATABLE :: lhs(:,:) ! left-hand side
REAL(KIND = 8), ALLOCATABLE :: rhs(:,:) ! right-hand side, x'y
REAL(KIND = 8), ALLOCATABLE :: temp_rhs(:) ! right-hand side one element
REAL(KIND = 8), ALLOCATABLE :: solutions(:,:) ! solutions
REAL(KIND = 8), ALLOCATABLE :: pre_sol(:) ! previous solution
REAL(KIND = 8), ALLOCATABLE :: sum_sol(:) ! sum solution INTEGER :: iteration ! iteration
REAL(KIND = 8) :: epsilon ! sum of squares between old and new solutions
INTEGER :: i,j,k ! loop
INTEGER :: status ! i/o status
CHARACTER(LEN = 40) :: error_msg ! error message
INTEGER, PARAMETER :: MAX_ITER = 50 ! maximum number of iteration
REAL(KIND = 8), PARAMETER :: criteria = 1.E-12 ! criteria for stopping no_of_ud = (no_of_trait)*(no_of_trait + 1) / 2 ! number of upper diagonal of variance-covariance
no_of_tc = 0 ! number of trait combination
DO i = 1, no_of_trait
no_of_tc = no_of_tc + 2 ** (i - 1)
END DO ALLOCATE(gvcv(no_of_ud)) ! genetic variance-covariance
ALLOCATE(gvcv_m(no_of_trait, no_of_trait)) ! genetic variance-covariance
ALLOCATE(rvcv(no_of_ud)) ! residual variance-covariance
ALLOCATE(temp_rvcv(no_of_ud)) ! temporary residual variance-covariance
ALLOCATE(rvcv_tc(no_of_tc, no_of_ud)) ! residual variance-covariance according to trait combination
ALLOCATE(rvcv_tc_m(no_of_tc, no_of_trait, no_of_trait)) ! residual variance-covariance matrix according to trait combination
ALLOCATE(trait_combi(no_of_ud)) ! trait combination
ALLOCATE(rhs(no_of_eq, no_of_trait * 2)) ! right hand side
ALLOCATE(loc_of_nz(no_of_trait,no_of_fixed + 1)) ! location of none-zero element
ALLOCATE(loc_of_nz_a(no_of_trait,no_of_fixed))
ALLOCATE(sum_sol(no_of_trait)) ! temporary solution
ALLOCATE(solutions(no_of_eq, no_of_trait)) ! solution
ALLOCATE(pre_sol(no_of_trait))
ALLOCATE(diag_ele(no_of_trait,no_of_trait))
ALLOCATE(inv_diag_ele(no_of_trait,no_of_trait))
ALLOCATE(temp_rhs(no_of_trait)) ! initialiaze solutions to zero
solutions = 0
! open file for data processing procedure
OPEN(UNIT = 31, FILE = 'ongoing_solve.dat', STATUS = 'REPLACE', ACTION = 'WRITE')
! open parameter file
OPEN(UNIT = 21, FILE = 'vcvtrt_3tr.par', STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg)
IF (status /= 0) THEN ! file open failed
WRITE (*,'(1X, A, A)') 'parameter file open failed -- error message : ', error_msg
STOP
END IF ! residual variance-covariance
READ(UNIT = 21, FMT = *) rvcv
write(31,*) 'residual variance-covariance = ', rvcv
! genetic variance-covariance
READ(UNIT = 21, FMT = *) gvcv
write(31,*) 'genetic variance-covariance = ', gvcv ! inverse of genetic variance-covariance
z = 0.0
CALL DJNVHF(gvcv,no_of_trait,iw,z,mr)
DO i = 1, no_of_trait
DO j = 1, no_of_trait
gvcv_m(i,j) = gvcv(IHMSSF(i,j,no_of_trait))
END DO
END DO
write(31,*) 'inverse of genetic variance-covariance', ((gvcv_m(i,j),i = 1, no_of_trait),j = 1, no_of_trait)
! trait combination
DO k = 1, no_of_tc
READ(UNIT = 21, FMT = *) trait_combi
write(31,*) 'trait combination ', k , '=', trait_combi
temp_rvcv = rvcv * trait_combi
write(31,*) 'trati combination residual vcv', temp_rvcv
z = 0.0
CALL DJNVHF(temp_rvcv,no_of_trait,iw,z,mr) ! inverse of upper diagonal matrix(residual)
write(31,*) 'inverse of trait combination rvcv', temp_rvcv
rvcv_tc(k,:) = temp_rvcv ! sotre the inverse of upper diagonal matrix(residual)
DO i = 1, no_of_trait
DO j = 1, no_of_trait
rvcv_tc_m(k,i,j) = temp_rvcv(IHMSSF(i,j,no_of_trait))
END DO
END DO
END DO
!open left-han side file
OPEN(UNIT = 11, FILE = 'sorted_lhs.dat', STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg)
IF (status /= 0) THEN ! file open failed
WRITE (*,'(1X, A, A)') 'Sorted lhs file open failed -- error message : ', error_msg
STOP
END IF ! count the lines of lhs
no_of_lhs = 0
DO
READ(UNIT = 11, FMT = *, IOSTAT = status)
IF (status /= 0) EXIT ! reach end of file
no_of_lhs = no_of_lhs + 1
END DO ALLOCATE(lhs(no_of_lhs, no_of_trait * (no_of_fixed + 1) + 4))
write(31,*) 'number of lines of left hand side = ', no_of_lhs ! store lhs to array
REWIND(11)
DO i = 1, no_of_lhs
READ(UNIT = 11, FMT = *, IOSTAT = status) lhs(i,:)
IF (status /= 0) EXIT ! reach end of file
END DO
CLOSE(11) write(31,*) 'left hand side, lines = ', no_of_lhs
DO i = 1, no_of_lhs
write(31,*) lhs(i,:)
END DO !open right-hand side file
OPEN(UNIT = 12, FILE = 'rhs.dat', STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg) IF (status /= 0) THEN ! file open failed
WRITE (*,'(1X, A, A)') 'RHS file open failed -- error message : ', error_msg
STOP
END IF
! read and store right-hand side
DO i = 1, no_of_eq
READ(UNIT = 12,FMT = *) rhs(i,:)
END DO
close(12) write(31,*) 'right hand side, number of equations = ', no_of_eq
DO i = 1, no_of_eq
write(31,*) i, rhs(i,:)
END DO ! start iteration
DO iteration = 1, MAX_ITER
pre_eq_no = 1
diag_ele = 0.0
temp_rhs = rhs(lhs(1,1),1:no_of_trait)
epsilon = 0.0 ! start reading and processing each line of lhs
DO i = 1, no_of_lhs cur_eq_no = lhs(i,1)
cur_trait_no = lhs(i,2)
data_type = lhs(i,2+no_of_trait*(no_of_fixed+1)+1)
tc_no = lhs(i,2+no_of_trait*(no_of_fixed+1)+2) ! if new equation, calculate solution
IF (pre_eq_no /= cur_eq_no) THEN
pre_sol = solutions(pre_eq_no, :)
CALL geninv(diag_ele, inv_diag_ele)
write(31,*) 'calculate solutions'
write(31,*) 'diag_ele = ', diag_ele
write(31,*) 'inv_diag_ele = ', inv_diag_ele
write(31,*) 'temp_rhs = ', temp_rhs
solutions(pre_eq_no, :) = MATMUL(inv_diag_ele, temp_rhs)
epsilon = epsilon + SUM((pre_sol - solutions(pre_eq_no, :))**2)
write(31,*) 'new sol =', MATMUL(inv_diag_ele, temp_rhs)
diag_ele = 0.0
temp_rhs = rhs(lhs(i,1), 1:no_of_trait) pre_eq_no = cur_eq_no
END IF ! adjust rhs and add diag_ele
SELECT CASE(data_type) CASE(1) ! fixed effect
DO k = 1, no_of_trait
loc_of_nz(k,:) = lhs(i,(k-1)*(no_of_fixed + 1)+ 3 : k * (no_of_fixed + 1) + 2)
END DO
sum_sol = 0
DO k = 1, no_of_trait
If (loc_of_nz(k,1) /= 0) THEN
sum_sol(k) = SUM(solutions(loc_of_nz(k,:),k))
END IF
END DO
temp_rhs(cur_trait_no) = temp_rhs(cur_trait_no) &
- SUM(rvcv_tc_m(tc_no,cur_trait_no,:) * sum_sol) &
+ rvcv_tc_m(tc_no,cur_trait_no,cur_trait_no) * solutions(cur_eq_no, cur_trait_no)
diag_ele(cur_trait_no,cur_trait_no) = diag_ele(cur_trait_no,cur_trait_no) + rvcv_tc_m(tc_no, cur_trait_no, cur_trait_no) CASE(2) ! animal as a fixed effect
DO k = 1, no_of_trait
loc_of_nz_a(k,:) = lhs(i,(k-1)*(no_of_fixed + 1)+ 3 : k * (no_of_fixed + 1) + 1)
END DO sum_sol = 0
DO k = 1, no_of_trait
If (loc_of_nz_a(k,1) /= 0) THEN
sum_sol(k) = SUM(solutions(loc_of_nz_a(k,:),k))
END IF
END DO temp_rhs(cur_trait_no) = temp_rhs(cur_trait_no) &
- SUM(rvcv_tc_m(tc_no,cur_trait_no,:) * sum_sol) diag_ele(cur_trait_no,:) = diag_ele(cur_trait_no,:) + rvcv_tc_m(tc_no, cur_trait_no,:)
CASE(1001) ! animal, both parents are known, animal is diagonal temp_rhs = temp_rhs - (-1.) * MATMUL(gvcv_m,solutions(lhs(i,4),:)) &
- (-1.) * MATMUL(gvcv_m,solutions(lhs(i,5),:)) diag_ele = diag_ele + 2.0 * gvcv_m CASE(1002) ! animal, both parents are known, sire or dam is diagonal
write(31,*) 'temp_rhs = ', temp_rhs
write(31,*) 'diag_ele = ', diag_ele
write(31,*) 'lhs(i,3) = ', lhs(i,3)
write(31,*) 'solutions(lhs(i,3) =', solutions(lhs(i,3),:)
temp_rhs = temp_rhs - (-1.0) * MATMUL(gvcv_m,solutions(lhs(i,3),:)) &
- ( 0.5) * MATMUL(gvcv_m,solutions(lhs(i,4),:)) diag_ele = diag_ele + ( 0.5) * gvcv_m
write(31,*) 'temp_rhs = ', temp_rhs
write(31,*) 'diag_ele = ', diag_ele CASE(1003) ! animal, one parent is known, animal is diagonal
temp_rhs = temp_rhs - (-2.0/3.0) * MATMUL(gvcv_m,solutions(lhs(i,4),:)) diag_ele = diag_ele + (4.0/3.0) * gvcv_m CASE(1004) ! animal effect, one parents is known, parent is diagonal
temp_rhs = temp_rhs - (-2.0/3.0) * MATMUL(gvcv_m,solutions(lhs(i,3),:)) diag_ele = diag_ele + (1.0/3.0) * gvcv_m CASE(1005) ! animal effect, no parents is known
diag_ele = diag_ele + 1.0 * gvcv_m
write(31,*) 'case 1005, diag_ele = ', diag_ele
CASE DEFAULT
WRITE(31,*) cur_eq_no, 'equation number has a ', data_type,' Invalid data_type. ', i,'th line of lhs'
END SELECT END DO
! end reading lhs ! calculate last solution
pre_sol = solutions(pre_eq_no, :)
CALL geninv(diag_ele, inv_diag_ele)
solutions(pre_eq_no, :) = MATMUL(inv_diag_ele, temp_rhs)
epsilon = (epsilon + SUM((pre_sol - solutions(pre_eq_no, :))**2))/(no_of_eq * no_of_trait) WRITE(31,*) iteration,'th iteration''s solutions'
DO k = 1, no_of_eq
WRITE(31,*) k, solutions(k,:)
END DO ! write iteration number and epsilon
write(*,*) 'iteration = ', iteration , ', epsilon = ', epsilon
IF (epsilon < criteria) THEN
EXIT
END IF END DO
! end iteration
CLOSE(31) ! open file for writing solution
OPEN(UNIT=13, FILE='sol.dat', STATUS='REPLACE', ACTION='WRITE', IOSTAT=status)
DO i = 1, no_of_eq
WRITE(13,*) i, solutions(i,:)
END DO
CLOSE(13)
END PROGRAM mt_dl_animal_model3_solve
위 프로그램을 mt_dl_animal_model3_solve.f95로 저장
컴파일 : gfortran inverse.f95 mt_dl_animal_model3_solve.f95 -o mt_dl_animal_model3_solve
위 프로그램을 구한 해 . sol.dat
1 8308.7891900281702 9578.4557094390275 9573.8483421794899
2 8004.3001912084628 0.0000000000000000 10498.868797644614
3 0.0000000000000000 7650.9025722133765 11399.268915019995
4 0.0000000000000000 11138.051660654915 8051.1897817375157
5 931.96551719210470 0.0000000000000000 0.0000000000000000
6 629.40119583558715 526.46396763173630 0.0000000000000000
7 -273.41821290427265 163.39076354750674 211.69614699415553
8 0.0000000000000000 0.0000000000000000 225.97203808115160
9 -183.66081924476862 -150.43353689705538 -101.28129508880323
10 -36.017903856112873 40.594179853707409 5.1718584319217555
11 292.90497160034624 146.45248598745229 128.14592527359096
12 -305.30774706542127 -172.18824770583947 -140.54869990905624
13 85.629022448683941 62.348888132333983 44.439260738691488
14 -171.25804230471451 -124.69777359537245 -88.878518917363479
15 24.805558117221523 51.471532479144905 24.805558107147725
16 146.45248559455848 73.226242835750497 64.072962490471284
17 -146.45248276807956 -73.226238957354951 -64.072958680859216
18 146.45248609657006 73.226243231169278 64.072962856736481
프로그램 컴파일 및 실행 화면
프로그램 파일
1247562147_dairy_3tr.dat
1247562147_pedi_3tr.dat
1247562147_vcvtrt_3tr.par
1247562147_mt_dl_animal_model3_setup.f95
1247562147_inverse.f95
1247562147_rhs.dat
1247562147_lhs.dat
1247562147_sorted_lhs.dat
1247562147_mt_dl_animal_model3_solve.f95
1247562147_sol.dat
모두 모아 한 파일
1247562147_3trait.zip
'Animal Breeding > Fortran program' 카테고리의 다른 글
| 5 형질 다형질 개체 모형 - same model but different level (0) | 2009.10.23 |
|---|---|
| 5형질 multple trait animal model - same model but different leve. (0) | 2009.10.09 |
| same model, different level, multiple trait animal model 3 - 2형질 (0) | 2009.07.13 |
| same model, different level, multiple trait animal model 1 (0) | 2009.06.16 |
| multiple trait animal model - solve (0) | 2009.06.06 |