主题:[讨论]求助!同样一段代码放在function中和主程序中结果不同
求助大家:
同样的一段代码,一个写入module的function,并在主程序中调用,一个直接放在主程序中,出来的结果却不一样
两者的差别主要是,后者为0的地方,前者为NaN或者极小、极大的值
怎样才能通过调用function得到正确的值呢?
谢谢...
提供代码:主要是fppb这个函数...
module trustregion
Implicit none
integer, parameter :: n = 300
integer, parameter :: p = 7, pn = 22, df = 3
real, parameter :: lambda = 2.0
integer :: id_begin(p), id_end(p), id(df)
data id_begin /2,5,8,11,14,17,20/
data id_end /4,7,10,13,16,19,22/
integer :: i
real :: x(n,pn), y(n)
real :: beta(pn) = (/(0.0,i = 1,pn)/)
real :: temp_beta(pn)
integer :: ite1 = 1, ite2 = 1
real :: tol = 1D-5
real :: delta0 = 0.5d0
integer :: itemax = 200
real :: fd(df), sd(df,df)
complex :: eigen_vector(df,df), eigen_value(df)
contains
!--------------------------------------------------------
!---------------- read data ----------------------------
subroutine loaddata()
integer :: i
open(100,file = "train_splice.txt")
do i = 1, n
read(100,*) y(i),x(i,:)
end do
end subroutine loaddata
!--------------------------------------------------------
!---------------- probability ----------------------------
function fprob(bt)
real :: bt(pn)
real :: fprob(n)
integer :: i
do i = 1,n
fprob(i) = exp(sum(x(i,:) * bt))/ (1 + exp(sum(x(i,:) * bt)))
end do
end function fprob
!--------------------------------------------------------
!---------------- beta first derivative----------------
function fpb(bt)
real :: bt(pn)
real :: fpb(pn)
integer :: j
do j = 1,pn
fpb(j) = - sum(x(:,j) * (y - fprob(bt)))
end do
end function fpb
!--------------------------------------------------------
!---------------- beta second derivative---------------
function fppb(bt)
real :: bt(pn), pr(n)
real :: fppb(pn,pn),tppb(pn,pn)
integer :: i,j,k
pr = fprob(bt)
do i = 1, n
do j = 1, pn
do k = 1, pn
tppb(j,k) = x(i,j) * x(i,k) * pr(i) * (1-pr(i))
end do
end do
fppb = fppb + tppb
end do
end function fppb
end module
!-----------------------------------------------
! main program
!-----------------------------------------------
program main
use trustregion
implicit none
real :: ppb1(pn,pn),ppb2(pn,pn)
real :: bt(pn), pr(n)
real :: tppb(pn,pn)
integer :: j,k
call loaddata()
ppb1 = fppb(beta)
write(*,*) "ppb1",ppb1
pr = fprob(bt)
do i = 1, n
do j = 1, pn
do k = 1, pn
tppb(j,k) = x(i,j) * x(i,k) * pr(i) * (1-pr(i))
end do
end do
ppb2 = ppb2 + tppb
end do
write(*,*) "ppb2",ppb2
end program main