5 형질 다형질 개체 모형 - same model but different level

5 형질 다형질 개체 모형 - same model but different level

 

지난 프로그램과 달라진 것은 입력 데이타의 구조 변경으로 인하여 데이터 읽는 부분만 바뀌었습니다.

 

자료(일부)

 

1302 12 20 10114 0 0 0 0 0 0 0 0 0 0 0 0
1296 11 20 8981 8 24 10418 5 28 9684 0 0 0 0 0 0
1291 10 21 7601 13 27 9042 10 30 9560 0 0 0 0 0 0
1288 11 22 7264 10 26 5098 6 30 7878 0 0 0 0 0 0
1287 8 21 8200 5 25 8181 2 28 8319 7 35 8853 0 0 0
1298 13 21 6693 15 26 5369 0 0 0 0 0 0 0 0 0
1332 9 17 8586 8 21 10698 9 26 8899 0 0 0 0 0 0
1327 14 19 5046 0 0 0 0 0 0 0 0 0 0 0 0
1325 6 17 8495 4 21 6888 2 25 8646 0 0 0 0 0 0
1268 9 24 8696 6 28 5308 4 32 8593 0 0 0 0 0 0

설명 : 개체, 1st 형질의 fixed1, fixed2, obs, 2nd 형질의fixed1, fixed2, obs, 5형질까지 반복

 

위 자료를 dairy_5tr.dat로 저장

 

 

혈통 자료(일부)

 

191 587 301
192 587 318
193 594 326
194 594 401
195 587 330
196 587 323
197 690 317
198 690 345
199 571 311
200 571 406

설명 : animal, sire, dam

 

위 자료를 pedi_5tr.dat로 저장

 

파라미터 자료

 

23000.0 14921.0 13885.0 12434.8 11213.9 32000.0 17007.9 14966.6 14403.0 31000.0 15614.8 14465.5 28000.0 14297.6 27000.0
12385.0 10947.8 10137.1 8924.0 8657.3 13715.0 12318.4 9617.2 9221.4 11759.0 9428.9 8641.4 9334.0 7790.7 9000.0
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
1 1 0 0 0 1 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 1 0 0 0 0 0
1 0 1 0 0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 1 1 0 0 1 0 0 0 0 0
1 1 1 0 0 1 1 0 0 1 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 1 0 0
1 0 0 1 0 0 0 0 0 0 0 0 1 0 0
0 0 0 0 0 1 0 1 0 0 0 0 1 0 0
1 1 0 1 0 1 0 1 0 0 0 0 1 0 0
0 0 0 0 0 0 0 0 0 1 1 0 1 0 0
1 0 1 1 0 0 0 0 0 1 1 0 1 0 0
0 0 0 0 0 1 1 1 0 1 1 0 1 0 0
1 1 1 1 0 1 1 1 0 1 1 0 1 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
1 0 0 0 1 0 0 0 0 0 0 0 0 0 1
0 0 0 0 0 1 0 0 1 0 0 0 0 0 1
1 1 0 0 1 1 0 0 1 0 0 0 0 0 1
0 0 0 0 0 0 0 0 0 1 0 1 0 0 1
1 0 1 0 1 0 0 0 0 1 0 1 0 0 1
0 0 0 0 0 1 1 0 1 1 0 1 0 0 1
1 1 1 0 1 1 1 0 1 1 0 1 0 0 1
0 0 0 0 0 0 0 0 0 0 0 0 1 1 1
1 0 0 1 1 0 0 0 0 0 0 0 1 1 1
0 0 0 0 0 1 0 1 1 0 0 0 1 1 1
1 1 0 1 1 1 0 1 1 0 0 0 1 1 1
0 0 0 0 0 0 0 0 0 1 1 1 1 1 1
1 0 1 1 1 0 0 0 0 1 1 1 1 1 1
0 0 0 0 0 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

1행 : residual variance-covariance

2행 : genetic variance-covariance

나머지 : trait combination

 

위 자료를 vcvtrt_5tr.par로 저장

 

 

LHS와RHS를 만드는 프로그램

 

PROGRAM mt_dl_animal_model3_5tr_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_5tr_setup
! purpose : read data file and write left-hand side and right-hand side
! Date - 2009.10.9.
! update
! 2009.10.23 - to change the data reading part because of the change of input data structure
! old structure : fixed, ..., fixed, cow, obs, ... obs
! new structrue : cow, 1st trait(fixed, .., fixed, obs), ..., 5th trait(fixed, ..., fixed, obs)

USE gi

IMPLICIT NONE

! data dictionary
INTEGER, PARAMETER :: no_of_trait = 5 ! number of trait
INTEGER, PARAMETER :: no_of_fixed = 2 ! number of fixed effects
INTEGER, PARAMETER :: no_of_eq = 1490 ! 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

REAL(KIND = 8), ALLOCATABLE :: temp_data(:) ! temporary animal and fixed effect
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_5tr.dat'
! get pedigree file name
pedi_filename = 'pedi_5tr.dat'
! get parameter file name
par_filename = 'vcvtrt_5tr.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_data(1 + (no_of_fixed + 1) * no_of_trait)) ! 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_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) ! 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_data
IF (status /= 0) EXIT ! end of file

! copy temp_data(animal, fixed and observation) to effects and observations
DO i = 1, no_of_trait
effects(i,1:no_of_fixed) = temp_data(2 + (i - 1) * (no_of_fixed + 1) : i * (no_of_fixed + 1))
effects(i,no_of_fixed + 1) = temp_data(1)
observations(i) = temp_data(1 + i * (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_5tr_setup

위 소스를 mt_dl_animal_model3_5tr_setup.f95로 저장

 

컴파일 방법 : gfortran inverse.f95 mt_dl_animal_model3_5tr_setup.f95 -o mt_dl_animal_model3_5tr_setup.exe

 

실행방법 : mt_dl_animal_model3_5tr_setup

 

 

정렬

 

정렬 방법 : sort -n -o sorted_lhs.dat lhs.dat

 

 

해를 구하는 프로그램

 

PROGRAM mt_dl_animal_model3_5tr_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 = 5 ! number of trait
INTEGER, PARAMETER :: no_of_fixed = 2 ! number of fixed effects
INTEGER, PARAMETER :: no_of_eq = 1490 ! 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 = 1000 ! 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_5tr.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_5tr_solve

 

위 소스를 mt_dl_animal_model3_5tr_solve.f95로 저장

 

컴파일 방법 : gfortran inverse.f95 mt_dl_animal_model3_5tr_solve.f95 -o mt_dl_animal_model3_5tr_solve.exe

 

실행방법 : mt_dl_animal_model3_5tr_solve

 

 

컴파일하고 실행하는 화면

 

 

관련 파일

1256295208_5tr_mtam.zip

다운로드

1256295208_5tr_mtam.zip