Single Trait Random Regression Model

 

Single Trait Random Regression Model

 

단형질 임의 회귀 모형

 

Data(R.A. Mrode, Linear Models for the Prediction of Animal Breeding Values, 2nd Edition, p138)

 

1 4 4 17.0
2 38 4 18.6
3 72 4 24.0
4 106 4 20.0
5 140 4 20.0
6 174 4 15.6
7 208 4 16.0
8 242 4 13.0
9 276 4 8.2
10 310 4 8.0
1 4 5 23.0
2 38 5 21.0
3 72 5 18.0
4 106 5 17.0
5 140 5 16.2
6 174 5 14.0
7 208 5 14.2
8 242 5 13.4
9 276 5 11.8
10 310 5 11.4
6 4 6 10.4
7 38 6 12.3
8 72 6 13.2
9 106 6 11.6
10 140 6 8.4
4 4 7 22.8
5 38 7 22.4
6 72 7 21.4
7 106 7 18.8
8 140 7 18.3
9 174 7 16.2
10 208 7 15.0
1 4 8 22.2
2 38 8 20.0
3 72 8 21.0
4 106 8 23.0
5 140 8 16.8
6 174 8 11.0
7 208 8 13.0
8 242 8 17.0
9 276 8 13.0
10 310 8 12.6
1st col : class variable - hert test day(HTD)2nd col : fixed regression - days in milk(DIM)3rd col : animal4th col : test dayfat yield(TDY)

 

위 자료를 strrm.dat로 저장

 

Pedigree

 

1 0 0
2 0 0
3 0 0
412
532
615
734
817
1st col : animal2nd col : sire3rd col : dam

 

위 혈통을 strrm.ped로 저장

 

LHS, RHS Preparation Program

 

PROGRAM strrm_setup ! program name
! : 'single trait random regression model' setup program
! prgrammer
! : Park Byoungho
! date
! : 2010.2.11.

IMPLICIT NONE

! data dictionary INTEGER, PARAMETER :: no_of_cv = 1 ! number of class variables
INTEGER, PARAMETER :: no_of_rc = 5 ! number of regression coefficient
INTEGER, PARAMETER :: no_of_effect = no_of_cv + no_of_rc + 2 ! number of effects(2 = animal effect, permanent environment effect)
INTEGER, PARAMETER :: order_of_ap = 3 ! order of animal polynomial

CHARACTER(LEN=256), PARAMETER :: data_filename = 'strrm.dat' ! data file name
CHARACTER(LEN=256), PARAMETER :: pedi_filename = 'strrm.ped' ! pedigree file name ! Legendre polynomial
REAL(KIND = 8), PARAMETER :: lambda(5,5) = &
RESHAPE((/SQRT(.5) , 0.0 , 0.0 , 0.0 , 0.0, &
0.0 , SQRT(1.5) , 0.0 , 0.0 , 0.0, &
-SQRT(2.5)*.5 , 0.0 , SQRT(2.5)*1.5 , 0.0 , 0.0, &
0.0 , -SQRT(3.5)*1.5, 0.0 , SQRT(3.5)*2.5, 0.0, &
SQRT(4.5)*3./8., 0.0 , -SQRT(4.5)*30./8., 0.0 , SQRT(4.5)*35./8./),(/5,5/)) !DO i = 1, 5
! WRITE(*,*) lambda(i,:)
!END DO

INTEGER, ALLOCATABLE :: input_ind(:) ! storage for independent variable when reading data
REAL(KIND = 8), ALLOCATABLE :: ind(:) ! storage for independent variable after making Legendre polynimal
REAL(KIND = 8) :: dep ! storage for dependent variable when reading data REAL(KIND = 8) :: dim_lp(500,5) ! Legendre polynomial of days in milk
INTEGER :: levels_of_cv ! levels of calss variable
INTEGER :: dim_max, dim_min, dim_range ! maximum, minimum and range of days in milk
REAL(KIND = 8) :: sut ! standardized unit of time
REAL(KIND = 8) :: sut_poly(5) ! polynomial of standardized unit of time
INTEGER :: no_of_animal ! number of animals
INTEGER :: no_of_eq ! number of equations REAL(KIND = 8), ALLOCATABLE :: rhs(:,:) ! right-hand side
INTEGER, ALLOCATABLE :: eq_no(:) ! storage for equation number REAL(KIND = 8) :: temp_ind ! temporary independent variable
INTEGER :: temp_eq_no INTEGER :: animal, sire, dam ! animal, sire, dam name for pedigree file INTEGER :: i,j ! for loop
INTEGER :: status ! I/O status
CHARACTER(LEN = 40) :: error_msg ! error message

