From efd64ab7e7bcd6c020cc53272208cbbe70721579 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?=
 <stephane.adjemian@univ-lemans.fr>
Date: Wed, 8 Feb 2012 13:52:23 +0100
Subject: [PATCH] Added mex file version for getPowerDeriv.m routine.

With this version of the mex file, the overhead cost of mex call is too important.
---
 matlab/dynare_config.m                       | 13 +++
 matlab/{ => getPowerDeriv}/getPowerDeriv.m   |  2 +-
 mex/build/getPowerDeriv.am                   |  5 ++
 mex/build/matlab/Makefile.am                 |  2 +-
 mex/build/matlab/configure.ac                |  3 +-
 mex/build/matlab/get_power_deriv/Makefile.am |  2 +
 mex/build/octave/Makefile.am                 |  2 +-
 mex/build/octave/configure.ac                |  3 +-
 mex/build/octave/get_power_deriv/Makefile.am |  3 +
 mex/sources/Makefile.am                      |  3 +-
 mex/sources/get_power_deriv/getPowerDeriv.cc | 87 ++++++++++++++++++++
 11 files changed, 119 insertions(+), 6 deletions(-)
 rename matlab/{ => getPowerDeriv}/getPowerDeriv.m (99%)
 create mode 100644 mex/build/getPowerDeriv.am
 create mode 100644 mex/build/matlab/get_power_deriv/Makefile.am
 create mode 100644 mex/build/octave/get_power_deriv/Makefile.am
 create mode 100644 mex/sources/get_power_deriv/getPowerDeriv.cc

diff --git a/matlab/dynare_config.m b/matlab/dynare_config.m
index 97b18745da..c38c254707 100644
--- a/matlab/dynare_config.m
+++ b/matlab/dynare_config.m
@@ -267,4 +267,17 @@ if verbose
     disp(' ')
 end
 
+% Test if getPowerDeriv DLL is present
+if exist('getPowerDeriv', 'file') == 3
+    message = '[mex] ';
+else
+    message = '[no]  ';
+    addpath([dynareroot 'getPowerDeriv']);
+end
+if verbose
+    disp([ message 'getPowerDeriv routine.' ])
+    disp(' ')
+end
+
+
 cd(origin);
diff --git a/matlab/getPowerDeriv.m b/matlab/getPowerDeriv/getPowerDeriv.m
similarity index 99%
rename from matlab/getPowerDeriv.m
rename to matlab/getPowerDeriv/getPowerDeriv.m
index e87eabe1e0..ab6e6b6d54 100644
--- a/matlab/getPowerDeriv.m
+++ b/matlab/getPowerDeriv/getPowerDeriv.m
@@ -39,4 +39,4 @@ else
         p = p-1;
     end
 end
-end
+end
\ No newline at end of file
diff --git a/mex/build/getPowerDeriv.am b/mex/build/getPowerDeriv.am
new file mode 100644
index 0000000000..4f7f4fd832
--- /dev/null
+++ b/mex/build/getPowerDeriv.am
@@ -0,0 +1,5 @@
+vpath %.cc $(top_srcdir)/../../sources/get_power_deriv
+
+noinst_PROGRAMS = getPowerDeriv
+
+nodist_getPowerDeriv_SOURCES = getPowerDeriv.cc
diff --git a/mex/build/matlab/Makefile.am b/mex/build/matlab/Makefile.am
index 58dbd369f8..8038a7ed9a 100644
--- a/mex/build/matlab/Makefile.am
+++ b/mex/build/matlab/Makefile.am
@@ -2,7 +2,7 @@ ACLOCAL_AMFLAGS = -I ../../../m4
 
 # libdynare++ must come before gensylv, k_order_perturbation, dynare_simul_
 if DO_SOMETHING
-SUBDIRS = mjdgges kronecker bytecode libdynare++ gensylv k_order_perturbation dynare_simul_ estimation block_kalman_filter sobol
+SUBDIRS = mjdgges kronecker bytecode libdynare++ gensylv k_order_perturbation dynare_simul_ estimation block_kalman_filter sobol get_power_deriv
 
 if HAVE_GSL
 SUBDIRS += ms_sbvar
diff --git a/mex/build/matlab/configure.ac b/mex/build/matlab/configure.ac
index 4267a3158c..ab97d1970f 100644
--- a/mex/build/matlab/configure.ac
+++ b/mex/build/matlab/configure.ac
@@ -144,6 +144,7 @@ AC_CONFIG_FILES([Makefile
                  kalman_steady_state/Makefile
                  ms_sbvar/Makefile
                  block_kalman_filter/Makefile
-	         sobol/Makefile])
+	         sobol/Makefile
+		 get_power_deriv/Makefile])
 
 AC_OUTPUT
diff --git a/mex/build/matlab/get_power_deriv/Makefile.am b/mex/build/matlab/get_power_deriv/Makefile.am
new file mode 100644
index 0000000000..e04608481a
--- /dev/null
+++ b/mex/build/matlab/get_power_deriv/Makefile.am
@@ -0,0 +1,2 @@
+include ../mex.am
+include ../../getPowerDeriv.am
diff --git a/mex/build/octave/Makefile.am b/mex/build/octave/Makefile.am
index a6d164c901..e11aa0423a 100644
--- a/mex/build/octave/Makefile.am
+++ b/mex/build/octave/Makefile.am
@@ -2,7 +2,7 @@ ACLOCAL_AMFLAGS = -I ../../../m4
 
 # libdynare++ must come before gensylv, k_order_perturbation, dynare_simul_
 if DO_SOMETHING
