From 3ea0baf21e1e6b4fb438dbb62c427e82adfff3d0 Mon Sep 17 00:00:00 2001 From: NormannR Date: Tue, 4 Oct 2022 15:25:32 +0200 Subject: [PATCH 1/3] Matlab logarithmic reduction: making norms and error codes consistent with cyclic reduction codes --- matlab/logarithmic_reduction.m | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/matlab/logarithmic_reduction.m b/matlab/logarithmic_reduction.m index b5a3292b9..3ab7d67fd 100644 --- a/matlab/logarithmic_reduction.m +++ b/matlab/logarithmic_reduction.m @@ -77,20 +77,27 @@ while cc>tol && kk<=maxit tmp1 = (eye(n)-tmp0*[tmp0(:,n+i);tmp0(:,i)])\[tmp0(:,i)*tmp0(:,i),tmp0(:,n+i)*tmp0(:,n+i)]; X1 = X0 + U0*tmp1(:,n+i); U1 = U0*tmp1(:,i); - cc = max(max(abs(X1-X0))); + cc = norm(X1-X0,1); + if isnan(cc) + info(1) = 412; + info(2) = -1.; + return + end X0 = X1; U0 = U1; tmp0 = tmp1; kk = kk+1; end if kk==maxit - disp(['logarithmic_reduction:: Convergence not achieved after ' int2str(maxit) ' iterations!']); - info = 1; + info(1) = 411; + info(2) = log(cc); end if nargin>5 && check - if max(max(abs(A*X1*X1 + B*X1 + C)))>tol - disp(['logarithmic_reduction:: Algotithm did not converge to the solution of the matrix quadratic equation!']); - info = 1; + res = norm(A*X1*X1 + B*X1 + C, 1); + if res>tol + info(1) = 413; + info(2) = log(res); + dprintf('The norm of the residual is %s whereas the tolerance criterion is %s', num2str(res), num2str(tol)); end end \ No newline at end of file From 855887b249d7c52eb8c9791f2a49fbdfd741a68b Mon Sep 17 00:00:00 2001 From: NormannR Date: Tue, 4 Oct 2022 22:34:07 +0200 Subject: [PATCH 2/3] Implements a logarithmic reduction Fortran routine and the associated test --- .../logarithmic_reduction.m | 0 mex/build/cycle_reduction.am | 7 +- mex/build/logarithmic_reduction.am | 19 ++ mex/build/matlab/Makefile.am | 2 +- mex/build/matlab/configure.ac | 4 +- .../matlab/logarithmic_reduction/Makefile.am | 2 + mex/build/octave/Makefile.am | 2 +- mex/build/octave/configure.ac | 3 +- .../octave/logarithmic_reduction/Makefile.am | 3 + mex/sources/Makefile.am | 3 +- mex/sources/blas_lapack.F08 | 87 ++++++- mex/sources/cycle_reduction/mexFunction.f08 | 59 +---- .../logarithmic_reduction/mexFunction.f08 | 213 ++++++++++++++++++ tests/Makefile.am | 8 +- tests/logarithmicreduction.m | 140 ++++++++++++ 15 files changed, 492 insertions(+), 60 deletions(-) rename matlab/{ => missing/mex/logarithmic_reduction}/logarithmic_reduction.m (100%) create mode 100644 mex/build/logarithmic_reduction.am create mode 100644 mex/build/matlab/logarithmic_reduction/Makefile.am create mode 100644 mex/build/octave/logarithmic_reduction/Makefile.am create mode 100644 mex/sources/logarithmic_reduction/mexFunction.f08 create mode 100644 tests/logarithmicreduction.m 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 Date: Tue, 4 Oct 2022 22:37:40 +0200 Subject: [PATCH 3/3] Fix cycle reduction: (i) making the norms consistent between cycle_reduction and its test; (ii) remove hard errors in cycle_reduction Fortran and Matlab routines --- matlab/missing/mex/cycle_reduction/cycle_reduction.m | 2 +- mex/sources/cycle_reduction/mexFunction.f08 | 2 +- tests/cyclereduction.m | 10 +++++----- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/matlab/missing/mex/cycle_reduction/cycle_reduction.m b/matlab/missing/mex/cycle_reduction/cycle_reduction.m index 18ba3637f..8f33c04c4 100644 --- a/matlab/missing/mex/cycle_reduction/cycle_reduction.m +++ b/matlab/missing/mex/cycle_reduction/cycle_reduction.m @@ -108,7 +108,7 @@ if (nargin == 5 && ~isempty(ch) ) if (res > cvg_tol) info(1) = 403 info(2) = log(res) - error(['The norm of the residual is ' num2str(res) ' whereas the tolerance criterion is ' num2str(cvg_tol)]); + dprintf('The norm of the residual is %s whereas the tolerance criterion is %s', num2str(res), num2str(cvg_tol)); end end diff --git a/mex/sources/cycle_reduction/mexFunction.f08 b/mex/sources/cycle_reduction/mexFunction.f08 index 1f9238fdd..cf39944f1 100644 --- a/mex/sources/cycle_reduction/mexFunction.f08 +++ b/mex/sources/cycle_reduction/mexFunction.f08 @@ -99,7 +99,7 @@ loop: do info(2) = log(residual) write (cvg_tol_str,"(es8.2)") cvg_tol write (residual_str,"(es8.2)") residual - call mexErrMsgTxt("The norm of the residual is "& + call mexPrintf("The norm of the residual is "& &// trim(residual_str) // & &", whereas the tolerance criterion is " & // trim(cvg_tol_str) // "." ) diff --git a/tests/cyclereduction.m b/tests/cyclereduction.m index 211e591cc..c8c09db68 100644 --- a/tests/cyclereduction.m +++ b/tests/cyclereduction.m @@ -54,8 +54,8 @@ tElapsed1 = 0.; try tic; [X1,info] = cycle_reduction_matlab(C,B,A,cvg_tol,[0.]); tElapsed1 = toc; disp(['Elapsed time for the Matlab cycle reduction algorithm is: ' num2str(tElapsed1) ' (n=' int2str(n) ').']) - R = C+B*X1+A*X1*X1; - if (max(abs(R(:))) > cvg_tol) + R = norm(C+B*X1+A*X1*X1,1); + if (R > cvg_tol) testFailed = testFailed+1; if debug dprintf('Matlab cycle_reduction solution is wrong') @@ -74,8 +74,8 @@ tElapsed2 = 0.; try tic; [X2,info] = cycle_reduction(C,B,A,cvg_tol,[0.]); tElapsed2 = toc; disp(['Elapsed time for the Fortran cycle reduction algorithm is: ' num2str(tElapsed2) ' (n=' int2str(n) ').']) - R = C+B*X2+A*X2*X2; - if (max(abs(R(:))) > cvg_tol) + R = norm(C+B*X2+A*X2*X2,1); + if (R > cvg_tol) testFailed = testFailed+1; if debug dprintf('Fortran cycle_reduction solution is wrong') @@ -90,7 +90,7 @@ end % 3. Compare solutions of the Fortran and Matlab routines NumberOfTests = NumberOfTests+1; -if (max(abs(X1(:) - X2(:))) > cvg_tol) +if (norm(X1 - X2, 1) > cvg_tol) testFailed = testFailed+1; if debug dprintf('Fortran and Matlab cycle reduction solutions differ');