same model, different level, multiple trait animal model 1

 

same model, different level, multiple trait animal model 1 - setup

 

same model, ie same fixed and random effect

 

but different level

 

 

사용한 자료

 

8 0 0 1 6 0 9570
9 1 5 1 6 9219 10530
10 2 6 4 7 10025 11490
11 2 6 3 7 7230 8600
12 1 7 0 0 8121 0
1st col : cowid2nd col : age 1 - fixed3rd col : contemporary group - fixed4th col : age 2 - fixed5th col : contemporary group - fixedmilk1 : observation 1milk2 : observation 2

 

위 데이터를 dairy1.dat로 저장

 

 

혈통 자료

 

8 13 16
9 14 16
10 15 17
11 13 16
12 14 17
13 0 0
14 0 0
15 0 0
16 0 0
17 0 0
animal, dam, sire

 

위 데이터를 pedi1.dat로 저장

 

 

파라미터 파일

 

30000 1150 31000
8000 4000 8300
1 0 0
0 0 1
1 1 1
1st row : residual variance - covariance 2nd row : genetic variance - covariance others : trait combination

 

위 데이터를 vcvtrt.par로 저장

 

 

프로그램

 

PROGRAM mt_dl_animal_model1_setup
! program name - multiple trait animal model setup
! same model but different fixed level
! programmer - Park Byoungho
! usage - mt_animal_model2_setup < parameter_file
! purpose : read data file and write left-hand side and right-hand side
! Date - 2009.6.15.
! update - none

USE gi

IMPLICIT NONE

! data dictionary
INTEGER, PARAMETER :: no_of_trait = 2 ! number of trait
INTEGER, PARAMETER :: no_of_fixed = 2 ! number of fixed effects
INTEGER, PARAMETER :: no_of_eq = 17 ! 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
INTEGER, ALLOCATABLE :: trait_combi(:) ! trait combination
INTEGER :: tc_no ! trait combination number

INTEGER :: cowid, age1, cg1, age2, cg2 ! animal and fixed effect
REAL(KIND = 8) :: milk1, milk2, r1, r2 ! data of trait for each line

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 ! do loop
REAL(KIND = 8) :: 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 = 'dairy1.dat'
! get pedigree file name
pedi_filename = 'pedi1.dat'
! get parameter file name
par_filename = 'vcvtrt.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))
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

xpy = 0
nobs = 0

 

! open file for data processing procedure
OPEN(UNIT = 31, FILE = 'ongoing_setup.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)
WRITE(31,*) 'inverse of trait combination rvcv', temp_rvcv
rvcv_tc(k,:) = temp_rvcv
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) cowid, age1, cg1, age2, cg2, milk1, milk2
IF (status /= 0) EXIT ! end of file

! find trait combination number
tc_no = 0
IF (milk1 /= 0.0) tc_no = tc_no + 1
IF (milk2 /= 0.0) tc_no = tc_no + 2 r1 = rvcv_tc(tc_no,1) * milk1 + rvcv_tc(tc_no,2) * milk2
r2 = rvcv_tc(tc_no,2) * milk1 + rvcv_tc(tc_no,3) * milk2 IF (tc_no == 1) THEN ! make LHS
WRITE(UNIT = 12, FMT = *) age1, cg1, cowid, zero, zero, tc_no, '1'
WRITE(UNIT = 12, FMT = *) cg1, age1, cowid, zero, zero, tc_no, '1'
WRITE(UNIT = 12, FMT = *) cowid, age1, cg1, zero, zero, tc_no, '1'

! make RHS, X'y
xpy(age1,1) = xpy(age1,1) + r1
xpy(cg1,1) = xpy(cg1,1) + r1
xpy(cowid,1) = xpy(cowid,1) + r1 ! count observations
nobs(age1,1) = nobs(age1,1) + 1
nobs(cg1,1) = nobs(cg1,1) + 1
nobs(cowid,1) = nobs(cowid,1) + 1
ELSE IF (tc_no == 2) THEN ! make LHS
WRITE(UNIT = 12, FMT = *) age2, cg2, cowid, zero, zero, tc_no, '2'
WRITE(UNIT = 12, FMT = *) cg2, age2, cowid, zero, zero, tc_no, '2'
WRITE(UNIT = 12, FMT = *) cowid, age2, cg2, zero, zero, tc_no, '2'

