【发布时间】:2015-08-10 00:55:15
【问题描述】:
我是 Fortran 的初学者,用它来解决经济学中的数学模型。 我写了一个运行时崩溃的程序(如下)。 Windows 显示该错误消息,但似乎没有任何反应。我不知道可能出了什么问题,我已经对其进行了调试,并且它没有在构建日志上显示任何错误消息。 我正在使用 gfortran。
!Value function iteration program with matrices, based on the matlab similar code
program valuefuncmat
implicit none
integer,parameter::dp=selected_real_kind(15) !double precision
integer,parameter::lower=0,upper=50 !lower and upper bounds on capital grid (lower = b in the model)
real(dp),parameter::hex1=.01,hex2=.1 !steps on capital grid, finer on low levels of capital
integer::i,j !loop counters
integer,parameter::sze1=((15-lower)/hex1)+1,sze2=( ((upper-15)/hex2)+1)
integer,parameter::sze=sze1+sze2
real(dp)::alinha1(sze1),alinha2(sze2),agrid(sze)
11 format(i5)
12 format(f35.10)
interface !interface of the functions used
function kron(A,B)
implicit none
integer,parameter::dp=selected_real_kind(15) !double precision
real(dp),dimension(:,:), intent(in)::A,B !dummy arguments
real(dp),dimension(size(A,1)*size(B,1),size(A,2)*size(B,2))::kron
end function kron
function repvec(A,n)
implicit none
integer,parameter::dp=selected_real_kind(15) !double precision, not essencial in the code
real(dp),dimension(:), intent(in)::A
integer::n
real(dp),dimension(n,size(A))::repvec
end function repvec
function repmat(A,n)
implicit none
integer,parameter::dp=selected_real_kind(15) !double precision, not essencial in the code
real(dp),dimension(:,:), intent(in)::A
integer::n
real(dp),dimension(size(A)*n,size(A))::repmat
end function repmat
function kronvec(A,V)
implicit none
integer,parameter::dp=selected_real_kind(15) !double precision, not essencial in the code
real(dp),dimension(:,:),intent(in)::A
integer,dimension(:),intent(in)::V
real(dp),dimension(size(A,1)*size(V),size(A,2))::kronvec
end function kronvec
function kronvec2(A,V)
implicit none
integer,parameter::dp=selected_real_kind(15) !double precision, not essencial in the code
integer,dimension(:,:),intent(in)::A
real(dp),dimension(:),intent(in)::V
real(dp),dimension(size(A,1)*size(V),size(A,2))::kronvec2
end function kronvec2
end interface
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!Markov Matrix and Productivity Grid!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(dp), dimension(3,3), parameter::P=reshape((/ .91,.08,.01,.13,.82,.05,.005,.07,.925 /),(/3,3/),(/.0/),(/2,1/)) !markov matrix
real(dp), dimension(3), parameter::theta_grid=(/ .1,.5,1.0 /) !the grid of productivity
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(dp),parameter::A=1.0,alpha=.25,beta=.98 !Parameters of Functional Forms
real(dp),parameter::r=0.05,w=1.0,T=1.0 !Parameters for testing the value function iteration
integer,parameter::y=1
real(dp),dimension(size(theta_grid),size(agrid))::vIN=.0,vOUT=1.0,apol
real(dp)::transf(size(theta_grid)*size(agrid),size(agrid)),amatrix(size(agrid),size(agrid))
real(dp)::C(size(agrid)*size(theta_grid),size(agrid))
real(dp)::Rzao(size(agrid)*size(theta_grid),size(agrid)),tvOUT(size(theta_grid)*size(agrid))
integer::indexa(size(theta_grid)*size(agrid)),um(size(agrid))=1,onesm(size(theta_grid))=1,umzao(size(agrid),size(agrid))=1
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!initializing the capital grid!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do i=1,sze1
alinha1(i)=lower+(i-1)*hex1
end do
do i=1,sze2
alinha2(i)=alinha1((sze1))+hex2 +(i-1)*hex2
end do
agrid(:sze1)=alinha1
agrid(sze1+1:sze)=alinha2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!Beginning the actual value function iteration algorithm!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
transf(1:y*size(agrid),:)=T !defining the types that will receive the transfer
Amatrix=transpose(repvec((1+r)*agrid,size(agrid)))-repvec(agrid,size(agrid)) !every combination of a and a' on square matrix of rank equal to the size of agrid
C=repmat(Amatrix,size(theta_grid))+w*kronvec2(umzao,theta_grid)+transf
where (C<=0) C=0.00000000000000001
Rzao=log(C)
do while (maxval(abs(vIN-vOUT))>0.001)
vIN=vOUT
tvOUT=maxval(Rzao+beta*matmul(kronvec(P,um),vIN),dim=2)
indexa=maxloc(Rzao+beta*matmul(kronvec(P,um),vIN),dim=2)
vOUT=reshape(tvOUT,(/size(theta_grid),size(agrid)/),order=(/2,1/))
end do
do i=1,size(vOUT,1)
write(*,'(20G12.4)'),vOUT(i,:)
end do
end program
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!Kronecker Product Routine!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function kron(A,B)
implicit none
integer,parameter::dp=selected_real_kind(15) !double precision, not essencial in the code
real(dp),dimension(:,:), intent(in)::A,B !dummy arguments
real(dp),dimension(size(A,1)*size(B,1),size(A,2)*size(B,2))::kron !output matrix of the kronecker product
integer::i,j !loop counters
do i=1,size(A,1)
do j=1,size(A,2)
kron(1+size(B,1)*(i-1):size(B,1)+size(B,1)*(i-1),1+size(B,2)*(j-1):size(B,2)+size(B,2)*(j-1))=A(i,j)*B
end do
end do
end function kron
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!Function that calculates a kronecker product between a matrix and a vector!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Vector must be integer!!!!!!!!!!!!!!!!!!!!!!!!!!!
function kronvec(A,V)
implicit none
integer,parameter::dp=selected_real_kind(15) !double precision, not essencial in the code
real(dp),dimension(:,:),intent(in)::A
integer,dimension(:),intent(in)::V
real(dp),dimension(size(A,1)*size(V),size(A,2))::kronvec
integer::i
do i=1,size(V)
kronvec(1+size(A,1)*(i-1):size(A,1)+size(A,1)*(i-1),:)=A*V(i)
end do
end function kronvec
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!Function that calculates a kronecker product between a matrix and a vector!!!!
!!!!!!!!!!!!!!!Matrix must be integer and vector real!!!!!!!!!!!!!!!!!!!!!!!!!!!
function kronvec2(A,V)
implicit none
integer,parameter::dp=selected_real_kind(15) !double precision, not essencial in the code
integer,dimension(:,:),intent(in)::A
real(dp),dimension(:),intent(in)::V
real(dp),dimension(size(A,1)*size(V),size(A,2))::kronvec2
integer::i
do i=1,size(V)
kronvec2(1+size(A,1)*(i-1):size(A,1)+size(A,1)*(i-1),:)=A*V(i)
end do
end function kronvec2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!Function that repeats a vector n times verticaly!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function repvec(A,n)
implicit none
integer,parameter::dp=selected_real_kind(15) !double precision, not essencial in the code
real(dp),dimension(:), intent(in)::A
integer::n
real(dp),dimension(n,size(A))::repvec
integer::i
do i=1,n
repvec(i,:)=A
end do
end function repvec
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!Function that repeats a matrix n times verticaly!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function repmat(A,n)
implicit none
integer,parameter::dp=selected_real_kind(15) !double precision, not essencial in the code
real(dp),dimension(:,:), intent(in)::A
integer::n
real(dp),dimension(size(A)*n,size(A))::repmat
integer::i
do i=1,size(A)
repmat(1+size(A)*(i-1):size(A)+size(A)*(i-1),:)=A
end do
end function repmat
【问题讨论】:
-
请参阅stackoverflow.com/help/mcve,了解如何让您的问题变得更好。
-
尝试用
gfortran -fcheck=bounds yourcode.f90编译。错误消息似乎来自repmat()例程,其中size(A)用于尺寸。但这似乎是size(A,1)或size(A,2)的错字。还要注意repmat()中do-loop的上限,应该是n而不是size(A)(根据例程的注释和repvec()比较)。 -
您必须添加所有相关信息(至少是错误消息!)并在此处发布之前尝试尽可能减少程序。编译调试时使用
gfortran -g -fbacktrace -fcheck=all -Wall。
标签: crash fortran fortran90 gfortran