ALLOCATE(input_ind(no_of_cv + 2)) ! independent variable
ALLOCATE(ind( no_of_cv + no_of_rc + 2)) ! independent variable after making Legendre polynimal
ALLOCATE(eq_no(no_of_cv + no_of_rc + 2)) ! storage for equation number ! open data file
OPEN(UNIT = 10, 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 ! look for the maximum and minimum of days in milk
! and levels of class variable(herd-test day)
! read first line of data file
READ(UNIT = 10, FMT = *, IOSTAT = status) input_ind, dep

levels_of_cv = input_ind(1)
dim_max = input_ind(2)
dim_min = input_ind(2)
no_of_animal = input_ind(3)

REWIND(10) ! read each line of data file
DO
READ(UNIT = 10, FMT = *, IOSTAT = status) input_ind, dep
IF (status /= 0) EXIT ! end of file IF (input_ind(1) > levels_of_cv) THEN
levels_of_cv = input_ind(1)
END IF IF (input_ind(2) < dim_min) THEN
dim_min = input_ind(2)
END IF IF (input_ind(2) > dim_max) THEN
dim_max = input_ind(2)
END IF IF (input_ind(3) > no_of_animal) THEN
no_of_animal = input_ind(3)
END IF END DO REWIND(10) no_of_eq = levels_of_cv + no_of_rc + no_of_animal * 2
ALLOCATE(rhs(no_of_eq, order_of_ap)) ! right-hand side

! initialize
rhs = 0.0 ! calculate the range of days in milk
dim_range = dim_max - dim_min

WRITE(*,*) levels_of_cv , ' = Levels of class variable '
WRITE(*,*) dim_min , ' = Minimum value of days in milk'
WRITE(*,*) dim_max , ' = Maximum value of days in milk'
WRITE(*,*) dim_range , ' = Range of days in milk '
WRITE(*,*) no_of_animal , ' = Number of animals ' ! open file for fixed_regression_setup result
OPEN(UNIT = 15, FILE = 'strrm_setup.rst', STATUS = 'REPLACE', ACTION = 'WRITE') ! write to fixed_regression_setup.rst
WRITE(15,*) levels_of_cv , ' = Levels of class variable '
WRITE(15,*) dim_min , ' = Minimum value of days in milk'
WRITE(15,*) dim_max , ' = Maximum value of days in milk'
WRITE(15,*) dim_range , ' = Range of days in milk '
WRITE(15,*) no_of_animal , ' = Number of animals ' CLOSE(15) ! close

!initialize
dim_lp = 0.0

! prepare Legendre polynomial of standardized unit of time(days in milk)
DO i = dim_min, dim_max
sut = -1. + 2. * (i - dim_min) / dim_range ! standardized unit of time
sut_poly(1) = 1.0 ! polynomial of standardized unit of time
sut_poly(2) = sut
sut_poly(3) = sut_poly(2) * sut
sut_poly(4) = sut_poly(3) * sut
sut_poly(5) = sut_poly(4) * sut
dim_lp(i,:) = MATMUL(sut_poly, lambda)
END DO


! left hand side
! open data file for writing left-hand side
OPEN(UNIT = 20, FILE = 'lhs.dat', STATUS = 'REPLACE', ACTION = 'WRITE')
! read each line
DO
READ(UNIT = 10, FMT = *, IOSTAT = status) input_ind, dep
IF (status /= 0) EXIT ! end of file ! assign equation number : class variable
eq_no(1) = input_ind(1) ! assign equation number : regression coefficient
DO i = 2, no_of_cv + no_of_rc
eq_no(i) = levels_of_cv + i - 1
END DO

! assign equation number : animal effect, permanent environmental effect
eq_no(no_of_cv + no_of_rc + 1) = levels_of_cv + no_of_rc + input_ind(3)
eq_no(no_of_cv + no_of_rc + 2) = levels_of_cv + no_of_rc + no_of_animal + input_ind(3) ! assign value : class variable
ind(1) = 1 ! Legendre polynomial of standardized unit of time(days in milk)
ind(2:6) = dim_lp(input_ind(2),:)

! assign value : animal effect, permanent environmental effect
ind(7:8) = 1.0

! for fixed_effect ( X'X, Q'Q, Z'Z)
DO i = 1, no_of_cv + no_of_rc + 2
! swap equation number
temp_eq_no = eq_no(i)
eq_no(i) = eq_no(1)
eq_no(1) = temp_eq_no

! swap independent variable
temp_ind = ind(i)
ind(i) = ind(1)
ind(1) = temp_ind

IF (eq_no(1) <= levels_of_cv + no_of_rc) THEN
WRITE(20,*) 1, input_ind(2), (eq_no(j),ind(1)*ind(j), j=1, no_of_cv + no_of_rc + 2)

rhs(eq_no(1),1) = rhs(eq_no(1),1) + ind(1) * dep
ELSE
WRITE(20,*) 2, input_ind(2), (eq_no(j),ind(1)*ind(j), j=1, no_of_cv + no_of_rc + 2)

rhs(eq_no(1),:) = rhs(eq_no(1),:) + dim_lp(INT(input_ind(2)),1:3) * dep
END IF

!rhs(eq_no(1)) = rhs(eq_no(1)) + ind(1) * dep
END DO END DO ! for processing variance ratio of permanent environmental effect
DO i = levels_of_cv + no_of_rc + no_of_animal + 1, no_of_eq
WRITE(20,*) 2001, 0, i, (0, j = 3, 1 + no_of_effect * 2)
END DO CLOSE(10) ! close data file

! open pedigree file
OPEN(UNIT = 30, 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 = 30, FMT = *, IOSTAT = status) animal, sire, dam
IF (status /= 0) EXIT ! end of file ! find location of equations
animal = animal + levels_of_cv + no_of_rc
IF (sire /= 0) sire = sire + levels_of_cv + no_of_rc ! if sire is 0, sire is unknown
IF (dam /= 0) dam = dam + levels_of_cv + no_of_rc ! if dam is 0, dam is unknown ! write animal effect
IF ( sire /= 0 .AND. dam /= 0) THEN ! both parents are known
WRITE(UNIT = 20, FMT = *) 1001, 0, animal, sire, dam, (0, i = 4, 2 + no_of_effect * 2)
WRITE(UNIT = 20, FMT = *) 1002, 0, sire, animal, dam, (0, i = 4, 2 + no_of_effect * 2)
WRITE(UNIT = 20, FMT = *) 1002, 0, dam, animal, sire, (0, i = 4, 2 + no_of_effect * 2)
ELSE IF (sire /= 0 .AND. dam == 0) THEN ! only sire is known
WRITE(UNIT = 20, FMT = *) 1003, 0, animal, sire, dam, (0, i = 4, 2 + no_of_effect * 2)
WRITE(UNIT = 20, FMT = *) 1004, 0, sire, animal, dam, (0, i = 4, 2 + no_of_effect * 2)
ELSE IF (sire == 0 .AND. dam /= 0) THEN ! only dam is known
WRITE(UNIT = 20, FMT = *) 1003, 0, animal, dam, sire, (0, i = 4, 2 + no_of_effect * 2)
WRITE(UNIT = 20, FMT = *) 1004, 0, dam, animal, sire, (0, i = 4, 2 + no_of_effect * 2)
ELSE ! no parents are known
WRITE(UNIT = 20, FMT = *) 1005, 0, animal, sire, dam, (0, i = 4, 2 + no_of_effect * 2)
END IF
END DO ! end of reading data CLOSE(30) ! close animal pedigree file CLOSE(20) ! close left hand side

! open file for right hand side
OPEN(UNIT = 40, FILE = 'rhs.dat', STATUS = 'REPLACE', ACTION = 'WRITE') ! write right-hand side
DO i = 1, no_of_eq
WRITE(40, *) rhs(i,:)
END DO CLOSE(40) ! close END PROGRAM strrm_setup

 

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

 

Soving Program

 

PROGRAM strrm_solve ! program name
! : iteration on data of single trait random regression model
! : solve program
! prgrammer
! : Park Byoungho
! date
! : 2010.2.11.

USE gi IMPLICIT NONE

! data dictionary INTEGER, PARAMETER :: no_of_cv = 1 ! number of class variables
INTEGER, PARAMETER :: no_of_rc = 5 ! number of regression coefficient
INTEGER, PARAMETER :: no_of_fe = no_of_cv + no_of_rc ! number of fixed effect: class variable + no of fixed regression coefficient
INTEGER, PARAMETER :: no_of_effect = no_of_cv + no_of_rc + 2 ! number of effects

REAL(KIND = 8), PARAMETER :: animal_vr(3,3) = & ! variance ratio of animal genetic variance and residual variance
RESHAPE((/ 3.297, 0.594, -1.381, &
0.594, 0.921, -0.289, &
-1.381, -0.289, 1.005/),(/3,3/))

REAL(KIND = 8), PARAMETER :: pe_vr(3,3) = & ! variance ratio of permanent environmental efeect variance and residual variance
RESHAPE((/ 6.872, -0.254, -1.101, &
-0.254, 3.171, 0.167, &
-1.101, 0.167, 2.457/),(/3,3/)) REAL(KIND = 8), PARAMETER :: residual_vr = 3.710 REAL(KIND = 8) :: animal_vr_i(3,3), pe_vr_i(3,3) ! (inverse of animal and permanent environment variance-covariance) * residual variance INTEGER, PARAMETER :: levels_of_cv = 10 ! levels of each calss variables
INTEGER, PARAMETER :: no_of_animal = 8
INTEGER, PARAMETER :: no_of_eq = levels_of_cv + no_of_rc + no_of_animal * 2 INTEGER :: no_of_lhs ! number of left hand side lines
REAL(KIND = 8), ALLOCATABLE :: lhs(:,:), rhs(:,:) ! left-hand side, right-hand side, x'y
REAL(KIND = 8), ALLOCATABLE :: solutions(:,:) ! solutions
INTEGER :: data_type ! fixed, animal, permanent environment
INTEGER :: dim ! days in milk INTEGER :: pre_eq_no ! previous equation number
INTEGER :: cur_eq_no ! current equation number REAL(KIND = 8) :: diag_ele(3,3) ! diagonal element
REAL(KIND = 8) :: diag_ele_i(3,3) ! diagonal element
REAL(KIND = 8) :: temp_rhs(3) ! temporary rihgt hand side
REAL(KIND = 8) :: pre_sol(3) ! previous solution ! Legendre polynomial
REAL(KIND = 8), PARAMETER :: lambda(3,3) = &
RESHAPE((/SQRT(.5) , 0.0 , 0.0 , &
0.0 , SQRT(1.5) , 0.0 , &
-SQRT(2.5)*.5 , 0.0 , SQRT(2.5)*1.5 /),(/3,3/)) INTEGER :: dim_min, dim_max, dim_range ! minimun, maximum and range of days in milk
REAL(KIND = 8) :: sut ! standardized unit of time
REAL(KIND = 8) :: sut_poly(3) ! polynomial of standardized unit of time REAL(KIND = 8) :: dim_lp(500,3) ! Legendre polynomial of days in milk
REAL(KIND = 8) :: dim_lp_m(500,3,3) ! outer product of Legendre polynomial of days in milk

INTEGER :: iteration ! iteration
REAL(KIND = 8) :: epsilon ! sum of squares of difference between old and new solutions
INTEGER :: i,j ! 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 ! calculate the inverse of animal and permanent environment variance-covariance
CALL geninv(animal_vr, animal_vr_i)
animal_vr_i = animal_vr_i * residual_vr

CALL geninv(pe_vr, pe_vr_i)
pe_vr_i = pe_vr_i * residual_vr

! allocate rhs, solution
ALLOCATE(rhs(no_of_eq, 3))
ALLOCATE(solutions(no_of_eq, 3))

! initialize
solutions = 0.0 !open left-hand side file
OPEN(UNIT = 10, 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 = 10, FMT = *, IOSTAT = status)
IF (status /= 0) EXIT ! reach end of file
no_of_lhs = no_of_lhs + 1
END DO ! decide array size of lhs
ALLOCATE(lhs(no_of_lhs, 2 + (no_of_cv + no_of_rc + 2) * 2))

!write(*,*) 'number of lines of left hand side = ', no_of_lhs ! store lhs to array
REWIND(10)
DO i = 1, no_of_lhs
READ(UNIT = 10, FMT = *, IOSTAT = status) lhs(i,:)
IF (status /= 0) EXIT ! reach end of file
END DO
CLOSE(10) !DO i = 1, no_of_lhs
! write(*,*) i, lhs(i,:)
!END DO !open right-hand side file
OPEN(UNIT = 20, 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 = 20,FMT = *) rhs(i,:)
END DO
close(20)! write(*,*) 'right hand side ...'
! DO i = 1, levels_of_cv + no_of_rc
! write(*,*) i, rhs(i)
! END DO ! prepare Legendre polynomial of standardized unit of time(days in milk)
! open setup result file
OPEN(UNIT = 25, FILE = 'strrm_setup.rst', STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg) IF (status /= 0) THEN ! file open failed
WRITE (*,'(1X, A, A)') 'Setup result file open failed -- error message : ', error_msg
STOP
END IF

READ(UNIT = 25, FMT ='(/ I12 / I12 / I12)', IOSTAT = status) dim_min, dim_max, dim_range

!initialize
dim_lp = 0.0
dim_lp_m = 0.0 ! prepare Legendre polynomial of standardized unit of time(days in milk)
DO i = dim_min, dim_max
sut = -1. + 2. * (i - dim_min) / dim_range ! standardized unit of time
sut_poly(1) = 1.0 ! polynomial of standardized unit of time
sut_poly(2) = sut
sut_poly(3) = sut_poly(2) * sut
dim_lp(i,:) = MATMUL(sut_poly, lambda)

DO j = 1, 3
dim_lp_m(i,j,:) = dim_lp(i,j) * dim_lp(i,:)
END DO
END DO ! solve MME
DO iteration = 1, MAX_ITER

pre_eq_no = 1
diag_ele = 0.0
temp_rhs = rhs(1,:)
epsilon = 0.0 ! start reading and processing each line of lhs
DO i = 1, no_of_lhs data_type = INT(lhs(i,1))
dim = INT(lhs(i,2))
cur_eq_no = INT(lhs(i,3)) IF(pre_eq_no /= cur_eq_no) THEN
pre_sol = solutions(pre_eq_no,:) ! store old solution
CALL geninv(diag_ele, diag_ele_i)
solutions(pre_eq_no,:) = MATMUL(diag_ele_i, 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
diag_ele = 0.0
temp_rhs = rhs(cur_eq_no,:)
END IF ! process off-diagnoal
SELECT CASE(data_type)
CASE(1) ! fixed effect
! adjust right-hand side
! temp_rhs - independend variable * solutions
DO j = 2, no_of_fe
temp_rhs(1) = temp_rhs(1) - lhs(i, 2*j+2) * solutions(INT(lhs(i,2*j+1)),1)
END DO temp_rhs(1) = temp_rhs(1) - lhs(i, 16) * SUM(dim_lp(dim,:) * solutions(INT(lhs(i,15)),:)) &
- lhs(i, 18) * SUM(dim_lp(dim,:) * solutions(INT(lhs(i,17)),:)) ! add diagonal element
diag_ele(1,1) = diag_ele(1,1) + lhs(i,4)
pre_eq_no = cur_eq_no CASE(2) ! fixed part of animal effect and permanent environmental effect
! adjust right-hand side
DO j = 2, (no_of_effect - 1)
temp_rhs = temp_rhs - dim_lp(dim,:) * lhs(i,2*j+2) * solutions(INT(lhs(i,2*j+1)),1)
END DO

temp_rhs = temp_rhs - MATMUL(dim_lp_m(dim,:,:) , solutions(INT(lhs(i,17)),:)) ! add diagonal element
diag_ele = diag_ele + dim_lp_m(dim,:,:)
pre_eq_no = cur_eq_no
CASE(1001) ! animal effect, both parents are known, animal is diagonal
temp_rhs = temp_rhs - (-0.5) * 2.0 * MATMUL(animal_vr_i , solutions(INT(lhs(i,4)),:))
temp_rhs = temp_rhs - (-0.5) * 2.0 * MATMUL(animal_vr_i , solutions(INT(lhs(i,5)),:))
diag_ele = diag_ele + 2.0 * animal_vr_i
pre_eq_no = cur_eq_no
CASE(1002) ! animal effect, both parents are known, sire or dam is diagonal
temp_rhs = temp_rhs - (-0.5) * 2.0 * MATMUL(animal_vr_i , solutions(INT(lhs(i,4)),:))
temp_rhs = temp_rhs - 0.25 * 2.0 * MATMUL(animal_vr_i , solutions(INT(lhs(i,5)),:))
diag_ele = diag_ele + 0.25 * 2.0 * animal_vr_i
pre_eq_no = cur_eq_no
CASE(1003) ! animal effect, one parent is known, animal is diagonal
temp_rhs = temp_rhs - (-0.5) * (4.0/3.0) * MATMUL(animal_vr_i , solutions(INT(lhs(i,3)),:))
diag_ele = diag_ele + (4.0/3.0) * animal_vr_i
pre_eq_no = cur_eq_no
CASE(1004) ! animal effect, one parents is known, parent is diagonal
temp_rhs = temp_rhs - (-0.5) * (4.0/3.0) * MATMUL(animal_vr_i , solutions(INT(lhs(i,3)),:))
diag_ele = diag_ele + 0.25 * (4.0/3.0) * animal_vr_i
pre_eq_no = cur_eq_no
CASE(1005) ! animal effect, no parents is known
diag_ele = diag_ele + animal_vr_i
pre_eq_no = cur_eq_no
CASE(2001) ! permanent environmental effect
diag_ele = diag_ele + pe_vr_i
pre_eq_no = cur_eq_no
END SELECT END DO
! end reading lhs ! calculate last solution
pre_sol = solutions(pre_eq_no,:) ! store old solution
CALL geninv(diag_ele, diag_ele_i)
solutions(pre_eq_no,:) = MATMUL(diag_ele_i, 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,'th iteration''s solutions'
! DO i = 1, no_of_rc
! WRITE(*,*) i, solutions(i)
! END DO epsilon = epsilon / (no_of_eq) ! write iteration number and epsilon
write(*,*) 'iteration = ', iteration , ', epsilon = ', epsilon
IF (epsilon < criteria) THEN
EXIT
END IF END DO
! end iteration ! open file for writing solution
OPEN(UNIT=30, FILE='sol.dat', STATUS='REPLACE', ACTION='WRITE', IOSTAT=status)
DO i = 1, no_of_eq
WRITE(30,*) i, solutions(i,:)
END DO
CLOSE(30)

 

END PROGRAM strrm_solve

 

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

 

Inverse Program

 

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 strrm_setup.f95 -o strrm_setup.exe
strrm_setup.exe
sort -nk 3 lhs.dat -o sorted_lhs.dat
gfortran inverse.f95 strrm_solve.f95 -o strrm_solve.exe
strrm_solve.exe
위 문장을strrim_run.bat로 저장

 

프로그램 실행 화면

 

 

 

파일

1265875683_strrm.zip

다운로드

1265875683_strrm.zip