From c9f771973df86419f4a0aa0f6dbeb944c60d84ac Mon Sep 17 00:00:00 2001 From: Michel Juillard Date: Thu, 27 Aug 2015 10:00:27 +0200 Subject: [PATCH] expand extended preprocessor + first implementation of Petsc interface --- others/cpp/tests/Makefile_petsc | 25 ++ others/cpp/tests/example1.mod | 20 +- others/cpp/tests/snes_1.cc | 685 ++++++++++++++++++++++++++++++++ preprocessor/DynamicModel.cc | 102 ++++- preprocessor/DynamicModel.hh | 2 + preprocessor/ExprNode.cc | 6 + preprocessor/ExprNode.hh | 2 + preprocessor/ModFile.cc | 2 +- 8 files changed, 832 insertions(+), 12 deletions(-) create mode 100644 others/cpp/tests/Makefile_petsc create mode 100644 others/cpp/tests/snes_1.cc diff --git a/others/cpp/tests/Makefile_petsc b/others/cpp/tests/Makefile_petsc new file mode 100644 index 000000000..87b829ce5 --- /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 1249be4ba..ec43f0871 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 000000000..a3ffb0c0b --- /dev/null +++ b/others/cpp/tests/snes_1.cc @@ -0,0 +1,685 @@ +#include +#include "../dynare_cpp_driver.hh" +#include +#include +#include +#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 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 -pc_type + */ + 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; itolerance; + 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 1a6f954b3..97616691a 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 " << endl; + + mDynamicModelFile << "#include " << 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 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 row_ptr(equations.size()); + fill(row_ptr.begin(),row_ptr.end(),0.0); + int k = 0; + for(vector::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::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 a9655220e..61b758e22 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 b9c440552..d2b779a44 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 0500ff6d5..394bf586a 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 a3de63a7e..1e8fb25b2 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);