! make RHS, X'y
xpy(age2,2) = xpy(age2,2) + r2
xpy(cg2,2) = xpy(cg2,2) + r2
xpy(cowid,2) = xpy(cowid,2) + r2

! count observations
nobs(age2,2) = nobs(age2,2) + 1
nobs(cg2,2) = nobs(cg2,2) + 1
nobs(cowid,2) = nobs(cowid,2) + 1

ELSE IF (tc_no == 3) THEN ! make LHS
WRITE(UNIT = 12, FMT = *) age1, age2, cg1, cg2, cowid, tc_no, '3'
WRITE(UNIT = 12, FMT = *) age2, age1, cg1, cg2, cowid, tc_no, '4'
WRITE(UNIT = 12, FMT = *) cg1, age1, age2, cg2, cowid, tc_no, '5'
WRITE(UNIT = 12, FMT = *) cg2, age1, age2, cg1, cowid, tc_no, '6'
WRITE(UNIT = 12, FMT = *) cowid, age1, age2, cg1, cg2, tc_no, '7'

! make RHS, X'y
xpy(age1,1) = xpy(age1,1) + r1
xpy(age2,2) = xpy(age2,2) + r2
xpy(cg1,1) = xpy(cg1,1) + r1
xpy(cg2,2) = xpy(cg2,2) + r2
xpy(cowid,1) = xpy(cowid,1) + r1
xpy(cowid,2) = xpy(cowid,2) + r2 ! count observations
nobs(age1,1) = nobs(age1,1) + 1
nobs(age2,2) = nobs(age2,2) + 1
nobs(cg1,1) = nobs(cg1,1) + 1
nobs(cg2,2) = nobs(cg2,2) + 1
nobs(cowid,1) = nobs(cowid,1) + 1
nobs(cowid,2) = nobs(cowid,2) + 1

END IF

END DO ! end of reading data 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
IF ( sire /= 0 .AND. dam /= 0) THEN
WRITE(UNIT = 12, FMT = *) animal, sire, dam, zero, zero, zero, '1001'
WRITE(UNIT = 12, FMT = *) sire, animal, dam, zero, zero, zero, '1002'
WRITE(UNIT = 12, FMT = *) dam, animal, sire, zero, zero, zero, '1002'
ELSE IF (sire /= 0 .AND. dam == 0) THEN
WRITE(UNIT = 12, FMT = *) animal, sire, dam, zero, zero, zero, '1003'
WRITE(UNIT = 12, FMT = *) sire, animal, dam, zero, zero, zero, '1004'
ELSE IF (sire == 0 .AND. dam /= 0) THEN
WRITE(UNIT = 12, FMT = *) animal, dam, sire, zero, zero, zero, '1003'
WRITE(UNIT = 12, FMT = *) dam, animal, sire, zero, zero, zero, '1004'
ELSE
WRITE(UNIT = 12, FMT = *) animal, sire, dam, zero, zero, zero, '1005'
END IF
END DO ! end of reading data 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 sideEND PROGRAM mt_dl_animal_model1_setup
위 프로그램을 mt_dl_animal_model1_setup.f95로 저장

 

컴파일 : gfortran inverse.f95 mt_dl_animal_model1_setup.f95 -o mt_dl_animal_model1_setup

 

위 프로그램을 컴파일할 때 module 필요

 

다음은 모듈 프로그램

 

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로 저장

 

프로그램으로 구한 LHS

 

