Implements a logarithmic reduction Fortran routine and the associated test
parent
3ea0baf21e
commit
855887b249
|
@ -1,17 +1,18 @@
|
||||||
mex_PROGRAMS = cycle_reduction
|
mex_PROGRAMS = cycle_reduction
|
||||||
|
|
||||||
nodist_cycle_reduction_SOURCES = \
|
nodist_cycle_reduction_SOURCES = \
|
||||||
mexFunction.f08 \
|
|
||||||
matlab_mex.F08 \
|
matlab_mex.F08 \
|
||||||
blas_lapack.F08
|
blas_lapack.F08 \
|
||||||
|
mexFunction.f08
|
||||||
|
|
||||||
BUILT_SOURCES = $(nodist_cycle_reduction_SOURCES)
|
BUILT_SOURCES = $(nodist_cycle_reduction_SOURCES)
|
||||||
CLEANFILES = $(nodist_cycle_reduction_SOURCES)
|
CLEANFILES = $(nodist_cycle_reduction_SOURCES)
|
||||||
|
|
||||||
mexFunction.o: matlab_mex.mod blas.mod lapack.mod
|
mexFunction.o: matlab_mex.mod blas.mod lapack.mod
|
||||||
|
|
||||||
|
mexFunction.mod: mexFunction.o
|
||||||
|
|
||||||
%.f08: $(top_srcdir)/../../sources/cycle_reduction/%.f08
|
%.f08: $(top_srcdir)/../../sources/cycle_reduction/%.f08
|
||||||
$(LN_S) -f $< $@
|
$(LN_S) -f $< $@
|
||||||
|
|
||||||
%.F08: $(top_srcdir)/../../sources/cycle_reduction/%.F08
|
%.F08: $(top_srcdir)/../../sources/cycle_reduction/%.F08
|
||||||
$(LN_S) -f $< $@
|
$(LN_S) -f $< $@
|
||||||
|
|
|
@ -0,0 +1,19 @@
|
||||||
|
mex_PROGRAMS = logarithmic_reduction
|
||||||
|
|
||||||
|
nodist_logarithmic_reduction_SOURCES = \
|
||||||
|
matlab_mex.F08 \
|
||||||
|
blas_lapack.F08 \
|
||||||
|
mexFunction.f08
|
||||||
|
|
||||||
|
BUILT_SOURCES = $(nodist_logarithmic_reduction_SOURCES)
|
||||||
|
CLEANFILES = $(nodist_logarithmic_reduction_SOURCES)
|
||||||
|
|
||||||
|
mexFunction.o: matlab_mex.mod blas.mod lapack.mod
|
||||||
|
|
||||||
|
mexFunction.mod: mexFunction.o
|
||||||
|
|
||||||
|
|
||||||
|
%.f08: $(top_srcdir)/../../sources/logarithmic_reduction/%.f08
|
||||||
|
$(LN_S) -f $< $@
|
||||||
|
%.F08: $(top_srcdir)/../../sources/logarithmic_reduction/%.F08
|
||||||
|
$(LN_S) -f $< $@
|
|
@ -1,6 +1,6 @@
|
||||||
ACLOCAL_AMFLAGS = -I ../../../m4
|
ACLOCAL_AMFLAGS = -I ../../../m4
|
||||||
|
|
||||||
SUBDIRS = mjdgges kronecker bytecode block_kalman_filter sobol perfect_foresight_problem num_procs block_trust_region disclyap_fast libkordersim local_state_space_iterations folded_to_unfolded_dr k_order_simul k_order_mean cycle_reduction
|
SUBDIRS = mjdgges kronecker bytecode block_kalman_filter sobol perfect_foresight_problem num_procs block_trust_region disclyap_fast libkordersim local_state_space_iterations folded_to_unfolded_dr k_order_simul k_order_mean cycle_reduction logarithmic_reduction
|
||||||
|
|
||||||
# libdynare++ must come before gensylv and k_order_perturbation
|
# libdynare++ must come before gensylv and k_order_perturbation
|
||||||
if ENABLE_MEX_DYNAREPLUSPLUS
|
if ENABLE_MEX_DYNAREPLUSPLUS
|
||||||
|
|
|
@ -170,6 +170,8 @@ AC_CONFIG_FILES([Makefile
|
||||||
num_procs/Makefile
|
num_procs/Makefile
|
||||||
block_trust_region/Makefile
|
block_trust_region/Makefile
|
||||||
disclyap_fast/Makefile
|
disclyap_fast/Makefile
|
||||||
cycle_reduction/Makefile])
|
cycle_reduction/Makefile
|
||||||
|
logarithmic_reduction/Makefile])
|
||||||
|
|
||||||
|
|
||||||
AC_OUTPUT
|
AC_OUTPUT
|
||||||
|
|
|
@ -0,0 +1,2 @@
|
||||||
|
include ../mex.am
|
||||||
|
include ../../logarithmic_reduction.am
|
|
@ -1,6 +1,6 @@
|
||||||
ACLOCAL_AMFLAGS = -I ../../../m4
|
ACLOCAL_AMFLAGS = -I ../../../m4
|
||||||
|
|
||||||
SUBDIRS = mjdgges kronecker bytecode block_kalman_filter sobol perfect_foresight_problem num_procs block_trust_region disclyap_fast libkordersim local_state_space_iterations folded_to_unfolded_dr k_order_simul k_order_mean cycle_reduction
|
SUBDIRS = mjdgges kronecker bytecode block_kalman_filter sobol perfect_foresight_problem num_procs block_trust_region disclyap_fast libkordersim local_state_space_iterations folded_to_unfolded_dr k_order_simul k_order_mean cycle_reduction logarithmic_reduction
|
||||||
|
|
||||||
# libdynare++ must come before gensylv and k_order_perturbation
|
# libdynare++ must come before gensylv and k_order_perturbation
|
||||||
if ENABLE_MEX_DYNAREPLUSPLUS
|
if ENABLE_MEX_DYNAREPLUSPLUS
|
||||||
|
|
|
@ -173,6 +173,7 @@ AC_CONFIG_FILES([Makefile
|
||||||
num_procs/Makefile
|
num_procs/Makefile
|
||||||
block_trust_region/Makefile
|
block_trust_region/Makefile
|
||||||
disclyap_fast/Makefile
|
disclyap_fast/Makefile
|
||||||
cycle_reduction/Makefile])
|
cycle_reduction/Makefile
|
||||||
|
logarithmic_reduction/Makefile])
|
||||||
|
|
||||||
AC_OUTPUT
|
AC_OUTPUT
|
||||||
|
|
|
@ -0,0 +1,3 @@
|
||||||
|
EXEEXT = .mex
|
||||||
|
include ../mex.am
|
||||||
|
include ../../logarithmic_reduction.am
|
|
@ -26,7 +26,8 @@ EXTRA_DIST = \
|
||||||
num_procs \
|
num_procs \
|
||||||
block_trust_region \
|
block_trust_region \
|
||||||
disclyap_fast \
|
disclyap_fast \
|
||||||
cycle_reduction
|
cycle_reduction \
|
||||||
|
logarithmic_reduction
|
||||||
|
|
||||||
clean-local:
|
clean-local:
|
||||||
rm -rf `find mex/sources -name *.o`
|
rm -rf `find mex/sources -name *.o`
|
||||||
|
|
|
@ -48,6 +48,57 @@ module blas
|
||||||
real(real64), dimension(*), intent(inout) :: y
|
real(real64), dimension(*), intent(inout) :: y
|
||||||
end subroutine dgemv
|
end subroutine dgemv
|
||||||
end interface
|
end interface
|
||||||
|
|
||||||
|
contains
|
||||||
|
|
||||||
|
! Updating a matrix using blas DGEMM
|
||||||
|
! C <- alpha*op(A)*op(B) + beta*C
|
||||||
|
subroutine matmul_add(opA, opB, alpha, A, B, beta, C)
|
||||||
|
character, intent(in) :: opA, opB
|
||||||
|
real(real64), dimension(:,:), intent(in) :: A, B
|
||||||
|
real(real64), dimension(:,:), intent(inout) :: C
|
||||||
|
real(real64), intent(in) :: alpha, beta
|
||||||
|
integer(blint) :: m, n, k, lda, ldb
|
||||||
|
if (opA == "N") then
|
||||||
|
m = int(size(A,1), blint)
|
||||||
|
k = int(size(A,2), blint)
|
||||||
|
lda = m
|
||||||
|
else
|
||||||
|
m = int(size(A,2), blint)
|
||||||
|
k = int(size(A,1), blint)
|
||||||
|
lda = k
|
||||||
|
end if
|
||||||
|
if (opB == "N") then
|
||||||
|
n = int(size(B,2), blint)
|
||||||
|
ldb = k
|
||||||
|
else
|
||||||
|
n = int(size(B,1), blint)
|
||||||
|
ldb = n
|
||||||
|
end if
|
||||||
|
#ifdef DEBUG
|
||||||
|
if ( (opA /= "N") .and. (opA /= "T") .and. (opA /= "C") ) then
|
||||||
|
print *, "opA must be either N, T or C"
|
||||||
|
end if
|
||||||
|
if ( (opB /= "N") .and. (opB /= "T") .and. (opB /= "C") ) then
|
||||||
|
print *, "opB must be either N, T or C"
|
||||||
|
end if
|
||||||
|
if (((opA == "N") .and. (opB == "N") .and. (size(A,2) /= size(B,1))) .or.&
|
||||||
|
&((opA == "N") .and. (opB /= "N") .and. (size(A,2) /= size(B,2))) .or.&
|
||||||
|
&((opA /= "N") .and. (opB == "N") .and. (size(A,1) /= size(B,1))) .or.&
|
||||||
|
&((opA /= "N") .and. (opB /= "N") .and. (size(A,1) /= size(B,2)))) &
|
||||||
|
then
|
||||||
|
print *, "Inconsistent number of columns of op(A) and number of rows &
|
||||||
|
&of op(B)"
|
||||||
|
end if
|
||||||
|
if (m /= size(C,1)) then
|
||||||
|
print *, "Inconsistent number of rows of op(A) and number of rows &
|
||||||
|
&of C"
|
||||||
|
end if
|
||||||
|
#endif
|
||||||
|
call dgemm(opA, opB, m, n, k, alpha, A, lda, B, &
|
||||||
|
ldb, beta, C, m)
|
||||||
|
end subroutine
|
||||||
|
|
||||||
end module blas
|
end module blas
|
||||||
|
|
||||||
module lapack
|
module lapack
|
||||||
|
@ -104,5 +155,39 @@ module lapack
|
||||||
real(real64) :: dlange
|
real(real64) :: dlange
|
||||||
end function dlange
|
end function dlange
|
||||||
end interface
|
end interface
|
||||||
|
|
||||||
|
contains
|
||||||
|
|
||||||
|
! Computing inv(A)*B using lapack DGESV
|
||||||
|
! B <- inv(A)*B
|
||||||
|
subroutine left_divide(A, B, ipiv, info)
|
||||||
|
real(real64), dimension(:,:), intent(inout) :: A
|
||||||
|
real(real64), dimension(:,:), intent(inout) :: B
|
||||||
|
integer(blint), dimension(:), intent(inout) :: ipiv
|
||||||
|
integer(blint), intent(inout) :: info
|
||||||
|
integer(blint) :: m, n
|
||||||
|
#ifdef DEBUG
|
||||||
|
if (size(ipiv) /= size(A,1)) then
|
||||||
|
print *, "Inconsistent number of rows of A and size of ipiv"
|
||||||
|
end if
|
||||||
|
if (size(A,2) /= size(B,1)) then
|
||||||
|
print *, "Inconsistent number of columns of A and number of rows of B"
|
||||||
|
end if
|
||||||
|
#endif
|
||||||
|
m = int(size(A,1), blint)
|
||||||
|
n = int(size(B,2), blint)
|
||||||
|
call dgesv(m, n, A, m, ipiv, B, m, info)
|
||||||
|
end subroutine left_divide
|
||||||
|
|
||||||
|
! Matlab-like norm function calling lapack DLANGE routine
|
||||||
|
real(real64) function norm(A, c)
|
||||||
|
real(real64), dimension(:,:), intent(in) :: A
|
||||||
|
character, intent(in) :: c
|
||||||
|
real(real64), dimension(:), allocatable :: work
|
||||||
|
integer(blint) :: nrow
|
||||||
|
nrow = int(size(A,1),blint)
|
||||||
|
norm = dlange(c, nrow, int(size(A,2), blint), &
|
||||||
|
&A, nrow, work)
|
||||||
|
end function norm
|
||||||
|
|
||||||
end module lapack
|
end module lapack
|
||||||
|
|
|
@ -18,51 +18,13 @@
|
||||||
! Implements the cyclic reduction algorithm described in
|
! Implements the cyclic reduction algorithm described in
|
||||||
! D.A. Bini, G. Latouche, B. Meini (2002), "Solving matrix polynomial equations arising in queueing problems", Linear Algebra and its Applications 340, pp. 222-244
|
! D.A. Bini, G. Latouche, B. Meini (2002), "Solving matrix polynomial equations arising in queueing problems", Linear Algebra and its Applications 340, pp. 222-244
|
||||||
! D.A. Bini, B. Meini (1996), "On the solution of a nonlinear matrix equation arising in queueing problems", SIAM J. Matrix Anal. Appl. 17, pp. 906-926.
|
! D.A. Bini, B. Meini (1996), "On the solution of a nonlinear matrix equation arising in queueing problems", SIAM J. Matrix Anal. Appl. 17, pp. 906-926.
|
||||||
|
module c_reduction
|
||||||
module reduction
|
|
||||||
use matlab_mex
|
|
||||||
use lapack
|
use lapack
|
||||||
use blas
|
use blas
|
||||||
|
use matlab_mex
|
||||||
implicit none
|
implicit none
|
||||||
|
|
||||||
contains
|
contains
|
||||||
! Matlab-like norm function calling lapack DLANGE routine
|
|
||||||
real(real64) function norm(A, c)
|
|
||||||
real(real64), dimension(:,:), intent(in) :: A
|
|
||||||
character, intent(in) :: c
|
|
||||||
real(real64), dimension(:), allocatable :: work
|
|
||||||
integer(blint) :: nrow
|
|
||||||
nrow = int(size(A,1),blint)
|
|
||||||
norm = dlange(c, nrow, int(size(A,2), blint), &
|
|
||||||
&A, nrow, work)
|
|
||||||
end function norm
|
|
||||||
|
|
||||||
! Updating a matrix using blas DGEMM
|
|
||||||
! C <- alpha*A*B + beta*C
|
|
||||||
subroutine updatemat(alpha, A, B, beta, C)
|
|
||||||
real(real64), dimension(:,:), intent(in) :: A, B
|
|
||||||
real(real64), dimension(:,:), intent(inout) :: C
|
|
||||||
real(real64), intent(in) :: alpha, beta
|
|
||||||
integer(blint) :: m, n, k
|
|
||||||
m = int(size(A,1), blint)
|
|
||||||
n = int(size(B,2), blint)
|
|
||||||
k = int(size(A,2), blint)
|
|
||||||
call dgemm("N", "N", m, n, k, alpha, A, m, B, &
|
|
||||||
k, beta, C, m)
|
|
||||||
end subroutine
|
|
||||||
|
|
||||||
! Computing inv(A)*B using lapack DGESV
|
|
||||||
! B <- inv(A)*B
|
|
||||||
subroutine invA_B(A, B, ipiv, info)
|
|
||||||
real(real64), dimension(:,:), intent(inout) :: A
|
|
||||||
real(real64), dimension(:,:), intent(inout) :: B
|
|
||||||
integer(blint), dimension(:), intent(inout) :: ipiv
|
|
||||||
integer(blint), intent(inout) :: info
|
|
||||||
integer(blint) :: m, n
|
|
||||||
m = int(size(A,1), blint)
|
|
||||||
n = int(size(B,2), blint)
|
|
||||||
call dgesv(m, n, A, m, ipiv, B, m, info)
|
|
||||||
end subroutine invA_B
|
|
||||||
|
|
||||||
! Cycle reduction algorithm
|
! Cycle reduction algorithm
|
||||||
subroutine cycle_reduction(A0, A1, A2, X, cvg_tol, check, info)
|
subroutine cycle_reduction(A0, A1, A2, X, cvg_tol, check, info)
|
||||||
|
@ -93,9 +55,9 @@ loop: do
|
||||||
! Computing [A0;A2]*(A1\[A0 A2])
|
! Computing [A0;A2]*(A1\[A0 A2])
|
||||||
A1_tmp = A1i
|
A1_tmp = A1i
|
||||||
invA1_A02 = A02
|
invA1_A02 = A02
|
||||||
call invA_B(A1_tmp, invA1_A02, ipiv, info_inv)
|
call left_divide(A1_tmp, invA1_A02, ipiv, info_inv)
|
||||||
call updatemat(1._real64, A02(:,1:n), invA1_A02, 0._real64, Q0)
|
call matmul_add("N", "N", 1._real64, A02(:,1:n), invA1_A02, 0._real64, Q0)
|
||||||
call updatemat(1._real64, A02(:,n+1:dn), invA1_A02, 0._real64, Q2)
|
call matmul_add("N", "N", 1._real64, A02(:,n+1:dn), invA1_A02, 0._real64, Q2)
|
||||||
! Updating A02, A1 and Ahat
|
! Updating A02, A1 and Ahat
|
||||||
A1i = A1i - Q0(:,n+1:dn) - Q2(:,1:n)
|
A1i = A1i - Q0(:,n+1:dn) - Q2(:,1:n)
|
||||||
A02(:,1:n) = -Q0(:,1:n)
|
A02(:,1:n) = -Q0(:,1:n)
|
||||||
|
@ -121,16 +83,16 @@ loop: do
|
||||||
|
|
||||||
! Computing X = -Ahat\A0
|
! Computing X = -Ahat\A0
|
||||||
X = -A0
|
X = -A0
|
||||||
call invA_B(Ahat, X, ipiv, info_inv)
|
call left_divide(Ahat, X, ipiv, info_inv)
|
||||||
|
|
||||||
! Checking residuals if necessary
|
! Checking residuals if necessary
|
||||||
if (check) then
|
if (check) then
|
||||||
! A1_tmp <- A2*X + A1
|
! A1_tmp <- A2*X + A1
|
||||||
A1_tmp = A1
|
A1_tmp = A1
|
||||||
call updatemat(1._real64, A2, X, 1._real64, A1_tmp)
|
call matmul_add("N", "N", 1._real64, A2, X, 1._real64, A1_tmp)
|
||||||
! A0_tmp <- A1_tmp*X + A0
|
! A0_tmp <- A1_tmp*X + A0
|
||||||
A0_tmp = A0
|
A0_tmp = A0
|
||||||
call updatemat(1._real64, A1_tmp, X, 1._real64, A0_tmp)
|
call matmul_add("N", "N", 1._real64, A1_tmp, X, 1._real64, A0_tmp)
|
||||||
residual = norm(A0_tmp, "1")
|
residual = norm(A0_tmp, "1")
|
||||||
if (residual>cvg_tol) then
|
if (residual>cvg_tol) then
|
||||||
info(1) = 403._c_double
|
info(1) = 403._c_double
|
||||||
|
@ -145,11 +107,10 @@ loop: do
|
||||||
end if
|
end if
|
||||||
end subroutine cycle_reduction
|
end subroutine cycle_reduction
|
||||||
|
|
||||||
end module reduction
|
end module c_reduction
|
||||||
|
|
||||||
subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
|
subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
|
||||||
use reduction
|
use c_reduction
|
||||||
use matlab_mex
|
|
||||||
implicit none
|
implicit none
|
||||||
|
|
||||||
type(c_ptr), dimension(*), intent(in), target :: prhs
|
type(c_ptr), dimension(*), intent(in), target :: prhs
|
||||||
|
|
|
@ -0,0 +1,213 @@
|
||||||
|
! Copyright © 2022 Dynare Team
|
||||||
|
!
|
||||||
|
! This file is part of Dynare.
|
||||||
|
!
|
||||||
|
! Dynare is free software: you can redistribute it and/or modify
|
||||||
|
! it under the terms of the GNU General Public License as published by
|
||||||
|
! the Free Software Foundation, either version 3 of the License, or
|
||||||
|
! (at your option) any later version.
|
||||||
|
!
|
||||||
|
! Dynare is distributed in the hope that it will be useful,
|
||||||
|
! but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||||
|
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||||
|
! GNU General Public License for more details.
|
||||||
|
!
|
||||||
|
! You should have received a copy of the GNU General Public License
|
||||||
|
! along with Dynare. If not, see <https://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
|
! Implements the logarithmic reduction algorithm described in
|
||||||
|
! D.A. Bini, G. Latouche, B. Meini (2002), "Solving matrix polynomial equations arising in queueing problems", Linear Algebra and its Applications 340, pp. 222-244
|
||||||
|
|
||||||
|
module l_reduction
|
||||||
|
use lapack
|
||||||
|
use blas
|
||||||
|
use matlab_mex
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
contains
|
||||||
|
|
||||||
|
! Logarithmic reduction algorithm
|
||||||
|
subroutine logarithmic_reduction(A0, A1, A2, G, cvg_tol, check, max_it, info)
|
||||||
|
real(real64), dimension(:,:), intent(in) :: A0, A1, A2
|
||||||
|
real(real64), dimension(:,:), intent(inout) :: G
|
||||||
|
real(real64), intent(in) :: cvg_tol
|
||||||
|
logical, intent(in) :: check
|
||||||
|
integer, intent(in) :: max_it
|
||||||
|
real(c_double), dimension(2), intent(inout) :: info
|
||||||
|
|
||||||
|
real(real64), dimension(:,:), allocatable :: tmp_inv, tmp, D02, Id, P, &
|
||||||
|
&A1_tmp, A0_tmp, G_new, P_new
|
||||||
|
integer :: it, n, dn, i, j
|
||||||
|
integer(blint) :: info_inv
|
||||||
|
integer(blint), dimension(:), allocatable :: ipiv
|
||||||
|
real(real64) :: crit
|
||||||
|
character(kind=c_char, len=10) :: cvg_tol_str, residual_str
|
||||||
|
|
||||||
|
n = size(A0,1)
|
||||||
|
dn = 2*n
|
||||||
|
allocate(D02(n,dn), ipiv(n), Id(n,n), tmp(n,dn), G_new(n,n), P_new(n,n))
|
||||||
|
info = [0._c_double,0._c_double]
|
||||||
|
! Set the identity matrix
|
||||||
|
do j=1,n
|
||||||
|
do i=1,n
|
||||||
|
if (i == j) then
|
||||||
|
Id(i,j) = 1.
|
||||||
|
else
|
||||||
|
Id(i,j) = 0.
|
||||||
|
end if
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
|
||||||
|
! Initialization: D02_0
|
||||||
|
D02(:,1:n) = A0
|
||||||
|
D02(:,n+1:dn) = A2
|
||||||
|
tmp_inv = -A1
|
||||||
|
call left_divide(tmp_inv, D02, ipiv, info_inv)
|
||||||
|
G = D02(:,1:n)
|
||||||
|
P = D02(:,n+1:dn)
|
||||||
|
it = 0
|
||||||
|
loop: do
|
||||||
|
! Computing E1_(i+1)
|
||||||
|
tmp_inv = Id
|
||||||
|
call matmul_add("N", "N", -1._real64, D02(:,1:n), D02(:,n+1:dn), 1._real64, tmp_inv)
|
||||||
|
call matmul_add("N", "N", -1._real64, D02(:,n+1:dn), D02(:,1:n), 1._real64, tmp_inv)
|
||||||
|
! Computing matrices [D0^2 D2^2], id est -E0_(i+1) and -E2_(i+1)
|
||||||
|
call matmul_add("N", "N", 1._real64, D02(:,1:n), D02(:,1:n), 0._real64, tmp(:,1:n))
|
||||||
|
call matmul_add("N", "N", 1._real64, D02(:,n+1:dn), D02(:,n+1:dn), 0._real64, tmp(:,n+1:dn))
|
||||||
|
! Computing D0_(i+1) and D2_(i+1)
|
||||||
|
call left_divide(tmp_inv, tmp, ipiv, info_inv)
|
||||||
|
! Computing G_(i+1) = G_i + P_i*D0_(i+1)
|
||||||
|
G_new = G
|
||||||
|
call matmul_add("N", "N", 1._real64, P, tmp(:,1:n), 1._real64, G_new)
|
||||||
|
! Computing P_(i+1)
|
||||||
|
call matmul_add("N", "N", 1._real64, P, tmp(:,n+1:dn), 0._real64, P_new)
|
||||||
|
crit = norm(G_new - G,"1")
|
||||||
|
! Checking for stopping conditions
|
||||||
|
if (crit < cvg_tol) then
|
||||||
|
exit loop
|
||||||
|
elseif (it == max_it) then
|
||||||
|
info(1) = 411._c_double
|
||||||
|
info(2) = real(log(crit), c_double)
|
||||||
|
exit loop
|
||||||
|
elseif (isnan(crit) .or. (info_inv /= 0_blint)) then
|
||||||
|
info(1) = 412._c_double
|
||||||
|
info(2) = -1._c_double
|
||||||
|
exit loop
|
||||||
|
end if
|
||||||
|
it = it + 1
|
||||||
|
G = G_new
|
||||||
|
P = P_new
|
||||||
|
D02 = tmp
|
||||||
|
end do loop
|
||||||
|
|
||||||
|
! Checking residuals if necessary
|
||||||
|
if (check) then
|
||||||
|
! A1_tmp <- A2*X + A1
|
||||||
|
A1_tmp = A1
|
||||||
|
call matmul_add("N", "N", 1._real64, A2, G, 1._real64, A1_tmp)
|
||||||
|
! A0_tmp <- A1_tmp*X + A0
|
||||||
|
A0_tmp = A0
|
||||||
|
call matmul_add("N", "N", 1._real64, A1_tmp, G, 1._real64, A0_tmp)
|
||||||
|
crit = norm(A0_tmp, "1")
|
||||||
|
if (crit>cvg_tol) then
|
||||||
|
info(1) = 413._c_double
|
||||||
|
info(2) = real(log(crit), c_double)
|
||||||
|
write (cvg_tol_str,"(es8.2)") cvg_tol
|
||||||
|
write (residual_str,"(es8.2)") crit
|
||||||
|
call mexPrintf("The norm of the residual is "&
|
||||||
|
&// trim(residual_str) // &
|
||||||
|
&", whereas the tolerance criterion is " &
|
||||||
|
// trim(cvg_tol_str) // "." )
|
||||||
|
end if
|
||||||
|
end if
|
||||||
|
end subroutine logarithmic_reduction
|
||||||
|
|
||||||
|
end module l_reduction
|
||||||
|
|
||||||
|
subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
|
||||||
|
use l_reduction
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
type(c_ptr), dimension(*), intent(in), target :: prhs
|
||||||
|
type(c_ptr), dimension(*), intent(out) :: plhs
|
||||||
|
integer(c_int), intent(in), value :: nlhs, nrhs
|
||||||
|
|
||||||
|
integer :: i, n, max_it
|
||||||
|
character(kind=c_char, len=2) :: num2str
|
||||||
|
real(real64) :: cvg_tol
|
||||||
|
real(real64), dimension(2) :: info
|
||||||
|
real(real64), dimension(:,:), pointer, contiguous :: A0, A1, A2, X
|
||||||
|
logical :: check
|
||||||
|
|
||||||
|
! 0. Checking the consistency and validity of input arguments
|
||||||
|
if (nrhs < 5) then
|
||||||
|
call mexErrMsgTxt("Must have at least 5 inputs")
|
||||||
|
end if
|
||||||
|
if (nrhs > 6) then
|
||||||
|
call mexErrMsgTxt("Too many input arguments")
|
||||||
|
end if
|
||||||
|
if (nlhs > 2) then
|
||||||
|
call mexErrMsgTxt("Too many output arguments")
|
||||||
|
end if
|
||||||
|
|
||||||
|
do i=1,3
|
||||||
|
if (.not. (c_associated(prhs(i)) .and. mxIsDouble(prhs(i)) .and. &
|
||||||
|
(.not. mxIsComplex(prhs(i))) .and. (.not. mxIsSparse(prhs(i))))) then
|
||||||
|
write (num2str,"(i2)") i
|
||||||
|
call mexErrMsgTxt("Argument " // trim(num2str) // " should be a real dense matrix")
|
||||||
|
end if
|
||||||
|
end do
|
||||||
|
do i=4,5
|
||||||
|
if (.not. (c_associated(prhs(i)) .and. mxIsScalar(prhs(i)) .and. &
|
||||||
|
mxIsNumeric(prhs(i)))) then
|
||||||
|
write (num2str,"(i2)") i
|
||||||
|
call mexErrMsgTxt("Argument " // trim(num2str) // " should be a numeric scalar")
|
||||||
|
end if
|
||||||
|
end do
|
||||||
|
|
||||||
|
check = .false.
|
||||||
|
if (nrhs == 6) then
|
||||||
|
if (.not. (c_associated(prhs(6)))) then
|
||||||
|
call mexErrMsgTxt("Argument 6 should be a Matlab object")
|
||||||
|
else
|
||||||
|
if (.not. (mxIsEmpty(prhs(6)))) then
|
||||||
|
check = .true.
|
||||||
|
end if
|
||||||
|
end if
|
||||||
|
end if
|
||||||
|
|
||||||
|
n = int(mxGetM(prhs(1))) ! Order of the considered matrices
|
||||||
|
if ((n /= mxGetN(prhs(1))) & ! Number of columns of A0
|
||||||
|
&.or. (n /= mxGetM(prhs(2))) & ! Number of lines of A1
|
||||||
|
&.or. (n /= mxGetN(prhs(2))) & ! Number of columns of A1
|
||||||
|
&.or. (n /= mxGetM(prhs(3))) & ! Number of lines of A2
|
||||||
|
&.or. (n /= mxGetN(prhs(3))) & ! Number of columns of A2
|
||||||
|
) then
|
||||||
|
call mexErrMsgTxt("Input dimension mismatch")
|
||||||
|
end if
|
||||||
|
|
||||||
|
! 1. Storing the relevant information in Fortran format
|
||||||
|
A2(1:n,1:n) => mxGetPr(prhs(1))
|
||||||
|
A1(1:n,1:n) => mxGetPr(prhs(2))
|
||||||
|
A0(1:n,1:n) => mxGetPr(prhs(3))
|
||||||
|
cvg_tol = mxGetScalar(prhs(4))
|
||||||
|
max_it = int(mxGetScalar(prhs(5)))
|
||||||
|
info = [0._c_double,0._c_double]
|
||||||
|
|
||||||
|
plhs(1) = mxCreateDoubleMatrix(int(n, mwSize), int(n, mwSize), mxREAL)
|
||||||
|
X(1:n,1:n) => mxGetPr(plhs(1))
|
||||||
|
|
||||||
|
! 2. Calling the Logarithmic Reduction algorithm
|
||||||
|
call logarithmic_reduction(A0, A1, A2, X, cvg_tol, check, max_it, info)
|
||||||
|
|
||||||
|
! 3. Editing the information output if necessary
|
||||||
|
if (nlhs == 2) then
|
||||||
|
if (info(1) == 0._c_double) then
|
||||||
|
plhs(2) = mxCreateDoubleScalar(0._c_double)
|
||||||
|
else
|
||||||
|
plhs(2) = mxCreateDoubleMatrix(1_mwSize, 2_mwSize, mxREAL)
|
||||||
|
mxGetPr(plhs(2)) = info
|
||||||
|
end if
|
||||||
|
end if
|
||||||
|
|
||||||
|
end subroutine mexFunction
|
|
@ -1224,7 +1224,9 @@ M_TRS_FILES += run_block_byte_tests_matlab.m.trs \
|
||||||
test_aggregate_routine_1_2_3.m.trs \
|
test_aggregate_routine_1_2_3.m.trs \
|
||||||
run_kronecker_tests.m.trs \
|
run_kronecker_tests.m.trs \
|
||||||
nonlinearsolvers.m.trs \
|
nonlinearsolvers.m.trs \
|
||||||
cyclereduction.m.trs
|
cyclereduction.m.trs \
|
||||||
|
logarithmicreduction.m.trs
|
||||||
|
|
||||||
|
|
||||||
M_XFAIL_TRS_FILES = $(patsubst %.mod, %.m.trs, $(XFAIL_MODFILES))
|
M_XFAIL_TRS_FILES = $(patsubst %.mod, %.m.trs, $(XFAIL_MODFILES))
|
||||||
|
|
||||||
|
@ -1237,7 +1239,9 @@ O_TRS_FILES += run_block_byte_tests_octave.o.trs \
|
||||||
histval_initval_file_unit_tests.o.trs \
|
histval_initval_file_unit_tests.o.trs \
|
||||||
run_kronecker_tests.o.trs \
|
run_kronecker_tests.o.trs \
|
||||||
nonlinearsolvers.o.trs \
|
nonlinearsolvers.o.trs \
|
||||||
cyclereduction.o.trs
|
cyclereduction.o.trs \
|
||||||
|
logarithmicreduction.o.trs
|
||||||
|
|
||||||
|
|
||||||
O_XFAIL_TRS_FILES = $(patsubst %.mod, %.o.trs, $(XFAIL_MODFILES))
|
O_XFAIL_TRS_FILES = $(patsubst %.mod, %.o.trs, $(XFAIL_MODFILES))
|
||||||
|
|
||||||
|
|
|
@ -0,0 +1,140 @@
|
||||||
|
debug = false;
|
||||||
|
|
||||||
|
if debug
|
||||||
|
[top_test_dir, ~, ~] = fileparts(mfilename('fullpath'));
|
||||||
|
else
|
||||||
|
top_test_dir = getenv('TOP_TEST_DIR');
|
||||||
|
end
|
||||||
|
|
||||||
|
addpath([top_test_dir filesep '..' filesep 'matlab']);
|
||||||
|
|
||||||
|
if ~debug
|
||||||
|
% Test Dynare Version
|
||||||
|
if ~strcmp(dynare_version(), getenv('DYNARE_VERSION'))
|
||||||
|
error('Incorrect version of Dynare is being tested')
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
dynare_config;
|
||||||
|
|
||||||
|
NumberOfTests = 0;
|
||||||
|
testFailed = 0;
|
||||||
|
|
||||||
|
if ~debug
|
||||||
|
skipline()
|
||||||
|
disp('*** TESTING: logarithmicreduction.m ***');
|
||||||
|
end
|
||||||
|
|
||||||
|
matlab_cr_path = [top_test_dir filesep '..' filesep 'matlab' filesep 'missing' filesep 'mex' filesep 'logarithmic_reduction'];
|
||||||
|
addpath(matlab_cr_path);
|
||||||
|
logarithmic_reduction_matlab = @logarithmic_reduction;
|
||||||
|
rmpath(matlab_cr_path);
|
||||||
|
|
||||||
|
if isoctave
|
||||||
|
addpath([top_test_dir filesep '..' filesep 'mex' filesep 'octave']);
|
||||||
|
else
|
||||||
|
addpath([top_test_dir filesep '..' filesep 'mex' filesep 'matlab']);
|
||||||
|
end
|
||||||
|
|
||||||
|
t0 = clock;
|
||||||
|
|
||||||
|
% Set the dimension of the problem to be solved.
|
||||||
|
n = 500;
|
||||||
|
% Set the equation to be solved
|
||||||
|
A = eye(n);
|
||||||
|
B = diag(30*ones(n,1)); B(1,1) = 20; B(end,end) = 20; B = B - diag(10*ones(n-1,1),-1); B = B - diag(10*ones(n-1,1),1);
|
||||||
|
C = diag(15*ones(n,1)); C = C - diag(5*ones(n-1,1),-1); C = C - diag(5*ones(n-1,1),1);
|
||||||
|
cvg_tol = 1e-12;
|
||||||
|
X1 = zeros(n,n);
|
||||||
|
X2 = zeros(n,n);
|
||||||
|
|
||||||
|
% 1. Solve the equation with the Matlab logarithmic reduction algorithm
|
||||||
|
NumberOfTests = NumberOfTests+1;
|
||||||
|
tElapsed1 = 0.;
|
||||||
|
try
|
||||||
|
tic; [X1,info] = logarithmic_reduction_matlab(A,B,C,cvg_tol,300,[0.]); tElapsed1 = toc;
|
||||||
|
disp(['Elapsed time for the Matlab logarithmic reduction algorithm is: ' num2str(tElapsed1) ' (n=' int2str(n) ').'])
|
||||||
|
R = norm(C+B*X1+A*X1*X1,1);
|
||||||
|
if (R > cvg_tol)
|
||||||
|
testFailed = testFailed+1;
|
||||||
|
if debug
|
||||||
|
dprintf('Matlab logarithmic_reduction solution is wrong')
|
||||||
|
end
|
||||||
|
end
|
||||||
|
catch
|
||||||
|
testFailed = testFailed+1;
|
||||||
|
if debug
|
||||||
|
dprintf('Matlab logarithmic_reduction failed')
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
% 2. Solve the equation with the Fortran logarithmic reduction algorithm
|
||||||
|
NumberOfTests = NumberOfTests+1;
|
||||||
|
tElapsed2 = 0.;
|
||||||
|
try
|
||||||
|
tic; [X2,info] = logarithmic_reduction(A,B,C,cvg_tol,300,[0.]); tElapsed2 = toc;
|
||||||
|
disp(['Elapsed time for the Fortran logarithmic reduction algorithm is: ' num2str(tElapsed2) ' (n=' int2str(n) ').'])
|
||||||
|
R = norm(C+B*X2+A*X2*X2,1);
|
||||||
|
if (R > cvg_tol)
|
||||||
|
testFailed = testFailed+1;
|
||||||
|
if debug
|
||||||
|
dprintf('Fortran logarithmic_reduction solution is wrong')
|
||||||
|
end
|
||||||
|
end
|
||||||
|
catch
|
||||||
|
testFailed = testFailed+1;
|
||||||
|
if debug
|
||||||
|
dprintf('Fortran logarithmic_reduction failed')
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
% 3. Compare solutions of the Fortran and Matlab routines
|
||||||
|
NumberOfTests = NumberOfTests+1;
|
||||||
|
if (norm(X1 - X2, 1) > cvg_tol)
|
||||||
|
testFailed = testFailed+1;
|
||||||
|
if debug
|
||||||
|
dprintf('Fortran and Matlab logarithmic reduction solutions differ');
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
% Compare the Fortran and Matlab execution time
|
||||||
|
if debug
|
||||||
|
if tElapsed1<tElapsed2
|
||||||
|
skipline()
|
||||||
|
dprintf('Matlab logarithmic reduction is %5.2f times faster than its Fortran counterpart.', tElapsed2/tElapsed1)
|
||||||
|
skipline()
|
||||||
|
else
|
||||||
|
skipline()
|
||||||
|
dprintf('Fortran logarithmic reduction is %5.2f times faster than its Matlab counterpart.', tElapsed1/tElapsed2)
|
||||||
|
skipline()
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
t1 = clock;
|
||||||
|
|
||||||
|
if ~debug
|
||||||
|
cd(getenv('TOP_TEST_DIR'));
|
||||||
|
else
|
||||||
|
dprintf('FAILED tests: %i', testFailed)
|
||||||
|
end
|
||||||
|
|
||||||
|
if isoctave
|
||||||
|
fid = fopen('logarithmicreduction.o.trs', 'w+');
|
||||||
|
else
|
||||||
|
fid = fopen('logarithmicreduction.m.trs', 'w+');
|
||||||
|
end
|
||||||
|
if testFailed
|
||||||
|
fprintf(fid,':test-result: FAIL\n');
|
||||||
|
else
|
||||||
|
fprintf(fid,':test-result: PASS\n');
|
||||||
|
end
|
||||||
|
fprintf(fid,':number-tests: %i\n', NumberOfTests);
|
||||||
|
fprintf(fid,':number-failed-tests: %i\n', testFailed);
|
||||||
|
fprintf(fid,':list-of-passed-tests: logarithmicreduction.m\n');
|
||||||
|
fprintf(fid,':elapsed-time: %f\n', etime(t1, t0));
|
||||||
|
fclose(fid);
|
||||||
|
|
||||||
|
if ~debug
|
||||||
|
exit;
|
||||||
|
end
|
||||||
|
|
Loading…
Reference in New Issue