diff --git a/matlab/logarithmic_reduction.m b/matlab/missing/mex/logarithmic_reduction/logarithmic_reduction.m
similarity index 100%
rename from matlab/logarithmic_reduction.m
rename to matlab/missing/mex/logarithmic_reduction/logarithmic_reduction.m
diff --git a/mex/build/cycle_reduction.am b/mex/build/cycle_reduction.am
index a5aa334c9..ad774a29f 100644
--- a/mex/build/cycle_reduction.am
+++ b/mex/build/cycle_reduction.am
@@ -1,17 +1,18 @@
mex_PROGRAMS = cycle_reduction
nodist_cycle_reduction_SOURCES = \
- mexFunction.f08 \
matlab_mex.F08 \
- blas_lapack.F08
+ blas_lapack.F08 \
+ mexFunction.f08
BUILT_SOURCES = $(nodist_cycle_reduction_SOURCES)
CLEANFILES = $(nodist_cycle_reduction_SOURCES)
mexFunction.o: matlab_mex.mod blas.mod lapack.mod
+mexFunction.mod: mexFunction.o
+
%.f08: $(top_srcdir)/../../sources/cycle_reduction/%.f08
$(LN_S) -f $< $@
-
%.F08: $(top_srcdir)/../../sources/cycle_reduction/%.F08
$(LN_S) -f $< $@
diff --git a/mex/build/logarithmic_reduction.am b/mex/build/logarithmic_reduction.am
new file mode 100644
index 000000000..dc79612f9
--- /dev/null
+++ b/mex/build/logarithmic_reduction.am
@@ -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 $< $@
diff --git a/mex/build/matlab/Makefile.am b/mex/build/matlab/Makefile.am
index e6617b39c..775b5dac7 100644
--- a/mex/build/matlab/Makefile.am
+++ b/mex/build/matlab/Makefile.am
@@ -1,6 +1,6 @@
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
if ENABLE_MEX_DYNAREPLUSPLUS
diff --git a/mex/build/matlab/configure.ac b/mex/build/matlab/configure.ac
index 97f5f3368..3cd63e872 100644
--- a/mex/build/matlab/configure.ac
+++ b/mex/build/matlab/configure.ac
@@ -170,6 +170,8 @@ AC_CONFIG_FILES([Makefile
num_procs/Makefile
block_trust_region/Makefile
disclyap_fast/Makefile
- cycle_reduction/Makefile])
+ cycle_reduction/Makefile
+ logarithmic_reduction/Makefile])
+
AC_OUTPUT
diff --git a/mex/build/matlab/logarithmic_reduction/Makefile.am b/mex/build/matlab/logarithmic_reduction/Makefile.am
new file mode 100644
index 000000000..fdc1ac2f0
--- /dev/null
+++ b/mex/build/matlab/logarithmic_reduction/Makefile.am
@@ -0,0 +1,2 @@
+include ../mex.am
+include ../../logarithmic_reduction.am
diff --git a/mex/build/octave/Makefile.am b/mex/build/octave/Makefile.am
index 7057795ed..4369fc38f 100644
--- a/mex/build/octave/Makefile.am
+++ b/mex/build/octave/Makefile.am
@@ -1,6 +1,6 @@
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
if ENABLE_MEX_DYNAREPLUSPLUS
diff --git a/mex/build/octave/configure.ac b/mex/build/octave/configure.ac
index 9806965a3..cf746d2a7 100644
--- a/mex/build/octave/configure.ac
+++ b/mex/build/octave/configure.ac
@@ -173,6 +173,7 @@ AC_CONFIG_FILES([Makefile
num_procs/Makefile
block_trust_region/Makefile
disclyap_fast/Makefile
- cycle_reduction/Makefile])
+ cycle_reduction/Makefile
+ logarithmic_reduction/Makefile])
AC_OUTPUT
diff --git a/mex/build/octave/logarithmic_reduction/Makefile.am b/mex/build/octave/logarithmic_reduction/Makefile.am
new file mode 100644
index 000000000..c2e5519d5
--- /dev/null
+++ b/mex/build/octave/logarithmic_reduction/Makefile.am
@@ -0,0 +1,3 @@
+EXEEXT = .mex
+include ../mex.am
+include ../../logarithmic_reduction.am
diff --git a/mex/sources/Makefile.am b/mex/sources/Makefile.am
index 9ebd3edf8..3c7335a5a 100644
--- a/mex/sources/Makefile.am
+++ b/mex/sources/Makefile.am
@@ -26,7 +26,8 @@ EXTRA_DIST = \
num_procs \
block_trust_region \
disclyap_fast \
- cycle_reduction
+ cycle_reduction \
+ logarithmic_reduction
clean-local:
rm -rf `find mex/sources -name *.o`
diff --git a/mex/sources/blas_lapack.F08 b/mex/sources/blas_lapack.F08
index a597ca500..6b981f3e2 100644
--- a/mex/sources/blas_lapack.F08
+++ b/mex/sources/blas_lapack.F08
@@ -48,6 +48,57 @@ module blas
real(real64), dimension(*), intent(inout) :: y
end subroutine dgemv
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
module lapack
@@ -104,5 +155,39 @@ module lapack
real(real64) :: dlange
end function dlange
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
diff --git a/mex/sources/cycle_reduction/mexFunction.f08 b/mex/sources/cycle_reduction/mexFunction.f08
index e0e5fe8a3..1f9238fdd 100644
--- a/mex/sources/cycle_reduction/mexFunction.f08
+++ b/mex/sources/cycle_reduction/mexFunction.f08
@@ -18,51 +18,13 @@
! 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, 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 reduction
- use matlab_mex
+module c_reduction
use lapack
use blas
+ use matlab_mex
implicit none
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
subroutine cycle_reduction(A0, A1, A2, X, cvg_tol, check, info)
@@ -93,9 +55,9 @@ loop: do
! Computing [A0;A2]*(A1\[A0 A2])
A1_tmp = A1i
invA1_A02 = A02
- call invA_B(A1_tmp, invA1_A02, ipiv, info_inv)
- call updatemat(1._real64, A02(:,1:n), invA1_A02, 0._real64, Q0)
- call updatemat(1._real64, A02(:,n+1:dn), invA1_A02, 0._real64, Q2)
+ call left_divide(A1_tmp, invA1_A02, ipiv, info_inv)
+ call matmul_add("N", "N", 1._real64, A02(:,1:n), invA1_A02, 0._real64, Q0)
+ call matmul_add("N", "N", 1._real64, A02(:,n+1:dn), invA1_A02, 0._real64, Q2)
! Updating A02, A1 and Ahat
A1i = A1i - Q0(:,n+1:dn) - Q2(:,1:n)
A02(:,1:n) = -Q0(:,1:n)
@@ -121,16 +83,16 @@ loop: do
! Computing X = -Ahat\A0
X = -A0
- call invA_B(Ahat, X, ipiv, info_inv)
+ call left_divide(Ahat, X, ipiv, info_inv)
! Checking residuals if necessary
if (check) then
! A1_tmp <- A2*X + 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 = 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")
if (residual>cvg_tol) then
info(1) = 403._c_double
@@ -145,11 +107,10 @@ loop: do
end if
end subroutine cycle_reduction
-end module reduction
+end module c_reduction
subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
- use reduction
- use matlab_mex
+ use c_reduction
implicit none
type(c_ptr), dimension(*), intent(in), target :: prhs
diff --git a/mex/sources/logarithmic_reduction/mexFunction.f08 b/mex/sources/logarithmic_reduction/mexFunction.f08
new file mode 100644
index 000000000..349f4a35e
--- /dev/null
+++ b/mex/sources/logarithmic_reduction/mexFunction.f08
@@ -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 .
+
+! 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
\ No newline at end of file
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 81f3ec588..7741b9a82 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -1224,7 +1224,9 @@ M_TRS_FILES += run_block_byte_tests_matlab.m.trs \
test_aggregate_routine_1_2_3.m.trs \
run_kronecker_tests.m.trs \
nonlinearsolvers.m.trs \
- cyclereduction.m.trs
+ cyclereduction.m.trs \
+ logarithmicreduction.m.trs
+
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 \
run_kronecker_tests.o.trs \
nonlinearsolvers.o.trs \
- cyclereduction.o.trs
+ cyclereduction.o.trs \
+ logarithmicreduction.o.trs
+
O_XFAIL_TRS_FILES = $(patsubst %.mod, %.o.trs, $(XFAIL_MODFILES))
diff --git a/tests/logarithmicreduction.m b/tests/logarithmicreduction.m
new file mode 100644
index 000000000..d6e765181
--- /dev/null
+++ b/tests/logarithmicreduction.m
@@ -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