multiple trait sire model - Gauss seidel solver

multiple trait sire model - Gauss seidel solver

 

지난번에 lhs와 rhs를 setup하였다.

 

이번에는 lhs와 rhs를 이용하여 solution을 구하여 본다.

 

lhs정렬 sort 명령어

 

sort -n -o sorted_lhs.dat lhs.dat

 

정렬한 sorted_lhs.dat

 

1 4 7 7 1
1 5 7 7 1
1 6 8 7 1
2 4 7 7 1
2 6 8 7 1
2 6 8 7 1
3 4 7 7 1
3 5 7 7 1
3 5 8 3 1
3 6 8 3 1
4 1 7 7 1
4 2 7 7 1
4 3 7 7 1
5 1 7 7 1
5 3 7 7 1
5 3 8 3 1
6 1 8 7 1
6 2 8 7 1
6 2 8 7 1
6 3 8 3 1
7 0 0 0 2
7 1 4 7 1
7 1 5 7 1
7 2 4 7 1
7 3 4 7 1
7 3 5 7 1
8 0 0 0 2
8 1 6 7 1
8 2 6 7 1
8 2 6 7 1
8 3 5 3 1
8 3 6 3 1
읽어서 저장할 rhs.dat

 

1.9571227 1.7957478 0.23437987
2.3159108 2.0812473 0.23651524
2.5084927 3.2069154 0.17655139
2.4685802 2.1218426 0.25010207
1.1732187 2.1414437 0.17449327
3.1397271 2.8206241 0.22285117
3.0944958 3.3640823 0.42459533
3.6870303 3.7198284 0.22285117
genetic and residual variace-covariance uppder diagonal과 trati combination있는 파일

 

1.25 1.75 2.5 7.5 62.5 500
15 12 90 85 200 4500
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
위 파일을 mt_sire_model_par.dat 로 저장

 

 

 

프로그램 소스

 

PROGRAM mt_sire_model_solve

! program name - multiple trait sire_model_solve
! programmer - Park Byoungho
! usage - mt_sire_model_solve < mt_sire_model_solve.par
! purpose : read sorted_lhs.dat and rhs.dat which are made in sire_model_setup.exe and are sorted
! find a solution using Gauss-Seidel iteration
! Date - 2009.5.25
! update - none

USE gi

IMPLICIT NONE

! data dictionary
INTEGER :: no_of_effects ! number of effects
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
REAL :: value ! observation
INTEGER :: class_code ! to identify fixed and random : 1-fixed, 2-random
INTEGER, ALLOCATABLE :: level_of_effects(:) ! level of each fixed and sire effect
REAL, ALLOCATABLE :: diagonal(:,:) ! diagonal element of lhs
REAL, ALLOCATABLE :: i_diagonal(:,:) ! inverse of diagonal element of lhs
INTEGER :: no_of_eq ! number of equations
INTEGER :: loc_of_sire ! starting location of sire effect
INTEGER :: no_of_trait ! number of trait

CHARACTER(LEN = 256) :: vcv_filename ! variance-covariance parameter file name

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_m(:,:,:) ! 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

REAL, ALLOCATABLE :: lhs(:,:), rhs(:,:)! left-hand side and right-hand side, x'y
REAL, ALLOCATABLE :: temp_rhs(:) ! right-hand side one line
REAL, ALLOCATABLE :: solutions(:,:) ! solutions
REAL, ALLOCATABLE :: pre_sol(:) ! previous solution
INTEGER :: iteration ! iteration
REAL :: epsilon ! sum of square 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, PARAMETER :: criteria = 1.E-8 ! criteria for stopping
INTEGER :: no_of_lhs ! number of lhs lines
INTEGER :: data_type ! data type : 1-fixed 2-sire

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

! get the number of effectss
WRITE (*,*) "Enter the number of effects including fixed and sire effects: "
READ (*,*) no_of_effects

ALLOCATE(loc_of_off(no_of_effects - 1), level_of_effects(no_of_effects))