-SUBDIRS = mjdgges kronecker bytecode libdynare++ gensylv k_order_perturbation dynare_simul_ qzcomplex ordschur block_kalman_filter sobol
+SUBDIRS = mjdgges kronecker bytecode libdynare++ gensylv k_order_perturbation dynare_simul_ qzcomplex ordschur block_kalman_filter sobol get_power_deriv
 
 if HAVE_GSL
 SUBDIRS += ms_sbvar
diff --git a/mex/build/octave/configure.ac b/mex/build/octave/configure.ac
index daa13a6af0..15a92f5ec0 100644
--- a/mex/build/octave/configure.ac
+++ b/mex/build/octave/configure.ac
@@ -131,6 +131,7 @@ AC_CONFIG_FILES([Makefile
                  kalman_steady_state/Makefile
                  ms_sbvar/Makefile
                  block_kalman_filter/Makefile
-		 sobol/Makefile])
+		 sobol/Makefile
+		 get_power_deriv/Makefile])
 
 AC_OUTPUT
diff --git a/mex/build/octave/get_power_deriv/Makefile.am b/mex/build/octave/get_power_deriv/Makefile.am
new file mode 100644
index 0000000000..1eaa678234
--- /dev/null
+++ b/mex/build/octave/get_power_deriv/Makefile.am
@@ -0,0 +1,3 @@
+EXEEXT = .mex
+include ../mex.am
+include ../../getPowerDeriv.am
diff --git a/mex/sources/Makefile.am b/mex/sources/Makefile.am
index 7444143f90..963a6e0e1d 100644
--- a/mex/sources/Makefile.am
+++ b/mex/sources/Makefile.am
@@ -14,7 +14,8 @@ EXTRA_DIST = \
 	kalman_steady_state \
 	ms-sbvar \
 	block_kalman_filter \
-	sobol
+	sobol \
+	get_power_deriv
 
 clean-local:
 	rm -rf `find mex/sources -name *.o`
diff --git a/mex/sources/get_power_deriv/getPowerDeriv.cc b/mex/sources/get_power_deriv/getPowerDeriv.cc
new file mode 100644
index 0000000000..a1634b46a4
--- /dev/null
+++ b/mex/sources/get_power_deriv/getPowerDeriv.cc
@@ -0,0 +1,87 @@
+/*
+** mex file for matlab getPowerDeriv.m.
+**
+** Copyright (C) 2012 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 <http://www.gnu.org/licenses/>.
+**
+** AUTHOR(S): stephane DOT adjemian AT univ DASH lemans DOT fr
+**/
+
+
+#include <math.h>
+#include <dynmex.h>
+#include "mex.h"
+
+#ifdef _MSC_VER
+# define nearbyint(x) (fabs((x)-floor(x)) < fabs((x)-ceil(x)) ? floor(x) : ceil(x))
+#endif
+
+double getPowerDeriv(double *x, double *p, int k)
+{
+  if ( fabs(*x) < 1e-12 && *p > 0 && k >= *p && fabs(*p-nearbyint(*p)) < 1e-12 )
+    return 0.0;
+  else
+    {
+      //int i = 0;
+      double dxp = pow(*x, *p-k);
+      for (int i=0; i<k; i++)
+        dxp *= (*p)--;
+      return dxp;
+    }
+}
+
+void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{
+  /*
+  ** INPUTS:
+  ** prhs[0]    [double]    scalar, x.
+  ** prhs[1]    [double]    scalar, p.
+  ** prhs[2]    [integer]   scalar, k.
+  **
+  ** OUTPUTS:
+  ** plhs[0]    [double]    scalar, The k-th derivative of x^p.
+  ** plhs[1]    [double]    scalar, info variable (equal to zero if call to getPowerDeriv is successfull).
+  */
+  if ( !( nrhs==3)  )
+    {
+      DYN_MEX_FUNC_ERR_MSG_TXT("getPowerDeriv:: Three input arguments are required (x,p and k)!");
+    }
+  if (nlhs>2)
+    {
+      DYN_MEX_FUNC_ERR_MSG_TXT("getPowerDeriv:: I only return one or two output arguments!");
+    }
+  if ( !(mxGetN(prhs[0])==1 && (mxGetM(prhs[0]))==1) )
+    {
+      DYN_MEX_FUNC_ERR_MSG_TXT("getPowerDeriv:: First argument, x, must be a double scalar!");
+    }
+  if ( !(mxGetN(prhs[1])==1 && (mxGetM(prhs[1]))==1) )
+    {
+      DYN_MEX_FUNC_ERR_MSG_TXT("getPowerDeriv:: Second argument, p, must be a double scalar!");
+    }
+  if ( !(mxGetN(prhs[2])==1 && (mxGetM(prhs[2]))==1) )
+    {
+      DYN_MEX_FUNC_ERR_MSG_TXT("getPowerDeriv:: Third argument, k, must be an integer scalar!");
+    }
+  double *x = mxGetPr(prhs[0]);
+  double *p = mxGetPr(prhs[1]);
+  int k = (int) mxGetScalar(prhs[2]);
+  plhs[0] = mxCreateDoubleScalar(getPowerDeriv(x,p,k));
+  if (nlhs>1)
+    {
+      plhs[1] = mxCreateDoubleScalar(0.0);
+    }
+}
-- 
GitLab