1
프로젝트에 대해 여러 개의 Fortran 90 코드를 R로 가져 오려고합니다. 그들은 처음에는 mex (Fortran 루틴의 matlab 통합) 유형 컴파일을 염두에두고 작성되었습니다. 이 인 코드 중 하나가 어떻게 생겼는지 :R CMD SHLIB를 사용하여 Mexfile 컴파일
# include <fintrf.h>
subroutine mexFunction(nlhs, plhs, nrhs, prhs)
!--------------------------------------------------------------
! MEX file for VFI3FCN routine
!
! log M_{t,t+1} = log \beta + gamma (x_t - x_{t+1})
! gamma = gamA + gamB (x_t - xbar)
!
!--------------------------------------------------------------
implicit none
mwPointer plhs(*), prhs(*)
integer nlhs, nrhs
mwPointer mxGetM, mxGetPr, mxCreateDoubleMatrix
mwPointer nk, nkp, nz, nx, nh
mwSize col_hxz, col_hz, col_xz
! check for proper number of arguments.
if(nrhs .ne. 31) then
call mexErrMsgTxt('31 input variables required.')
elseif(nlhs .ne. 4) then
call mexErrMsgTxt('4 output variables required.')
endif
! get the size of the input array.
nk = mxGetM(prhs(5))
nx = mxGetM(prhs(7))
nz = mxGetM(prhs(11))
nh = mxGetM(prhs(14))
nkp = mxGetM(prhs(16))
col_hxz = nx*nz*nh
col_xz = nx*nz
col_hz = nz*nh
! create matrix for the return arguments.
plhs(1) = mxCreateDoubleMatrix(nk, col_hxz, 0)
plhs(2) = mxCreateDoubleMatrix(nk, col_hxz, 0)
plhs(3) = mxCreateDoubleMatrix(nk, col_hxz, 0)
plhs(4) = mxCreateDoubleMatrix(nk, col_hxz, 0)
call vfi3fcnIEccB(%val(mxGetPr(plhs(1))), nkp)
return
end
subroutine vfi3fcnIEccB(optK, V, I, div, & ! output variables
alp1, alp2, alp3, V0, k, nk, x, xbar, nx, Qx, z, nz, Qz, h, nh, kp, &
alpha, beta, delta, f, gamA, gamB, gP, gN, istar, kmin, kmtrx, ksubm, hmtrx, xmtrx, zmtrx, &
nkp, col_hxz, col_xz, col_hz)
use ifwin
implicit none
! specify input and output variables
integer, intent(in) :: nk, nkp, nx, nz, nh, col_hxz, col_xz, col_hz
real*8, intent(out) :: V(nk, col_hxz), optK(nk, col_hxz), I(nk, col_hxz), div(nk, col_hxz)
real*8, intent(in) :: V0(nk, col_hxz), k(nk), kp(nkp), x(nx), z(nz), Qx(nx, nx), Qz(nz, nz), h(nh)
real*8, intent(in) :: alp1, alp2, alp3, xbar, kmin, alpha, gP, gN, beta, delta, gamA, gamB, f, istar
real*8, intent(in) :: kmtrx(nk, col_hxz), ksubm(nk, col_hz), zmtrx(nk, col_hxz), xmtrx(nk, col_hxz), hmtrx(nk, col_hxz)
! specify intermediate variables
real*8 :: Res(nk, col_hxz), Obj(nk, col_hxz), optKold(nk, col_hxz), Vold(nk, col_hxz), tmpEMV(nkp, col_hz), tmpI(nkp), &
tmpObj(nkp, col_hz), tmpA(nk, col_hxz), tmpQ(nx*nh, nh), detM(nx), stoM(nx), g(nkp), tmpInd(nh, nz)
real*8 :: Qh(nh, nh, nx), Qxh(nx*nh, nx*nh), Qzxh(col_hxz, col_hxz)
real*8 :: hp, d(nh), errK, errV, T1, lapse
integer :: ix, ih, iter, optJ(col_hz), ik, iz, ind(nh, col_xz), subindex(nx, col_hz)
logical*4 :: statConsole
! construct the transition matrix for kh --- there are nx number of these transition matrix: 3-d
Qh = 0.0
do ix = 1, nx
do ih = 1, nh
! compute the predicted next period kh
hp = alp1 + alp2*h(ih) + alp3*(x(ix) - xbar)
! construct transition probability vector
d = abs(h - hp) + 1D-32
Qh(:, ih, ix) = (1/d)/sum(1/d)
end do
end do
! construct the compound transition matrix over (z x h) space
! compound the (x h) space
Qxh = 0.0
do ix = 1, nx
call kron(tmpQ, Qx(:, ix), Qh(:, :, ix), nx, 1, nh, nh)
Qxh(:, (ix - 1)*nh + 1 : ix*nh) = tmpQ
end do
! compound the (z x h) space: h changes the faster, followed by x, and z changes the slowest
call kron(Qzxh, Qz, Qxh, nz, nz, nx*nh, nx*nh)
! available funds for the firm
Res = dexp(xmtrx + zmtrx + hmtrx)*(kmtrx**alpha) + (1 - delta)*kmtrx - f
! initializing
Obj = 0.0
optK = 0.0
optKold = optK + 1.0
Vold = V0
! Some Intermediate Variables Used in Stochastic Discount Factor
detM = beta*dexp((gamA - gamB*xbar)*x + gamB*x**2)
stoM = -(gamA - gamB*xbar + gamB*x)
! Intermediate index vector to facilitate submatrix extracting
ind = reshape((/1 : col_hxz : 1/), (/nh, col_xz/))
do ix = 1, nx
tmpInd = ind(:, ix : col_xz : nx)
do iz = 1, nz
subindex(ix, (iz - 1)*nh + 1 : iz*nh) = tmpInd(:, iz)
end do
end do
! start iterations
errK = 1.0
errV = 1.0
iter = 0
T1 = secnds(0.0)
do
if (errV <= 1D-3 .AND. errK <= 1D-8) then
exit
else
iter = iter + 1
do ix = 1, nx
! next period value function by linear interpolation: nkp by nz*nh matrix
call interp1(tmpEMV, k, detM(ix)*(matmul(dexp(stoM(ix)*xmtrx)*Vold, Qzxh(:, subindex(ix, :)))) - ksubm, kp, &
nk, nkp, col_hz)
! maximize the right-hand size of Bellman equation on EACH grid point of capital stock
do ik = 1, nk
! with istar tmpI is no longer investment but a linear transformation of that
tmpI = (kp - (1.0 - delta)*k(ik))/k(ik) - istar
where (tmpI >= 0.0)
g = gP
elsewhere
g = gN
end where
tmpObj = tmpEMV - spread((g/2.0)*(tmpI**2)*k(ik), 2, col_hz)
! direct discrete maximization
Obj(ik, subindex(ix, :)) = maxval(tmpObj, 1)
optJ = maxloc(tmpObj, 1)
optK(ik, subindex(ix, :)) = kp(optJ)
end do
end do
! update value function and impose limited liability condition
V = max(Res + Obj, 1D-16)
! convergence criterion
errK = maxval(abs(optK - optKold))
errV = maxval(abs(V - Vold))
! revise Initial Guess
Vold = V
optKold = optK
! visual
if (modulo(iter, 50) == 0) then
lapse = secnds(T1)
statConsole = AllocConsole()
print "(a, f10.7, a, f10.7, a, f8.1, a)", " errV:", errV, " errK:", errK, " Time:", lapse, "s"
end if
end if
end do
! visual check on errors
lapse = secnds(T1)
statConsole = AllocConsole()
print "(a, f10.7, a, f10.7, a, f8.1, a)", " errV:", errV, " errK:", errK, " Time:", lapse, "s"
! optimal investment and dividend
I = optK - (1.0 - delta)*kmtrx
tmpA = I/kmtrx - istar
where (tmpA >= 0)
div = Res - optK - (gP/2.0)*(tmpA**2)*kmtrx
elsewhere
div = Res - optK - (gN/2.0)*(tmpA**2)*kmtrx
end where
return
end
subroutine interp1(v, x, y, u, m, n, col)
!-------------------------------------------------------------------------------------------------------
! Linear interpolation routine similar to interp1 with 'linear' as method parameter in Matlab
!
! OUTPUT:
! v - function values on non-grid points (n by col matrix)
!
! INPUT:
! x - grid (m by one vector)
! y - function defined on the grid x (m by col matrix)
! u - non-grid points on which y(x) is to be interpolated (n by one vector)
! m - length of x and y vectors
! n - length of u and v vectors
! col - number of columns of v and y matrices
!
! Four ways to pass array arguments:
! 1. Use explicit-shape arrays and pass the dimension as an argument(most efficient)
! 2. Use assumed-shape arrays and use interface to call external subroutine
! 3. Use assumed-shape arrays and make subroutine internal by using "contains"
! 4. Use assumed-shape arrays and put interface in a module then use module
!
! This subroutine is equavilent to the following matlab call
! v = interp1(x, y, u, 'linear', 'extrap') with x (m by 1), y (m by col), u (n by 1), and v (n by col)
!------------------------------------------------------------------------------------------------------
implicit none
integer :: m, n, col, i, j
real*8, intent(out) :: v(n, col)
real*8, intent(in) :: x(m), y(m, col), u(n)
real*8 :: prob
do i = 1, n
if (u(i) < x(1)) then
! extrapolation to the left
v(i, :) = y(1, :) - (y(2, :) - y(1, :)) * ((x(1) - u(i))/(x(2) - x(1)))
else if (u(i) > x(m)) then
! extrapolation to the right
v(i, :) = y(m, :) + (y(m, :) - y(m-1, :)) * ((u(i) - x(m))/(x(m) - x(m-1)))
else
! interpolation
! find the j such that x(j) <= u(i) < x(j+1)
call bisection(x, u(i), m, j)
prob = (u(i) - x(j))/(x(j+1) - x(j))
v(i, :) = y(j, :)*(1 - prob) + y(j+1, :)*prob
end if
end do
end subroutine interp1
subroutine bisection(list, element, m, k)
!--------------------------------------------------------------------------------
! find index k in list such that (list(k) <= element < list(k+1)
!--------------------------------------------------------------------------------
implicit none
integer*4 :: m, k, first, last, half
real*8 :: list(m), element
first = 1
last = m
do
if (first == (last-1)) exit
half = (first + last)/2
if (element < list(half)) then
! discard second half
last = half
else
! discard first half
first = half
end if
end do
k = first
end subroutine bisection
subroutine kron(K, A, B, rowA, colA, rowB, colB)
!--------------------------------------------------------------------------------
! Perform K = kron(A, B); translated directly from kron.m
!
! OUTPUT:
! K -- rowA*rowB by colA*colB matrix
!
! INPUT:
! A -- rowA by colA matrix
! B -- rowB by colB matrix
! rowA, colA, rowB, colB -- integers containing shape information
!--------------------------------------------------------------------------------
implicit none
integer, intent(in) :: rowA, colA, rowB, colB
real*8, intent(in) :: A(rowA, colA), B(rowB, colB)
real*8, intent(out) :: K(rowA*rowB, colA*colB)
integer :: t1(rowA*rowB), t2(colA*colB), i, ia(rowA*rowB), ja(colA*colB), ib(rowA*rowB), jb(colA*colB)
t1 = (/ (i, i = 0, (rowA*rowB - 1)) /)
ia = int(t1/rowB) + 1
ib = mod(t1, rowB) + 1
t2 = (/ (i, i = 0, (colA*colB - 1)) /)
ja = int(t2/colB) + 1
jb = mod(t2, colB) + 1
K = A(ia, ja)*B(ib, jb)
end subroutine kron
내 초기 계획은 mexFunction 서브 루틴을 제거하고 R CMD SHLIB 명령을 사용하여 주요 포트란 서브 루틴을 컴파일했지만 나는이 곳을 찾을 알고하지 Rtools 컴파일러로 실행 ifwin 라이브러리는 인텔 Intell Fortran 컴파일러 폴더에 있습니다.
그래서 내 첫 번째 질문은 다음과 같습니다
1) ifwin 라이브러리와 내가 포함 할 필요가 다른 라이브러리를 찾을 수 rtools에게하는 방법이 있나요? 또는 컴파일러가 필요한 라이브러리를 찾고 컴파일하도록 종속 라이브러리를 R CMD SHLIB 명령에 포함시키는 방법이 있습니까?
2) 대답이 '아니오'인 경우, Matlab의 컴파일 된 버전을 R에서 어떻게 사용할 수 있습니까? mex Zhang_4.f90 명령을 사용하여 오류없이 MATLAB에서 파일을 컴파일 할 수 있습니다.
3) Visual Studio 2015에서 환경을 설정하여 인텔 컴파일러를 사용하는 R에서 Fortran 서브 루틴을 컴파일 할 수 있습니까? 나는 mexFunction 서브 루틴을 꺼내 나머지 코드를 컴파일하려고
, 내가받을 다음과 같은 오류 : 다음 사용하지 않는 코드를 다시 작성 다른 방법이라고 생각하지 않습니다
D:\SS_Programming\Fortran>R CMD SHLIB Zhang_4.f90
c:/Rtools/mingw_64/bin/gfortran -O2 -mtune=core2 -c Zhang_4.f90 -o
Zhang_4.o
Zhang_4.f90:6.4:
use ifwin
1
Fatal Error: Can't open module file 'ifwin.mod' for reading at (1): No
such file or directory
make: *** [Zhang_4.o] Error 1
Warning message:
running command 'make -f "C:/PROGRA~1/R/R-34~1.2/etc/x64/Makeconf" -f
"C:/PROGRA~1/R/R-34~1.2/share/make/winshlib.mk"
SHLIB_LDFLAGS='$(SHLIB_FCLDFLAGS)' SHLIB_LD='$(SHLIB_FCLD)'
SHLIB="Zhang_4.dll" SHLIB_LIBADD='$(FCLIBS)' WIN=64 TCLBIN=64
OBJECTS="Zhang_4.o"' had status 2
코드 ('mexFunction')는 명시 적으로'mex' 심볼을 호출합니다. R에서 어떻게 작동할까요? 또는 해당 부분을 삭제 하시겠습니까? 인텔 컴파일러는 IFWIN을 찾을 위치를 알아야합니다. –
블라디미르 답장을 보내 주셔서 감사합니다. 내가 말했듯이, mexFunction 서브 루틴없이 R 컴파일러로 코드 컴파일을 시도했다. R 컴파일러가 ifwin 라이브러리를 어디서 찾을 지 모르기 때문에 rtools 컴파일러를 가리키는 방법을 모르므로 작동하지 않습니다. 저는 지난 3 일 동안 행운이없는 온라인 솔루션을 찾고있었습니다. –
시도한 정확한 명령과 정확한 오류 메시지를 게시하십시오. [ask] 및 [mcve]를 참조하십시오. –