1 6 8 0 0 2 2
6 1 8 0 0 2 2
8 1 6 0 0 2 2
1 1 5 6 9 3 3
1 1 5 6 9 3 4
5 1 1 6 9 3 5
6 1 1 5 9 3 6
9 1 1 5 6 3 7
2 4 6 7 10 3 3
4 2 6 7 10 3 4
6 2 4 7 10 3 5
7 2 4 6 10 3 6
10 2 4 6 7 3 7
2 3 6 7 11 3 3
3 2 6 7 11 3 4
6 2 3 7 11 3 5
7 2 3 6 11 3 6
11 2 3 6 7 3 7
1 7 12 0 0 1 1
7 1 12 0 0 1 1
12 1 7 0 0 1 1
8 13 16 0 0 0 1001
13 8 16 0 0 0 1002
16 8 13 0 0 0 1002
9 14 16 0 0 0 1001
14 9 16 0 0 0 1002
16 9 14 0 0 0 1002
10 15 17 0 0 0 1001
15 10 17 0 0 0 1002
17 10 15 0 0 0 1002
11 13 16 0 0 0 1001
13 11 16 0 0 0 1002
16 11 13 0 0 0 1002
12 14 17 0 0 0 1001
14 12 17 0 0 0 1002
17 12 14 0 0 0 1002
13 0 0 0 0 0 1005
14 0 0 0 0 0 1005
15 0 0 0 0 0 1005
16 0 0 0 0 0 1005
17 0 0 0 0 0 1005

 

 

프로그램으로 구한 RHS

 

0.56539810563947124 0.63745474769401966 2 2
0.55110789267533677 0.0000000000000000 2 0
0.0000000000000000 0.26886136468257282 0 1
0.0000000000000000 0.35875882639560025 0 1
0.29469810563947119 0.0000000000000000 1 0
0.55110789267533677 0.63745474769401966 2 2
0.27070000000000000 0.62762019107817313 1 2
0.0000000000000000 0.30870967741935484 0 1
0.29469810563947119 0.32874507027466482 1 1
0.32041424498816867 0.35875882639560025 1 1
0.23069364768716805 0.26886136468257282 1 1
0.27070000000000000 0.0000000000000000 1 0
0.0000000000000000 0.0000000000000000 0 0
0.0000000000000000 0.0000000000000000 0 0
0.0000000000000000 0.0000000000000000 0 0
0.0000000000000000 0.0000000000000000 0 0
0.0000000000000000 0.0000000000000000 0 0

 

LHS의 정렬

 

sort -n -o sorted_lhs.dat lhs.dat

 

정렬된 sorted_lhs.dat

 

1 1 5 6 9 3 3
1 1 5 6 9 3 4
1 6 8 0 0 2 2
1 7 12 0 0 1 1
2 3 6 7 11 3 3
2 4 6 7 10 3 3
3 2 6 7 11 3 4
4 2 6 7 10 3 4
5 1 1 6 9 3 5
6 1 1 5 9 3 6
6 1 8 0 0 2 2
6 2 3 7 11 3 5
6 2 4 7 10 3 5
7 1 12 0 0 1 1
7 2 3 6 11 3 6
7 2 4 6 10 3 6
8 1 6 0 0 2 2
8 13 16 0 0 0 1001
9 1 1 5 6 3 7
9 14 16 0 0 0 1001
10 2 4 6 7 3 7
10 15 17 0 0 0 1001
11 2 3 6 7 3 7
11 13 16 0 0 0 1001
12 1 7 0 0 1 1
12 14 17 0 0 0 1001
13 0 0 0 0 0 1005
13 8 16 0 0 0 1002
13 11 16 0 0 0 1002
14 0 0 0 0 0 1005
14 9 16 0 0 0 1002
14 12 17 0 0 0 1002
15 0 0 0 0 0 1005
15 10 17 0 0 0 1002
16 0 0 0 0 0 1005
16 8 13 0 0 0 1002
16 9 14 0 0 0 1002
16 11 13 0 0 0 1002
17 0 0 0 0 0 1005
17 10 15 0 0 0 1002
17 12 14 0 0 0 1002

 

정렬된 sorted_lhs.dat와 rhs.dat를 이용하여 해를 구하는 프로그램

 

