expand extended preprocessor + first implementation of Petsc interface
Michel Juillard
parent
154980ec8c
commit
c9f771973d
|
|
@ -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++
|
||||
|
|
@ -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;
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
}
|
||||
|
|
@ -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
|
||||
{
|
||||
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
|
|
|
|||
Loading…
Reference in New Issue