From 079a6483618589a8caa9f2b49ee95b0e4ed3692c Mon Sep 17 00:00:00 2001
From: Houtan Bastani <houtan@dynare.org>
Date: Fri, 14 Nov 2014 11:16:37 +0100
Subject: [PATCH] dmm: switch to mixed c/fortran for octave compilation
---
mex/build/dmm.am | 1 +
mex/sources/dmm/design.f90 | 68 +++++++++++--------------------
mex/sources/dmm/designInternal.cc | 46 +++++++++++++++++++++
mex/sources/dmm/mexinterface.f90 | 26 ++++++++----
mex/sources/dmm/mexprint.f90 | 9 ++--
5 files changed, 91 insertions(+), 59 deletions(-)
create mode 100644 mex/sources/dmm/designInternal.cc
diff --git a/mex/build/dmm.am b/mex/build/dmm.am
index 0838e408f1..82a9cac7c1 100644
--- a/mex/build/dmm.am
+++ b/mex/build/dmm.am
@@ -114,6 +114,7 @@ nodist_dmm_SOURCES = \
$(DMMDIR)/randlib/spofa.f \
$(DMMMEXSRCDIR)/design.f90 \
$(DMMMEXSRCDIR)/mexprint.f90 \
+ $(DMMMEXSRCDIR)/designInternal.cc \
$(DMMMEXSRCDIR)/dmm.cc
CLEANFILES = *.mod
diff --git a/mex/sources/dmm/design.f90 b/mex/sources/dmm/design.f90
index 8b42659a64..95732d2854 100644
--- a/mex/sources/dmm/design.f90
+++ b/mex/sources/dmm/design.f90
@@ -20,66 +20,44 @@
#if defined(MATLAB_MEX_FILE)
#include "fintrf.h"
#endif
+
+#ifndef MWPOINTER
+#define MWPOINTER integer(4)
+#endif
+
+#ifndef mwPointer
+#define mwPointer MWPOINTER
+#endif
+
SUBROUTINE DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
USE MEXINTERFACE
IMPLICIT NONE
+
! INPUT
- INTEGER ny,nz,nx,nu,ns(6),nt
- DOUBLE PRECISION theta(nt), nsd(6)
+ INTEGER(C_INT), INTENT(IN) :: ny,nz,nx,nu,nt
+ INTEGER, DIMENSION(6), INTENT(IN) :: ns
+ DOUBLE PRECISION, DIMENSION(nt), INTENT(IN) :: theta
! OUTPUT
- DOUBLE PRECISION c(ny,max(1,nz),ns(1)),H(ny,nx,ns(2))
- DOUBLE PRECISION G(ny,nu,ns(3)),a(nx,ns(4))
- DOUBLE PRECISION F(nx,nx,ns(5)),R(nx,nu,ns(6))
+ DOUBLE PRECISION, DIMENSION(ny,max(1,nz),ns(1)), INTENT(OUT) :: c
+ DOUBLE PRECISION, DIMENSION(ny,nx,ns(2)), INTENT(OUT) :: H
+ DOUBLE PRECISION, DIMENSION(ny,nu,ns(3)), INTENT(OUT) :: G
+ DOUBLE PRECISION, DIMENSION(nx,ns(4)), INTENT(OUT) :: a
+ DOUBLE PRECISION, DIMENSION(nx,nx,ns(5)), INTENT(OUT) :: F
+ DOUBLE PRECISION, DIMENSION(nx,nu,ns(6)), INTENT(OUT) :: R
! COMMON
CHARACTER*200 mfile
COMMON /M/ mfile
! LOCAL
- mwPointer INPUT(6), OUTPUT(6)
- INTEGER STATUS, I
- CHARACTER(len=200) :: toprint
- integer*4, PARAMETER :: mxREAL = 0
-
- ! Matlab mex/mx functions
- INTEGER mexCallMATLAB
- mwPointer mxGetPr, mxCreateDoubleScalar, mxCreateDoubleMatrix
- EXTERNAL mxGetPr, mxCreateDoubleScalar, mxCreateDoubleMatrix, mxCopyPtrToReal8, mexCallMATLAB
+ INTEGER I
+ INTEGER(C_INT), DIMENSION(6) :: nsC
DO I=1,6
- nsd(I) = ns(I)
+ nsC(I) = ns(I)
END DO
+ CALL designInternal(ny,nz,nx,nu,nsC,nt,theta,TRIM(mfile)//C_NULL_CHAR,c,H,G,a,F,R)
- ! Set input values
- INPUT(1) = mxCreateDoubleScalar(ny*1.d0)
- INPUT(2) = mxCreateDoubleScalar(nz*1.d0)
- INPUT(3) = mxCreateDoubleScalar(nx*1.d0)
- INPUT(4) = mxCreateDoubleScalar(nu*1.d0)
- INPUT(5) = mxCreateDoubleMatrix(1, 6, mxREAL)
- CALL mxCopyReal8ToPtr(nsd, mxGetPr(INPUT(5)), 6)
- INPUT(6) = mxCreateDoubleMatrix(1, nt, mxREAL)
- CALL mxCopyReal8ToPtr(theta, mxGetPr(INPUT(6)), nt)
-
- ! Call design .m function
- STATUS = mexCallMATLAB(6, OUTPUT, 6, INPUT, mfile)
- IF (STATUS .ne. 0) THEN
- WRITE(toprint, '(A,A,A,I2,A)') '\nCall to ',mfile,' failed with code ',STATUS,'\n'
- CALL mexErrMsgTxt(toprint)
- ENDIF
-
- ! Copy Matlab output into Fortran
- CALL mxCopyPtrToReal8(mxGetPr(OUTPUT(1)), c, ny*max(1,nz)*ns(1))
- CALL mxCopyPtrToReal8(mxGetPr(OUTPUT(2)), H, ny*nx*ns(2))
- CALL mxCopyPtrToReal8(mxGetPr(OUTPUT(3)), G, ny*nu*ns(3))
- CALL mxCopyPtrToReal8(mxGetPr(OUTPUT(4)), a, nx*ns(4))
- CALL mxCopyPtrToReal8(mxGetPr(OUTPUT(5)), F, nx*nx*ns(5))
- CALL mxCopyPtrToReal8(mxGetPr(OUTPUT(6)), R, nx*nu*ns(6))
-
- ! Free memory
- DO I = 1,6
- CALL mxDestroyArray(INPUT(I))
- CALL mxDestroyArray(OUTPUT(I))
- END DO
END SUBROUTINE DESIGN
#endif
diff --git a/mex/sources/dmm/designInternal.cc b/mex/sources/dmm/designInternal.cc
new file mode 100644
index 0000000000..e47b8fa0f2
--- /dev/null
+++ b/mex/sources/dmm/designInternal.cc
@@ -0,0 +1,46 @@
+#include <dynmex.h>
+#include <algorithm>
+using namespace std;
+
+extern "C"
+void designInternal(int ny, int nz, int nx, int nu, int ns[6], int nt,
+ double *theta, char *mfile,
+ double *c, double *H, double *G, double *a, double *F, double *R)
+{
+ mxArray *lhs[6], *rhs[6];
+ nz = max(1, nz);
+ const mwSize dimsC[] = {ny, nz, ns[0]};
+ const mwSize dimsH[] = {ny, nx, ns[1]};
+ const mwSize dimsG[] = {ny, nu, ns[2]};
+ const mwSize dimsA[] = {nx, ns[3]};
+ const mwSize dimsF[] = {nx, nx, ns[4]};
+ const mwSize dimsR[] = {nx, nu, ns[5]};
+
+ lhs[0] = mxCreateNumericArray(3, dimsC, mxDOUBLE_CLASS, mxREAL);
+ lhs[1] = mxCreateNumericArray(3, dimsH, mxDOUBLE_CLASS, mxREAL);
+ lhs[2] = mxCreateNumericArray(3, dimsG, mxDOUBLE_CLASS, mxREAL);
+ lhs[3] = mxCreateNumericArray(2, dimsA, mxDOUBLE_CLASS, mxREAL);
+ lhs[4] = mxCreateNumericArray(3, dimsF, mxDOUBLE_CLASS, mxREAL);
+ lhs[5] = mxCreateNumericArray(3, dimsR, mxDOUBLE_CLASS, mxREAL);
+
+ rhs[0] = mxCreateDoubleScalar(ny);
+ rhs[1] = mxCreateDoubleScalar(nz);
+ rhs[2] = mxCreateDoubleScalar(nx);
+ rhs[3] = mxCreateDoubleScalar(nu);
+ rhs[4] = mxCreateDoubleMatrix(1, 6, mxREAL);
+ rhs[5] = mxCreateDoubleMatrix(1, nt, mxREAL);
+
+ double *passns = mxGetPr(rhs[4]);
+ for (int i = 0; i < 6; i++)
+ passns[i] = (double) ns[i];
+ memcpy(mxGetPr(rhs[5]), theta, sizeof(double) * nt);
+
+ mexCallMATLAB(6, lhs, 6, rhs, mfile);
+
+ memcpy(c, mxGetPr(lhs[0]), sizeof(double) * mxGetNumberOfElements(lhs[0]));
+ memcpy(H, mxGetPr(lhs[1]), sizeof(double) * mxGetNumberOfElements(lhs[1]));
+ memcpy(G, mxGetPr(lhs[2]), sizeof(double) * mxGetNumberOfElements(lhs[2]));
+ memcpy(a, mxGetPr(lhs[3]), sizeof(double) * mxGetNumberOfElements(lhs[3]));
+ memcpy(F, mxGetPr(lhs[4]), sizeof(double) * mxGetNumberOfElements(lhs[4]));
+ memcpy(R, mxGetPr(lhs[5]), sizeof(double) * mxGetNumberOfElements(lhs[5]));
+}
diff --git a/mex/sources/dmm/mexinterface.f90 b/mex/sources/dmm/mexinterface.f90
index 365b80dbda..08cb7e92c0 100644
--- a/mex/sources/dmm/mexinterface.f90
+++ b/mex/sources/dmm/mexinterface.f90
@@ -23,17 +23,25 @@ MODULE MEXINTERFACE
INTERFACE
- INTEGER(4) FUNCTION mexPrintf(toprint) BIND(C, NAME="mexPrintf")
- USE ISO_C_BINDING, ONLY : C_CHAR
- IMPLICIT NONE
- CHARACTER(KIND=C_CHAR), INTENT(IN) :: toprint(*)
- END FUNCTION mexPrintf
+ SUBROUTINE mexPrintf(string) BIND(C, NAME="mexPrintf")
+ USE ISO_C_BINDING, ONLY: C_CHAR
+ CHARACTER(KIND=C_CHAR) :: string(*)
+ END SUBROUTINE mexPrintf
- SUBROUTINE mexErrMsgTxt(toprint) BIND(C, NAME="mexErrMsgTxt")
- USE ISO_C_BINDING, ONLY : C_CHAR
- IMPLICIT NONE
- CHARACTER(KIND=C_CHAR), INTENT(IN) :: toprint(*)
+ SUBROUTINE mexErrMsgTxt(string) BIND(C, NAME="mexErrMsgTxt")
+ USE ISO_C_BINDING, ONLY: C_CHAR
+ CHARACTER(KIND=C_CHAR) :: string(*)
END SUBROUTINE mexErrMsgTxt
+ SUBROUTINE designInternal(ny,nz,nx,nu,ns,nt,theta,mfile,c,H,G,a,F,R) BIND(C, NAME="designInternal")
+ USE, INTRINSIC :: ISO_C_BINDING
+ IMPLICIT NONE
+ INTEGER(C_INT), VALUE, INTENT(IN) :: ny, nz, nx, nu, nt
+ INTEGER(C_INT), DIMENSION(*), INTENT(IN) :: ns
+ REAL(C_DOUBLE), DIMENSION(*), INTENT(IN) :: theta
+ CHARACTER(KIND=C_CHAR), DIMENSION(*), INTENT(IN) :: mfile
+ REAL(C_DOUBLE), DIMENSION(*), INTENT(OUT) :: c, H, G, a, F, R
+ END SUBROUTINE designInternal
+
END INTERFACE
END MODULE MEXINTERFACE
diff --git a/mex/sources/dmm/mexprint.f90 b/mex/sources/dmm/mexprint.f90
index 869cd3f558..03e920694f 100644
--- a/mex/sources/dmm/mexprint.f90
+++ b/mex/sources/dmm/mexprint.f90
@@ -17,13 +17,12 @@
! along with Dynare. If not, see <http://www.gnu.org/licenses/>.
!
-INTEGER(4) FUNCTION mexPrint(toprint)
+SUBROUTINE mexPrint(string)
USE MEXINTERFACE
- USE ISO_C_BINDING
IMPLICIT NONE
- CHARACTER(len=200), INTENT(IN) :: toprint
+ CHARACTER(len=200), INTENT(IN) :: string
- mexPrint = mexPrintf(TRIM(toprint)//NEW_LINE('A')//c_null_char)
+ CALL mexPrintf(TRIM(string)//NEW_LINE('A')//c_null_char)
+END SUBROUTINE mexPrint
-END FUNCTION mexPrint
--
GitLab