! get the number of equations
WRITE (*,*) "Enter the level of each fixed and sire effect : "
READ (*,*) level_of_effects

! get the number of traits
WRITE (*,*) "Enter the number of traits : "
READ (*,*) no_of_trait

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

no_of_eq = SUM(level_of_effects) ! total number of equations
loc_of_sire = SUM(level_of_effects(1:no_of_effects-1)) + 1 ! (sum of level of fixed effects) + 1
ALLOCATE(rhs(no_of_eq ,no_of_trait), 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_m(no_of_tc,no_of_trait,no_of_trait))
ALLOCATE(trait_combi(no_of_ud))

! initialiaze solutions to zero
solutions = 0

! get the genetic and residual variance-covariance file name
WRITE (*,*) "Enter the genetic and residual variance-covariance file name "
READ (*,*) vcv_filename

! open parameter file
OPEN(UNIT = 21, FILE = vcv_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

! genetic variance-covariance
READ(UNIT = 21, FMT = *) 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

! residual variance-covariance
READ(UNIT = 21, FMT = *) rvcv

! read trait combination and inverse rvcv
DO k = 1, no_of_tc
READ(UNIT = 21, FMT = *) trait_combi
temp_rvcv = rvcv * trait_combi
z = 0.0
CALL DJNVHF(temp_rvcv,no_of_trait,iw,z,mr)
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-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_effects + 2))

! 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)

!open right-hand side file
OPEN(UNIT = 11, 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 = 11,FMT = *, IOSTAT = status) rhs(i,:)
IF (status /= 0) EXIT ! reach end of file
END DO
close(11)
! iteration
DO iteration = 1, MAX_ITER

pre_eq_no = 1
diagonal = 0.0
temp_rhs = rhs(1,:)
epsilon = 0.0

! lhs
DO i = 1, no_of_lhs
cur_eq_no = lhs(i,1)
loc_of_off = lhs(i,2:no_of_effects) ! location of off diagonal
tc_no = lhs(i,no_of_effects + 1) ! trati combination number
data_type = lhs(i,no_of_effects + 2)

! nrew 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

diagonal = 0.0
temp_rhs = rhs(cur_eq_no,:)
END IF

! process off-diagnoal
SELECT CASE(data_type)
CASE(1) ! fixed effect
temp_rhs = temp_rhs - MATMUL(rvcv_tc_m(tc_no,:,:), SUM(solutions(loc_of_off,:),1))
diagonal = diagonal + rvcv_tc_m(tc_no,:,:)
pre_eq_no = cur_eq_no
CASE(2) ! sire effect
diagonal = diagonal + gvcv_m
pre_eq_no = cur_eq_no
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

! 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 = *) solutions(i,:)
END DO
CLOSE(13)

END PROGRAM mt_sire_model_solve

 

 

위 프로그램을 mt_sire_model_solve.f95로 저장

 

 

inverse 하고, upper diagonal 을 matrix로 바꾸는 function

 

MODULE gi CONTAINS

SUBROUTINE geninv(x,gix) REAL, DIMENSION(:,:), INTENT(IN) :: x ! square matrix to inverse
REAL, 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 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로 저장

 

위 프로그램을 이용하여 나온 solution

 

24.048454 74.278763 527.86261
25.404404 82.552605 567.17712
24.160355 85.626450 572.11475
3.7658474 5.9781690 36.903538
-4.4663043 -4.8003001 -2.8413122
0.52534527 -0.88338321 -25.695047
2.97499169E-02 -0.13078016 -1.9550242
-2.97498535E-02 0.13078132 1.9550378

 

컴파일 및 프로그램 실행 화면

 

 

1243315321_mt_sire_model_solve.f95

다운로드

1243315321_mt_sire_model_solve.f95

1243315321_inverse.f95

다운로드

1243315321_inverse.f95

1243315321_mt_sire_model_solve.par

다운로드

1243315321_mt_sire_model_solve.par