diff --git a/others/cpp/tests/Makefile_petsc b/others/cpp/tests/Makefile_petsc
new file mode 100644
index 0000000000000000000000000000000000000000..87b829ce5ee0c327490615f186cf43950b352aba
--- /dev/null
+++ b/others/cpp/tests/Makefile_petsc
@@ -0,0 +1,25 @@
+include ${PETSC_DIR}/lib/petsc/conf/variables
+include ${PETSC_DIR}/lib/petsc/conf/rules
+
+CFLAGS = ${PETSC_CC_INCLUDES} -I. -I..
+CPPFLAGS = ${PETSC_CC_INCLUDES} -I. -I..
+
+all: main_tao_1
+
+main_tao_1: main_tao_1.o pf_jacobian.o example1.o example1_steadystate.o example1_first_derivatives.o example1_residuals.o
+ ${CLINKER} -o main_tao_1 main_tao_1.o example1.o example1_steadystate.o example1_residuals.o example1_first_derivatives.o ${PETSC_LIB} -lstdc++ -lm
+
+snes_1: snes_1.o example1.o example1_steadystate.o example1_first_derivatives.o example1_residuals.o
+ ${CLINKER} -o snes_1 snes_1.o example1.o example1_steadystate.o example1_residuals.o example1_first_derivatives.o ${PETSC_LIB} -lstdc++ -lm
+
+main_tao_1.o: main_tao_1.cc main_tao.h ../dynare_cpp_driver.hh
+snes_1.o: snes_1.cc main_tao.h ../dynare_cpp_driver.hh
+pf_jacobian.o: pf_jacobian.cc main_tao.h
+example1.o: example1.cc ../dynare_cpp_driver.hh
+example1_steadystate.o: example1_steadystate.cc
+example1_residuals.o: example1_residuals.c
+example1_first_derivatives.o: example1_first_derivatives.c
+example1.cc: example1.mod
+ /home/michel/dynare/git/master/matlab/preprocessor64/dynare_m example1.mod output=first language=C++
+example1_residuals.c: example1.mod
+ /home/michel/dynare/git/master/matlab/preprocessor64/dynare_m example1.mod output=first language=C++
diff --git a/others/cpp/tests/example1.mod b/others/cpp/tests/example1.mod
index 1249be4ba08f8acca8905dfd3e1ba7e1a87b7a38..ec43f08718f5b1971a15af3a13cc7688e128a942 100644
--- a/others/cpp/tests/example1.mod
+++ b/others/cpp/tests/example1.mod
@@ -1,5 +1,5 @@
// Example 1 from Collard's guide to Dynare
-var y, c, k, a, h, b;
+var h, c, y, k, a, b;
varexo e, u;
parameters beta, rho, alpha, delta, theta, psi, tau;
@@ -15,11 +15,11 @@ theta = 2.95;
phi = 0.1;
model;
-c*theta*h^(1+psi)=(1-alpha)*y;
-k = beta*(((exp(b)*c)/(exp(b(+1))*c(+1)))
- *(exp(b(+1))*alpha*y(+1)+(1-delta)*k));
-y = exp(a)*(k(-1)^alpha)*(h^(1-alpha));
-k = exp(b)*(y-c)+(1-delta)*k(-1);
+c*theta*exp(h)^(1+psi)=(1-alpha)*y;
+exp(k) = beta*(((exp(b)*c)/(exp(b(+1))*c(+1)))
+ *(exp(b(+1))*alpha*y(+1)+(1-delta)*exp(k)));
+y = exp(a)*(exp(k(-1))^alpha)*(exp(h)^(1-alpha));
+exp(k) = exp(b)*(y-c)+(1-delta)*exp(k(-1));
a = rho*a(-1)+tau*b(-1) + e;
b = tau*a(-1)+rho*b(-1) + u;
end;
@@ -30,10 +30,10 @@ b = 0;
k_h = ((1-beta*(1-delta))/(beta*alpha))^(1/(alpha-1));
y_h = k_h^alpha;
c_h = y_h - delta*k_h;
-h = (y_h*(1-alpha)/(c_h*theta))^(1/(1+psi));
-k = k_h*h;
-c = c_h*h;
-y = y_h*h;
+h = log((y_h*(1-alpha)/(c_h*theta))^(1/(1+psi)));
+k = log(k_h*exp(h));
+c = c_h*exp(h);
+y = y_h*exp(h);
end;
shocks;
diff --git a/others/cpp/tests/snes_1.cc b/others/cpp/tests/snes_1.cc
new file mode 100644
index 0000000000000000000000000000000000000000..a3ffb0c0bb495457d7d961756c429fc3ef727ff5
--- /dev/null
+++ b/others/cpp/tests/snes_1.cc
@@ -0,0 +1,685 @@
+#include <iostream>
+#include "../dynare_cpp_driver.hh"
+#include <petscdm.h>
+#include <petscdmda.h>
+#include <petscsnes.h>
+#include "main_tao.h"
+
+DynareInfo *preprocessorOutput(void);
+PetscErrorCode FormFunction(SNES snes,Vec y,Vec f,void *ctx);
+PetscErrorCode FormJacobian(SNES snes,Vec y,Mat J, Mat B,void *ctx);
+PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void*);
+PetscErrorCode PreCheck(SNESLineSearch,Vec,Vec,PetscBool*,void*);
+PetscErrorCode PostCheck(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);
+PetscErrorCode PostSetSubKSP(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);
+PetscErrorCode MatrixFreePreconditioner(PC,Vec,Vec);
+
+static char help[] = "Testing";
+
+/*
+ User-defined context for monitoring
+*/
+typedef struct {
+ PetscViewer viewer;
+} MonitorCtx;
+
+/*
+ User-defined context for checking candidate iterates that are
+ determined by line search methods
+*/
+typedef struct {
+ Vec last_step; /* previous iterate */
+ PetscReal tolerance; /* tolerance for changes between successive iterates */
+ AppCtx *user;
+} StepCheckCtx;
+
+typedef struct {
+ PetscInt its0; /* num of prevous outer KSP iterations */
+} SetSubKSPCtx;
+
+#undef __FUNCT__
+#define __FUNCT__ "main"
+int main(int argc, char **argv)
+{
+ DynareInfo model_info;
+
+ int endo_nbr = model_info.get_endo_nbr();
+ int exo_nbr = model_info.get_exo_nbr();
+ double *params = model_info.get_params_data();
+ // Steady state
+ double *steady_state = new double[endo_nbr];
+ int info;
+ steadystate(NULL,params, steady_state, &info);
+ vector<int> NNZDerivatives = model_info.get_NNZDerivatives();
+ AppCtx user; /* user-defined work context */
+ user.exogenous = new double[exo_nbr];
+ user.params = params;
+ user.steady_state = steady_state;
+ user.periods = 40;
+ user.first_col = endo_nbr;
+ user.endo_nbr = endo_nbr;
+ user.exo_nbr = exo_nbr;
+ user.row_ptr = new int[endo_nbr+1];
+ user.nnz = NNZDerivatives[0];
+ user.col_ptr = new int[NNZDerivatives[0]];
+ user.val_ptr = new double[NNZDerivatives[0]];
+ user.initial_values = new double[user.endo_nbr];
+ user.terminal_values = new double[user.endo_nbr];
+ for (int i=0; i < user.endo_nbr; ++i)
+ {
+ user.initial_values[i] = user.steady_state[i];
+ user.terminal_values[i] = user.steady_state[i];
+ }
+ /* Initialize PETSc */
+ PetscInitialize(&argc, &argv, (char *)0, help);
+ PetscErrorCode ierr;
+ /* Get number of processors */
+ PetscInt size;
+ ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
+ /* Set periods a multiple of processor nbr */
+ user.periods = size*ceil((double)user.periods/size);
+ user.nb_row_x = user.periods + 1;
+ user.X = new double[user.nb_row_x*user.exo_nbr];
+ for(double* px=user.X; px < user.X+user.nb_row_x*user.exo_nbr; px++) *px = 0.0;
+ user.X[1] = 0.01;
+ user.X[1+user.nb_row_x] = 0.01;
+ std::cout << size << " " << user.periods << " " << std::endl;
+ PetscInt N = user.periods*user.endo_nbr;
+ /* Create DMA */
+ DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_GHOSTED,N,1,user.endo_nbr,NULL,&user.da);
+
+ /* Allocate vector and Jacobian matrix */\
+ Vec Y, R;
+ DMCreateGlobalVector(user.da,&Y);
+ VecDuplicate(Y,&R);
+
+ Mat J ;
+ ierr = MatCreate(PETSC_COMM_WORLD,&J);CHKERRQ(ierr);
+ ierr = MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,N,N);CHKERRQ(ierr);
+ ierr = MatSetFromOptions(J);CHKERRQ(ierr);
+ ierr = MatSetUp(J);CHKERRQ(ierr);
+
+ /*
+ Get local grid boundaries (for 1-dimensional DMDA):
+ xs, xm - starting grid index, width of local grid (no ghost points)
+ */
+ PetscInt xs, xm, xsg, xmg;
+ DMDAGetCorners(user.da,&xs,NULL,NULL,&xm,NULL,NULL);
+ std::cout << xs << " " << xm << std::endl;
+ DMDAGetGhostCorners(user.da,&xsg,NULL,NULL,&xmg,NULL,NULL);
+ std::cout << "ghost " << xsg << " " << xmg << std::endl;
+ /*
+ Get pointers to vector data
+ */
+ PetscScalar *YY;
+ DMDAVecGetArray(user.da,Y,&YY);
+
+ /*
+ Compute local vector entries
+ */
+ for (int i=xs; i<xs+xm; i += user.endo_nbr)
+ for (int j=0; j < user.endo_nbr; j++)
+ YY[i+j] = steady_state[j];
+ /*
+ Restore vectors
+ */
+ DMDAVecRestoreArray(user.da,Y,&YY);
+
+ SNES snes;
+ ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr);
+
+ // SNES snes;
+ SNESLineSearch linesearch; /* SNESLineSearch context */
+ MonitorCtx monP; /* monitoring context */
+ StepCheckCtx checkP; /* step-checking context */
+ SetSubKSPCtx checkP1;
+ PetscBool pre_check,post_check,post_setsubksp; /* flag indicating whether we're checking candidate iterates */
+ PetscReal abstol,rtol,stol,norm;
+ PetscInt its,maxit,maxf;
+ PetscBool flg;
+
+ /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ Create nonlinear solver context
+ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr);
+
+ /*
+ Set function evaluation routine and vector. Whenever the nonlinear
+ solver needs to compute the nonlinear function, it will call this
+ routine.
+ - Note that the final routine argument is the user-defined
+ context that provides application-specific data for the
+ function evaluation routine.
+ */
+ ierr = SNESSetFunction(snes,R,FormFunction,&user);CHKERRQ(ierr);
+
+ /*
+ Set Jacobian matrix data structure and default Jacobian evaluation
+ routine. Whenever the nonlinear solver needs to compute the
+ Jacobian matrix, it will call this routine.
+ - Note that the final routine argument is the user-defined
+ context that provides application-specific data for the
+ Jacobian evaluation routine.
+ */
+ ierr = SNESSetJacobian(snes,J,J,FormJacobian,&user);CHKERRQ(ierr);
+
+ /*
+ Optional allow user provided preconditioner
+ */
+ ierr = PetscOptionsHasName(NULL,"-user_precond",&flg);CHKERRQ(ierr);
+ if (flg) {
+ KSP ksp;
+ PC pc;
+ ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr);
+ ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
+ ierr = PCSetType(pc,PCSHELL);CHKERRQ(ierr);
+ ierr = PCShellSetApply(pc,MatrixFreePreconditioner);CHKERRQ(ierr);
+ }
+
+ /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ Customize nonlinear solver; set runtime options
+ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ /*
+ Set an optional user-defined monitoring routine
+ */
+ ierr = PetscViewerDrawOpen(PETSC_COMM_WORLD,0,0,0,0,400,400,&monP.viewer);CHKERRQ(ierr);
+ ierr = SNESMonitorSet(snes,Monitor,&monP,0);CHKERRQ(ierr);
+
+ /*
+ Set names for some vectors to facilitate monitoring (optional)
+ */
+ ierr = PetscObjectSetName((PetscObject)Y,"Approximate Solution");CHKERRQ(ierr);
+ // ierr = PetscObjectSetName((PetscObject)U,"Exact Solution");CHKERRQ(ierr);
+
+ /*
+ Set SNES/KSP/KSP/PC runtime options, e.g.,
+ -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
+ */
+ ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
+
+ /*
+ Set an optional user-defined routine to check the validity of candidate
+ iterates that are determined by line search methods
+ */
+ ierr = SNESGetLineSearch(snes, &linesearch);CHKERRQ(ierr);
+ ierr = PetscOptionsHasName(NULL,"-post_check_iterates",&post_check);CHKERRQ(ierr);
+
+ if (post_check) {
+ ierr = PetscPrintf(PETSC_COMM_WORLD,"Activating post step checking routine\n");CHKERRQ(ierr);
+ ierr = SNESLineSearchSetPostCheck(linesearch,PostCheck,&checkP);CHKERRQ(ierr);
+ ierr = VecDuplicate(Y,&(checkP.last_step));CHKERRQ(ierr);
+
+ checkP.tolerance = 1.0;
+ checkP.user = &user;
+
+ ierr = PetscOptionsGetReal(NULL,"-check_tol",&checkP.tolerance,NULL);CHKERRQ(ierr);
+ }
+
+ ierr = PetscOptionsHasName(NULL,"-post_setsubksp",&post_setsubksp);CHKERRQ(ierr);
+ if (post_setsubksp) {
+ ierr = PetscPrintf(PETSC_COMM_WORLD,"Activating post setsubksp\n");CHKERRQ(ierr);
+ ierr = SNESLineSearchSetPostCheck(linesearch,PostSetSubKSP,&checkP1);CHKERRQ(ierr);
+ }
+
+ ierr = PetscOptionsHasName(NULL,"-pre_check_iterates",&pre_check);CHKERRQ(ierr);
+ if (pre_check) {
+ ierr = PetscPrintf(PETSC_COMM_WORLD,"Activating pre step checking routine\n");CHKERRQ(ierr);
+ ierr = SNESLineSearchSetPreCheck(linesearch,PreCheck,&checkP);CHKERRQ(ierr);
+ }
+
+ /*
+ Print parameters used for convergence testing (optional) ... just
+ to demonstrate this routine; this information is also printed with
+ the option -snes_view
+ */
+ ierr = SNESGetTolerances(snes,&abstol,&rtol,&stol,&maxit,&maxf);CHKERRQ(ierr);
+ ierr = PetscPrintf(PETSC_COMM_WORLD,"atol=%g, rtol=%g, stol=%g, maxit=%D, maxf=%D\n",(double)abstol,(double)rtol,(double)stol,maxit,maxf);CHKERRQ(ierr);
+
+ /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ Evaluate initial guess; then solve nonlinear system
+ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ /*
+ Note: The user should initialize the vector, x, with the initial guess
+ for the nonlinear solver prior to calling SNESSolve(). In particular,
+ to employ an initial guess of zero, the user should explicitly set
+ this vector to zero by calling VecSet().
+ */
+ ierr = SNESSolve(snes,NULL,Y);CHKERRQ(ierr);
+ ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr);
+ ierr = PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);CHKERRQ(ierr);
+
+ /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ Check solution and clean up
+ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ /*
+ Free work space. All PETSc objects should be destroyed when they
+ are no longer needed.
+ */
+ ierr = PetscViewerDestroy(&monP.viewer);CHKERRQ(ierr);
+ if (post_check) {ierr = VecDestroy(&checkP.last_step);CHKERRQ(ierr);}
+ ierr = SNESDestroy(&snes);CHKERRQ(ierr);
+ ierr = DMDestroy(&user.da);CHKERRQ(ierr);
+
+
+ ierr = MatDestroy(&J);CHKERRQ(ierr);
+ ierr = VecDestroy(&Y);CHKERRQ(ierr);
+ ierr = VecDestroy(&R);CHKERRQ(ierr);
+ ierr = PetscFinalize();CHKERRQ(ierr);
+
+ PetscFunctionReturn(0);
+}
+
+PetscErrorCode FormFunction(SNES snes,Vec y,Vec f,void *ctx)
+{
+ AppCtx *user = (AppCtx*) ctx;
+ DM da = user->da;
+ PetscScalar *yy,*ff;
+ PetscInt M,ys,ym;
+ Vec ylocal;
+
+ DMGetLocalVector(da,&ylocal);
+ /*
+ Scatter ghost points to local vector, using the 2-step process
+ DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
+ By placing code between these two statements, computations can
+ be done while messages are in transition.
+ */
+ DMGlobalToLocalBegin(da,y,INSERT_VALUES,ylocal);
+ DMGlobalToLocalEnd(da,y,INSERT_VALUES,ylocal);
+
+ /*
+ Get pointers to vector data.
+ - The vector xlocal includes ghost point; the vectors x and f do
+ NOT include ghost points.
+ - Using DMDAVecGetArray() allows accessing the values using global ordering
+ */
+ DMDAVecGetArray(da,ylocal,&yy);
+ DMDAVecGetArray(da,f,&ff);
+
+ /*
+ Get local grid boundaries (for 1-dimensional DMDA):
+ ys, ym - starting grid index, width of local grid (no ghost points)
+ */
+ DMDAGetCorners(da,&ys,NULL,NULL,&ym,NULL,NULL);
+ DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
+
+ /*
+ Set function values for boundary points; define local interior grid point range:
+ xsi - starting interior grid index
+ xei - ending interior grid index
+ */
+ if (ys == 0) /* left boundary */
+ {
+ PetscReal *y1 = new PetscReal[3*user->endo_nbr];
+ for (int i=0; i < user->endo_nbr; ++i) y1[i] = user->initial_values[i];
+ for (int i=0; i < 2*user->endo_nbr; ++i) y1[i+user->endo_nbr] = yy[i];
+ Residuals(y1, user->X, user->nb_row_x, user->params, user->steady_state, 1, ff);
+ ys += user->endo_nbr;
+ ym -= user->endo_nbr;
+ }
+
+ /*
+ Compute function over locally owned part of the grid (interior points only)
+ */
+ while ( (ym >= user->endo_nbr) && (ys + 2*user->endo_nbr <= M) )
+ {
+ int it = ys/user->endo_nbr + 2;
+ Residuals(yy+ys-user->endo_nbr, user->X, user->nb_row_x, user->params, user->steady_state, it, ff+ys);
+ ys += user->endo_nbr;
+ ym -= user->endo_nbr;
+ }
+
+
+ if ( (ym >= user->endo_nbr) && (ys + 2*user->endo_nbr >= M) )
+ {
+ int it = ys/user->endo_nbr + 1;
+ PetscReal *y1 = new PetscReal[3*user->endo_nbr];
+ for (int i=0; i < 2*user->endo_nbr; ++i) y1[i] = yy[ys+i-user->endo_nbr];
+ for (int i=0; i < user->endo_nbr; ++i) y1[i+2*user->endo_nbr] = user->terminal_values[i];
+ Residuals(y1, user->X, user->nb_row_x, user->params, user->steady_state, it, ff+ys);
+ }
+
+ /*
+ Restore vectors
+ */
+ DMDAVecRestoreArray(da,ylocal,&yy);
+ DMDAVecRestoreArray(da,f,&ff);
+ DMRestoreLocalVector(da,&ylocal);
+ return(0);
+}
+
+PetscErrorCode FormJacobian(SNES snes,Vec y,Mat J, Mat B,void *ctx)
+{
+ AppCtx *user = (AppCtx*) ctx;
+ DM da = user->da;
+ PetscScalar *yy;
+ PetscInt M,ys,ym,ierr;
+ Vec ylocal;
+
+ DMGetLocalVector(da,&ylocal);
+ /*
+ Scatter ghost points to local vector, using the 2-step process
+ DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
+ By placing code between these two statements, computations can
+ be done while messages are in transition.
+ */
+ DMGlobalToLocalBegin(da,y,INSERT_VALUES,ylocal);
+ DMGlobalToLocalEnd(da,y,INSERT_VALUES,ylocal);
+
+ /*
+ Get pointers to vector data.
+ - The vector xlocal includes ghost point; the vectors x and f do
+ NOT include ghost points.
+ - Using DMDAVecGetArray() allows accessing the values using global ordering
+ */
+ DMDAVecGetArray(da,ylocal,&yy);
+
+ /*
+ Get local grid boundaries (for 1-dimensional DMDA):
+ ys, ym - starting grid index, width of local grid (no ghost points)
+ */
+ DMDAGetCorners(da,&ys,NULL,NULL,&ym,NULL,NULL);
+ DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
+
+ /*
+ Set function values for boundary points; define local interior grid point range:
+ xsi - starting interior grid index
+ xei - ending interior grid index
+ */
+ int row = 0;
+ if (ys == 0) /* left boundary */
+ {
+ PetscReal *y1 = new PetscReal[3*user->endo_nbr];
+ for (int i=0; i < user->endo_nbr; ++i) y1[i] = user->initial_values[i];
+ for (int i=0; i < 2*user->endo_nbr; ++i) y1[i+user->endo_nbr] = yy[i];
+ FirstDerivatives(y1, user->X, user->nb_row_x, user->params, user->steady_state, 1, NULL, user->row_ptr, user->col_ptr, user->val_ptr);
+ for (int* r=user->row_ptr; r < user->row_ptr+user->endo_nbr; r++)
+ {
+ int first_col = 0;
+ int ncol = 0;
+ int *pc = user->col_ptr + *r;
+ while(*(pc) < user->endo_nbr)
+ {
+ ++first_col;
+ ++pc;
+ }
+ while(pc < ((user->col_ptr)+*(r+1)))
+ {
+ if (*pc < 3*(user->endo_nbr))
+ {
+ ++ncol;
+ *pc -= user->endo_nbr;
+ }
+ ++pc;
+ }
+ ierr = MatSetValues(J,1,&row,ncol,user->col_ptr + *r + first_col,user->val_ptr + *r + first_col,INSERT_VALUES);
+ CHKERRQ(ierr);
+ ++row;
+ }
+ ys += user->endo_nbr;
+ ym -= user->endo_nbr;
+ }
+
+ /*
+ Compute function over locally owned part of the grid (interior points only)
+ */
+ while ( (ym >= user->endo_nbr) && (ys + 2*user->endo_nbr <= M) )
+ {
+ int it = ys/user->endo_nbr + 2;
+ FirstDerivatives(yy+ys-user->endo_nbr, user->X, user->nb_row_x, user->params, user->steady_state, it, NULL, user->row_ptr, user->col_ptr, user->val_ptr);
+ for (int* r=user->row_ptr; r < user->row_ptr+user->endo_nbr; r++)
+ {
+ int ncol = 0;
+ for(int *pc = user->col_ptr + *r; pc < user->col_ptr + *(r+1); ++pc)
+ if(*pc < 3*user->endo_nbr)
+ {
+ *pc += ys - user->endo_nbr;
+ ++ncol;
+ }
+ ierr = MatSetValues(J,1,&row,ncol,user->col_ptr + *r,user->val_ptr + *r,INSERT_VALUES);
+ CHKERRQ(ierr);
+ ++row;
+ }
+ ys += user->endo_nbr;
+ ym -= user->endo_nbr;
+ }
+
+
+ if ( (ym >= user->endo_nbr) && (ys + 2*user->endo_nbr >= M) )
+ {
+ int it = ys/user->endo_nbr + 1;
+ PetscReal *y1 = new PetscReal[3*user->endo_nbr];
+ for (int i=0; i < 2*user->endo_nbr; ++i) y1[i] = yy[ys+i-user->endo_nbr];
+ for (int i=0; i < user->endo_nbr; ++i) y1[i+2*user->endo_nbr] = user->terminal_values[i];
+ FirstDerivatives(y1, user->X, user->nb_row_x, user->params, user->steady_state, it, NULL, user->row_ptr, user->col_ptr, user->val_ptr);
+ for (int* r=user->row_ptr; r < user->row_ptr+user->endo_nbr; r++)
+ {
+ int *pc = user->col_ptr + *r;
+ int ncol = 0;
+ while((*pc < 2*user->endo_nbr) && (pc < user->col_ptr + *(r+1)))
+ {
+ ++ncol;
+ *pc += ys - user->endo_nbr;
+ ++pc;
+ }
+ ierr = MatSetValues(J,1,&row,ncol,user->col_ptr + *r,user->val_ptr + *r,INSERT_VALUES);
+ CHKERRQ(ierr);
+ ++row;
+ }
+ }
+
+ /*
+ Restore vectors
+ */
+ ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+ DMDAVecRestoreArray(da,ylocal,&yy);
+ DMRestoreLocalVector(da,&ylocal);
+ ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+ return(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "Monitor"
+/*
+ Monitor - Optional user-defined monitoring routine that views the
+ current iterate with an x-window plot. Set by SNESMonitorSet().
+
+ Input Parameters:
+ snes - the SNES context
+ its - iteration number
+ norm - 2-norm function value (may be estimated)
+ ctx - optional user-defined context for private data for the
+ monitor routine, as set by SNESMonitorSet()
+
+ Note:
+ See the manpage for PetscViewerDrawOpen() for useful runtime options,
+ such as -nox to deactivate all x-window output.
+ */
+PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
+{
+ PetscErrorCode ierr;
+ MonitorCtx *monP = (MonitorCtx*) ctx;
+ Vec x;
+
+ PetscFunctionBeginUser;
+ ierr = PetscPrintf(PETSC_COMM_WORLD,"iter = %D,SNES Function norm %g\n",its,(double)fnorm);CHKERRQ(ierr);
+ ierr = SNESGetSolution(snes,&x);CHKERRQ(ierr);
+ ierr = VecView(x,monP->viewer);CHKERRQ(ierr);
+ PetscFunctionReturn(0);
+}
+
+/* ------------------------------------------------------------------- */
+#undef __FUNCT__
+#define __FUNCT__ "PreCheck"
+/*
+ PreCheck - Optional user-defined routine that checks the validity of
+ candidate steps of a line search method. Set by SNESLineSearchSetPreCheck().
+
+ Input Parameters:
+ snes - the SNES context
+ xcurrent - current solution
+ y - search direction and length
+
+ Output Parameters:
+ y - proposed step (search direction and length) (possibly changed)
+ changed_y - tells if the step has changed or not
+ */
+PetscErrorCode PreCheck(SNESLineSearch linesearch,Vec xcurrent,Vec y, PetscBool *changed_y, void * ctx)
+{
+ PetscFunctionBeginUser;
+ *changed_y = PETSC_FALSE;
+ PetscFunctionReturn(0);
+}
+
+/* ------------------------------------------------------------------- */
+#undef __FUNCT__
+#define __FUNCT__ "PostCheck"
+/*
+ PostCheck - Optional user-defined routine that checks the validity of
+ candidate steps of a line search method. Set by SNESLineSearchSetPostCheck().
+
+ Input Parameters:
+ snes - the SNES context
+ ctx - optional user-defined context for private data for the
+ monitor routine, as set by SNESLineSearchSetPostCheck()
+ xcurrent - current solution
+ y - search direction and length
+ x - the new candidate iterate
+
+ Output Parameters:
+ y - proposed step (search direction and length) (possibly changed)
+ x - current iterate (possibly modified)
+
+ */
+PetscErrorCode PostCheck(SNESLineSearch linesearch,Vec xcurrent,Vec y,Vec x,PetscBool *changed_y,PetscBool *changed_x, void * ctx)
+{
+ PetscErrorCode ierr;
+ PetscInt i,iter,xs,xm;
+ StepCheckCtx *check;
+ AppCtx *user;
+ PetscScalar *xa,*xa_last,tmp;
+ PetscReal rdiff;
+ DM da;
+ SNES snes;
+
+ PetscFunctionBeginUser;
+ *changed_x = PETSC_FALSE;
+ *changed_y = PETSC_FALSE;
+
+ ierr = SNESLineSearchGetSNES(linesearch, &snes);CHKERRQ(ierr);
+ check = (StepCheckCtx*)ctx;
+ user = check->user;
+ ierr = SNESGetIterationNumber(snes,&iter);CHKERRQ(ierr);
+ ierr = SNESLineSearchGetPreCheck(linesearch, NULL, (void**)&check);CHKERRQ(ierr);
+
+ /* iteration 1 indicates we are working on the second iteration */
+ if (iter > 0) {
+ da = user->da;
+ ierr = PetscPrintf(PETSC_COMM_WORLD,"Checking candidate step at iteration %D with tolerance %g\n",iter,(double)check->tolerance);CHKERRQ(ierr);
+
+ /* Access local array data */
+ ierr = DMDAVecGetArray(da,check->last_step,&xa_last);CHKERRQ(ierr);
+ ierr = DMDAVecGetArray(da,x,&xa);CHKERRQ(ierr);
+ ierr = DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr);
+
+ /*
+ If we fail the user-defined check for validity of the candidate iterate,
+ then modify the iterate as we like. (Note that the particular modification
+ below is intended simply to demonstrate how to manipulate this data, not
+ as a meaningful or appropriate choice.)
+ */
+ for (i=xs; i<xs+xm; i++) {
+ if (!PetscAbsScalar(xa[i])) rdiff = 2*check->tolerance;
+ else rdiff = PetscAbsScalar((xa[i] - xa_last[i])/xa[i]);
+ if (rdiff > check->tolerance) {
+ tmp = xa[i];
+ xa[i] = .5*(xa[i] + xa_last[i]);
+ *changed_x = PETSC_TRUE;
+ ierr = PetscPrintf(PETSC_COMM_WORLD," Altering entry %D: x=%g, x_last=%g, diff=%g, x_new=%g\n",
+ i,(double)PetscAbsScalar(tmp),(double)PetscAbsScalar(xa_last[i]),(double)rdiff,(double)PetscAbsScalar(xa[i]));CHKERRQ(ierr);
+ }
+ }
+ ierr = DMDAVecRestoreArray(da,check->last_step,&xa_last);CHKERRQ(ierr);
+ ierr = DMDAVecRestoreArray(da,x,&xa);CHKERRQ(ierr);
+ }
+ ierr = VecCopy(x,check->last_step);CHKERRQ(ierr);
+ PetscFunctionReturn(0);
+}
+
+
+/* ------------------------------------------------------------------- */
+#undef __FUNCT__
+#define __FUNCT__ "PostSetSubKSP"
+/*
+ PostSetSubKSP - Optional user-defined routine that reset SubKSP options when hierarchical bjacobi PC is used
+ e.g,
+ mpiexec -n 8 ./ex3 -nox -n 10000 -ksp_type fgmres -pc_type bjacobi -pc_bjacobi_blocks 4 -sub_ksp_type gmres -sub_ksp_max_it 3 -post_setsubksp -sub_ksp_rtol 1.e-16
+ Set by SNESLineSearchSetPostCheck().
+
+ Input Parameters:
+ linesearch - the LineSearch context
+ xcurrent - current solution
+ y - search direction and length
+ x - the new candidate iterate
+
+ Output Parameters:
+ y - proposed step (search direction and length) (possibly changed)
+ x - current iterate (possibly modified)
+
+ */
+PetscErrorCode PostSetSubKSP(SNESLineSearch linesearch,Vec xcurrent,Vec y,Vec x,PetscBool *changed_y,PetscBool *changed_x, void * ctx)
+{
+ PetscErrorCode ierr;
+ SetSubKSPCtx *check;
+ PetscInt iter,its,sub_its,maxit;
+ KSP ksp,sub_ksp,*sub_ksps;
+ PC pc;
+ PetscReal ksp_ratio;
+ SNES snes;
+
+ PetscFunctionBeginUser;
+ ierr = SNESLineSearchGetSNES(linesearch, &snes);CHKERRQ(ierr);
+ check = (SetSubKSPCtx*)ctx;
+ ierr = SNESGetIterationNumber(snes,&iter);CHKERRQ(ierr);
+ ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr);
+ ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
+ ierr = PCBJacobiGetSubKSP(pc,NULL,NULL,&sub_ksps);CHKERRQ(ierr);
+ sub_ksp = sub_ksps[0];
+ ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr); /* outer KSP iteration number */
+ ierr = KSPGetIterationNumber(sub_ksp,&sub_its);CHKERRQ(ierr); /* inner KSP iteration number */
+
+ if (iter) {
+ ierr = PetscPrintf(PETSC_COMM_WORLD," ...PostCheck snes iteration %D, ksp_it %d %d, subksp_it %d\n",iter,check->its0,its,sub_its);CHKERRQ(ierr);
+ ksp_ratio = ((PetscReal)(its))/check->its0;
+ maxit = (PetscInt)(ksp_ratio*sub_its + 0.5);
+ if (maxit < 2) maxit = 2;
+ ierr = KSPSetTolerances(sub_ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,maxit);CHKERRQ(ierr);
+ ierr = PetscPrintf(PETSC_COMM_WORLD," ...ksp_ratio %g, new maxit %d\n\n",ksp_ratio,maxit);CHKERRQ(ierr);
+ }
+ check->its0 = its; /* save current outer KSP iteration number */
+ PetscFunctionReturn(0);
+}
+
+/* ------------------------------------------------------------------- */
+/*
+ MatrixFreePreconditioner - This routine demonstrates the use of a
+ user-provided preconditioner. This code implements just the null
+ preconditioner, which of course is not recommended for general use.
+
+ Input Parameters:
++ pc - preconditioner
+- x - input vector
+
+ Output Parameter:
+. y - preconditioned vector
+*/
+PetscErrorCode MatrixFreePreconditioner(PC pc,Vec x,Vec y)
+{
+ PetscErrorCode ierr;
+ ierr = VecCopy(x,y);CHKERRQ(ierr);
+ return 0;
+}
diff --git a/preprocessor/DynamicModel.cc b/preprocessor/DynamicModel.cc
index 1a6f954b3db23095dd993c4add3a3078f17bc03c..97616691ae29afe443a1aa5e3304065a6f36cc97 100644
--- a/preprocessor/DynamicModel.cc
+++ b/preprocessor/DynamicModel.cc
@@ -4467,7 +4467,7 @@ DynamicModel::writeResidualsC(const string &basename, bool cuda) const
deriv_node_temp_terms_t tef_terms;
ostringstream model_output; // Used for storing model equations
- writeModelEquations(model_output, oCDynamicModel);
+ writeModelEquations(model_output, oCDynamic2Model);
mDynamicModelFile << " double lhs, rhs;" << endl
<< endl
@@ -4540,6 +4540,106 @@ DynamicModel::writeFirstDerivativesC(const string &basename, bool cuda) const
// using compressed sparse row format (CSR)
void
+DynamicModel::writeFirstDerivativesC_csr(const string &basename, bool cuda) const
+{
+ string filename = basename + "_first_derivatives.c";
+ ofstream mDynamicModelFile, mDynamicMexFile;
+
+ mDynamicModelFile.open(filename.c_str(), ios::out | ios::binary);
+ if (!mDynamicModelFile.is_open())
+ {
+ cerr << "Error: Can't open file " << filename << " for writing" << endl;
+ exit(EXIT_FAILURE);
+ }
+ mDynamicModelFile << "/*" << endl
+ << " * " << filename << " : Computes first order derivatives of the model for Dynare" << endl
+ << " *" << endl
+ << " * Warning : this file is generated automatically by Dynare" << endl
+ << " * from model " << basename << "(.mod)" << endl
+ << " */" << endl
+ << "#include <math.h>" << endl;
+
+ mDynamicModelFile << "#include <stdlib.h>" << endl;
+
+ mDynamicModelFile << "#define max(a, b) (((a) > (b)) ? (a) : (b))" << endl
+ << "#define min(a, b) (((a) > (b)) ? (b) : (a))" << endl;
+
+ // Write function definition if oPowerDeriv is used
+ writePowerDerivCHeader(mDynamicModelFile);
+
+ mDynamicModelFile << "void FirstDerivatives(const double *y, double *x, int nb_row_x, double *params, double *steady_state, int it_, double *residual, int *row_ptr, int *col_ptr, double *value)" << endl
+ << "{" << endl;
+
+ int cols_nbr = 3*symbol_table.endo_nbr() + symbol_table.exo_nbr() + symbol_table.exo_det_nbr();
+ // this is always empty here, but needed by d1->writeOutput
+ deriv_node_temp_terms_t tef_terms;
+
+
+ // Indexing derivatives in column order
+ vector<derivative> D;
+ for (first_derivatives_t::const_iterator it = first_derivatives.begin();
+ it != first_derivatives.end(); it++)
+ {
+ int eq = it->first.first;
+ int dynvar = it->first.second;
+ int lag = getLagByDerivID(dynvar);
+ int symb = getSymbIDByDerivID(dynvar);
+ SymbolType type = getTypeByDerivID(dynvar);
+ int col_id;
+ switch(type)
+ {
+ case eEndogenous:
+ col_id = symb+(lag+1)*symbol_table.endo_nbr();
+ break;
+ case eExogenous:
+ col_id = symb+3*symbol_table.endo_nbr();
+ break;
+ case eExogenousDet:
+ col_id = symb+3*symbol_table.endo_nbr()+symbol_table.exo_nbr();
+ break;
+ default:
+ std::cerr << "This case shouldn't happen" << std::endl;
+ exit(1);
+ }
+ derivative deriv(col_id + eq*cols_nbr,col_id,eq,it->second);
+ D.push_back(deriv);
+ }
+ sort(D.begin(), D.end(), derivative_less_than() );
+
+ // writing sparse Jacobian
+ vector<int> row_ptr(equations.size());
+ fill(row_ptr.begin(),row_ptr.end(),0.0);
+ int k = 0;
+ for(vector<derivative>::const_iterator it = D.begin(); it != D.end(); ++it)
+ {
+ row_ptr[it->row_nbr]++;
+ mDynamicModelFile << "col_ptr[" << k << "] "
+ << "=" << it->col_nbr << ";" << endl;
+ mDynamicModelFile << "value[" << k << "] = ";
+ // oCstaticModel makes reference to the static variables
+ it->value->writeOutput(mDynamicModelFile, oCDynamic2Model, temporary_terms, tef_terms);
+ mDynamicModelFile << ";" << endl;
+ k++;
+ }
+
+ // row_ptr must point to the relative address of the first element of the row
+ int cumsum = 0;
+ mDynamicModelFile << "int row_ptr_data[" << row_ptr.size() + 1 << "] = { 0";
+ for (vector<int>::iterator it=row_ptr.begin(); it != row_ptr.end(); ++it)
+ {
+ cumsum += *it;
+ mDynamicModelFile << ", " << cumsum;
+ }
+ mDynamicModelFile << "};" << endl
+ << "int i;" << endl
+ << "for (i=0; i < " << row_ptr.size() + 1 << "; i++) row_ptr[i] = row_ptr_data[i];" << endl;
+ mDynamicModelFile << "}" << endl;
+
+ // writePowerDeriv(mDynamicModelFile, true);
+ mDynamicModelFile.close();
+
+}
+ void
DynamicModel::writeSecondDerivativesC_csr(const string &basename, bool cuda) const
{
diff --git a/preprocessor/DynamicModel.hh b/preprocessor/DynamicModel.hh
index a9655220e11aec9930429fb77025e1ee9ea11e5e..61b758e229c149c158ccbc7bf926c961285cb718 100644
--- a/preprocessor/DynamicModel.hh
+++ b/preprocessor/DynamicModel.hh
@@ -480,6 +480,8 @@ public:
void writeResidualsC(const string &basename, bool cuda) const;
//! Writes C file containing first order derivatives of model evaluated at steady state
void writeFirstDerivativesC(const string &basename, bool cuda) const;
+ //! Writes C file containing first order derivatives of model evaluated at steady state (conpressed sparse column)
+ void writeFirstDerivativesC_csr(const string &basename, bool cuda) const;
//! Writes C file containing second order derivatives of model evaluated at steady state (compressed sparse column)
void writeSecondDerivativesC_csr(const string &basename, bool cuda) const;
//! Writes C file containing third order derivatives of model evaluated at steady state (compressed sparse column)
diff --git a/preprocessor/ExprNode.cc b/preprocessor/ExprNode.cc
index b9c4405528ed4f596f0f159f00d71199d57e251e..d2b779a449caa423998f8c9efcc0d03d6eeb6f1c 100644
--- a/preprocessor/ExprNode.cc
+++ b/preprocessor/ExprNode.cc
@@ -659,6 +659,10 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
i = datatree.getDynJacobianCol(datatree.getDerivID(symb_id, lag)) + ARRAY_SUBSCRIPT_OFFSET(output_type);
output << "y" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
break;
+ case oCDynamic2Model:
+ i = symb_id + (lag+1)*datatree.symbol_table.endo_nbr() + ARRAY_SUBSCRIPT_OFFSET(output_type);
+ output << "y" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
+ break;
case oCStaticModel:
case oMatlabStaticModel:
case oMatlabStaticModelSparse:
@@ -710,6 +714,7 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
output << "x(it_, " << i << ")";
break;
case oCDynamicModel:
+ case oCDynamic2Model:
if (lag == 0)
output << "x[it_+" << i << "*nb_row_x]";
else if (lag > 0)
@@ -755,6 +760,7 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
output << "x(it_, " << i << ")";
break;
case oCDynamicModel:
+ case oCDynamic2Model:
if (lag == 0)
output << "x[it_+" << i << "*nb_row_x]";
else if (lag > 0)
diff --git a/preprocessor/ExprNode.hh b/preprocessor/ExprNode.hh
index 0500ff6d501edf45b2bf713884520c38fffc8f6d..394bf586a0eca6e7afb1adab54fb41cc16510d0d 100644
--- a/preprocessor/ExprNode.hh
+++ b/preprocessor/ExprNode.hh
@@ -65,6 +65,7 @@ enum ExprNodeOutputType
oMatlabStaticModelSparse, //!< Matlab code, static block decomposed model
oMatlabDynamicModelSparse, //!< Matlab code, dynamic block decomposed model
oCDynamicModel, //!< C code, dynamic model
+ oCDynamic2Model, //!< C code, dynamic model, alternative numbering of endogenous variables
oCStaticModel, //!< C code, static model
oMatlabOutsideModel, //!< Matlab code, outside model block (for example in initval)
oLatexStaticModel, //!< LaTeX code, static model
@@ -87,6 +88,7 @@ enum ExprNodeOutputType
|| (output_type) == oSteadyStateFile)
#define IS_C(output_type) ((output_type) == oCDynamicModel \
+ || (output_type) == oCDynamic2Model \
|| (output_type) == oCStaticModel \
|| (output_type) == oCDynamicSteadyStateOperator \
|| (output_type) == oCSteadyStateFile)
diff --git a/preprocessor/ModFile.cc b/preprocessor/ModFile.cc
index a3de63a7ed4ee3cfdbeb0a2b90a6bb4bceeb1a5d..1e8fb25b2813566dcfc47d0ef200716148f683d1 100644
--- a/preprocessor/ModFile.cc
+++ b/preprocessor/ModFile.cc
@@ -972,7 +972,7 @@ ModFile::writeExternalFilesCC(const string &basename, FileOutputType output) con
// dynamic_model.writeResidualsC(basename, cuda);
// dynamic_model.writeParamsDerivativesFileC(basename, cuda);
dynamic_model.writeResidualsC(basename, cuda);
- dynamic_model.writeFirstDerivativesC(basename, cuda);
+ dynamic_model.writeFirstDerivativesC_csr(basename, cuda);
if (output == second)
dynamic_model.writeSecondDerivativesC_csr(basename, cuda);