same model, different level, multiple trait animal model 1 - slove

 

PROGRAM mt_dl_animal_model1_solve

! program name - multiple trait, different level, animal_model1_solve
! programmer - Park Byoungho
! usage - mt_dl_animal_model1_solve
! purpose : read sorted_lhs.dat and rhs.dat which are made in mt_animal_model_setup.exe and are sorted
! find a solution using Gauss-Seidel iteration
! Date - 2009.7.10.
! update - none

USE gi

IMPLICIT NONE

! data dictionary
INTEGER, PARAMETER :: no_of_trait = 2 ! number of trait
INTEGER, PARAMETER :: no_of_fixed = 2 ! number of fixed effects
INTEGER, PARAMETER :: no_of_eq = 17 ! 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 :: pre_eq_no ! previous equation number
INTEGER :: cur_eq_no ! current equation number

INTEGER, ALLOCATABLE :: loc_of_off(:) ! location of off-diagonal in the equation
INTEGER :: class_code ! to identify fixed and random : 1-fixed, 2-random
INTEGER, ALLOCATABLE :: level_of_effects(:) ! level of each fixed and animal effect
REAL(KIND=8), ALLOCATABLE :: diagonal(:,:) ! diagonal element of lhs
REAL(KIND=8), ALLOCATABLE :: i_diagonal(:,:) ! inverse of diagonal element of lhs

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(:,:) ! genetic variance-covariance matrix
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 INTEGER, ALLOCATABLE :: trait_combi(:) ! trait combination
INTEGER :: no_of_ud ! number of upper digonal
INTEGER :: no_of_tc ! number of trait combination
INTEGER :: tc_no ! trait combination number

INTEGER, ALLOCATABLE :: lhs(:,:) ! left-hand side
INTEGER :: lhs2, lhs3, lhs4, lhs5
REAL(KIND = 8), ALLOCATABLE :: rhs(:,:)! left-hand side and right-hand side, x'y
REAL(KIND = 8), ALLOCATABLE :: temp_rhs(:) ! right-hand side one line
REAL(KIND = 8), ALLOCATABLE :: solutions(:,:) ! solutions
REAL(KIND = 8), ALLOCATABLE :: pre_sol(:) ! previous solution
INTEGER :: iteration ! iteration
REAL(KIND = 8) :: epsilon ! sum of square between old and new solutions
INTEGER :: i,j,k,l ! loop
INTEGER :: status ! i/o status
CHARACTER(LEN = 40) :: error_msg ! error message
INTEGER, PARAMETER :: MAX_ITER = 200 ! maximum number of iteration
REAL(KIND = 8), PARAMETER :: criteria = 1.E-8 ! criteria for stopping
INTEGER :: no_of_lhs ! number of lhs lines
INTEGER :: data_type ! data type : 1-fixed 2-6:animal

INTEGER, DIMENSION(99) :: iw ! for inverse
REAL(KIND = 8) :: z ! for inverse
INTEGER :: mr ! for inverse

ALLOCATE(gvcv_m(no_of_trait,no_of_trait))
ALLOCATE(diagonal(no_of_trait,no_of_trait))
ALLOCATE(i_diagonal(no_of_trait,no_of_trait))
ALLOCATE(temp_rhs(no_of_trait))
ALLOCATE(pre_sol(no_of_trait)) ! previous solution

ALLOCATE(rhs(no_of_eq ,no_of_trait * 2), solutions(no_of_eq ,no_of_trait))

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))
ALLOCATE(trait_combi(no_of_ud))

! initialiaze solutions to zero
solutions = 0

! get parameter file name
par_filename = 'vcvtrt.par'

! open file for data processing procedure
OPEN(UNIT = 31, FILE = 'ongoing_solve.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 ! 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 vcv =', gvcv_m

! 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)
WRITE(31,*) 'inverse of trait combination rvcv', temp_rvcv
rvcv_tc(k,:) = temp_rvcv
END DO

!open left-hand 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 + 3))

