dmm: switch to mixed c/fortran for octave compilation

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

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

mex/sources/dmm/designInternal.cc

+#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]));
+}

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

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