Sire Model - 역행렬을 이용하여 해를 구하기

 

Sire Model - 역행렬을 이용하여 해를 구하기

 

자료

 

1 1 1 28
1 2 1 19
1 3 2 25
2 1 1 30
2 3 2 21
2 3 2 30
3 2 1 21
3 1 1 27
3 3 2 25
3 2 2 19

 

위 자료를 sire_model.dat로 저장

 

자료 설명:1st column : 생년2nd column : age of dam3rd column : sire(random)4th column : observation

 

 

프로그램 소스

 

PROGRAM sire_model
! program name - sire_model
! programmer - Park Byoungho
! usage - sire_model < parameter_file
! purpose : read data file, construct left-hand side and right-hand side and find solution
! Date - 2009.4.15.
! update - none USE gi

IMPLICIT NONE

! data dictionary
CHARACTER(LEN = 40) :: filename ! data file name
INTEGER :: No_Of_Effects ! number of effects(including sire effect)
INTEGER, ALLOCATABLE :: no_of_levels(:) ! number of levels of each fixed and sire effect
INTEGER :: first_equation_of_sire ! first equaton number of sire effect

REAL :: var_ratio ! variance ratio
INTEGER, ALLOCATABLE :: data_line(:) ! fixed and sire effects for each line
REAL :: y ! observation for each line REAL, ALLOCATABLE :: lhs(:,:), gi_lhs(:,:), rhs(:), solution(:) ! left-hand side, right-hand side and solutions

INTEGER :: No_Of_Equations ! number of equation
INTEGER :: i, j ! do loop

INTEGER :: status ! I/O status
CHARACTER(LEN = 40) :: error_msg ! error message

! get data file name
WRITE (*,*) "Enter data file name : "
READ (*,*) filename ! get number of effects
WRITE (*,*) "Enter the number of effects(including fixed and sire effects) : "
READ (*,*) No_Of_Effects ! allocate dimensions of no_of_levels and fixed_data
ALLOCATE(no_of_levels(no_of_effects), data_line(no_of_effects))

WRITE (*,*) "Enter the number of levels for each fixed and sire effect : "
READ (*,*) no_of_levels

! get the variance ratio
WRITE (*,*) "Enter the variance ratio : "
READ (*,*) var_ratio ! calculate first equation number of sire and total number of equation
first_equation_of_sire = SUM(no_of_levels(1:no_of_effects-1)) + 1
no_of_equations = SUM(no_of_levels) ALLOCATE(lhs(no_of_equations,no_of_equations), &
gi_lhs(no_of_equations,no_of_equations), &
rhs(no_of_equations), &
solution(no_of_equations) &
) ! initialize
lhs = 0
rhs = 0
gi_lhs = 0
solution = 0

! open data file
OPEN(UNIT = 11, FILE = filename, STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg)

IF (status /= 0) THEN ! file open failed
WRITE (*,'(1X, A, A)') 'File open failed -- error message : ', error_msg
STOP
END IF DO ! read each line
READ(UNIT = 11, FMT = *, IOSTAT = status) data_line, y
IF (status /= 0) EXIT ! end of file ! convert fixed effect into equation number
DO i = 2, no_of_effects
data_line(i) = data_line(i) + SUM(no_of_levels(1:i - 1))
END DO
!add 1 to appropriate location of lhs(ex : 9 elements for 2 fixed and 1 sire effect)
lhs(data_line,data_line) = lhs(data_line,data_line) + 1 ! add observation to each fixed and sire location of rhs
rhs(data_line) = rhs(data_line) + y END DO ! add variance ratio to sire effect diagonal part of lhs
DO i = first_equation_of_sire, no_of_equations
lhs(i, i) = lhs(i, i) + var_ratio
END DO ! constraint
lhs(1,:)=0
lhs(:,1)=0
rhs(1)=0 CALL geninv(lhs, gi_lhs)

solution = MATMUL(gi_lhs, rhs)

WRITE(*,*) 'left-hand side ='
WRITE(*,'(8(F3.0, 1X))') ((lhs(i,j), j = 1, no_of_equations), i = 1, no_of_equations)

WRITE(*,*) 'right-hand side ='
WRITE(*,'(F6.1)') (rhs(i),i = 1, no_of_equations) WRITE(*,*) 'ginverse of lhs ='
WRITE(*,'(8(F6.4, 1X))') ((gi_lhs(i,j), j = 1, no_of_equations), i = 1, no_of_equations)

WRITE(*,*) 'solutions ='
WRITE(*,*) (solution(i),i = 1, no_of_equations)

END PROGRAM sire_model

 

위 소스를 sire_model_1.f95로 저장

 

 

다음은 역행렬을 구하기 위한 소스

 

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

 

 

실행을 편리하게 하기 위하여

 

sire_model.dat
3
3 3 2
10
위 소스를 sire_model_1.par로 저장

 

 

프로그램 컴파일 및 실행