! 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 = *, IOSTAT = status) rhs(i,:)
IF (status /= 0) EXIT ! reach end of file
END DO
CLOSE(12) WRITE(31,*) 'right hand side, number of equations = ', no_of_eq
DO i = 1, no_of_eq
WRITE(31,*) rhs(i,:)
END DO

! iteration
DO iteration = 1, MAX_ITER

pre_eq_no = 1
diagonal = 0.0
i_diagonal = 0.0
temp_rhs = rhs(1,1:no_of_trait)
epsilon = 0.0 !write(31,*) 'temp_rhs = ', temp_rhs ! process each lhs line
DO i = 1, no_of_lhs

cur_eq_no = lhs(i,1) lhs2 = lhs(i,2)
lhs3 = lhs(i,3)
lhs4 = lhs(i,4)
lhs5 = lhs(i,5) tc_no = lhs(i, 6) ! trati combination number
data_type = lhs(i, 7)

! new equation, or calculate solution
IF (pre_eq_no /= cur_eq_no) THEN
pre_sol = solutions(pre_eq_no,:) ! store old solution
CALL geninv(diagonal, i_diagonal)
solutions(pre_eq_no,:) = MATMUL(i_diagonal, temp_rhs) ! get a new solution
epsilon = epsilon + SUM((pre_sol - solutions(pre_eq_no,:)) * (pre_sol - solutions(pre_eq_no,:))) ! calculate sum of squares of difference between old solution and new solution

 

!write(31,*) 'diagonal = ', diagonal
!write(31,*) 'i_diagonal = ', i_diagonal
!write(31,*) 'new_sol = ', MATMUL(i_diagonal, temp_rhs)
!write(31,*) 'equation = ', cur_eq_no !WRITE(31,*) 'iteration ', iteration , 'solutions, equation = ', pre_eq_no
!DO j = 1, no_of_eq
! WRITE(UNIT = 31, FMT = *) j, solutions(j,:)
!END DO diagonal = 0.0
i_diagonal = 0.0
temp_rhs = rhs(cur_eq_no,1:no_of_trait) !write(31,*) 'temp_rhs = ', temp_rhs pre_eq_no = cur_eq_no
END IF

! process off-diagnoal
SELECT CASE(data_type)
CASE(1) ! fixed effect
diagonal(1,1) = diagonal(1,1) + rvcv_tc(tc_no,1)
temp_rhs(1) = temp_rhs(1) - rvcv_tc(tc_no,1)*solutions(lhs2,1) &
- rvcv_tc(tc_no,1)*solutions(lhs3,1) !write(31,*) 'temp_rhs = ', temp_rhs CASE(2) ! fixed effect
diagonal(2,2) = diagonal(2,2) + rvcv_tc(tc_no,3)
temp_rhs(2) = temp_rhs(2) - rvcv_tc(tc_no,3)*solutions(lhs2,2) &
- rvcv_tc(tc_no,3)*solutions(lhs3,2) !write(31,*) 'temp_rhs = ', temp_rhs

 

CASE(3) ! fixed effect
diagonal(1,1) = diagonal(1,1) + rvcv_tc(tc_no,1)
temp_rhs(1) = temp_rhs(1) - rvcv_tc(tc_no,2)*solutions(lhs2,2) &
- rvcv_tc(tc_no,1)*solutions(lhs3,1) - rvcv_tc(tc_no,2)*solutions(lhs4,2) &
- rvcv_tc(tc_no,1)*solutions(lhs5,1) - rvcv_tc(tc_no,2)*solutions(lhs5,2) !write(31,*) 'temp_rhs = ', temp_rhs

 

CASE(4) ! fixed effect
diagonal(2,2) = diagonal(2,2) + rvcv_tc(tc_no,3)
temp_rhs(2) = temp_rhs(2) - rvcv_tc(tc_no,2)*solutions(lhs2,1) &
- rvcv_tc(tc_no,2)*solutions(lhs3,1) - rvcv_tc(tc_no,3)*solutions(lhs4,2) &
- rvcv_tc(tc_no,2)*solutions(lhs5,1) - rvcv_tc(tc_no,3)*solutions(lhs5,2) !write(31,*) 'temp_rhs = ', temp_rhs

 

CASE(5) ! fixed effect
diagonal(1,1) = diagonal(1,1) + rvcv_tc(tc_no,1)
temp_rhs(1) = temp_rhs(1) - rvcv_tc(tc_no,1)*solutions(lhs2,1) - rvcv_tc(tc_no,2)*solutions(lhs3,2) &
- rvcv_tc(tc_no,2)*solutions(lhs4,2) &
- rvcv_tc(tc_no,1)*solutions(lhs5,1) - rvcv_tc(tc_no,2)*solutions(lhs5,2) !write(31,*) 'temp_rhs = ', temp_rhs

 

CASE(6) ! fixed effect
diagonal(2,2) = diagonal(2,2) + rvcv_tc(tc_no,3)
temp_rhs(2) = temp_rhs(2) - rvcv_tc(tc_no,2)*solutions(lhs2,1) - rvcv_tc(tc_no,3)*solutions(lhs3,2) &
- rvcv_tc(tc_no,2)*solutions(lhs4,1) &
- rvcv_tc(tc_no,2)*solutions(lhs5,1) - rvcv_tc(tc_no,3)*solutions(lhs5,2) !write(31,*) 'temp_rhs = ', temp_rhs

 

CASE(7) ! fixed effect
diagonal(1,1) = diagonal(1,1) + rvcv_tc(tc_no,1)
diagonal(1,2) = diagonal(1,2) + rvcv_tc(tc_no,2)
diagonal(2,1) = diagonal(2,1) + rvcv_tc(tc_no,2)
diagonal(2,2) = diagonal(2,2) + rvcv_tc(tc_no,3)
temp_rhs(1) = temp_rhs(1) - rvcv_tc(tc_no,1)*solutions(lhs2,1) - rvcv_tc(tc_no,2)*solutions(lhs3,2) &
- rvcv_tc(tc_no,1)*solutions(lhs4,1) - rvcv_tc(tc_no,2)*solutions(lhs5,2)
temp_rhs(2) = temp_rhs(2) - rvcv_tc(tc_no,2)*solutions(lhs2,1) - rvcv_tc(tc_no,3)*solutions(lhs3,2) &
- rvcv_tc(tc_no,2)*solutions(lhs4,1) - rvcv_tc(tc_no,3)*solutions(lhs5,2) !write(31,*) 'temp_rhs = ', temp_rhs
CASE(1001) ! animal effect, both parents are known, animal is diagonal
temp_rhs = temp_rhs - (-0.5) * 2.0 * MATMUL(gvcv_m, solutions(lhs2,:)) &
- (-0.5) * 2.0 * MATMUL(gvcv_m, solutions(lhs3,:))
diagonal = diagonal + 2.0 * gvcv_m !write(31,*) 'temp_rhs = ', temp_rhs CASE(1002) ! animal effect, both parents are known, sire or dam is diagonal
temp_rhs = temp_rhs - (-0.5) * 2.0 * MATMUL(gvcv_m, solutions(lhs2,:)) &
- 0.25 * 2.0 * MATMUL(gvcv_m, solutions(lhs3,:))
diagonal = diagonal + 0.25 * 2.0 * gvcv_m !write(31,*) 'temp_rhs = ', temp_rhs CASE(1003) ! animal effect, one parent is known, animal is diagonal
temp_rhs = temp_rhs - (-0.5) * (4.0/3.0) * MATMUL(gvcv_m, solutions(lhs2,:))
diagonal = diagonal + (4.0/3.0) * gvcv_m !write(31,*) 'temp_rhs = ', temp_rhs CASE(1004) ! animal effect, one parents is known, parent is diagonal
temp_rhs = temp_rhs - (-0.5) * (4.0/3.0) * MATMUL(gvcv_m, solutions(lhs2,:))
diagonal = diagonal + 0.25 * (4.0/3.0) * gvcv_m !write(31,*) 'temp_rhs = ', temp_rhs CASE(1005) ! animal effect, no parents is known
diagonal = diagonal + 1.0 * gvcv_m !write(31,*) 'temp_rhs = ', temp_rhs END SELECT
END DO ! one iteration ends - end of array

! calculate last solution
pre_sol = solutions(pre_eq_no,:) ! store old solution
CALL geninv(diagonal, i_diagonal)
solutions(pre_eq_no,:) = MATMUL(i_diagonal, temp_rhs) ! get a new solution
epsilon = epsilon + SUM((pre_sol - solutions(pre_eq_no,:)) * (pre_sol - solutions(pre_eq_no,:))) ! calculate sum of squares of difference between old solution and new solution

! writing temporary solutions
WRITE(31,*) 'iteration ', iteration , 'solutions, epsilon = ', epsilon/(no_of_eq * no_of_trait)
DO i = 1, no_of_eq
WRITE(UNIT = 31, FMT = *) i, solutions(i,:)
END DO ! write iteration number and epsilon
epsilon = epsilon / (no_of_eq * no_of_trait)
WRITE(*,*) 'iteration = ', iteration , ', epsilon = ', epsilon
IF (epsilon < criteria) THEN
EXIT
END IF

END DO ! iteration

! open file for writing solutions
OPEN(UNIT = 13, FILE = 'sol.dat', STATUS = 'REPLACE', ACTION = 'WRITE', IOSTAT = status, IOMSG = error_msg)
DO i = 1, no_of_eq
WRITE(UNIT = 13, FMT = *) i, solutions(i,:)
END DO
CLOSE(13) CLOSE(31)

END PROGRAM mt_dl_animal_model1_solve
위 프로그램을 mt_dl_animal_model1_solve.f95로 저장

 

컴파일 : gfortran inverse.f95 mt_dl_animal_model1_solve.f95 -o mt_dl_animal_model1_solve

 

위 프로그램을 이용하여 구한 해

 

1 8454.2180156649902 9927.6787965349304
2 8261.1165896839902 0.0000000000000000
3 0.0000000000000000 8800.0105249646367
4 0.0000000000000000 11287.159930923246
5 786.53640606703379 0.0000000000000000
6 372.58461329470367 177.24048154387867
7 -418.84716617122098 14.282410176315935
8 -183.66064678099795 -150.43326040524298
9 -36.017708571063032 40.594451419793934
10 292.90505035639950 146.45254115918314
11 -305.30757371828741 -172.18797068288069
12 85.629109864515328 62.348959719157151
13 -171.25801190825430 -124.69771265006048
14 24.805603893806175 51.471578446150843
15 146.45250567278080 73.226256012354938
16 -146.45228764359885 -73.225942496537854
17 146.45255271881069 73.226291824543708

 

컴파일 및 실행화면

 

 

파일

 

 

1247191813_dairy1.dat

다운로드

1247191813_dairy1.dat

1247191813_pedi1.dat

다운로드

1247191813_pedi1.dat

1247191813_vcvtrt.par

다운로드

1247191813_vcvtrt.par

1247191813_mt_dl_animal_model1_setup.f95

다운로드

1247191813_mt_dl_animal_model1_setup.f95

1247191813_mt_dl_animal_model1_solve.f95

다운로드

1247191813_mt_dl_animal_model1_solve.f95

1247191813_inverse.f95

다운로드

1247191813_inverse.f95

1247191813_lhs.dat

다운로드

1247191813_lhs.dat

1247191813_rhs.dat

다운로드

1247191813_rhs.dat

1247191813_sorted_lhs.dat

다운로드

1247191813_sorted_lhs.dat

1247191813_sol.dat

다운로드

1247191813_sol.dat

 

위 모든 파일을 한방에

1247191813_mt_dl_animal_model1.zip

다운로드

1247191813_mt_dl_animal_model1.zip

 

1247191186_dairy1.dat

다운로드